Cantera  3.1.0a1
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 
6 #include "cantera/oneD/OneDim.h"
9 #include "cantera/base/AnyMap.h"
10 
11 #include <fstream>
12 #include <ctime>
13 
14 using namespace std;
15 
16 namespace Cantera
17 {
18 
19 OneDim::OneDim()
20 {
21  m_newt = make_unique<MultiNewton>(1);
22 }
23 
24 OneDim::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 
36 OneDim::~OneDim()
37 {
38 }
39 
40 size_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 
50 std::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 
64 void 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 
87 MultiJac& OneDim::jacobian()
88 {
89  return *m_jac;
90 }
91 MultiNewton& OneDim::newton()
92 {
93  return *m_newt;
94 }
95 
96 void 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 
106 void 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 
123 void OneDim::saveStats()
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 
141 void OneDim::clearStats()
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 
154 void OneDim::resize()
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 
212 int 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 
224 void OneDim::evalSSJacobian(double* x, double* xnew)
225 {
226  double rdt_save = m_rdt;
227  m_jac_ok = false;
228  setSteadyMode();
229  eval(npos, x, xnew, 0.0, 0);
230  m_jac->eval(x, xnew, 0.0);
231  m_rdt = rdt_save;
232 }
233 
234 Domain1D* OneDim::pointDomain(size_t i)
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 
246 void 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 
278 double 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 
288 void 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 
307 void OneDim::setSteadyMode()
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 
324 void OneDim::init()
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 
336 double 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  }
370  if (m_time_step_callback) {
371  m_time_step_callback->eval(dt);
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 
403 void 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.
Definition: ctexceptions.h:66
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
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 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
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.