21 m_newt = make_unique<MultiNewton>(1);
24 OneDim::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) {
40 size_t OneDim::domainIndex(
const string& name)
42 for (
size_t n = 0; n < m_dom.size(); n++) {
43 if (domain(n).
id() == name) {
47 throw CanteraError(
"OneDim::domainIndex",
"no domain named >>"+name+
"<<");
50 std::tuple<string, size_t, string> OneDim::component(
size_t i) {
52 for (n = nDomains()-1; n !=
npos; n--) {
58 size_t offset = i - start(n);
64 void OneDim::addDomain(shared_ptr<Domain1D> d)
68 size_t n = m_dom.size();
70 m_dom.back()->append(d.get());
75 m_connect.push_back(d);
83 d->setContainer(
this, m_dom.size()-1);
96 void OneDim::setJacAge(
int ss_age,
int ts_age)
98 m_ss_jac_age = ss_age;
100 m_ts_jac_age = ts_age;
102 m_ts_jac_age = m_ss_jac_age;
106 void OneDim::writeStats(
int printTime)
109 writelog(
"\nStatistics:\n\n Grid Timesteps Functions Time Jacobians Time\n");
110 size_t n = m_gridpts.size();
111 for (
size_t i = 0; i < n; i++) {
113 writelog(
"{:5d} {:5d} {:6d} {:9.4f} {:5d} {:9.4f}\n",
114 m_gridpts[i], m_timeSteps[i], m_funcEvals[i], m_funcElapsed[i],
115 m_jacEvals[i], m_jacElapsed[i]);
117 writelog(
"{:5d} {:5d} {:6d} NA {:5d} NA\n",
118 m_gridpts[i], m_timeSteps[i], m_funcEvals[i], m_jacEvals[i]);
123 void OneDim::saveStats()
126 int nev = m_jac->nEvals();
127 if (nev > 0 && m_nevals > 0) {
128 m_gridpts.push_back(m_pts);
129 m_jacEvals.push_back(m_jac->nEvals());
130 m_jacElapsed.push_back(m_jac->elapsedTime());
131 m_funcEvals.push_back(m_nevals);
133 m_funcElapsed.push_back(m_evaltime);
135 m_timeSteps.push_back(m_nsteps);
141 void OneDim::clearStats()
145 m_jacElapsed.clear();
147 m_funcElapsed.clear();
154 void OneDim::resize()
164 for (
size_t i = 0; i < nDomains(); i++) {
165 const auto& d = m_dom[i];
167 size_t np = d->nPoints();
168 size_t nv = d->nComponents();
169 for (
size_t n = 0; n < np; n++) {
170 m_nvars.push_back(nv);
179 size_t bw1 = d->bandwidth();
181 bw1 = std::max<size_t>(2*d->nComponents(), 1) - 1;
183 m_bw = std::max(m_bw, bw1);
188 size_t bw2 = m_dom[i-1]->bandwidth();
190 bw2 = m_dom[i-1]->nComponents();
192 bw2 += d->nComponents() - 1;
193 m_bw = std::max(m_bw, bw2);
195 m_size = d->loc() + d->size();
198 m_state->resize(size());
200 m_newt->resize(size());
201 m_mask.resize(size());
204 m_jac = make_unique<MultiJac>(*
this);
207 for (
size_t i = 0; i < nDomains(); i++) {
208 m_dom[i]->setJac(m_jac.get());
215 eval(
npos, x, xnew, 0.0, 0);
216 m_jac->eval(x, xnew, 0.0);
217 m_jac->updateTransient(m_rdt, m_mask.data());
221 return m_newt->solve(x, xnew, *
this, *m_jac, loglevel);
224 void OneDim::evalSSJacobian(
double* x,
double* xnew)
226 double rdt_save = m_rdt;
229 eval(
npos, x, xnew, 0.0, 0);
230 m_jac->eval(x, xnew, 0.0);
246 void OneDim::eval(
size_t j,
double* x,
double* r,
double rdt,
int count)
248 clock_t t0 = clock();
250 m_interrupt->eval(m_nevals);
252 fill(r, r + m_size, 0.0);
254 fill(m_mask.begin(), m_mask.end(), 0);
261 for (
const auto& d : m_bulk) {
262 d->eval(j, x, r, m_mask.data(), rdt);
266 for (
const auto& d : m_connect) {
267 d->eval(j, x, r, m_mask.data(), rdt);
272 clock_t t1 = clock();
273 m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
278 double OneDim::ssnorm(
double* x,
double* r)
280 eval(
npos, x, r, 0.0, 0);
282 for (
size_t i = 0; i < m_size; i++) {
283 ss = std::max(fabs(r[i]),ss);
288 void OneDim::initTimeInteg(
double dt,
double* x)
290 double rdt_old = m_rdt;
295 if (fabs(rdt_old - m_rdt) >
Tiny) {
296 m_jac->updateTransient(m_rdt, m_mask.data());
307 void OneDim::setSteadyMode()
314 m_jac->updateTransient(m_rdt, m_mask.data());
336 double OneDim::timeStep(
int nsteps,
double dt,
double* x,
double* r,
int loglevel)
339 newton().setOptions(m_ts_jac_age);
341 debuglog(
"\n\n step size (s) log10(ss) \n", loglevel);
342 debuglog(
"===============================\n", loglevel);
345 int successiveFailures = 0;
349 double ss = ssnorm(x, r);
350 writelog(
" {:>4d} {:10.4g} {:10.4g}", n, dt, log10(ss));
357 int m =
solve(x, r, loglevel-1);
362 successiveFailures = 0;
366 copy(r, r + m_size, x);
370 if (m_time_step_callback) {
371 m_time_step_callback->eval(dt);
373 dt = std::min(dt, m_tmax);
374 if (m_nsteps >= m_nsteps_max) {
376 "Took maximum number of timesteps allowed ({}) without "
377 "reaching steady-state solution.", m_nsteps_max);
380 successiveFailures++;
383 debuglog(
"...failure.\n", loglevel);
384 if (successiveFailures > 2) {
387 successiveFailures = 0;
392 "Time integration failed.");
403 void OneDim::resetBadValues(
double* x)
405 for (
auto dom : m_dom) {
406 dom->resetBadValues(x);
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.
Domain1D * left() const
Return a pointer to the left neighbor.
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 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,.
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 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.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.