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