21 m_newt = make_unique<MultiNewton>(1);
24OneDim::OneDim(vector<shared_ptr<Domain1D>>& domains)
27 m_newt = make_unique<MultiNewton>(1);
28 m_state = make_shared<vector<double>>();
29 for (
auto& dom : domains) {
36OneDim::OneDim(vector<Domain1D*> domains)
39 "To be removed after Cantera 3.0; superseded by "
40 "OneDim::OneDim(vector<shared_ptr<Domain1D>>).");
43 m_newt = make_unique<MultiNewton>(1);
44 m_state = make_shared<vector<double>>();
45 for (
size_t i = 0; i < domains.size(); i++) {
56size_t OneDim::domainIndex(
const string& name)
58 for (
size_t n = 0; n <
m_dom.size(); n++) {
59 if (
domain(n).
id() == name) {
63 throw CanteraError(
"OneDim::domainIndex",
"no domain named >>"+name+
"<<");
84 size_t n =
m_dom.size();
86 m_dom.back()->append(d.get());
91 m_sharedConnect.push_back(d);
94 m_sharedBulk.push_back(d);
99 m_sharedDom.push_back(d);
100 m_dom.push_back(d.get());
102 d->setContainer(
this,
m_dom.size()-1);
109 "To be removed after Cantera 3.0; superseded by "
110 "OneDim::addDomain(shared_ptr<Domain1D>).");
114 size_t n =
m_dom.size();
116 m_dom.back()->append(d);
142void OneDim::setJacAge(
int ss_age,
int ts_age)
144 m_ss_jac_age = ss_age;
146 m_ts_jac_age = ts_age;
148 m_ts_jac_age = m_ss_jac_age;
155 writelog(
"\nStatistics:\n\n Grid Timesteps Functions Time Jacobians Time\n");
156 size_t n = m_gridpts.size();
157 for (
size_t i = 0; i < n; i++) {
159 writelog(
"{:5d} {:5d} {:6d} {:9.4f} {:5d} {:9.4f}\n",
160 m_gridpts[i],
m_timeSteps[i], m_funcEvals[i], m_funcElapsed[i],
161 m_jacEvals[i], m_jacElapsed[i]);
163 writelog(
"{:5d} {:5d} {:6d} NA {:5d} NA\n",
164 m_gridpts[i],
m_timeSteps[i], m_funcEvals[i], m_jacEvals[i]);
172 int nev =
m_jac->nEvals();
173 if (nev > 0 && m_nevals > 0) {
174 m_gridpts.push_back(m_pts);
175 m_jacEvals.push_back(
m_jac->nEvals());
176 m_jacElapsed.push_back(
m_jac->elapsedTime());
177 m_funcEvals.push_back(m_nevals);
179 m_funcElapsed.push_back(m_evaltime);
191 m_jacElapsed.clear();
193 m_funcElapsed.clear();
210 for (
size_t i = 0; i <
nDomains(); i++) {
215 for (
size_t n = 0; n < np; n++) {
216 m_nvars.push_back(nv);
234 size_t bw2 =
m_dom[i-1]->bandwidth();
236 bw2 =
m_dom[i-1]->nComponents();
247 m_mask.resize(
size());
250 m_jac = make_unique<MultiJac>(*
this);
253 for (
size_t i = 0; i <
nDomains(); i++) {
262 m_jac->eval(x, xnew, 0.0);
267 return m_newt->solve(x, xnew, *
this, *
m_jac, loglevel);
270void OneDim::evalSSJacobian(
double* x,
double* xnew)
272 double rdt_save =
m_rdt;
276 m_jac->eval(x, xnew, 0.0);
292void OneDim::eval(
size_t j,
double* x,
double* r,
double rdt,
int count)
294 clock_t t0 = clock();
300 fill(m_mask.begin(), m_mask.end(), 0);
307 for (
const auto& d :
m_bulk) {
308 d->eval(j, x, r, m_mask.data(),
rdt);
313 d->eval(j, x, r, m_mask.data(),
rdt);
318 clock_t t1 = clock();
319 m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
328 for (
size_t i = 0; i <
m_size; i++) {
329 ss = std::max(fabs(r[i]),ss);
336 double rdt_old =
m_rdt;
387 debuglog(
"\n\n step size (s) log10(ss) \n", loglevel);
388 debuglog(
"===============================\n", loglevel);
391 int successiveFailures = 0;
396 writelog(
" {:>4d} {:10.4g} {:10.4g}", n, dt, log10(ss));
403 int m =
solve(x, r, loglevel-1);
408 successiveFailures = 0;
419 dt = std::min(dt,
m_tmax);
422 "Took maximum number of timesteps allowed ({}) without "
426 successiveFailures++;
429 debuglog(
"...failure.\n", loglevel);
430 if (successiveFailures > 2) {
433 successiveFailures = 0;
438 "Time integration failed.");
449void OneDim::resetBadValues(
double* x)
451 for (
auto dom :
m_dom) {
452 dom->resetBadValues(x);
459 "To be removed after Cantera 3.0; unused.");
461 for (
size_t i = 0; i <
m_dom.size(); i++) {
A map of string keys to values whose type can vary at runtime.
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
size_t nComponents() const
Number of components at each grid point.
size_t bandwidth()
Set the Jacobian bandwidth for this domain.
Domain1D * left() const
Return a pointer to the left neighbor.
void setContainer(OneDim *c, size_t index)
Specify the container object for this domain, and the position of this domain in the list.
size_t nPoints() const
Number of grid points in this domain.
Domain1D * right() const
Return a pointer to the right neighbor.
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
virtual void init()
Initialize.
void initTimeInteg(double dt, const double *x0)
Prepare to do time stepping with time step dt.
void setData(shared_ptr< vector< double > > &data)
Set shared data pointer.
void setSteadyMode()
Prepare to solve the steady-state problem.
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
virtual double eval(double t) const
Evaluate the function.
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Newton iterator for multi-domain, one-dimensional problems.
void setOptions(int maxJacAge=5)
Set options.
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()
void init()
Initialize all domains.
vector< Domain1D * > m_bulk
size_t m_size
solution vector size
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
AnyMap serialize(const double *soln) const
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
unique_ptr< MultiNewton > m_newt
Newton iterator.
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.
void addDomain(shared_ptr< Domain1D > d)
Add a domain. Domains are added left-to-right.
double rdt() const
Reciprocal of the time step.
void initTimeInteg(double dt, double *x)
Prepare for time stepping beginning with solution x and timestep dt.
size_t nDomains() const
Number of domains.
Domain1D * right()
Pointer to right-most domain (last added).
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...
double m_rdt
reciprocal of time step
shared_ptr< vector< double > > m_state
Solution vector.
Func1 * m_interrupt
Function called at the start of every call to eval.
vector< Domain1D * > m_connect
unique_ptr< MultiJac > m_jac
Jacobian evaluator.
size_t m_bw
Jacobian bandwidth.
double m_tfactor
factor time step is multiplied by if time stepping fails ( < 1 )
vector< Domain1D * > m_dom
bool m_jac_ok
if true, Jacobian is current
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
double m_tmin
minimum timestep size
double m_tmax
maximum timestep size
void writeStats(int printTime=1)
Write statistics about the number of iterations and Jacobians at each grid level.
void clearStats()
Clear saved statistics.
void setSteadyMode()
Prepare to solve the steady-state problem.
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
vector< int > m_timeSteps
Number of time steps taken in each call to solve() (for example, for each successive grid refinement)
Domain1D & domain(size_t i) const
Return a reference to domain i.
Domain1D * left()
Pointer to left-most domain (first added).
MultiNewton & newton()
Return a reference to the Newton iterator.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
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.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const double Tiny
Small number to compare differences of mole fractions against.
offset
Offsets of solution components in the 1D solution array.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.