Cantera  2.4.0
Domain1D.h
Go to the documentation of this file.
1  //! @file Domain1D.h
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
6 #ifndef CT_DOMAIN1D_H
7 #define CT_DOMAIN1D_H
8 
11 #include "cantera/base/global.h"
12 #include "refine.h"
13 
14 namespace Cantera
15 {
16 
17 // domain types
18 const int cFlowType = 50;
19 const int cFreeFlow = 51;
20 const int cAxisymmetricStagnationFlow = 52;
21 const int cConnectorType = 100;
22 const int cSurfType = 102;
23 const int cInletType = 104;
24 const int cSymmType = 105;
25 const int cOutletType = 106;
26 const int cEmptyType = 107;
27 const int cOutletResType = 108;
28 const int cPorousType = 109;
29 
30 class MultiJac;
31 class OneDim;
32 class XML_Node;
33 
34 /**
35  * Base class for one-dimensional domains.
36  * @ingroup onedim
37  */
38 class Domain1D
39 {
40 public:
41  /**
42  * Constructor.
43  * @param nv Number of variables at each grid point.
44  * @param points Number of grid points.
45  * @param time (unused)
46  */
47  Domain1D(size_t nv=1, size_t points=1, double time=0.0);
48 
49  virtual ~Domain1D() {}
50  Domain1D(const Domain1D&) = delete;
51  Domain1D& operator=(const Domain1D&) = delete;
52 
53  //! Domain type flag.
54  int domainType() {
55  return m_type;
56  }
57 
58  //! The left-to-right location of this domain.
59  size_t domainIndex() {
60  return m_index;
61  }
62 
63  //! True if the domain is a connector domain.
64  bool isConnector() {
65  return (m_type >= cConnectorType);
66  }
67 
68  //! The container holding this domain.
69  const OneDim& container() const {
70  return *m_container;
71  }
72 
73  //! Specify the container object for this domain, and the position of this
74  //! domain in the list.
75  void setContainer(OneDim* c, size_t index) {
76  m_container = c;
77  m_index = index;
78  }
79 
80  //! Set the Jacobian bandwidth. See the discussion of method bandwidth().
81  void setBandwidth(int bw = -1) {
82  m_bw = bw;
83  }
84 
85  //! Set the Jacobian bandwidth for this domain.
86  /**
87  * When class OneDim computes the bandwidth of the overall multi-domain
88  * problem (in OneDim::resize()), it calls this method for the bandwidth
89  * of each domain. If setBandwidth has not been called, then a negative
90  * bandwidth is returned, in which case OneDim assumes that this domain is
91  * dense -- that is, at each point, all components depend on the value of
92  * all other components at that point. In this case, the bandwidth is bw =
93  * 2*nComponents() - 1. However, if this domain contains some components
94  * that are uncoupled from other components at the same point, then this
95  * default bandwidth may greatly overestimate the true bandwidth, with a
96  * substantial penalty in performance. For such domains, use method
97  * setBandwidth to specify the bandwidth before passing this domain to the
98  * Sim1D or OneDim constructor.
99  */
100  size_t bandwidth() {
101  return m_bw;
102  }
103 
104  /*!
105  * Initialize. This method is called by OneDim::init() for each domain once
106  * at the beginning of a simulation. Base class method does nothing, but may
107  * be overloaded.
108  */
109  virtual void init() { }
110 
111  virtual void setInitialState(doublereal* xlocal = 0) {}
112  virtual void setState(size_t point, const doublereal* state, doublereal* x) {}
113 
114  /*!
115  * When called, this function should reset "bad" values in the state vector
116  * such as negative species concentrations. This function may be called
117  * after a failed solution attempt.
118  */
119  virtual void resetBadValues(double* xg) {}
120 
121  /*!
122  * Resize the domain to have nv components and np grid points. This method
123  * is virtual so that subclasses can perform other actions required to
124  * resize the domain.
125  */
126  virtual void resize(size_t nv, size_t np);
127 
128  //! Return a reference to the grid refiner.
130  return *m_refiner;
131  }
132 
133  //! Number of components at each grid point.
134  size_t nComponents() const {
135  return m_nv;
136  }
137 
138  //! Check that the specified component index is in range.
139  //! Throws an exception if n is greater than nComponents()-1
140  void checkComponentIndex(size_t n) const {
141  if (n >= m_nv) {
142  throw IndexError("checkComponentIndex", "points", n, m_nv-1);
143  }
144  }
145 
146  //! Check that an array size is at least nComponents().
147  //! Throws an exception if nn is less than nComponents(). Used before calls
148  //! which take an array pointer.
149  void checkComponentArraySize(size_t nn) const {
150  if (m_nv > nn) {
151  throw ArraySizeError("checkComponentArraySize", nn, m_nv);
152  }
153  }
154 
155  //! Number of grid points in this domain.
156  size_t nPoints() const {
157  return m_points;
158  }
159 
160  //! Check that the specified point index is in range.
161  //! Throws an exception if n is greater than nPoints()-1
162  void checkPointIndex(size_t n) const {
163  if (n >= m_points) {
164  throw IndexError("checkPointIndex", "points", n, m_points-1);
165  }
166  }
167 
168  //! Check that an array size is at least nPoints().
169  //! Throws an exception if nn is less than nPoints(). Used before calls
170  //! which take an array pointer.
171  void checkPointArraySize(size_t nn) const {
172  if (m_points > nn) {
173  throw ArraySizeError("checkPointArraySize", nn, m_points);
174  }
175  }
176 
177  //! Name of the nth component. May be overloaded.
178  virtual std::string componentName(size_t n) const;
179 
180  void setComponentName(size_t n, const std::string& name) {
181  m_name[n] = name;
182  }
183 
184  //! index of component with name \a name.
185  size_t componentIndex(const std::string& name) const;
186 
187  void setBounds(size_t n, doublereal lower, doublereal upper) {
188  m_min[n] = lower;
189  m_max[n] = upper;
190  }
191 
192  //! Set tolerances for time-stepping mode
193  /*!
194  * @param rtol Relative tolerance
195  * @param atol Absolute tolerance
196  * @param n component index these tolerances apply to. If set to -1 (the
197  * default), these tolerances will be applied to all solution
198  * components.
199  */
200  void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos);
201 
202  //! Set tolerances for steady-state mode
203  /*!
204  * @param rtol Relative tolerance
205  * @param atol Absolute tolerance
206  * @param n component index these tolerances apply to. If set to -1 (the
207  * default), these tolerances will be applied to all solution
208  * components.
209  */
210  void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos);
211 
212  //! Relative tolerance of the nth component.
213  doublereal rtol(size_t n) {
214  return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]);
215  }
216 
217  //! Absolute tolerance of the nth component.
218  doublereal atol(size_t n) {
219  return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]);
220  }
221 
222  //! Steady relative tolerance of the nth component
223  double steady_rtol(size_t n) {
224  return m_rtol_ss[n];
225  }
226 
227  //! Steady absolute tolerance of the nth component
228  double steady_atol(size_t n) {
229  return m_atol_ss[n];
230  }
231 
232  //! Transient relative tolerance of the nth component
233  double transient_rtol(size_t n) {
234  return m_rtol_ts[n];
235  }
236 
237  //! Transient absolute tolerance of the nth component
238  double transient_atol(size_t n) {
239  return m_atol_ts[n];
240  }
241 
242  //! Upper bound on the nth component.
243  doublereal upperBound(size_t n) const {
244  return m_max[n];
245  }
246 
247  //! Lower bound on the nth component
248  doublereal lowerBound(size_t n) const {
249  return m_min[n];
250  }
251 
252  //! Prepare to do time stepping with time step dt
253  /*!
254  * Copy the internally-stored solution at the last time step to array x0.
255  */
256  void initTimeInteg(doublereal dt, const doublereal* x0) {
257  std::copy(x0 + loc(), x0 + loc() + size(), m_slast.begin());
258  m_rdt = 1.0/dt;
259  }
260 
261  //! Prepare to solve the steady-state problem
262  /*!
263  * Set the internally-stored reciprocal of the time step to 0.0
264  */
265  void setSteadyMode() {
266  m_rdt = 0.0;
267  }
268 
269  //! True if in steady-state mode
270  bool steady() {
271  return (m_rdt == 0.0);
272  }
273 
274  //! True if not in steady-state mode
275  bool transient() {
276  return (m_rdt != 0.0);
277  }
278 
279  /*!
280  * Set this if something has changed in the governing
281  * equations (e.g. the value of a constant has been changed,
282  * so that the last-computed Jacobian is no longer valid.
283  */
284  void needJacUpdate();
285 
286  //! Evaluate the residual function at point j. If j == npos,
287  //! evaluate the residual function at all points.
288  /*!
289  * This function must be implemented in classes derived from Domain1D.
290  *
291  * @param j Grid point at which to update the residual
292  * @param[in] x State vector
293  * @param[out] r residual vector
294  * @param[out] mask Boolean mask indicating whether each solution
295  * component has a time derivative (1) or not (0).
296  * @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-
297  * state.)
298  */
299  virtual void eval(size_t j, doublereal* x, doublereal* r,
300  integer* mask, doublereal rdt=0.0) {
301  throw NotImplementedError("Domain1D::eval");
302  }
303 
304  size_t index(size_t n, size_t j) const {
305  return m_nv*j + n;
306  }
307  doublereal value(const doublereal* x, size_t n, size_t j) const {
308  return x[index(n,j)];
309  }
310 
311  virtual void setJac(MultiJac* jac) {}
312 
313  //! Save the current solution for this domain into an XML_Node
314  /*!
315  * Base class version of the general domain1D save function. Derived classes
316  * should call the base class method in addition to saving their own data.
317  *
318  * @param o XML_Node to save the solution to.
319  * @param sol Current value of the solution vector. The object will pick
320  * out which part of the solution vector pertains to this
321  * object.
322  * @return XML_Node created to represent this domain
323  */
324  virtual XML_Node& save(XML_Node& o, const doublereal* const sol);
325 
326  //! Restore the solution for this domain from an XML_Node
327  /*!
328  * Base class version of the general Domain1D restore function. Derived
329  * classes should call the base class method in addition to restoring
330  * their own data.
331  *
332  * @param dom XML_Node for this domain
333  * @param soln Current value of the solution vector, local to this object.
334  * @param loglevel 0 to suppress all output; 1 to show warnings; 2 for
335  * verbose output
336  */
337  virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
338 
339  size_t size() const {
340  return m_nv*m_points;
341  }
342 
343  /**
344  * Find the index of the first grid point in this domain, and
345  * the start of its variables in the global solution vector.
346  */
347  void locate();
348 
349  /**
350  * Location of the start of the local solution vector in the global
351  * solution vector,
352  */
353  virtual size_t loc(size_t j = 0) const {
354  return m_iloc;
355  }
356 
357  /**
358  * The index of the first (i.e., left-most) grid point belonging to this
359  * domain.
360  */
361  size_t firstPoint() const {
362  return m_jstart;
363  }
364 
365  /**
366  * The index of the last (i.e., right-most) grid point belonging to this
367  * domain.
368  */
369  size_t lastPoint() const {
370  return m_jstart + m_points - 1;
371  }
372 
373  /**
374  * Set the left neighbor to domain 'left.' Method 'locate' is called to
375  * update the global positions of this domain and all those to its right.
376  */
378  m_left = left;
379  locate();
380  }
381 
382  //! Set the right neighbor to domain 'right.'
384  m_right = right;
385  }
386 
387  //! Append domain 'right' to this one, and update all links.
389  linkRight(right);
390  right->linkLeft(this);
391  }
392 
393  //! Return a pointer to the left neighbor.
394  Domain1D* left() const {
395  return m_left;
396  }
397 
398  //! Return a pointer to the right neighbor.
399  Domain1D* right() const {
400  return m_right;
401  }
402 
403  //! Value of component n at point j in the previous solution.
404  double prevSoln(size_t n, size_t j) const {
405  return m_slast[m_nv*j + n];
406  }
407 
408  //! Specify an identifying tag for this domain.
409  void setID(const std::string& s) {
410  m_id = s;
411  }
412 
413  std::string id() const {
414  if (m_id != "") {
415  return m_id;
416  } else {
417  return fmt::format("domain {}", m_index);
418  }
419  }
420 
421  virtual void showSolution_s(std::ostream& s, const doublereal* x) {}
422 
423  //! Print the solution.
424  virtual void showSolution(const doublereal* x);
425 
426  doublereal z(size_t jlocal) const {
427  return m_z[jlocal];
428  }
429  doublereal zmin() const {
430  return m_z[0];
431  }
432  doublereal zmax() const {
433  return m_z[m_points - 1];
434  }
435 
436  void setProfile(const std::string& name, double* values, double* soln);
437 
438  vector_fp& grid() {
439  return m_z;
440  }
441  const vector_fp& grid() const {
442  return m_z;
443  }
444  doublereal grid(size_t point) const {
445  return m_z[point];
446  }
447 
448  //! called to set up initial grid, and after grid refinement
449  virtual void setupGrid(size_t n, const doublereal* z);
450 
451  /**
452  * Writes some or all initial solution values into the global solution
453  * array, beginning at the location pointed to by x. This method is called
454  * by the Sim1D constructor, and allows default values or ones that have
455  * been set locally prior to installing this domain into the container to be
456  * written to the global solution vector.
457  */
458  virtual void _getInitialSoln(doublereal* x);
459 
460  //! Initial value of solution component \a n at grid point \a j.
461  virtual doublereal initialValue(size_t n, size_t j);
462 
463  /**
464  * In some cases, a domain may need to set parameters that depend on the
465  * initial solution estimate. In such cases, the parameters may be set in
466  * method _finalize. This method is called just before the Newton solver is
467  * called, and the x array is guaranteed to be the local solution vector for
468  * this domain that will be used as the initial guess. If no such parameters
469  * need to be set, then method _finalize does not need to be overloaded.
470  */
471  virtual void _finalize(const doublereal* x) {}
472 
473  /**
474  * In some cases, for computational efficiency some properties (e.g.
475  * transport coefficients) may not be updated during Jacobian evaluations.
476  * Set this to `true` to force these properties to be udpated even while
477  * calculating Jacobian elements.
478  */
479  void forceFullUpdate(bool update) {
480  m_force_full_update = update;
481  }
482 
483 protected:
484  doublereal m_rdt;
485  size_t m_nv;
486  size_t m_points;
487  vector_fp m_slast;
488  vector_fp m_max;
489  vector_fp m_min;
490  vector_fp m_rtol_ss, m_rtol_ts;
491  vector_fp m_atol_ss, m_atol_ts;
492  vector_fp m_z;
493  OneDim* m_container;
494  size_t m_index;
495  int m_type;
496 
497  //! Starting location within the solution vector for unknowns that
498  //! correspond to this domain
499  /*!
500  * Remember there may be multiple domains associated with this problem
501  */
502  size_t m_iloc;
503 
504  size_t m_jstart;
505 
506  Domain1D* m_left, *m_right;
507 
508  //! Identity tag for the domain
509  std::string m_id;
510  std::string m_desc;
511  std::unique_ptr<Refiner> m_refiner;
512  std::vector<std::string> m_name;
513  int m_bw;
514  bool m_force_full_update;
515 };
516 }
517 
518 #endif
Domain1D(size_t nv=1, size_t points=1, double time=0.0)
Constructor.
Definition: Domain1D.cpp:17
Container class for multiple-domain 1D problems.
Definition: OneDim.h:25
Array size error.
Definition: ctexceptions.h:128
const OneDim & container() const
The container holding this domain.
Definition: Domain1D.h:69
virtual void _getInitialSoln(doublereal *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
Definition: Domain1D.cpp:243
double steady_rtol(size_t n)
Steady relative tolerance of the nth component.
Definition: Domain1D.h:223
void locate()
Find the index of the first grid point in this domain, and the start of its variables in the global s...
Definition: Domain1D.cpp:164
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition: Domain1D.h:265
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:187
size_t m_iloc
Starting location within the solution vector for unknowns that correspond to this domain...
Definition: Domain1D.h:502
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
Domain1D * right() const
Return a pointer to the right neighbor.
Definition: Domain1D.h:399
size_t firstPoint() const
The index of the first (i.e., left-most) grid point belonging to this domain.
Definition: Domain1D.h:361
virtual void eval(size_t j, doublereal *x, doublereal *r, integer *mask, doublereal rdt=0.0)
Evaluate the residual function at point j.
Definition: Domain1D.h:299
void forceFullUpdate(bool update)
In some cases, for computational efficiency some properties (e.g.
Definition: Domain1D.h:479
void checkPointIndex(size_t n) const
Check that the specified point index is in range.
Definition: Domain1D.h:162
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
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: Domain1D.cpp:123
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
virtual XML_Node & save(XML_Node &o, const doublereal *const sol)
Save the current solution for this domain into an XML_Node.
Definition: Domain1D.cpp:110
void linkLeft(Domain1D *left)
Set the left neighbor to domain &#39;left.
Definition: Domain1D.h:377
virtual void resetBadValues(double *xg)
Definition: Domain1D.h:119
void setBandwidth(int bw=-1)
Set the Jacobian bandwidth. See the discussion of method bandwidth().
Definition: Domain1D.h:81
Base class for one-dimensional domains.
Definition: Domain1D.h:38
size_t lastPoint() const
The index of the last (i.e., right-most) grid point belonging to this domain.
Definition: Domain1D.h:369
std::string m_id
Identity tag for the domain.
Definition: Domain1D.h:509
void checkComponentArraySize(size_t nn) const
Check that an array size is at least nComponents().
Definition: Domain1D.h:149
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.cpp:33
size_t bandwidth()
Set the Jacobian bandwidth for this domain.
Definition: Domain1D.h:100
void setID(const std::string &s)
Specify an identifying tag for this domain.
Definition: Domain1D.h:409
Domain1D * left() const
Return a pointer to the left neighbor.
Definition: Domain1D.h:394
virtual void _finalize(const doublereal *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
Definition: Domain1D.h:471
virtual void showSolution(const doublereal *x)
Print the solution.
Definition: Domain1D.cpp:194
doublereal lowerBound(size_t n) const
Lower bound on the nth component.
Definition: Domain1D.h:248
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: Domain1D.cpp:56
Refiner & refiner()
Return a reference to the grid refiner.
Definition: Domain1D.h:129
size_t componentIndex(const std::string &name) const
index of component with name name.
Definition: Domain1D.cpp:65
int domainType()
Domain type flag.
Definition: Domain1D.h:54
size_t domainIndex()
The left-to-right location of this domain.
Definition: Domain1D.h:59
void checkPointArraySize(size_t nn) const
Check that an array size is at least nPoints().
Definition: Domain1D.h:171
double prevSoln(size_t n, size_t j) const
Value of component n at point j in the previous solution.
Definition: Domain1D.h:404
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:134
void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for steady-state mode.
Definition: Domain1D.cpp:89
double transient_atol(size_t n)
Transient absolute tolerance of the nth component.
Definition: Domain1D.h:238
void append(Domain1D *right)
Append domain &#39;right&#39; to this one, and update all links.
Definition: Domain1D.h:388
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
bool steady()
True if in steady-state mode.
Definition: Domain1D.h:270
void checkComponentIndex(size_t n) const
Check that the specified component index is in range.
Definition: Domain1D.h:140
double steady_atol(size_t n)
Steady absolute tolerance of the nth component.
Definition: Domain1D.h:228
Contains declarations for string manipulation functions within Cantera.
double transient_rtol(size_t n)
Transient relative tolerance of the nth component.
Definition: Domain1D.h:233
void linkRight(Domain1D *right)
Set the right neighbor to domain &#39;right.&#39;.
Definition: Domain1D.h:383
doublereal atol(size_t n)
Absolute tolerance of the nth component.
Definition: Domain1D.h:218
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
Definition: refine.h:16
doublereal rtol(size_t n)
Relative tolerance of the nth component.
Definition: Domain1D.h:213
An array index is out of range.
Definition: ctexceptions.h:158
void needJacUpdate()
Definition: Domain1D.cpp:102
bool isConnector()
True if the domain is a connector domain.
Definition: Domain1D.h:64
doublereal upperBound(size_t n) const
Upper bound on the nth component.
Definition: Domain1D.h:243
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:156
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition: Domain1D.cpp:76
virtual void setupGrid(size_t n, const doublereal *z)
called to set up initial grid, and after grid refinement
Definition: Domain1D.cpp:184
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
virtual void init()
Definition: Domain1D.h:109
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:353
virtual doublereal initialValue(size_t n, size_t j)
Initial value of solution component n at grid point j.
Definition: Domain1D.cpp:252
void initTimeInteg(doublereal dt, const doublereal *x0)
Prepare to do time stepping with time step dt.
Definition: Domain1D.h:256