12#include "cantera/oneD/refine.h"
32 for (
size_t n = 0; n <
nDomains(); n++) {
44 for (
size_t dom=0; dom<
nDomains(); dom++) {
47 for (
size_t comp=0; comp<ncomp; comp++) {
59 "Index out of bounds: {} > {}", iloc,
m_state->size());
60 (*m_state)[iloc] =
value;
63double Sim1D::value(
size_t dom,
size_t comp,
size_t localPoint)
const
67 "Index out of bounds: {} > {}", iloc,
m_state->size());
75 "Index out of bounds: {} > {}", iloc,
m_state->size());
80 const vector<double>& pos,
const vector<double>& values)
82 if (pos.front() != 0.0 || pos.back() != 1.0) {
84 "`pos` vector must span the range [0, 1]. Got a vector spanning "
85 "[{}, {}] instead.", pos.front(), pos.back());
90 for (
size_t n = 0; n < d.
nPoints(); n++) {
92 double frac = (zpt - z0)/(z1 - z0);
98void Sim1D::save(
const string& fname,
const string& name,
const string& desc,
99 bool overwrite,
int compression,
const string& basis)
101 size_t dot = fname.find_last_of(
".");
103 if (extension ==
"csv") {
104 for (
auto dom :
m_dom) {
105 auto arr = dom->asArray(
m_state->data() + dom->loc());
106 if (dom->size() > 1) {
107 arr->writeEntry(fname, overwrite, basis);
115 "Species basis '{}' not implemented for HDF5 or YAML output.", basis);
117 if (extension ==
"h5" || extension ==
"hdf" || extension ==
"hdf5") {
119 for (
auto dom :
m_dom) {
120 auto arr = dom->asArray(
m_state->data() + dom->loc());
121 arr->writeEntry(fname, name, dom->id(), overwrite, compression);
125 if (extension ==
"yaml" || extension ==
"yml") {
128 if (std::ifstream(fname).good()) {
133 for (
auto dom :
m_dom) {
134 auto arr = dom->asArray(
m_state->data() + dom->loc());
135 arr->writeEntry(data, name, dom->id(), overwrite);
139 std::ofstream out(fname);
140 out << data.toYamlString();
144 throw CanteraError(
"Sim1D::save",
"Unsupported file format '{}'.", extension);
148 const string& desc,
bool overwrite,
int compression)
150 vector<double> res(
m_state->size(), -999);
154 vector<double> backup(*
m_state);
156 save(fname, name, desc, overwrite, compression);
163AnyMap legacyH5(shared_ptr<SolutionArray> arr,
const AnyMap& header={})
165 auto meta = arr->meta();
168 map<string, string> meta_pairs = {
169 {
"type",
"Domain1D_type"},
171 {
"emissivity-left",
"emissivity_left"},
172 {
"emissivity-right",
"emissivity_right"},
174 for (
const auto& [newName, oldName] : meta_pairs) {
175 if (meta.hasKey(oldName)) {
176 out[newName] = meta[oldName];
180 map<string, string> tol_pairs = {
181 {
"transient-abstol",
"transient_abstol"},
182 {
"steady-abstol",
"steady_abstol"},
183 {
"transient-reltol",
"transient_reltol"},
184 {
"steady-reltol",
"steady_reltol"},
186 for (
const auto& [newName, oldName] : tol_pairs) {
187 if (meta.hasKey(oldName)) {
188 out[
"tolerances"][newName] = meta[oldName];
192 if (meta.hasKey(
"phase")) {
193 out[
"phase"][
"name"] = meta[
"phase"][
"name"];
194 out[
"phase"][
"source"] = meta[
"phase"][
"source"];
197 if (arr->size() <= 1) {
201 map<string, string> header_pairs = {
202 {
"transport-model",
"transport_model"},
203 {
"radiation-enabled",
"radiation_enabled"},
204 {
"energy-enabled",
"energy_enabled"},
205 {
"Soret-enabled",
"soret_enabled"},
207 for (
const auto& [newName, oldName] : header_pairs) {
208 if (header.hasKey(oldName)) {
209 out[newName] = header[oldName];
213 map<string, string> refiner_pairs = {
219 {
"max-points",
"max_grid_points"},
221 for (
const auto& [newName, oldName] : refiner_pairs) {
222 if (header.hasKey(oldName)) {
223 out[
"refine-criteria"][newName] = header[oldName];
227 if (header.hasKey(
"fixed_temperature")) {
228 double temp = header.
getDouble(
"fixed_temperature", -1.);
229 auto profile = arr->getComponent(
"T").as<vector<double>>();
231 while (profile[ix] <= temp && ix < arr->size()) {
235 auto grid = arr->getComponent(
"grid").as<vector<double>>();
236 out[
"fixed-point"][
"location"] = grid[ix - 1];
237 out[
"fixed-point"][
"temperature"] = temp;
248 size_t dot = fname.find_last_of(
".");
250 if (extension ==
"xml") {
252 "Restoring from XML is no longer supported.");
255 if (extension ==
"h5" || extension ==
"hdf" || extension ==
"hdf5") {
256 map<string, shared_ptr<SolutionArray>> arrs;
259 for (
auto dom :
m_dom) {
262 arr->readEntry(fname, name, dom->id());
265 "Encountered exception when reading entry '{}' from '{}':\n{}",
268 dom->resize(dom->nComponents(), arr->size());
269 if (!header.
hasKey(
"generator")) {
270 arr->meta() = legacyH5(arr, header);
272 arrs[dom->id()] = arr;
276 for (
auto dom :
m_dom) {
278 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
281 "Encountered exception when restoring domain '{}' from HDF:\n{}",
286 }
else if (extension ==
"yaml" || extension ==
"yml") {
288 map<string, shared_ptr<SolutionArray>> arrs;
291 for (
auto dom :
m_dom) {
294 arr->readEntry(root, name, dom->id());
297 "Encountered exception when reading entry '{}' from '{}':\n{}",
300 dom->resize(dom->nComponents(), arr->size());
301 arrs[dom->id()] = arr;
305 for (
auto dom :
m_dom) {
307 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
310 "Encountered exception when restoring domain '{}' from YAML:\n{}",
317 "Unknown file extension '{}'; supported extensions include "
318 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
326 for (
size_t n = 0; n < np; n++) {
333 for (
size_t n = 0; n <
nDomains(); n++) {
336 +
" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
345 throw CanteraError(
"Sim1D::restoreTimeSteppingSolution",
346 "No successful time steps taken on this grid.");
355 "No successful steady state solution");
358 for (
size_t n = 0; n <
nDomains(); n++) {
366 for (
size_t n = 0; n <
nDomains(); n++) {
373 for (
size_t n = 0; n <
nDomains(); n++) {
382 for (
size_t i = 0; i < n; i++) {
393 }
else if (m > -10) {
397 "ERROR: OneDim::solve returned m = {}", m);
409 int attempt_counter = 0;
410 const int max_history = 10;
415 while (new_points > 0) {
421 writeline(
'.', 78,
true,
true);
425 attempt_counter = (attempt_counter % max_history) + 1;
430 debuglog(
"\nAttempt Newton solution of steady-state problem.", loglevel);
435 writelog(
"\nNewton steady-state solve succeeded.\n\n");
437 for (
size_t mm = 1; mm <
nDomains(); mm+=2) {
449 loglevel, attempt_counter);
453 debuglog(
"\nNewton steady-state solve failed.\n", loglevel);
455 loglevel, attempt_counter);
458 writelog(
"\nAttempt {} timesteps.", nsteps);
468 writelog(
"\nFinal timestep info: dt= {:<10.4g} log(ss)= {:<10.4g}\n", dt,
477 dt = std::min(dt,
m_tmax);
481 writeline(
'.', 78,
true,
true);
488 new_points =
refine(loglevel);
494 writeDebugInfo(
"Regridding",
"After regridding", loglevel, attempt_counter);
496 debuglog(
"grid refinement disabled.\n", loglevel);
500 if (new_points < 0) {
504 for (
auto dom :
m_dom) {
515 vector<double> znew, xnew;
516 vector<size_t> dsize;
521 for (
size_t n = 0; n <
nDomains(); n++) {
540 size_t nstart = znew.size();
541 for (
size_t m = 0; m < npnow; m++) {
544 znew.push_back(d.
z(m));
547 for (
size_t i = 0; i < comp; i++) {
548 xnew.push_back(
value(n, i, m));
556 double zmid = 0.5*(d.
z(m) + d.
z(m+1));
557 znew.push_back(zmid);
562 for (
size_t i = 0; i < comp; i++) {
563 double xmid = 0.5*(
value(n, i, m) +
value(n, i, m+1));
564 xnew.push_back(xmid);
570 writelog(
"refine: discarding point at {}\n", d.
z(m));
574 dsize.push_back(znew.size() - nstart);
581 size_t gridstart = 0, gridsize;
582 for (
size_t n = 0; n <
nDomains(); n++) {
588 gridstart += gridsize;
595 return added || -discarded;
600 std::filesystem::remove(
"debug_sim1d.yaml");
604 int loglevel,
int attempt_counter)
608 file_header = fmt::format(
"solution_{}_{}", attempt_counter, header_suffix);
609 save(
"debug_sim1d.yaml", file_header, message,
true);
612 file_header = fmt::format(
"residual_{}_{}", attempt_counter, header_suffix);
613 saveResidual(
"debug_sim1d.yaml", file_header, message,
true);
620 vector<double> znew, xnew;
622 double z1 = 0.0, z2 = 0.0;
623 vector<size_t> dsize;
625 for (
size_t n = 0; n <
nDomains(); n++) {
628 size_t mfixed =
npos;
633 size_t nstart = znew.size();
634 if (d_free && d_free->
isFree()) {
635 for (
size_t m = 0; m < npnow - 1; m++) {
636 bool fixedpt =
false;
640 double thresh = min(1., 1.e-1 * (t2 - t1));
643 if (fabs(t - t1) <= thresh) {
646 }
else if (fabs(t2 - t) <= thresh) {
649 }
else if ((t1 < t) && (t < t2)) {
651 zfixed = (z1 - z2) / (t1 - t2) * (t - t2) + z2;
664 for (
size_t m = 0; m < npnow; m++) {
666 znew.push_back(d.
z(m));
669 for (
size_t i = 0; i < comp; i++) {
670 xnew.push_back(
value(n, i, m));
674 znew.push_back(zfixed);
676 double interp_factor = (zfixed - z2) / (z1 - z2);
679 for (
size_t i = 0; i < comp; i++) {
680 double xmid = interp_factor*(
value(n, i, m) -
value(n, i, m+1)) +
value(n,i,m+1);
681 xnew.push_back(xmid);
685 dsize.push_back(znew.size() - nstart);
691 size_t gridstart = 0;
692 for (
size_t n = 0; n <
nDomains(); n++) {
694 size_t gridsize = dsize[n];
696 gridstart += gridsize;
709 double t_fixed = std::numeric_limits<double>::quiet_NaN();
710 for (
size_t n = 0; n <
nDomains(); n++) {
722 double z_fixed = std::numeric_limits<double>::quiet_NaN();
723 for (
size_t n = 0; n <
nDomains(); n++) {
735 bool two_point_domain_found =
false;
736 for (
size_t n = 0; n <
nDomains(); n++) {
751 two_point_domain_found =
true;
753 double current_val, next_val;
754 for (
size_t m = 0; m < np-1; m++) {
757 if ((current_val - temperature) * (next_val - temperature) < 0.0) {
761 if (std::abs(current_val - temperature) <
762 std::abs(next_val - temperature)) {
774 if (!two_point_domain_found) {
776 "No domain with two-point control enabled was found.");
779 "No control point with temperature {} was able to be found in the"
780 "flame's temperature range.", temperature);
786 bool two_point_domain_found =
false;
787 for (
size_t n = 0; n <
nDomains(); n++) {
802 two_point_domain_found =
true;
804 double current_val, next_val;
805 for (
size_t m = np-1; m > 0; m--) {
808 if ((current_val - temperature) * (next_val - temperature) < 0.0) {
812 if (std::abs(current_val - temperature) <
813 std::abs(next_val - temperature)) {
825 if (!two_point_domain_found) {
827 "No domain with two-point control enabled was found.");
830 "No control point with temperature {} was able to be found in the"
831 "flame's temperature range.", temperature);
837 double slope,
double curve,
double prune)
843 for (
size_t n = 0; n <
nDomains(); n++) {
857 "Must specify domain to get criteria from");
867 for (
size_t n = 0; n <
nDomains(); n++) {
880 for (
size_t n = 0; n <
nDomains(); n++) {
895 warn_deprecated(
"Sim1D::jacobian",
"To be removed after Cantera 3.2.");
906 for (
auto& D :
m_dom) {
907 D->forceFullUpdate(
true);
910 for (
auto& D :
m_dom) {
911 D->forceFullUpdate(
false);
914 auto multijac = dynamic_pointer_cast<MultiJac>(
m_jac);
916 throw CanteraError(
"Sim1D::solveAdjoint",
"Banded (MultiJac) required");
921 for (
size_t i = 0; i <
size(); i++) {
922 size_t j1 = (i > bw) ? i - bw : 0;
923 size_t j2 = (i + bw >=
size()) ?
size() - 1: i + bw;
924 for (
size_t j = j1; j <= j2; j++) {
925 Jt(j,i) = multijac->
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.
vector< double > & grid()
Access the array of grid coordinates [m].
double zmin() const
Get the coordinate [m] of the first (leftmost) grid point in this domain.
size_t nPoints() const
Number of grid points in this domain.
virtual string domainType() const
Domain type flag.
double z(size_t jlocal) const
Get the coordinate [m] of the point with local index jlocal
Refiner & refiner()
Return a reference to the grid refiner.
virtual void show(const double *x)
Print the solution.
virtual string componentName(size_t n) const
Name of component n. May be overloaded.
string type() const
String indicating the domain implemented.
double zmax() const
Get the coordinate [m] of the last (rightmost) grid point in this domain.
virtual void _getInitialSoln(double *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
size_t index(size_t n, size_t j) const
Returns the index of the solution vector, which corresponds to component n at grid point j.
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 setupGrid(size_t n, const double *z)
called to set up initial grid, and after grid refinement
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
void setLeftControlPointTemperature(double temperature)
Sets the temperature of the left control point.
void setLeftControlPointCoordinate(double z_left)
Sets the coordinate of the left control point.
void setRightControlPointCoordinate(double z_right)
Sets the coordinate of the right control point.
bool twoPointControlEnabled() const
Returns the status of the two-point control.
void setRightControlPointTemperature(double temperature)
Sets the temperature of the right control point.
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.
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< SystemJacobian > m_jac
Jacobian evaluator.
shared_ptr< vector< double > > m_state
Solution vector.
vector< int > m_mask
Transient mask. See transientMask().
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
void evalSSJacobian(double *x, double *rsd)
Evaluate the steady-state Jacobian, accessible via jacobian()
double m_tmax
maximum timestep size
int m_ss_jac_age
Maximum age of the Jacobian in steady-state mode.
void setSteadyMode()
Prepare to solve the steady-state problem.
Domain1D & domain(size_t i) const
Return a reference to domain i.
vector< shared_ptr< Domain1D > > m_dom
All domains comprising the system.
MultiNewton & newton()
Return a reference to the Newton iterator.
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
bool newPointNeeded(size_t j)
Returns true if a new grid point is needed to the right of grid index j.
size_t maxPoints() const
Returns the maximum number of points allowed in the domain.
int nNewPoints()
Returns the number of new grid points that were needed.
void show()
Displays the results of the grid refinement analysis.
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.
int analyze(size_t n, const double *z, const double *x)
Determine locations in the grid that need additional grid points and update the internal state of the...
bool keepPoint(size_t j)
Returns true if the grid point at index j should be kept.
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
void getInitialSoln()
Get the initial value of the system state from each domain in the simulation.
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.
void saveResidual(const string &fname, const string &name, const string &desc, bool overwrite=false, int compression=0)
Save the residual of the current solution to a container file.
vector< double > m_xnew
a work array used to hold the residual or the new solution
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.
void writeDebugInfo(const string &header_suffix, const string &message, int loglevel, int attempt_counter)
Write solver debugging information to a YAML file based on the specified log level.
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 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.
void clearDebugFile()
Deletes a debug_sim1d.yaml file if it exists.
void setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component for all domains.
void solve(int loglevel=0, bool refine_grid=true)
Performs the hybrid Newton steady/time-stepping solution.
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 evalSSJacobian()
Evaluate the Jacobian in steady-state mode.
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
AnyMap restore(const string &fname, const string &name)
Retrieve data and settings from a previously saved simulation.
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.
void setTimeStep(double stepsize, size_t n, const int *tsteps)
Set the number of time steps to try when the steady Newton solver is unsuccessful.
void setRightControlPoint(double temperature)
Set the right control point location using the specified temperature.
double value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
double workValue(size_t dom, size_t comp, size_t localPoint) const
Get an entry in the work vector, which may contain either a new system state or the current residual ...
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
void save(const string &fname, const string &name, const string &desc, bool overwrite=false, int compression=0, const string &basis="")
Save current simulation data to a container file or CSV format.
void setLeftControlPoint(double temperature)
Set the left control point location using the specified temperature.
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.
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"
@ c_offset_T
temperature [kelvin]
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.