21 Sim1D::Sim1D(vector<Domain1D*>& domains) :
28 for (
size_t n = 0; n <
nDomains(); n++) {
39 for (
size_t dom=0; dom<
nDomains(); dom++) {
42 for (
size_t comp=0; comp<ncomp; comp++) {
50 void Sim1D::setValue(
size_t dom,
size_t comp,
size_t localPoint, doublereal value)
54 "Index out of bounds: {} > {}", iloc,
m_x.size());
58 doublereal
Sim1D::value(
size_t dom,
size_t comp,
size_t localPoint)
const 62 "Index out of bounds: {} > {}", iloc,
m_x.size());
66 doublereal Sim1D::workValue(
size_t dom,
size_t comp,
size_t localPoint)
const 70 "Index out of bounds: {} > {}", iloc,
m_x.size());
77 if (pos.front() != 0.0 || pos.back() != 1.0) {
79 "`pos` vector must span the range [0, 1]. Got a vector spanning " 80 "[{}, {}] instead.", pos.front(), pos.back());
83 doublereal z0 = d.zmin();
84 doublereal z1 = d.zmax();
85 for (
size_t n = 0; n < d.
nPoints(); n++) {
87 double frac = (zpt - z0)/(z1 - z0);
93 void Sim1D::save(
const std::string& fname,
const std::string&
id,
94 const std::string& desc,
int loglevel)
96 OneDim::save(fname,
id, desc,
m_x.data(), loglevel);
99 void Sim1D::saveResidual(
const std::string& fname,
const std::string&
id,
100 const std::string& desc,
int loglevel)
104 OneDim::save(fname,
id, desc, &res[0], loglevel);
115 throw CanteraError(
"Sim1D::restore",
"No solution with id = "+
id);
120 throw CanteraError(
"Sim1D::restore",
"Solution does not contain the " 121 " correct number of domains. Found {} expected {}.\n",
124 for (
size_t m = 0; m <
nDomains(); m++) {
126 if (loglevel > 0 && xd[m]->attrib(
"id") != dom.id()) {
127 writelog(
"Warning: domain names do not match: '" +
128 (*xd[m])[
"id"] + +
"' and '" + dom.id() +
"'\n");
134 for (
size_t m = 0; m <
nDomains(); m++) {
143 for (
size_t n = 0; n < np; n++) {
148 void Sim1D::showSolution(ostream& s)
150 for (
size_t n = 0; n <
nDomains(); n++) {
157 void Sim1D::showSolution()
159 for (
size_t n = 0; n <
nDomains(); n++) {
162 +
" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
171 throw CanteraError(
"Sim1D::restoreTimeSteppingSolution",
172 "No successful time steps taken on this grid.");
181 "No successful steady state solution");
184 for (
size_t n = 0; n <
nDomains(); n++) {
190 void Sim1D::getInitialSoln()
192 for (
size_t n = 0; n <
nDomains(); n++) {
199 for (
size_t n = 0; n <
nDomains(); n++) {
204 void Sim1D::setTimeStep(
double stepsize,
size_t n,
const int* tsteps)
208 for (
size_t i = 0; i < n; i++) {
219 }
else if (m > -10) {
223 "ERROR: OneDim::solve returned m = {}", m);
227 void Sim1D::solve(
int loglevel,
bool refine_grid)
232 int soln_number = -1;
235 while (new_points > 0) {
241 writeline(
'.', 78,
true,
true);
247 debuglog(
"Attempt Newton solution of steady-state problem...", loglevel);
254 for (
size_t mm = 1; mm <
nDomains(); mm+=2) {
267 save(
"debug_sim1d.xml",
"debug",
268 "After successful Newton solve");
271 saveResidual(
"debug_sim1d.xml",
"residual",
272 "After successful Newton solve");
279 save(
"debug_sim1d.xml",
"debug",
280 "After unsuccessful Newton solve");
283 saveResidual(
"debug_sim1d.xml",
"residual",
284 "After unsuccessful Newton solve");
287 writelog(
"Take {} timesteps ", nsteps);
293 save(
"debug_sim1d.xml",
"debug",
"After timestepping");
296 saveResidual(
"debug_sim1d.xml",
"residual",
297 "After timestepping");
301 writelog(
" {:10.4g} {:10.4g}\n", dt,
310 dt = std::min(dt,
m_tmax);
314 writeline(
'.', 78,
true,
true);
321 new_points =
refine(loglevel);
327 if (new_points && loglevel > 6) {
328 save(
"debug_sim1d.xml",
"debug",
"After regridding");
330 if (new_points && loglevel > 7) {
331 saveResidual(
"debug_sim1d.xml",
"residual",
335 debuglog(
"grid refinement disabled.\n", loglevel);
343 int ianalyze, np = 0;
345 std::vector<size_t> dsize;
350 for (
size_t n = 0; n <
nDomains(); n++) {
358 ianalyze = r.analyze(d.grid().size(), d.grid().data(), &
m_x[
start(n)]);
367 np += r.nNewPoints();
372 size_t nstart = znew.size();
373 for (
size_t m = 0; m < npnow; m++) {
374 if (r.keepPoint(m)) {
376 znew.push_back(d.grid(m));
379 for (
size_t i = 0; i < comp; i++) {
380 xnew.push_back(
value(n, i, m));
386 if (r.newPointNeeded(m) && m + 1 < npnow) {
388 double zmid = 0.5*(d.grid(m) + d.grid(m+1));
389 znew.push_back(zmid);
394 for (
size_t i = 0; i < comp; i++) {
395 double xmid = 0.5*(
value(n, i, m) +
value(n, i, m+1));
396 xnew.push_back(xmid);
401 writelog(
"refine: discarding point at {}\n", d.grid(m));
405 dsize.push_back(znew.size() - nstart);
412 size_t gridstart = 0, gridsize;
413 for (
size_t n = 0; n <
nDomains(); n++) {
417 gridstart += gridsize;
432 doublereal z1 = 0.0, z2 = 0.0, t1,t2;
434 std::vector<size_t> dsize;
436 for (
size_t n = 0; n <
nDomains(); n++) {
445 size_t nstart = znew.size();
446 if (d_free && d_free->
domainType() == cFreeFlow) {
447 for (
size_t m = 0; m < npnow-1; m++) {
448 if (
value(n,2,m) == t) {
454 }
else if ((
value(n,2,m)<t) && (
value(n,2,m+1)>t)) {
461 zfixed = (z1-z2)/(t1-t2)*(t-t2)+z2;
471 for (
size_t m = 0; m < npnow; m++) {
473 znew.push_back(d.grid(m));
476 for (
size_t i = 0; i < comp; i++) {
477 xnew.push_back(
value(n, i, m));
479 if (m==m1 && addnewpt) {
481 znew.push_back(zfixed);
483 double interp_factor = (zfixed-z2) / (z1-z2);
486 for (
size_t i = 0; i < comp; i++) {
487 double xmid = interp_factor*(
value(n, i, m) -
value(n, i, m+1)) +
value(n,i,m+1);
488 xnew.push_back(xmid);
492 dsize.push_back(znew.size() - nstart);
498 size_t gridstart = 0;
499 for (
size_t n = 0; n <
nDomains(); n++) {
501 size_t gridsize = dsize[n];
503 gridstart += gridsize;
515 doublereal slope, doublereal curve, doublereal prune)
521 for (
size_t n = 0; n <
nDomains(); n++) {
535 "Must specify domain to get criteria from");
545 for (
size_t n = 0; n <
nDomains(); n++) {
558 for (
size_t n = 0; n <
nDomains(); n++) {
576 void Sim1D::evalSSJacobian()
578 OneDim::evalSSJacobian(
m_x.data(),
m_xnew.data());
583 for (
auto& D : m_dom) {
584 D->forceFullUpdate(
true);
587 for (
auto& D : m_dom) {
588 D->forceFullUpdate(
false);
594 for (
size_t i = 0; i <
size(); i++) {
595 size_t j1 = (i > bw) ? i - bw : 0;
596 size_t j2 = (i + bw >=
size()) ?
size() - 1: i + bw;
597 for (
size_t j = j1; j <= j2; j++) {
598 Jt(j,i) =
m_jac->value(i,j);
Container class for multiple-domain 1D problems.
std::vector< XML_Node * > getChildren(const std::string &name) const
Get a vector of pointers to XML_Node containing all of the children of the current node which match t...
std::vector< vector_fp > m_grid_last_ss
the grids for each domain after the last successful steady-state solve (stored before grid refinement...
void showSolution(std::ostream &s)
Print to stream s the current solution for all domains.
virtual void _getInitialSoln(doublereal *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
void restore(const std::string &fname, const std::string &id, int loglevel=2)
Initialize the solution with a previously-saved solution.
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
vector_fp m_x
the solution vector
vector_fp m_xlast_ts
the solution vector after the last successful timestepping
virtual void resize()
Call after one or more grids has changed size, e.g. after being refined.
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
doublereal ssnorm(doublereal *x, doublereal *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x. ...
void eval(size_t j, double *x, double *r, doublereal rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
vector_fp m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement) ...
size_t size() const
Total solution vector length;.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator.
std::unique_ptr< MultiJac > m_jac
Jacobian evaluator.
const size_t npos
index returned by functions to indicate "no position"
void setValue(size_t dom, size_t comp, size_t localPoint, doublereal value)
Set a single value in the solution vector.
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
void setRefineCriteria(int dom=-1, doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
int solve(doublereal *x0, doublereal *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
Class XML_Node is a tree-based representation of the contents of an XML file.
void setInitialGuess(const std::string &component, vector_fp &locs, vector_fp &vals)
Set initial guess for one component for all domains.
virtual doublereal eval(doublereal t) const
Evaluate the function.
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
size_t maxPoints() const
Returns the maximum number of points allowed in the domain.
size_t loc(size_t jg)
Location in the solution vector of the first component of global point jg.
vector_fp getRefineCriteria(int dom)
Get the grid refinement criteria.
MultiNewton & newton()
Return a reference to the Newton iterator.
vector_int m_steps
array of number of steps to take before re-attempting the steady-state solution
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
size_t nDomains() const
Number of domains.
Domain1D & domain(size_t i) const
Return a reference to domain i.
Base class for one-dimensional domains.
virtual void resize(size_t nv, size_t np)
int setFixedTemperature(doublereal t)
Add node for fixed temperature point of freely propagating flame.
void setCriteria(doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
void setMaxPoints(int npmax)
Set the maximum number of points allowed in the domain.
void finalize()
Calls method _finalize in each domain.
virtual void _finalize(const doublereal *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
size_t bandwidth() const
Jacobian bandwidth.
Classes providing support for XML data files.
virtual void resize()
Call after one or more grids has changed size, e.g. after being refined.
virtual void showSolution(const doublereal *x)
Print the solution.
std::tuple< std::string, size_t, std::string > component(size_t i)
Return the domain, local point index, and component name for the i-th component of the global solutio...
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
Base class for exceptions thrown by Cantera classes.
doublereal m_tstep
timestep
void build(const std::string &filename)
Populate the XML tree from an input file.
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
void setFlatProfile(size_t dom, size_t comp, doublereal v)
Set component 'comp' of domain 'dom' to value 'v' at all points.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Refiner & refiner()
Return a reference to the grid refiner.
int domainType()
Domain type flag.
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
int intValue(const std::string &val)
Translate a string into one integer value.
vector_fp getCriteria()
Get the grid refinement criteria.
vector_fp m_xnew
a work array used to hold the residual or the new solution
size_t nComponents() const
Number of components at each grid point.
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
double m_tfixed
Temperature at the point used to fix the flame location.
doublereal value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
int refine(int loglevel=0)
Refine the grid in all domains.
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
void setProfile(size_t dom, size_t comp, const vector_fp &pos, const vector_fp &values)
Specify a profile for one component of one domain.
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
size_t nPoints() const
Number of grid points in this domain.
Namespace for the Cantera kernel.
double m_zfixed
Location of the point where temperature is fixed.
Header for a file containing miscellaneous numerical functions.
int m_nsteps
Number of time steps taken in the current call to solve()
void setGridMin(int dom, double gridmin)
Set the minimum grid spacing in the specified domain(s).
doublereal m_tmax
maximum timestep size
void setOptions(int maxJacAge=5)
Set options.
virtual void setupGrid(size_t n, const doublereal *z)
called to set up initial grid, and after grid refinement
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
XML_Node * findID(const std::string &id, const int depth=100) const
This routine carries out a recursive search for an XML node based on the XML element attribute...
doublereal linearInterp(doublereal x, const vector_fp &xpts, const vector_fp &fpts)
Linearly interpolate a function defined on a discrete grid.
A class for banded matrices, involving matrix inversion processes.