Cantera  3.0.0
Loading...
Searching...
No Matches
OneDim.cpp
Go to the documentation of this file.
1//! @file OneDim.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
10
11#include <fstream>
12#include <ctime>
13
14using namespace std;
15
16namespace Cantera
17{
18
19OneDim::OneDim()
20{
21 m_newt = make_unique<MultiNewton>(1);
22}
23
24OneDim::OneDim(vector<shared_ptr<Domain1D>>& domains)
25{
26 // create a Newton iterator, and add each domain.
27 m_newt = make_unique<MultiNewton>(1);
28 m_state = make_shared<vector<double>>();
29 for (auto& dom : domains) {
30 addDomain(dom);
31 }
32 init();
33 resize();
34}
35
36OneDim::OneDim(vector<Domain1D*> domains)
37{
38 warn_deprecated("OneDim::OneDim(vector<Domain1D*>)",
39 "To be removed after Cantera 3.0; superseded by "
40 "OneDim::OneDim(vector<shared_ptr<Domain1D>>).");
41
42 // create a Newton iterator, and add each domain.
43 m_newt = make_unique<MultiNewton>(1);
44 m_state = make_shared<vector<double>>();
45 for (size_t i = 0; i < domains.size(); i++) {
46 addDomain(domains[i]);
47 }
48 init();
49 resize();
50}
51
52OneDim::~OneDim()
53{
54}
55
56size_t OneDim::domainIndex(const string& name)
57{
58 for (size_t n = 0; n < m_dom.size(); n++) {
59 if (domain(n).id() == name) {
60 return n;
61 }
62 }
63 throw CanteraError("OneDim::domainIndex","no domain named >>"+name+"<<");
64}
65
66std::tuple<string, size_t, string> OneDim::component(size_t i) {
67 size_t n;
68 for (n = nDomains()-1; n != npos; n--) {
69 if (i >= start(n)) {
70 break;
71 }
72 }
73 Domain1D& dom = domain(n);
74 size_t offset = i - start(n);
75 size_t pt = offset / dom.nComponents();
76 size_t comp = offset - pt*dom.nComponents();
77 return make_tuple(dom.id(), pt, dom.componentName(comp));
78}
79
80void OneDim::addDomain(shared_ptr<Domain1D> d)
81{
82 // if 'd' is not the first domain, link it to the last domain
83 // added (the rightmost one)
84 size_t n = m_dom.size();
85 if (n > 0) {
86 m_dom.back()->append(d.get());
87 }
88
89 // every other domain is a connector
90 if (n % 2 == 0) {
91 m_sharedConnect.push_back(d);
92 m_connect.push_back(d.get());
93 } else {
94 m_sharedBulk.push_back(d);
95 m_bulk.push_back(d.get());
96 }
97
98 // add it also to the global domain list, and set its container and position
99 m_sharedDom.push_back(d);
100 m_dom.push_back(d.get());
101 d->setData(m_state);
102 d->setContainer(this, m_dom.size()-1);
103 resize();
104}
105
107{
108 warn_deprecated("OneDim::addDomain(Domain1D*)",
109 "To be removed after Cantera 3.0; superseded by "
110 "OneDim::addDomain(shared_ptr<Domain1D>).");
111
112 // if 'd' is not the first domain, link it to the last domain
113 // added (the rightmost one)
114 size_t n = m_dom.size();
115 if (n > 0) {
116 m_dom.back()->append(d);
117 }
118
119 // every other domain is a connector
120 if (n % 2 == 0) {
121 m_connect.push_back(d);
122 } else {
123 m_bulk.push_back(d);
124 }
125
126 // add it also to the global domain list, and set its container and position
127 m_dom.push_back(d);
128 d->setData(m_state);
129 d->setContainer(this, m_dom.size()-1);
130 resize();
131}
132
134{
135 return *m_jac;
136}
138{
139 return *m_newt;
140}
141
142void OneDim::setJacAge(int ss_age, int ts_age)
143{
144 m_ss_jac_age = ss_age;
145 if (ts_age > 0) {
146 m_ts_jac_age = ts_age;
147 } else {
148 m_ts_jac_age = m_ss_jac_age;
149 }
150}
151
152void OneDim::writeStats(int printTime)
153{
154 saveStats();
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++) {
158 if (printTime) {
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]);
162 } else {
163 writelog("{:5d} {:5d} {:6d} NA {:5d} NA\n",
164 m_gridpts[i], m_timeSteps[i], m_funcEvals[i], m_jacEvals[i]);
165 }
166 }
167}
168
170{
171 if (m_jac) {
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);
178 m_nevals = 0;
179 m_funcElapsed.push_back(m_evaltime);
180 m_evaltime = 0.0;
181 m_timeSteps.push_back(m_nsteps);
182 m_nsteps = 0;
183 }
184 }
185}
186
188{
189 m_gridpts.clear();
190 m_jacEvals.clear();
191 m_jacElapsed.clear();
192 m_funcEvals.clear();
193 m_funcElapsed.clear();
194 m_timeSteps.clear();
195 m_nevals = 0;
196 m_evaltime = 0.0;
197 m_nsteps = 0;
198}
199
201{
202 m_bw = 0;
203 m_nvars.clear();
204 m_loc.clear();
205 size_t lc = 0;
206
207 // save the statistics for the last grid
208 saveStats();
209 m_pts = 0;
210 for (size_t i = 0; i < nDomains(); i++) {
211 Domain1D* d = m_dom[i];
212
213 size_t np = d->nPoints();
214 size_t nv = d->nComponents();
215 for (size_t n = 0; n < np; n++) {
216 m_nvars.push_back(nv);
217 m_loc.push_back(lc);
218 lc += nv;
219 m_pts++;
220 }
221
222 // update the Jacobian bandwidth
223
224 // bandwidth of the local block
225 size_t bw1 = d->bandwidth();
226 if (bw1 == npos) {
227 bw1 = std::max<size_t>(2*d->nComponents(), 1) - 1;
228 }
229 m_bw = std::max(m_bw, bw1);
230
231 // bandwidth of the block coupling the first point of this
232 // domain to the last point of the previous domain
233 if (i > 0) {
234 size_t bw2 = m_dom[i-1]->bandwidth();
235 if (bw2 == npos) {
236 bw2 = m_dom[i-1]->nComponents();
237 }
238 bw2 += d->nComponents() - 1;
239 m_bw = std::max(m_bw, bw2);
240 }
241 m_size = d->loc() + d->size();
242 }
243
244 m_state->resize(size());
245
246 m_newt->resize(size());
247 m_mask.resize(size());
248
249 // delete the current Jacobian evaluator and create a new one
250 m_jac = make_unique<MultiJac>(*this);
251 m_jac_ok = false;
252
253 for (size_t i = 0; i < nDomains(); i++) {
254 m_dom[i]->setJac(m_jac.get());
255 }
256}
257
258int OneDim::solve(double* x, double* xnew, int loglevel)
259{
260 if (!m_jac_ok) {
261 eval(npos, x, xnew, 0.0, 0);
262 m_jac->eval(x, xnew, 0.0);
263 m_jac->updateTransient(m_rdt, m_mask.data());
264 m_jac_ok = true;
265 }
266
267 return m_newt->solve(x, xnew, *this, *m_jac, loglevel);
268}
269
270void OneDim::evalSSJacobian(double* x, double* xnew)
271{
272 double rdt_save = m_rdt;
273 m_jac_ok = false;
275 eval(npos, x, xnew, 0.0, 0);
276 m_jac->eval(x, xnew, 0.0);
277 m_rdt = rdt_save;
278}
279
281{
282 Domain1D* d = right();
283 while (d) {
284 if (d->loc() <= i) {
285 return d;
286 }
287 d = d->left();
288 }
289 return 0;
290}
291
292void OneDim::eval(size_t j, double* x, double* r, double rdt, int count)
293{
294 clock_t t0 = clock();
295 if (m_interrupt) {
296 m_interrupt->eval(m_nevals);
297 }
298 fill(r, r + m_size, 0.0);
299 if (j == npos) {
300 fill(m_mask.begin(), m_mask.end(), 0);
301 }
302 if (rdt < 0.0) {
303 rdt = m_rdt;
304 }
305
306 // iterate over the bulk domains first
307 for (const auto& d : m_bulk) {
308 d->eval(j, x, r, m_mask.data(), rdt);
309 }
310
311 // then over the connector domains
312 for (const auto& d : m_connect) {
313 d->eval(j, x, r, m_mask.data(), rdt);
314 }
315
316 // increment counter and time
317 if (count) {
318 clock_t t1 = clock();
319 m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
320 m_nevals++;
321 }
322}
323
324double OneDim::ssnorm(double* x, double* r)
325{
326 eval(npos, x, r, 0.0, 0);
327 double ss = 0.0;
328 for (size_t i = 0; i < m_size; i++) {
329 ss = std::max(fabs(r[i]),ss);
330 }
331 return ss;
332}
333
334void OneDim::initTimeInteg(double dt, double* x)
335{
336 double rdt_old = m_rdt;
337 m_rdt = 1.0/dt;
338
339 // if the stepsize has changed, then update the transient part of the
340 // Jacobian
341 if (fabs(rdt_old - m_rdt) > Tiny) {
342 m_jac->updateTransient(m_rdt, m_mask.data());
343 }
344
345 // iterate over all domains, preparing each one to begin time stepping
346 Domain1D* d = left();
347 while (d) {
348 d->initTimeInteg(dt, x);
349 d = d->right();
350 }
351}
352
354{
355 if (m_rdt == 0) {
356 return;
357 }
358
359 m_rdt = 0.0;
360 m_jac->updateTransient(m_rdt, m_mask.data());
361
362 // iterate over all domains, preparing them for steady-state solution
363 Domain1D* d = left();
364 while (d) {
365 d->setSteadyMode();
366 d = d->right();
367 }
368}
369
371{
372 if (!m_init) {
373 Domain1D* d = left();
374 while (d) {
375 d->init();
376 d = d->right();
377 }
378 }
379 m_init = true;
380}
381
382double OneDim::timeStep(int nsteps, double dt, double* x, double* r, int loglevel)
383{
384 // set the Jacobian age parameter to the transient value
385 newton().setOptions(m_ts_jac_age);
386
387 debuglog("\n\n step size (s) log10(ss) \n", loglevel);
388 debuglog("===============================\n", loglevel);
389
390 int n = 0;
391 int successiveFailures = 0;
392
393 while (n < nsteps) {
394 if (loglevel > 0) {
395 double ss = ssnorm(x, r);
396 writelog(" {:>4d} {:10.4g} {:10.4g}", n, dt, log10(ss));
397 }
398
399 // set up for time stepping with stepsize dt
400 initTimeInteg(dt,x);
401
402 // solve the transient problem
403 int m = solve(x, r, loglevel-1);
404
405 // successful time step. Copy the new solution in r to
406 // the current solution in x.
407 if (m >= 0) {
408 successiveFailures = 0;
409 m_nsteps++;
410 n += 1;
411 debuglog("\n", loglevel);
412 copy(r, r + m_size, x);
413 if (m == 100) {
414 dt *= 1.5;
415 }
418 }
419 dt = std::min(dt, m_tmax);
420 if (m_nsteps >= m_nsteps_max) {
421 throw CanteraError("OneDim::timeStep",
422 "Took maximum number of timesteps allowed ({}) without "
423 "reaching steady-state solution.", m_nsteps_max);
424 }
425 } else {
426 successiveFailures++;
427 // No solution could be found with this time step.
428 // Decrease the stepsize and try again.
429 debuglog("...failure.\n", loglevel);
430 if (successiveFailures > 2) {
431 //debuglog("Resetting negative species concentrations.\n", loglevel);
432 resetBadValues(x);
433 successiveFailures = 0;
434 } else {
435 dt *= m_tfactor;
436 if (dt < m_tmin) {
437 throw CanteraError("OneDim::timeStep",
438 "Time integration failed.");
439 }
440 }
441 }
442 }
443
444 // return the value of the last stepsize, which may be smaller
445 // than the initial stepsize
446 return dt;
447}
448
449void OneDim::resetBadValues(double* x)
450{
451 for (auto dom : m_dom) {
452 dom->resetBadValues(x);
453 }
454}
455
456AnyMap OneDim::serialize(const double* soln) const
457{
458 warn_deprecated("OneDim::serialize",
459 "To be removed after Cantera 3.0; unused.");
460 AnyMap state;
461 for (size_t i = 0; i < m_dom.size(); i++) {
462 state[m_dom[i]->id()] = m_dom[i]->serialize(soln + start(i));
463 }
464 return state;
465}
466
467}
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
Definition Domain1D.h:41
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:160
size_t bandwidth()
Set the Jacobian bandwidth for this domain.
Definition Domain1D.h:126
Domain1D * left() const
Return a pointer to the left neighbor.
Definition Domain1D.h:465
void setContainer(OneDim *c, size_t index)
Specify the container object for this domain, and the position of this domain in the list.
Definition Domain1D.h:101
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:182
Domain1D * right() const
Return a pointer to the right neighbor.
Definition Domain1D.h:470
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition Domain1D.cpp:56
virtual void init()
Initialize.
Definition Domain1D.h:135
void initTimeInteg(double dt, const double *x0)
Prepare to do time stepping with time step dt.
Definition Domain1D.h:282
void setData(shared_ptr< vector< double > > &data)
Set shared data pointer.
Definition Domain1D.h:563
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition Domain1D.h:291
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
Definition Domain1D.h:418
virtual double eval(double t) const
Evaluate the function.
Definition Func1.cpp:61
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Definition MultiJac.h:24
Newton iterator for multi-domain, one-dimensional problems.
Definition MultiNewton.h:24
void setOptions(int maxJacAge=5)
Set options.
Definition MultiNewton.h:68
int solve(double *x0, double *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
Definition OneDim.cpp:258
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition OneDim.h:95
int m_nsteps
Number of time steps taken in the current call to solve()
Definition OneDim.h:377
void init()
Initialize all domains.
Definition OneDim.cpp:370
vector< Domain1D * > m_bulk
Definition OneDim.h:358
size_t m_size
solution vector size
Definition OneDim.h:350
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
Definition OneDim.h:380
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
Definition OneDim.cpp:200
AnyMap serialize(const double *soln) const
Definition OneDim.cpp:456
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
Definition OneDim.cpp:169
unique_ptr< MultiNewton > m_newt
Newton iterator.
Definition OneDim.h:345
size_t size() const
Total solution vector length;.
Definition OneDim.h:106
void eval(size_t j, double *x, double *r, double rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition OneDim.cpp:292
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition OneDim.cpp:324
void addDomain(shared_ptr< Domain1D > d)
Add a domain. Domains are added left-to-right.
Definition OneDim.cpp:80
double rdt() const
Reciprocal of the time step.
Definition OneDim.h:160
void initTimeInteg(double dt, double *x)
Prepare for time stepping beginning with solution x and timestep dt.
Definition OneDim.cpp:334
size_t nDomains() const
Number of domains.
Definition OneDim.h:64
Domain1D * right()
Pointer to right-most domain (last added).
Definition OneDim.h:116
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...
Definition OneDim.cpp:66
double m_rdt
reciprocal of time step
Definition OneDim.h:346
shared_ptr< vector< double > > m_state
Solution vector.
Definition OneDim.h:342
Func1 * m_interrupt
Function called at the start of every call to eval.
Definition OneDim.h:371
vector< Domain1D * > m_connect
Definition OneDim.h:357
unique_ptr< MultiJac > m_jac
Jacobian evaluator.
Definition OneDim.h:344
size_t m_bw
Jacobian bandwidth.
Definition OneDim.h:349
double m_tfactor
factor time step is multiplied by if time stepping fails ( < 1 )
Definition OneDim.h:340
vector< Domain1D * > m_dom
Definition OneDim.h:356
bool m_jac_ok
if true, Jacobian is current
Definition OneDim.h:347
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
Definition OneDim.cpp:382
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
Definition OneDim.cpp:280
double m_tmin
minimum timestep size
Definition OneDim.h:336
double m_tmax
maximum timestep size
Definition OneDim.h:337
void writeStats(int printTime=1)
Write statistics about the number of iterations and Jacobians at each grid level.
Definition OneDim.cpp:152
void clearStats()
Clear saved statistics.
Definition OneDim.cpp:187
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition OneDim.cpp:353
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
Definition OneDim.h:374
vector< int > m_timeSteps
Number of time steps taken in each call to solve() (for example, for each successive grid refinement)
Definition OneDim.h:394
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition OneDim.h:69
Domain1D * left()
Pointer to left-most domain (first added).
Definition OneDim.h:111
MultiNewton & newton()
Return a reference to the Newton iterator.
Definition OneDim.cpp:137
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
Definition OneDim.cpp:133
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:158
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
offset
Offsets of solution components in the 1D solution array.
Definition StFlow.h:24
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926