12#include "cantera/oneD/refine.h"
32 for (
size_t n = 0; n <
nDomains(); n++) {
46 "To be removed after Cantera 3.0; superseded by "
47 "Sim1D::Sim1D(vector<shared_ptr<Domain1D>>&).");
52 for (
size_t n = 0; n <
nDomains(); n++) {
64 for (
size_t dom=0; dom<
nDomains(); dom++) {
67 for (
size_t comp=0; comp<ncomp; comp++) {
79 "Index out of bounds: {} > {}", iloc,
m_state->size());
80 (*m_state)[iloc] =
value;
83double Sim1D::value(
size_t dom,
size_t comp,
size_t localPoint)
const
87 "Index out of bounds: {} > {}", iloc,
m_state->size());
91double Sim1D::workValue(
size_t dom,
size_t comp,
size_t localPoint)
const
95 "Index out of bounds: {} > {}", iloc,
m_state->size());
100 const vector<double>& pos,
const vector<double>& values)
102 if (pos.front() != 0.0 || pos.back() != 1.0) {
104 "`pos` vector must span the range [0, 1]. Got a vector spanning "
105 "[{}, {}] instead.", pos.front(), pos.back());
108 double z0 = d.zmin();
109 double z1 = d.zmax();
110 for (
size_t n = 0; n < d.
nPoints(); n++) {
112 double frac = (zpt - z0)/(z1 - z0);
119 const string& desc,
int loglevel)
122 "To be removed after Cantera 3.0; use version without loglevel instead.");
123 save(fname,
id, desc,
true);
125 writelog(
"Solution saved to file '{}' as entry '{}'.\n", fname,
id);
129void Sim1D::save(
const string& fname,
const string& name,
const string& desc,
130 bool overwrite,
int compression,
const string& basis)
132 size_t dot = fname.find_last_of(
".");
134 if (extension ==
"csv") {
135 for (
auto dom :
m_dom) {
136 auto arr = dom->asArray(
m_state->data() + dom->loc());
137 if (dom->size() > 1) {
138 arr->writeEntry(fname, overwrite, basis);
146 "Species basis '{}' not implemented for HDF5 or YAML output.", basis);
148 if (extension ==
"h5" || extension ==
"hdf" || extension ==
"hdf5") {
150 for (
auto dom :
m_dom) {
151 auto arr = dom->asArray(
m_state->data() + dom->loc());
152 arr->writeEntry(fname, name, dom->id(), overwrite, compression);
156 if (extension ==
"yaml" || extension ==
"yml") {
159 if (std::ifstream(fname).good()) {
164 for (
auto dom :
m_dom) {
165 auto arr = dom->asArray(
m_state->data() + dom->loc());
166 arr->writeEntry(data, name, dom->id(), overwrite);
170 std::ofstream out(fname);
171 out << data.toYamlString();
175 throw CanteraError(
"Sim1D::save",
"Unsupported file format '{}'.", extension);
179 const string& desc,
int loglevel)
182 "To be removed after Cantera 3.0; use version without loglevel instead.");
185 writelog(
"Solution saved to file '{}' as entry '{}'.\n", fname,
id);
190 const string& desc,
bool overwrite,
int compression)
192 vector<double> res(
m_state->size(), -999);
196 vector<double> backup(*
m_state);
198 save(fname, name, desc, overwrite, compression);
205AnyMap legacyH5(shared_ptr<SolutionArray> arr,
const AnyMap& header={})
207 auto meta = arr->meta();
210 map<string, string> meta_pairs = {
211 {
"type",
"Domain1D_type"},
213 {
"emissivity-left",
"emissivity_left"},
214 {
"emissivity-right",
"emissivity_right"},
216 for (
const auto& [newName, oldName] : meta_pairs) {
217 if (meta.hasKey(oldName)) {
218 out[newName] = meta[oldName];
222 map<string, string> tol_pairs = {
223 {
"transient-abstol",
"transient_abstol"},
224 {
"steady-abstol",
"steady_abstol"},
225 {
"transient-reltol",
"transient_reltol"},
226 {
"steady-reltol",
"steady_reltol"},
228 for (
const auto& [newName, oldName] : tol_pairs) {
229 if (meta.hasKey(oldName)) {
230 out[
"tolerances"][newName] = meta[oldName];
234 if (meta.hasKey(
"phase")) {
235 out[
"phase"][
"name"] = meta[
"phase"][
"name"];
236 out[
"phase"][
"source"] = meta[
"phase"][
"source"];
239 if (arr->size() <= 1) {
243 map<string, string> header_pairs = {
244 {
"transport-model",
"transport_model"},
245 {
"radiation-enabled",
"radiation_enabled"},
246 {
"energy-enabled",
"energy_enabled"},
247 {
"Soret-enabled",
"soret_enabled"},
248 {
"species-enabled",
"species_enabled"},
250 for (
const auto& [newName, oldName] : header_pairs) {
251 if (header.hasKey(oldName)) {
252 out[newName] = header[oldName];
256 map<string, string> refiner_pairs = {
262 {
"max-points",
"max_grid_points"},
264 for (
const auto& [newName, oldName] : refiner_pairs) {
265 if (header.hasKey(oldName)) {
266 out[
"refine-criteria"][newName] = header[oldName];
270 if (header.hasKey(
"fixed_temperature")) {
271 double temp = header.
getDouble(
"fixed_temperature", -1.);
272 auto profile = arr->getComponent(
"T").as<vector<double>>();
274 while (profile[ix] <= temp && ix < arr->size()) {
278 auto grid = arr->getComponent(
"grid").as<vector<double>>();
279 out[
"fixed-point"][
"location"] = grid[ix - 1];
280 out[
"fixed-point"][
"temperature"] = temp;
292 "To be removed after Cantera 3.0; use version without loglevel instead.");
298 size_t dot = fname.find_last_of(
".");
300 if (extension ==
"xml") {
302 "Restoring from XML is no longer supported.");
305 if (extension ==
"h5" || extension ==
"hdf" || extension ==
"hdf5") {
306 map<string, shared_ptr<SolutionArray>> arrs;
309 for (
auto dom :
m_dom) {
312 arr->readEntry(fname, name, dom->id());
315 "Encountered exception when reading entry '{}' from '{}':\n{}",
318 dom->resize(dom->nComponents(), arr->size());
319 if (!header.
hasKey(
"generator")) {
320 arr->meta() = legacyH5(arr, header);
322 arrs[dom->id()] = arr;
326 for (
auto dom :
m_dom) {
328 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
331 "Encountered exception when restoring domain '{}' from HDF:\n{}",
336 }
else if (extension ==
"yaml" || extension ==
"yml") {
338 map<string, shared_ptr<SolutionArray>> arrs;
341 for (
auto dom :
m_dom) {
344 arr->readEntry(root, name, dom->id());
347 "Encountered exception when reading entry '{}' from '{}':\n{}",
350 dom->resize(dom->nComponents(), arr->size());
351 arrs[dom->id()] = arr;
355 for (
auto dom :
m_dom) {
357 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
360 "Encountered exception when restoring domain '{}' from YAML:\n{}",
367 "Unknown file extension '{}'; supported extensions include "
368 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
376 for (
size_t n = 0; n < np; n++) {
384 "To be removed after Cantera 3.0; replaced by 'show'.");
391 "To be removed after Cantera 3.0; replaced by 'show'.");
397 for (
size_t n = 0; n <
nDomains(); n++) {
406 for (
size_t n = 0; n <
nDomains(); n++) {
409 +
" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
418 throw CanteraError(
"Sim1D::restoreTimeSteppingSolution",
419 "No successful time steps taken on this grid.");
428 "No successful steady state solution");
431 for (
size_t n = 0; n <
nDomains(); n++) {
437void Sim1D::getInitialSoln()
439 for (
size_t n = 0; n <
nDomains(); n++) {
446 for (
size_t n = 0; n <
nDomains(); n++) {
451void Sim1D::setTimeStep(
double stepsize,
size_t n,
const int* tsteps)
455 for (
size_t i = 0; i < n; i++) {
466 }
else if (m > -10) {
470 "ERROR: OneDim::solve returned m = {}", m);
474void Sim1D::solve(
int loglevel,
bool refine_grid)
481 while (new_points > 0) {
487 writeline(
'.', 78,
true,
true);
493 debuglog(
"Attempt Newton solution of steady-state problem...", loglevel);
500 for (
size_t mm = 1; mm <
nDomains(); mm+=2) {
513 save(
"debug_sim1d.yaml",
"debug",
514 "After successful Newton solve");
518 "After successful Newton solve");
524 save(
"debug_sim1d.yaml",
"debug",
525 "After unsuccessful Newton solve");
529 "After unsuccessful Newton solve");
532 writelog(
"Take {} timesteps ", nsteps);
537 save(
"debug_sim1d.yaml",
"debug",
"After timestepping");
541 "After timestepping");
545 writelog(
" {:10.4g} {:10.4g}\n", dt,
554 dt = std::min(dt,
m_tmax);
558 writeline(
'.', 78,
true,
true);
565 new_points =
refine(loglevel);
571 if (new_points && loglevel > 6) {
572 save(
"debug_sim1d.yaml",
"debug",
"After regridding");
574 if (new_points && loglevel > 7) {
579 debuglog(
"grid refinement disabled.\n", loglevel);
587 int ianalyze, np = 0;
588 vector<double> znew, xnew;
589 vector<size_t> dsize;
594 for (
size_t n = 0; n <
nDomains(); n++) {
602 ianalyze = r.analyze(d.grid().size(), d.grid().data(),
612 np += r.nNewPoints();
617 size_t nstart = znew.size();
618 for (
size_t m = 0; m < npnow; m++) {
619 if (r.keepPoint(m)) {
621 znew.push_back(d.grid(m));
624 for (
size_t i = 0; i < comp; i++) {
625 xnew.push_back(
value(n, i, m));
631 if (r.newPointNeeded(m) && m + 1 < npnow) {
633 double zmid = 0.5*(d.grid(m) + d.grid(m+1));
634 znew.push_back(zmid);
639 for (
size_t i = 0; i < comp; i++) {
640 double xmid = 0.5*(
value(n, i, m) +
value(n, i, m+1));
641 xnew.push_back(xmid);
646 writelog(
"refine: discarding point at {}\n", d.grid(m));
650 dsize.push_back(znew.size() - nstart);
657 size_t gridstart = 0, gridsize;
658 for (
size_t n = 0; n <
nDomains(); n++) {
664 gridstart += gridsize;
677 vector<double> znew, xnew;
679 double z1 = 0.0, z2 = 0.0;
680 vector<size_t> dsize;
682 for (
size_t n = 0; n <
nDomains(); n++) {
685 size_t mfixed =
npos;
690 size_t nstart = znew.size();
691 if (d_free && d_free->
isFree()) {
692 for (
size_t m = 0; m < npnow - 1; m++) {
693 bool fixedpt =
false;
697 double thresh = min(1., 1.e-1 * (t2 - t1));
700 if (fabs(t - t1) <= thresh) {
703 }
else if (fabs(t2 - t) <= thresh) {
706 }
else if ((t1 < t) && (t < t2)) {
708 zfixed = (z1 - z2) / (t1 - t2) * (t - t2) + z2;
721 for (
size_t m = 0; m < npnow; m++) {
723 znew.push_back(d.grid(m));
726 for (
size_t i = 0; i < comp; i++) {
727 xnew.push_back(
value(n, i, m));
731 znew.push_back(zfixed);
733 double interp_factor = (zfixed - z2) / (z1 - z2);
736 for (
size_t i = 0; i < comp; i++) {
737 double xmid = interp_factor*(
value(n, i, m) -
value(n, i, m+1)) +
value(n,i,m+1);
738 xnew.push_back(xmid);
742 dsize.push_back(znew.size() - nstart);
748 size_t gridstart = 0;
749 for (
size_t n = 0; n <
nDomains(); n++) {
751 size_t gridsize = dsize[n];
753 gridstart += gridsize;
766 double t_fixed = std::numeric_limits<double>::quiet_NaN();
767 for (
size_t n = 0; n <
nDomains(); n++) {
779 double z_fixed = std::numeric_limits<double>::quiet_NaN();
780 for (
size_t n = 0; n <
nDomains(); n++) {
791 double slope,
double curve,
double prune)
797 for (
size_t n = 0; n <
nDomains(); n++) {
811 "Must specify domain to get criteria from");
821 for (
size_t n = 0; n <
nDomains(); n++) {
834 for (
size_t n = 0; n <
nDomains(); n++) {
852void Sim1D::evalSSJacobian()
859 for (
auto& D :
m_dom) {
860 D->forceFullUpdate(
true);
863 for (
auto& D :
m_dom) {
864 D->forceFullUpdate(
false);
870 for (
size_t i = 0; i <
size(); i++) {
871 size_t j1 = (i > bw) ? i - bw : 0;
872 size_t j2 = (i + bw >=
size()) ?
size() - 1: i + bw;
873 for (
size_t j = j1; j <= j2; j++) {
874 Jt(j,i) =
m_jac->value(i,j);
A map of string keys to values whose type can vary at runtime.
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
static void clearCachedFile(const string &filename)
Remove the specified file from the input cache if it is present.
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
A class for banded matrices, involving matrix inversion processes.
int solve(const double *const b, double *const x)
Solve the matrix problem Ax = b.
double & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Base class for one-dimensional domains.
size_t nComponents() const
Number of components at each grid point.
virtual void _finalize(const double *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
virtual string type() const
String indicating the domain implemented.
size_t nPoints() const
Number of grid points in this domain.
Refiner & refiner()
Return a reference to the grid refiner.
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
virtual void _getInitialSoln(double *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
virtual void show(std::ostream &s, const double *x)
Print the solution.
virtual void setupGrid(size_t n, const double *z)
called to set up initial grid, and after grid refinement
virtual double eval(double t) const
Evaluate the function.
void setOptions(int maxJacAge=5)
Set options.
Container class for multiple-domain 1D problems.
int solve(double *x0, double *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
int m_nsteps
Number of time steps taken in the current call to solve()
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
size_t size() const
Total solution vector length;.
void eval(size_t j, double *x, double *r, double rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
size_t nDomains() const
Number of domains.
std::tuple< string, size_t, string > component(size_t i)
Return the domain, local point index, and component name for the i-th component of the global solutio...
size_t bandwidth() const
Jacobian bandwidth.
shared_ptr< vector< double > > m_state
Solution vector.
unique_ptr< MultiJac > m_jac
Jacobian evaluator.
vector< Domain1D * > m_dom
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
double m_tmax
maximum timestep size
void setSteadyMode()
Prepare to solve the steady-state problem.
Domain1D & domain(size_t i) const
Return a reference to domain i.
MultiNewton & newton()
Return a reference to the Newton iterator.
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
size_t maxPoints() const
Returns the maximum number of points allowed in the domain.
vector< double > getCriteria()
Get the grid refinement criteria.
void setMaxPoints(int npmax)
Set the maximum number of points allowed in the domain.
void setCriteria(double ratio=10.0, double slope=0.8, double curve=0.8, double prune=-0.1)
Set grid refinement criteria.
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
void resize() override
Call after one or more grids has changed size, for example after being refined.
vector< double > m_xnew
a work array used to hold the residual or the new solution
void saveResidual(const string &fname, const string &id, const string &desc, int loglevel)
Save the residual of the current solution to a container file.
void setProfile(size_t dom, size_t comp, const vector< double > &pos, const vector< double > &values)
Specify a profile for one component of one domain.
double fixedTemperatureLocation()
Return location of the point where temperature is fixed.
vector< vector< double > > m_grid_last_ss
the grids for each domain after the last successful steady-state solve (stored before grid refinement...
void finalize()
Calls method _finalize in each domain.
void setValue(size_t dom, size_t comp, size_t localPoint, double value)
Set a single value in the solution vector.
int refine(int loglevel=0)
Refine the grid in all domains.
void show()
Show logging information on current solution for all domains.
double fixedTemperature()
Return temperature at the point used to fix the flame location.
vector< double > m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement)
void save(const string &fname, const string &id, const string &desc, int loglevel)
Save the current solution to a container file.
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
int setFixedTemperature(double t)
Add node for fixed temperature point of freely propagating flame.
AnyMap restore(const string &fname, const string &id, int loglevel)
Initialize the solution with a previously-saved solution.
void setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component for all domains.
void showSolution()
Show logging information on current solution for all domains.
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
vector< int > m_steps
array of number of steps to take before re-attempting the steady-state solution
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
void setFlatProfile(size_t dom, size_t comp, double v)
Set component 'comp' of domain 'dom' to value 'v' at all points.
double value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
vector< double > getRefineCriteria(int dom)
Get the grid refinement criteria.
Sim1D()
Default constructor.
void setGridMin(int dom, double gridmin)
Set the minimum grid spacing in the specified domain(s).
void setRefineCriteria(int dom=-1, double ratio=10.0, double slope=0.8, double curve=0.8, double prune=-0.1)
Set grid refinement criteria.
vector< double > m_xlast_ts
the solution vector after the last successful timestepping
static AnyMap readHeader(const string &fname, const string &name)
Read header information from a HDF container file.
static void writeHeader(const string &fname, const string &name, const string &desc, bool overwrite=false)
Write header data to a HDF container file.
static shared_ptr< SolutionArray > create(const shared_ptr< Solution > &sol, int size=0, const AnyMap &meta={})
Instantiate a new SolutionArray reference.
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
double m_tfixed
Temperature at the point used to fix the flame location.
double m_zfixed
Location of the point where temperature is fixed.
bool isFree() const
Retrieve flag indicating whether flow is freely propagating.
Header for a file containing miscellaneous numerical functions.
string toLowerCopy(const string &input)
Convert to lower case.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
double linearInterp(double x, const vector< double > &xpts, const vector< double > &fpts)
Linearly interpolate a function defined on a discrete grid.
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Contains declarations for string manipulation functions within Cantera.