Cantera  2.4.0
OneDim.h
Go to the documentation of this file.
1 /**
2  * @file OneDim.h
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at http://www.cantera.org/license.txt for license and copyright information.
7 
8 #ifndef CT_ONEDIM_H
9 #define CT_ONEDIM_H
10 
11 #include "Domain1D.h"
12 #include "MultiJac.h"
13 
14 namespace Cantera
15 {
16 
17 class Func1;
18 class MultiNewton;
19 
20 /**
21  * Container class for multiple-domain 1D problems. Each domain is
22  * represented by an instance of Domain1D.
23  * @ingroup onedim
24  */
25 class OneDim
26 {
27 public:
28  OneDim();
29 
30  //! Construct a OneDim container for the domains in the list *domains*.
31  OneDim(std::vector<Domain1D*> domains);
32  virtual ~OneDim();
33  OneDim(const OneDim&) = delete;
34  OneDim& operator=(const OneDim&) = delete;
35 
36  /// Add a domain. Domains are added left-to-right.
37  void addDomain(Domain1D* d);
38 
39  //! Return a reference to the Jacobian evaluator.
40  MultiJac& jacobian();
41 
42  /// Return a reference to the Newton iterator.
44 
45  /**
46  * Solve F(x) = 0, where F(x) is the multi-domain residual function.
47  * @param x0 Starting estimate of solution.
48  * @param x1 Final solution satisfying F(x1) = 0.
49  * @param loglevel Controls amount of diagnostic output.
50  */
51  int solve(doublereal* x0, doublereal* x1, int loglevel);
52 
53  /// Number of domains.
54  size_t nDomains() const {
55  return m_dom.size();
56  }
57 
58  /// Return a reference to domain i.
59  Domain1D& domain(size_t i) const {
60  return *m_dom[i];
61  }
62 
63  size_t domainIndex(const std::string& name);
64 
65  //! Check that the specified domain index is in range.
66  //! Throws an exception if n is greater than nDomains()-1
67  void checkDomainIndex(size_t n) const {
68  if (n >= m_dom.size()) {
69  throw IndexError("checkDomainIndex", "domains", n, m_dom.size()-1);
70  }
71  }
72 
73  //! Check that an array size is at least nDomains().
74  //! Throws an exception if nn is less than nDomains(). Used before calls
75  //! which take an array pointer.
76  void checkDomainArraySize(size_t nn) const {
77  if (m_dom.size() > nn) {
78  throw ArraySizeError("checkDomainArraySize", nn, m_dom.size());
79  }
80  }
81 
82  /// The index of the start of domain i in the solution vector.
83  size_t start(size_t i) const {
84  return m_dom[i]->loc();
85  }
86 
87  /// Total solution vector length;
88  size_t size() const {
89  return m_size;
90  }
91 
92  /// Pointer to left-most domain (first added).
94  return m_dom[0];
95  }
96 
97  /// Pointer to right-most domain (last added).
99  return m_dom.back();
100  }
101 
102  /// Number of solution components at global point jg.
103  size_t nVars(size_t jg) {
104  return m_nvars[jg];
105  }
106 
107  //! Location in the solution vector of the first component of global point
108  //! jg.
109  size_t loc(size_t jg) {
110  return m_loc[jg];
111  }
112 
113  //! Return the domain, local point index, and component name for the i-th
114  //! component of the global solution vector
115  std::tuple<std::string, size_t, std::string> component(size_t i);
116 
117  /// Jacobian bandwidth.
118  size_t bandwidth() const {
119  return m_bw;
120  }
121 
122  /*!
123  * Initialize all domains. On the first call, this methods calls the init
124  * method of each domain, proceeding from left to right. Subsequent calls
125  * do nothing.
126  */
127  void init();
128 
129  /// Total number of points.
130  size_t points() {
131  return m_pts;
132  }
133 
134  /**
135  * Steady-state max norm (infinity norm) of the residual evaluated using
136  * solution x. On return, array r contains the steady-state residual
137  * values. Used only for diagnostic output.
138  */
139  doublereal ssnorm(doublereal* x, doublereal* r);
140 
141  /// Reciprocal of the time step.
142  doublereal rdt() const {
143  return m_rdt;
144  }
145 
146  //! Prepare for time stepping beginning with solution *x* and timestep *dt*.
147  void initTimeInteg(doublereal dt, doublereal* x);
148 
149  /// True if transient mode.
150  bool transient() const {
151  return (m_rdt != 0.0);
152  }
153 
154  /// True if steady mode.
155  bool steady() const {
156  return (m_rdt == 0.0);
157  }
158 
159  /*!
160  * Prepare to solve the steady-state problem. After invoking this method,
161  * subsequent calls to solve() will solve the steady-state problem. Sets
162  * the reciprocal of the time step to zero, and, if it was previously non-
163  * zero, signals that a new Jacobian will be needed.
164  */
165  void setSteadyMode();
166 
167  /**
168  * Evaluate the multi-domain residual function
169  *
170  * @param j if j != npos, only evaluate residual for points j-1, j,
171  * and j + 1; otherwise, evaluate at all grid points.
172  * @param x solution vector
173  * @param r on return, contains the residual vector
174  * @param rdt Reciprocal of the time step. if omitted, then
175  * the default value is used.
176  * @param count Set to zero to omit this call from the statistics
177  */
178  void eval(size_t j, double* x, double* r, doublereal rdt=-1.0,
179  int count = 1);
180 
181  //! Return a pointer to the domain global point *i* belongs to.
182  /*!
183  * The domains are scanned right-to-left, and the first one with starting
184  * location less or equal to i is returned.
185  */
186  Domain1D* pointDomain(size_t i);
187 
188  //! Call after one or more grids has changed size, e.g. after being refined.
189  virtual void resize();
190 
191  vector_int& transientMask() {
192  return m_mask;
193  }
194 
195  /*!
196  * Take time steps using Backward Euler.
197  *
198  * @param nsteps number of steps
199  * @param dt initial step size
200  * @param x current solution vector
201  * @param r solution vector after time stepping
202  * @param loglevel controls amount of printed diagnostics
203  * @returns size of last timestep taken
204  */
205  double timeStep(int nsteps, double dt, double* x,
206  double* r, int loglevel);
207 
208  void resetBadValues(double* x);
209 
210  //! Write statistics about the number of iterations and Jacobians at each
211  //! grid level
212  /*!
213  * @param printTime Boolean that indicates whether time should be printed
214  * out The default is true. It's turned off for test
215  * problems where we don't want to print any times
216  */
217  void writeStats(int printTime = 1);
218 
219  void save(const std::string& fname, std::string id,
220  const std::string& desc, doublereal* sol, int loglevel);
221 
222  // options
223  void setMinTimeStep(doublereal tmin) {
224  m_tmin = tmin;
225  }
226  void setMaxTimeStep(doublereal tmax) {
227  m_tmax = tmax;
228  }
229  void setTimeStepFactor(doublereal tfactor) {
230  m_tfactor = tfactor;
231  }
232 
233  //! Set the maximum number of timeteps allowed before successful
234  //! steady-state solve
235  void setMaxTimeStepCount(int nmax) {
236  m_nsteps_max = nmax;
237  }
238 
239  //! Return the maximum number of timeteps allowed before successful
240  //! steady-state solve
241  int maxTimeStepCount() const {
242  return m_nsteps_max;
243  }
244 
245  void setJacAge(int ss_age, int ts_age=-1);
246 
247  /**
248  * Save statistics on function and Jacobian evaluation, and reset the
249  * counters. Statistics are saved only if the number of Jacobian
250  * evaluations is greater than zero. The statistics saved are:
251  *
252  * - number of grid points
253  * - number of Jacobian evaluations
254  * - CPU time spent evaluating Jacobians
255  * - number of non-Jacobian function evaluations
256  * - CPU time spent evaluating functions
257  * - number of time steps
258  */
259  void saveStats();
260 
261  //! Clear saved statistics
262  void clearStats();
263 
264  //! Return total grid size in each call to solve()
265  const std::vector<size_t>& gridSizeStats() {
266  saveStats();
267  return m_gridpts;
268  }
269 
270  //! Return CPU time spent evaluating Jacobians in each call to solve()
272  saveStats();
273  return m_jacElapsed;
274  }
275 
276  //! Return CPU time spent on non-Jacobian function evaluations in each call
277  //! to solve()
279  saveStats();
280  return m_funcElapsed;
281  }
282 
283  //! Return number of Jacobian evaluations made in each call to solve()
285  saveStats();
286  return m_jacEvals;
287  }
288 
289  //! Return number of non-Jacobian function evaluations made in each call to
290  //! solve()
292  saveStats();
293  return m_funcEvals;
294  }
295 
296  //! Return number of time steps taken in each call to solve()
298  saveStats();
299  return m_timeSteps;
300  }
301 
302  //! Set a function that will be called every time #eval is called.
303  //! Can be used to provide keyboard interrupt support in the high-level
304  //! language interfaces.
305  void setInterrupt(Func1* interrupt) {
306  m_interrupt = interrupt;
307  }
308 
309  //! Set a function that will be called after each successful timestep. The
310  //! function will be called with the size of the timestep as the argument.
311  //! Intended to be used for observing solver progress for debugging
312  //! purposes.
313  void setTimeStepCallback(Func1* callback) {
314  m_time_step_callback = callback;
315  }
316 
317 protected:
318  void evalSSJacobian(doublereal* x, doublereal* xnew);
319 
320  doublereal m_tmin; //!< minimum timestep size
321  doublereal m_tmax; //!< maximum timestep size
322 
323  //! factor time step is multiplied by if time stepping fails ( < 1 )
324  doublereal m_tfactor;
325 
326  std::unique_ptr<MultiJac> m_jac; //!< Jacobian evaluator
327  std::unique_ptr<MultiNewton> m_newt; //!< Newton iterator
328  doublereal m_rdt; //!< reciprocal of time step
329  bool m_jac_ok; //!< if true, Jacobian is current
330 
331  size_t m_bw; //!< Jacobian bandwidth
332  size_t m_size; //!< solution vector size
333 
334  std::vector<Domain1D*> m_dom, m_connect, m_bulk;
335 
336  bool m_init;
337  std::vector<size_t> m_nvars;
338  std::vector<size_t> m_loc;
339  vector_int m_mask;
340  size_t m_pts;
341  doublereal m_solve_time;
342 
343  // options
344  int m_ss_jac_age, m_ts_jac_age;
345 
346  //! Function called at the start of every call to #eval.
348 
349  //! User-supplied function called after each successful timestep.
351 
352  //! Number of time steps taken in the current call to solve()
353  int m_nsteps;
354 
355  //! Maximum number of timesteps allowed per call to solve()
357 
358 private:
359  // statistics
360  int m_nevals;
361  doublereal m_evaltime;
362  std::vector<size_t> m_gridpts;
363  vector_int m_jacEvals;
364  vector_fp m_jacElapsed;
365  vector_int m_funcEvals;
366  vector_fp m_funcElapsed;
367 
368  //! Number of time steps taken in each call to solve() (e.g. for each
369  //! successive grid refinement)
371 };
372 
373 }
374 
375 #endif
Container class for multiple-domain 1D problems.
Definition: OneDim.h:25
Array size error.
Definition: ctexceptions.h:128
Domain1D * left()
Pointer to left-most domain (first added).
Definition: OneDim.h:93
void addDomain(Domain1D *d)
Add a domain. Domains are added left-to-right.
Definition: OneDim.cpp:79
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Definition: OneDim.cpp:348
vector_int m_timeSteps
Number of time steps taken in each call to solve() (e.g.
Definition: OneDim.h:370
const vector_int & jacobianCountStats()
Return number of Jacobian evaluations made in each call to solve()
Definition: OneDim.h:284
virtual void resize()
Call after one or more grids has changed size, e.g. after being refined.
Definition: OneDim.cpp:168
doublereal ssnorm(doublereal *x, doublereal *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x. ...
Definition: OneDim.cpp:290
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
void checkDomainIndex(size_t n) const
Check that the specified domain index is in range.
Definition: OneDim.h:67
size_t size() const
Total solution vector length;.
Definition: OneDim.h:88
MultiJac & jacobian()
Return a reference to the Jacobian evaluator.
Definition: OneDim.cpp:101
const vector_fp & jacobianTimeStats()
Return CPU time spent evaluating Jacobians in each call to solve()
Definition: OneDim.h:271
bool steady() const
True if steady mode.
Definition: OneDim.h:155
size_t m_bw
Jacobian bandwidth.
Definition: OneDim.h:331
std::unique_ptr< MultiJac > m_jac
Jacobian evaluator.
Definition: OneDim.h:326
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
size_t loc(size_t jg)
Location in the solution vector of the first component of global point jg.
Definition: OneDim.h:109
MultiNewton & newton()
Return a reference to the Newton iterator.
Definition: OneDim.cpp:105
size_t nDomains() const
Number of domains.
Definition: OneDim.h:54
const vector_int & evalCountStats()
Return number of non-Jacobian function evaluations made in each call to solve()
Definition: OneDim.h:291
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition: OneDim.h:59
Base class for one-dimensional domains.
Definition: Domain1D.h:38
void setMaxTimeStepCount(int nmax)
Set the maximum number of timeteps allowed before successful steady-state solve.
Definition: OneDim.h:235
size_t points()
Total number of points.
Definition: OneDim.h:130
std::unique_ptr< MultiNewton > m_newt
Newton iterator.
Definition: OneDim.h:327
void setTimeStepCallback(Func1 *callback)
Set a function that will be called after each successful timestep.
Definition: OneDim.h:313
const vector_int & timeStepStats()
Return number of time steps taken in each call to solve()
Definition: OneDim.h:297
size_t m_size
solution vector size
Definition: OneDim.h:332
size_t bandwidth() const
Jacobian bandwidth.
Definition: OneDim.h:118
void setSteadyMode()
Definition: OneDim.cpp:319
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
Base class for &#39;functor&#39; classes that evaluate a function of one variable.
Definition: Func1.h:41
doublereal m_tfactor
factor time step is multiplied by if time stepping fails ( < 1 )
Definition: OneDim.h:324
int maxTimeStepCount() const
Return the maximum number of timeteps allowed before successful steady-state solve.
Definition: OneDim.h:241
void checkDomainArraySize(size_t nn) const
Check that an array size is at least nDomains().
Definition: OneDim.h:76
const vector_fp & evalTimeStats()
Return CPU time spent on non-Jacobian function evaluations in each call to solve() ...
Definition: OneDim.h:278
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition: OneDim.h:83
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
Definition: OneDim.cpp:137
void clearStats()
Clear saved statistics.
Definition: OneDim.cpp:155
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
Definition: OneDim.h:356
void writeStats(int printTime=1)
Write statistics about the number of iterations and Jacobians at each grid level. ...
Definition: OneDim.cpp:120
doublereal m_rdt
reciprocal of time step
Definition: OneDim.h:328
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Definition: MultiJac.h:22
void initTimeInteg(doublereal dt, doublereal *x)
Prepare for time stepping beginning with solution x and timestep dt.
Definition: OneDim.cpp:300
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
const std::vector< size_t > & gridSizeStats()
Return total grid size in each call to solve()
Definition: OneDim.h:265
void setInterrupt(Func1 *interrupt)
Set a function that will be called every time eval is called.
Definition: OneDim.h:305
Func1 * m_interrupt
Function called at the start of every call to eval.
Definition: OneDim.h:347
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
Definition: OneDim.h:350
doublereal rdt() const
Reciprocal of the time step.
Definition: OneDim.h:142
Newton iterator for multi-domain, one-dimensional problems.
Definition: MultiNewton.h:19
An array index is out of range.
Definition: ctexceptions.h:158
Domain1D * right()
Pointer to right-most domain (last added).
Definition: OneDim.h:98
size_t nVars(size_t jg)
Number of solution components at global point jg.
Definition: OneDim.h:103
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
Definition: OneDim.cpp:246
int m_nsteps
Number of time steps taken in the current call to solve()
Definition: OneDim.h:353
doublereal m_tmax
maximum timestep size
Definition: OneDim.h:321
bool m_jac_ok
if true, Jacobian is current
Definition: OneDim.h:329
doublereal m_tmin
minimum timestep size
Definition: OneDim.h:320