Cantera  3.1.0a1
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()
37{
38}
39
40size_t OneDim::domainIndex(const string& name)
41{
42 for (size_t n = 0; n < m_dom.size(); n++) {
43 if (domain(n).id() == name) {
44 return n;
45 }
46 }
47 throw CanteraError("OneDim::domainIndex","no domain named >>"+name+"<<");
48}
49
50std::tuple<string, size_t, string> OneDim::component(size_t i) {
51 size_t n;
52 for (n = nDomains()-1; n != npos; n--) {
53 if (i >= start(n)) {
54 break;
55 }
56 }
57 Domain1D& dom = domain(n);
58 size_t offset = i - start(n);
59 size_t pt = offset / dom.nComponents();
60 size_t comp = offset - pt*dom.nComponents();
61 return make_tuple(dom.id(), pt, dom.componentName(comp));
62}
63
64void OneDim::addDomain(shared_ptr<Domain1D> d)
65{
66 // if 'd' is not the first domain, link it to the last domain
67 // added (the rightmost one)
68 size_t n = m_dom.size();
69 if (n > 0) {
70 m_dom.back()->append(d.get());
71 }
72
73 // every other domain is a connector
74 if (n % 2 == 0) {
75 m_connect.push_back(d);
76 } else {
77 m_bulk.push_back(d);
78 }
79
80 // add it also to the global domain list, and set its container and position
81 m_dom.push_back(d);
82 d->setData(m_state);
83 d->setContainer(this, m_dom.size()-1);
84 resize();
85}
86
88{
89 return *m_jac;
90}
92{
93 return *m_newt;
94}
95
96void OneDim::setJacAge(int ss_age, int ts_age)
97{
98 m_ss_jac_age = ss_age;
99 if (ts_age > 0) {
100 m_ts_jac_age = ts_age;
101 } else {
102 m_ts_jac_age = m_ss_jac_age;
103 }
104}
105
106void OneDim::writeStats(int printTime)
107{
108 saveStats();
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++) {
112 if (printTime) {
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]);
116 } else {
117 writelog("{:5d} {:5d} {:6d} NA {:5d} NA\n",
118 m_gridpts[i], m_timeSteps[i], m_funcEvals[i], m_jacEvals[i]);
119 }
120 }
121}
122
124{
125 if (m_jac) {
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);
132 m_nevals = 0;
133 m_funcElapsed.push_back(m_evaltime);
134 m_evaltime = 0.0;
135 m_timeSteps.push_back(m_nsteps);
136 m_nsteps = 0;
137 }
138 }
139}
140
142{
143 m_gridpts.clear();
144 m_jacEvals.clear();
145 m_jacElapsed.clear();
146 m_funcEvals.clear();
147 m_funcElapsed.clear();
148 m_timeSteps.clear();
149 m_nevals = 0;
150 m_evaltime = 0.0;
151 m_nsteps = 0;
152}
153
155{
156 m_bw = 0;
157 m_nvars.clear();
158 m_loc.clear();
159 size_t lc = 0;
160
161 // save the statistics for the last grid
162 saveStats();
163 m_pts = 0;
164 for (size_t i = 0; i < nDomains(); i++) {
165 const auto& d = m_dom[i];
166
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);
171 m_loc.push_back(lc);
172 lc += nv;
173 m_pts++;
174 }
175
176 // update the Jacobian bandwidth
177
178 // bandwidth of the local block
179 size_t bw1 = d->bandwidth();
180 if (bw1 == npos) {
181 bw1 = std::max<size_t>(2*d->nComponents(), 1) - 1;
182 }
183 m_bw = std::max(m_bw, bw1);
184
185 // bandwidth of the block coupling the first point of this
186 // domain to the last point of the previous domain
187 if (i > 0) {
188 size_t bw2 = m_dom[i-1]->bandwidth();
189 if (bw2 == npos) {
190 bw2 = m_dom[i-1]->nComponents();
191 }
192 bw2 += d->nComponents() - 1;
193 m_bw = std::max(m_bw, bw2);
194 }
195 m_size = d->loc() + d->size();
196 }
197
198 m_state->resize(size());
199
200 m_newt->resize(size());
201 m_mask.resize(size());
202
203 // delete the current Jacobian evaluator and create a new one
204 m_jac = make_unique<MultiJac>(*this);
205 m_jac_ok = false;
206
207 for (size_t i = 0; i < nDomains(); i++) {
208 m_dom[i]->setJac(m_jac.get());
209 }
210}
211
212int OneDim::solve(double* x, double* xnew, int loglevel)
213{
214 if (!m_jac_ok) {
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());
218 m_jac_ok = true;
219 }
220
221 return m_newt->solve(x, xnew, *this, *m_jac, loglevel);
222}
223
224void OneDim::evalSSJacobian(double* x, double* xnew)
225{
226 double rdt_save = m_rdt;
227 m_jac_ok = false;
229 eval(npos, x, xnew, 0.0, 0);
230 m_jac->eval(x, xnew, 0.0);
231 m_rdt = rdt_save;
232}
233
235{
236 Domain1D* d = right();
237 while (d) {
238 if (d->loc() <= i) {
239 return d;
240 }
241 d = d->left();
242 }
243 return 0;
244}
245
246void OneDim::eval(size_t j, double* x, double* r, double rdt, int count)
247{
248 clock_t t0 = clock();
249 if (m_interrupt) {
250 m_interrupt->eval(m_nevals);
251 }
252 fill(r, r + m_size, 0.0);
253 if (j == npos) {
254 fill(m_mask.begin(), m_mask.end(), 0);
255 }
256 if (rdt < 0.0) {
257 rdt = m_rdt;
258 }
259
260 // iterate over the bulk domains first
261 for (const auto& d : m_bulk) {
262 d->eval(j, x, r, m_mask.data(), rdt);
263 }
264
265 // then over the connector domains
266 for (const auto& d : m_connect) {
267 d->eval(j, x, r, m_mask.data(), rdt);
268 }
269
270 // increment counter and time
271 if (count) {
272 clock_t t1 = clock();
273 m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
274 m_nevals++;
275 }
276}
277
278double OneDim::ssnorm(double* x, double* r)
279{
280 eval(npos, x, r, 0.0, 0);
281 double ss = 0.0;
282 for (size_t i = 0; i < m_size; i++) {
283 ss = std::max(fabs(r[i]),ss);
284 }
285 return ss;
286}
287
288void OneDim::initTimeInteg(double dt, double* x)
289{
290 double rdt_old = m_rdt;
291 m_rdt = 1.0/dt;
292
293 // if the stepsize has changed, then update the transient part of the
294 // Jacobian
295 if (fabs(rdt_old - m_rdt) > Tiny) {
296 m_jac->updateTransient(m_rdt, m_mask.data());
297 }
298
299 // iterate over all domains, preparing each one to begin time stepping
300 Domain1D* d = left();
301 while (d) {
302 d->initTimeInteg(dt, x);
303 d = d->right();
304 }
305}
306
308{
309 if (m_rdt == 0) {
310 return;
311 }
312
313 m_rdt = 0.0;
314 m_jac->updateTransient(m_rdt, m_mask.data());
315
316 // iterate over all domains, preparing them for steady-state solution
317 Domain1D* d = left();
318 while (d) {
319 d->setSteadyMode();
320 d = d->right();
321 }
322}
323
325{
326 if (!m_init) {
327 Domain1D* d = left();
328 while (d) {
329 d->init();
330 d = d->right();
331 }
332 }
333 m_init = true;
334}
335
336double OneDim::timeStep(int nsteps, double dt, double* x, double* r, int loglevel)
337{
338 // set the Jacobian age parameter to the transient value
339 newton().setOptions(m_ts_jac_age);
340
341 debuglog("\n\n step size (s) log10(ss) \n", loglevel);
342 debuglog("===============================\n", loglevel);
343
344 int n = 0;
345 int successiveFailures = 0;
346
347 while (n < nsteps) {
348 if (loglevel > 0) {
349 double ss = ssnorm(x, r);
350 writelog(" {:>4d} {:10.4g} {:10.4g}", n, dt, log10(ss));
351 }
352
353 // set up for time stepping with stepsize dt
354 initTimeInteg(dt,x);
355
356 // solve the transient problem
357 int m = solve(x, r, loglevel-1);
358
359 // successful time step. Copy the new solution in r to
360 // the current solution in x.
361 if (m >= 0) {
362 successiveFailures = 0;
363 m_nsteps++;
364 n += 1;
365 debuglog("\n", loglevel);
366 copy(r, r + m_size, x);
367 if (m == 100) {
368 dt *= 1.5;
369 }
372 }
373 dt = std::min(dt, m_tmax);
374 if (m_nsteps >= m_nsteps_max) {
375 throw CanteraError("OneDim::timeStep",
376 "Took maximum number of timesteps allowed ({}) without "
377 "reaching steady-state solution.", m_nsteps_max);
378 }
379 } else {
380 successiveFailures++;
381 // No solution could be found with this time step.
382 // Decrease the stepsize and try again.
383 debuglog("...failure.\n", loglevel);
384 if (successiveFailures > 2) {
385 //debuglog("Resetting negative species concentrations.\n", loglevel);
386 resetBadValues(x);
387 successiveFailures = 0;
388 } else {
389 dt *= m_tfactor;
390 if (dt < m_tmin) {
391 throw CanteraError("OneDim::timeStep",
392 "Time integration failed.");
393 }
394 }
395 }
396 }
397
398 // return the value of the last stepsize, which may be smaller
399 // than the initial stepsize
400 return dt;
401}
402
403void OneDim::resetBadValues(double* x)
404{
405 for (auto dom : m_dom) {
406 dom->resetBadValues(x);
407 }
408}
409
410}
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
Definition Domain1D.h:28
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:145
Domain1D * left() const
Return a pointer to the left neighbor.
Definition Domain1D.h:431
Domain1D * right() const
Return a pointer to the right neighbor.
Definition Domain1D.h:436
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition Domain1D.cpp:49
virtual void init()
Initialize.
Definition Domain1D.h:120
void initTimeInteg(double dt, const double *x0)
Prepare to do time stepping with time step dt.
Definition Domain1D.h:267
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition Domain1D.h:276
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:384
virtual double eval(double t) const
Evaluate the function.
Definition Func1.cpp:28
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:212
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition OneDim.h:88
int m_nsteps
Number of time steps taken in the current call to solve()
Definition OneDim.h:363
void init()
Initialize all domains.
Definition OneDim.cpp:324
size_t m_size
solution vector size
Definition OneDim.h:340
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
Definition OneDim.h:366
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
Definition OneDim.cpp:154
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
Definition OneDim.cpp:123
unique_ptr< MultiNewton > m_newt
Newton iterator.
Definition OneDim.h:335
size_t size() const
Total solution vector length;.
Definition OneDim.h:99
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:246
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition OneDim.cpp:278
void addDomain(shared_ptr< Domain1D > d)
Add a domain. Domains are added left-to-right.
Definition OneDim.cpp:64
double rdt() const
Reciprocal of the time step.
Definition OneDim.h:153
void initTimeInteg(double dt, double *x)
Prepare for time stepping beginning with solution x and timestep dt.
Definition OneDim.cpp:288
size_t nDomains() const
Number of domains.
Definition OneDim.h:57
Domain1D * right()
Pointer to right-most domain (last added).
Definition OneDim.h:109
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:50
double m_rdt
reciprocal of time step
Definition OneDim.h:336
shared_ptr< vector< double > > m_state
Solution vector.
Definition OneDim.h:332
Func1 * m_interrupt
Function called at the start of every call to eval.
Definition OneDim.h:357
unique_ptr< MultiJac > m_jac
Jacobian evaluator.
Definition OneDim.h:334
size_t m_bw
Jacobian bandwidth.
Definition OneDim.h:339
double m_tfactor
factor time step is multiplied by if time stepping fails ( < 1 )
Definition OneDim.h:330
bool m_jac_ok
if true, Jacobian is current
Definition OneDim.h:337
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
Definition OneDim.cpp:336
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
Definition OneDim.cpp:234
double m_tmin
minimum timestep size
Definition OneDim.h:326
double m_tmax
maximum timestep size
Definition OneDim.h:327
void writeStats(int printTime=1)
Write statistics about the number of iterations and Jacobians at each grid level.
Definition OneDim.cpp:106
void clearStats()
Clear saved statistics.
Definition OneDim.cpp:141
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition OneDim.cpp:307
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
Definition OneDim.h:360
vector< int > m_timeSteps
Number of time steps taken in each call to solve() (for example, for each successive grid refinement)
Definition OneDim.h:380
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition OneDim.h:62
Domain1D * left()
Pointer to left-most domain (first added).
Definition OneDim.h:104
MultiNewton & newton()
Return a reference to the Newton iterator.
Definition OneDim.cpp:91
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
Definition OneDim.cpp:87
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:180
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