12 #include "cantera/oneD/refine.h"
25 Sim1D::Sim1D(vector<shared_ptr<Domain1D>>& domains) :
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;
63 double Sim1D::value(
size_t dom,
size_t comp,
size_t localPoint)
const
67 "Index out of bounds: {} > {}", iloc,
m_state->size());
71 double Sim1D::workValue(
size_t dom,
size_t comp,
size_t localPoint)
const
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);
98 void 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);
163 AnyMap 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"},
206 {
"species-enabled",
"species_enabled"},
208 for (
const auto& [newName, oldName] : header_pairs) {
209 if (header.hasKey(oldName)) {
210 out[newName] = header[oldName];
214 map<string, string> refiner_pairs = {
220 {
"max-points",
"max_grid_points"},
222 for (
const auto& [newName, oldName] : refiner_pairs) {
223 if (header.hasKey(oldName)) {
224 out[
"refine-criteria"][newName] = header[oldName];
228 if (header.hasKey(
"fixed_temperature")) {
229 double temp = header.
getDouble(
"fixed_temperature", -1.);
230 auto profile = arr->getComponent(
"T").as<vector<double>>();
232 while (profile[ix] <= temp && ix < arr->size()) {
236 auto grid = arr->getComponent(
"grid").as<vector<double>>();
237 out[
"fixed-point"][
"location"] = grid[ix - 1];
238 out[
"fixed-point"][
"temperature"] = temp;
249 size_t dot = fname.find_last_of(
".");
251 if (extension ==
"xml") {
253 "Restoring from XML is no longer supported.");
256 if (extension ==
"h5" || extension ==
"hdf" || extension ==
"hdf5") {
257 map<string, shared_ptr<SolutionArray>> arrs;
260 for (
auto dom : m_dom) {
263 arr->readEntry(fname, name, dom->id());
266 "Encountered exception when reading entry '{}' from '{}':\n{}",
269 dom->resize(dom->nComponents(), arr->size());
270 if (!header.
hasKey(
"generator")) {
271 arr->meta() = legacyH5(arr, header);
273 arrs[dom->id()] = arr;
277 for (
auto dom : m_dom) {
279 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
282 "Encountered exception when restoring domain '{}' from HDF:\n{}",
287 }
else if (extension ==
"yaml" || extension ==
"yml") {
289 map<string, shared_ptr<SolutionArray>> arrs;
292 for (
auto dom : m_dom) {
295 arr->readEntry(root, name, dom->id());
298 "Encountered exception when reading entry '{}' from '{}':\n{}",
301 dom->resize(dom->nComponents(), arr->size());
302 arrs[dom->id()] = arr;
306 for (
auto dom : m_dom) {
308 dom->fromArray(*arrs[dom->id()],
m_state->data() + dom->loc());
311 "Encountered exception when restoring domain '{}' from YAML:\n{}",
318 "Unknown file extension '{}'; supported extensions include "
319 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
327 for (
size_t n = 0; n < np; n++) {
334 for (
size_t n = 0; n <
nDomains(); n++) {
343 for (
size_t n = 0; n <
nDomains(); n++) {
346 +
" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
355 throw CanteraError(
"Sim1D::restoreTimeSteppingSolution",
356 "No successful time steps taken on this grid.");
365 "No successful steady state solution");
368 for (
size_t n = 0; n <
nDomains(); n++) {
374 void Sim1D::getInitialSoln()
376 for (
size_t n = 0; n <
nDomains(); n++) {
383 for (
size_t n = 0; n <
nDomains(); n++) {
388 void Sim1D::setTimeStep(
double stepsize,
size_t n,
const int* tsteps)
392 for (
size_t i = 0; i < n; i++) {
403 }
else if (m > -10) {
407 "ERROR: OneDim::solve returned m = {}", m);
411 void Sim1D::solve(
int loglevel,
bool refine_grid)
418 while (new_points > 0) {
424 writeline(
'.', 78,
true,
true);
430 debuglog(
"Attempt Newton solution of steady-state problem...", loglevel);
437 for (
size_t mm = 1; mm <
nDomains(); mm+=2) {
450 save(
"debug_sim1d.yaml",
"debug",
451 "After successful Newton solve");
455 "After successful Newton solve");
461 save(
"debug_sim1d.yaml",
"debug",
462 "After unsuccessful Newton solve");
466 "After unsuccessful Newton solve");
469 writelog(
"Take {} timesteps ", nsteps);
474 save(
"debug_sim1d.yaml",
"debug",
"After timestepping");
478 "After timestepping");
482 writelog(
" {:10.4g} {:10.4g}\n", dt,
491 dt = std::min(dt,
m_tmax);
495 writeline(
'.', 78,
true,
true);
502 new_points =
refine(loglevel);
508 if (new_points && loglevel > 6) {
509 save(
"debug_sim1d.yaml",
"debug",
"After regridding");
511 if (new_points && loglevel > 7) {
516 debuglog(
"grid refinement disabled.\n", loglevel);
524 int ianalyze, np = 0;
525 vector<double> znew, xnew;
526 vector<size_t> dsize;
531 for (
size_t n = 0; n <
nDomains(); n++) {
539 ianalyze = r.analyze(d.grid().size(), d.grid().data(),
549 np += r.nNewPoints();
554 size_t nstart = znew.size();
555 for (
size_t m = 0; m < npnow; m++) {
556 if (r.keepPoint(m)) {
558 znew.push_back(d.grid(m));
561 for (
size_t i = 0; i < comp; i++) {
562 xnew.push_back(
value(n, i, m));
568 if (r.newPointNeeded(m) && m + 1 < npnow) {
570 double zmid = 0.5*(d.grid(m) + d.grid(m+1));
571 znew.push_back(zmid);
576 for (
size_t i = 0; i < comp; i++) {
577 double xmid = 0.5*(
value(n, i, m) +
value(n, i, m+1));
578 xnew.push_back(xmid);
583 writelog(
"refine: discarding point at {}\n", d.grid(m));
587 dsize.push_back(znew.size() - nstart);
594 size_t gridstart = 0, gridsize;
595 for (
size_t n = 0; n <
nDomains(); n++) {
601 gridstart += gridsize;
614 vector<double> znew, xnew;
616 double z1 = 0.0, z2 = 0.0;
617 vector<size_t> dsize;
619 for (
size_t n = 0; n <
nDomains(); n++) {
622 size_t mfixed =
npos;
627 size_t nstart = znew.size();
628 if (d_free && d_free->
isFree()) {
629 for (
size_t m = 0; m < npnow - 1; m++) {
630 bool fixedpt =
false;
634 double thresh = min(1., 1.e-1 * (t2 - t1));
637 if (fabs(t - t1) <= thresh) {
640 }
else if (fabs(t2 - t) <= thresh) {
643 }
else if ((t1 < t) && (t < t2)) {
645 zfixed = (z1 - z2) / (t1 - t2) * (t - t2) + z2;
658 for (
size_t m = 0; m < npnow; m++) {
660 znew.push_back(d.grid(m));
663 for (
size_t i = 0; i < comp; i++) {
664 xnew.push_back(
value(n, i, m));
668 znew.push_back(zfixed);
670 double interp_factor = (zfixed - z2) / (z1 - z2);
673 for (
size_t i = 0; i < comp; i++) {
674 double xmid = interp_factor*(
value(n, i, m) -
value(n, i, m+1)) +
value(n,i,m+1);
675 xnew.push_back(xmid);
679 dsize.push_back(znew.size() - nstart);
685 size_t gridstart = 0;
686 for (
size_t n = 0; n <
nDomains(); n++) {
688 size_t gridsize = dsize[n];
690 gridstart += gridsize;
703 double t_fixed = std::numeric_limits<double>::quiet_NaN();
704 for (
size_t n = 0; n <
nDomains(); n++) {
716 double z_fixed = std::numeric_limits<double>::quiet_NaN();
717 for (
size_t n = 0; n <
nDomains(); n++) {
728 double slope,
double curve,
double prune)
734 for (
size_t n = 0; n <
nDomains(); n++) {
748 "Must specify domain to get criteria from");
758 for (
size_t n = 0; n <
nDomains(); n++) {
771 for (
size_t n = 0; n <
nDomains(); n++) {
789 void Sim1D::evalSSJacobian()
796 for (
auto& D : m_dom) {
797 D->forceFullUpdate(
true);
800 for (
auto& D : m_dom) {
801 D->forceFullUpdate(
false);
807 for (
size_t i = 0; i <
size(); i++) {
808 size_t j1 = (i > bw) ? i - bw : 0;
809 size_t j2 = (i + bw >=
size()) ?
size() - 1: i + bw;
810 for (
size_t j = j1; j <= j2; j++) {
811 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.
size_t nPoints() const
Number of grid points in this domain.
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Refiner & refiner()
Return a reference to the grid refiner.
string type() const
String indicating the domain implemented.
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.
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.
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.
vector< double > getCriteria()
Get the 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.
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.
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 setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component 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 .
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.
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.
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.
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"
Contains declarations for string manipulation functions within Cantera.