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