20 : m_tmin(1.0e-16), m_tmax(1e8), m_tfactor(0.5),
21 m_rdt(0.0), m_jac_ok(false),
23 m_init(false), m_pts(0), m_solve_time(0.0),
24 m_ss_jac_age(20), m_ts_jac_age(20),
25 m_interrupt(0), m_time_step_callback(0),
26 m_nsteps(0), m_nsteps_max(500),
27 m_nevals(0), m_evaltime(0.0)
29 m_newt.reset(
new MultiNewton(1));
32 OneDim::OneDim(vector<Domain1D*> domains) :
33 m_tmin(1.0e-16), m_tmax(1e8), m_tfactor(0.5),
34 m_rdt(0.0), m_jac_ok(false),
36 m_init(false), m_solve_time(0.0),
37 m_ss_jac_age(20), m_ts_jac_age(20),
38 m_interrupt(0), m_time_step_callback(0),
39 m_nsteps(0), m_nsteps_max(500),
40 m_nevals(0), m_evaltime(0.0)
44 for (
size_t i = 0; i < domains.size(); i++) {
55 size_t OneDim::domainIndex(
const std::string& name)
57 for (
size_t n = 0; n < m_dom.size(); n++) {
58 if (
domain(n).id() == name) {
62 throw CanteraError(
"OneDim::domainIndex",
"no domain named >>"+name+
"<<");
73 size_t offset = i -
start(n);
83 size_t n = m_dom.size();
85 m_dom.back()->append(d);
90 m_connect.push_back(d);
110 void OneDim::setJacAge(
int ss_age,
int ts_age)
112 m_ss_jac_age = ss_age;
114 m_ts_jac_age = ts_age;
116 m_ts_jac_age = m_ss_jac_age;
123 writelog(
"\nStatistics:\n\n Grid Timesteps Functions Time Jacobians Time\n");
124 size_t n = m_gridpts.size();
125 for (
size_t i = 0; i < n; i++) {
127 writelog(
"{:5d} {:5d} {:6d} {:9.4f} {:5d} {:9.4f}\n",
128 m_gridpts[i],
m_timeSteps[i], m_funcEvals[i], m_funcElapsed[i],
129 m_jacEvals[i], m_jacElapsed[i]);
131 writelog(
"{:5d} {:5d} {:6d} NA {:5d} NA\n",
132 m_gridpts[i],
m_timeSteps[i], m_funcEvals[i], m_jacEvals[i]);
140 int nev =
m_jac->nEvals();
141 if (nev > 0 && m_nevals > 0) {
142 m_gridpts.push_back(m_pts);
143 m_jacEvals.push_back(
m_jac->nEvals());
144 m_jacElapsed.push_back(
m_jac->elapsedTime());
145 m_funcEvals.push_back(m_nevals);
147 m_funcElapsed.push_back(m_evaltime);
159 m_jacElapsed.clear();
161 m_funcElapsed.clear();
178 for (
size_t i = 0; i <
nDomains(); i++) {
183 for (
size_t n = 0; n < np; n++) {
184 m_nvars.push_back(nv);
202 size_t bw2 = m_dom[i-1]->bandwidth();
204 bw2 = m_dom[i-1]->nComponents();
213 m_mask.resize(
size());
219 for (
size_t i = 0; i <
nDomains(); i++) {
220 m_dom[i]->setJac(
m_jac.get());
228 m_jac->eval(x, xnew, 0.0);
233 return m_newt->solve(x, xnew, *
this, *
m_jac, loglevel);
236 void OneDim::evalSSJacobian(doublereal* x, doublereal* xnew)
238 doublereal rdt_save =
m_rdt;
242 m_jac->eval(x, xnew, 0.0);
258 void OneDim::eval(
size_t j,
double* x,
double* r, doublereal rdt,
int count)
260 clock_t t0 = clock();
266 fill(m_mask.begin(), m_mask.end(), 0);
273 for (
const auto& d : m_bulk) {
274 d->eval(j, x, r, m_mask.data(),
rdt);
278 for (
const auto& d : m_connect) {
279 d->eval(j, x, r, m_mask.data(),
rdt);
284 clock_t t1 = clock();
285 m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
294 for (
size_t i = 0; i <
m_size; i++) {
295 ss = std::max(fabs(r[i]),ss);
302 doublereal rdt_old =
m_rdt;
349 doublereal* r,
int loglevel)
354 debuglog(
"\n\n step size (s) log10(ss) \n", loglevel);
355 debuglog(
"===============================\n", loglevel);
358 int successiveFailures = 0;
362 doublereal ss =
ssnorm(x, r);
363 writelog(
" {:>4d} {:10.4g} {:10.4g}", n, dt, log10(ss));
370 int m =
solve(x, r, loglevel-1);
375 successiveFailures = 0;
386 dt = std::min(dt,
m_tmax);
389 "Took maximum number of timesteps allowed ({}) without " 393 successiveFailures++;
396 debuglog(
"...failure.\n", loglevel);
397 if (successiveFailures > 2) {
400 successiveFailures = 0;
405 "Time integration failed.");
416 void OneDim::resetBadValues(
double* x)
418 for (
auto dom : m_dom) {
419 dom->resetBadValues(x);
423 void OneDim::save(
const std::string& fname, std::string
id,
424 const std::string& desc, doublereal* sol,
429 struct tm* newtime = localtime(&aclock);
431 XML_Node root(
"ctml");
434 root.build(fin, fname);
436 XML_Node* same_ID = root.findID(
id);
438 same_ID->parent()->removeChild(same_ID);
442 XML_Node& sim = root.addChild(
"simulation");
443 sim.addAttribute(
"id",
id);
444 addString(sim,
"timestamp",asctime(newtime));
449 Domain1D* d =
left();
456 throw CanteraError(
"OneDim::save",
"could not open file "+fname);
460 debuglog(
"Solution saved to file "+fname+
" as solution "+
id+
".\n", loglevel);
Domain1D * left()
Pointer to left-most domain (first added).
void addDomain(Domain1D *d)
Add a domain. Domains are added left-to-right.
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
vector_int m_timeSteps
Number of time steps taken in each call to solve() (e.g.
virtual void resize()
Call after one or more grids has changed size, e.g. after being refined.
void setSteadyMode()
Prepare to solve the steady-state problem.
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.
size_t size() const
Total solution vector length;.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator.
size_t m_bw
Jacobian bandwidth.
std::unique_ptr< MultiJac > m_jac
Jacobian evaluator.
const size_t npos
index returned by functions to indicate "no position"
Domain1D * right() const
Return a pointer to the right neighbor.
void addString(XML_Node &node, const std::string &titleString, const std::string &valueString, const std::string &typeString)
This function adds a child node with the name string with a string value to the current node...
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
int solve(doublereal *x0, doublereal *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
void setContainer(OneDim *c, size_t index)
Specify the container object for this domain, and the position of this domain in the list...
virtual doublereal eval(doublereal t) const
Evaluate the function.
MultiNewton & newton()
Return a reference to the Newton iterator.
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.
size_t bandwidth()
Set the Jacobian bandwidth for this domain.
std::unique_ptr< MultiNewton > m_newt
Newton iterator.
size_t m_size
solution vector size
Domain1D * left() const
Return a pointer to the left neighbor.
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...
doublereal m_tfactor
factor time step is multiplied by if time stepping fails ( < 1 )
Base class for exceptions thrown by Cantera classes.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
void clearStats()
Clear saved statistics.
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
void writeStats(int printTime=1)
Write statistics about the number of iterations and Jacobians at each grid level. ...
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
doublereal m_rdt
reciprocal of time step
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
size_t nComponents() const
Number of components at each grid point.
void initTimeInteg(doublereal dt, doublereal *x)
Prepare for time stepping beginning with solution x and timestep dt.
const doublereal Tiny
Small number to compare differences of mole fractions against.
Func1 * m_interrupt
Function called at the start of every call to eval.
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
doublereal rdt() const
Reciprocal of the time step.
Newton iterator for multi-domain, one-dimensional problems.
Domain1D * right()
Pointer to right-most domain (last added).
size_t nPoints() const
Number of grid points in this domain.
Namespace for the Cantera kernel.
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
int m_nsteps
Number of time steps taken in the current call to solve()
doublereal m_tmax
maximum timestep size
void setOptions(int maxJacAge=5)
Set options.
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
bool m_jac_ok
if true, Jacobian is current
void initTimeInteg(doublereal dt, const doublereal *x0)
Prepare to do time stepping with time step dt.
doublereal m_tmin
minimum timestep size