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