Cantera  2.0
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
16 /**
17  * Container class for multiple-domain 1D problems. Each domain is
18  * represented by an instance of Domain1D.
19  */
20 class OneDim
21 {
22
23 public:
24
25  // Default constructor.
26  OneDim();
27
28  // Constructor.
29  OneDim(std::vector<Domain1D*> domains);
30
31  /// Destructor.
32  virtual ~OneDim();
33
34  /// Add a domain.
35  void addDomain(Domain1D* d);
36
37  //! Return a reference to the Jacobian evaluator.
38  MultiJac& jacobian();
39
40  /// Return a reference to the Newton iterator.
42
43  /**
44  * Solve F(x) = 0, where F(x) is the multi-domain residual function.
45  * @param x0 Starting estimate of solution.
46  * @param x1 Final solution satisfying F(x1) = 0.
47  * @param loglevel Controls amount of diagnostic output.
48  */
49  int solve(doublereal* x0, doublereal* x1, int loglevel);
50
51  /// Number of domains.
52  size_t nDomains() const {
53  return m_nd;
54  }
55
56  /// Return a reference to domain i.
57  Domain1D& domain(size_t i) const {
58  return *m_dom[i];
59  }
60
61  size_t domainIndex(std::string name);
62
63  //! Check that the specified domain index is in range
64  //! Throws an exception if n is greater than nDomains()-1
65  void checkDomainIndex(size_t n) const {
66  if (n >= m_nd) {
67  throw IndexError("checkDomainIndex", "domains", n, m_nd-1);
68  }
69  }
70
71  //! Check that an array size is at least nDomains()
72  //! Throws an exception if nn is less than nDomains(). Used before calls
73  //! which take an array pointer.
74  void checkDomainArraySize(size_t nn) const {
75  if (m_nd > nn) {
76  throw ArraySizeError("checkDomainArraySize", nn, m_nd);
77  }
78  }
79
80  /// The index of the start of domain i in the solution vector.
81  size_t start(size_t i) const {
82  return m_dom[i]->loc();
83  }
84
85  /// Total solution vector length;
86  size_t size() const {
87  return m_size;
88  }
89
90  /// Pointer to left-most domain (first added).
92  return m_dom[0];
93  }
94
95  /// Pointer to right-most domain (last added).
97  return m_dom.back();
98  }
99
100  /// Number of solution components at global point jg.
101  size_t nVars(size_t jg) {
102  return m_nvars[jg];
103  }
104
105  /**
106  * Location in the solution vector of the first component of
107  * global point jg.
108  */
109  size_t loc(size_t jg) {
110  return m_loc[jg];
111  }
112
113  /// Jacobian bandwidth.
114  size_t bandwidth() const {
115  return m_bw;
116  }
117
118  /// Initialize.
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 of the residual evaluated using solution x.
128  * On return, array r contains the steady-state residual values.
129  */
130  doublereal ssnorm(doublereal* x, doublereal* r);
131
132  /// Reciprocal of the time step.
133  doublereal rdt() const {
134  return m_rdt;
135  }
136
137  /// Prepare for time stepping beginning with solution x.
138  void initTimeInteg(doublereal dt, doublereal* x);
139
140  /// True if transient mode.
141  bool transient() const {
142  return (m_rdt != 0.0);
143  }
144
145  /// True if steady mode.
146  bool steady() const {
147  return (m_rdt == 0.0);
148  }
149
150
151  /**
152  * Set steady mode. After invoking this method, subsequent
153  * calls to solve() will solve the steady-state problem.
154  */
156
157
158  /**
159  * Evaluate the multi-domain residual function
160  *
161  * @param j if j > 0, only evaluate residual for points j-1, j,
162  * and j + 1; otherwise, evaluate at all grid points.
163  * @param x solution vector
164  * @param r on return, contains the residual vector
165  * @param rdt Reciprocal of the time step. if omitted, then
166  * the default value is used.
167  * @param count Set to zero to omit this call from the statistics
168  */
169  void eval(size_t j, double* x, double* r, doublereal rdt=-1.0,
170  int count = 1);
171
172  /// Pointer to the domain global point i belongs to.
173  Domain1D* pointDomain(size_t i);
174
175  void resize();
176
177  //doublereal solveTime() { return m_solve_time; }
178
180  vector_int& transientMask() {
182  }
183
184  double timeStep(int nsteps, double dt, double* x,
185  double* r, int loglevel);
186
187  //! Write statistics about the number of iterations and Jacobians at each grid level
188  /*!
189  * @param printTime Boolean that indicates whether time should be printed out
190  * The default is true. It's turned off for test problems where
191  * we don't want to print any times
192  */
193  void writeStats(int printTime = 1);
194
195  void save(std::string fname, std::string id, std::string desc, doublereal* sol);
196
197  // options
198  void setMinTimeStep(doublereal tmin) {
199  m_tmin = tmin;
200  }
201  void setMaxTimeStep(doublereal tmax) {
202  m_tmax = tmax;
203  }
204  void setTimeStepFactor(doublereal tfactor) {
205  m_tfactor = tfactor;
206  }
207  void setJacAge(int ss_age, int ts_age=-1) {
208  m_ss_jac_age = ss_age;
209  if (ts_age > 0) {
210  m_ts_jac_age = ts_age;
211  } else {
212  m_ts_jac_age = m_ss_jac_age;
213  }
214  }
215  void saveStats();
216
217 protected:
218
219  void evalSSJacobian(doublereal* x, doublereal* xnew);
220
221  doublereal m_tmin; // minimum timestep size
222  doublereal m_tmax; // maximum timestep size
223  doublereal m_tfactor; // factor time step is multiplied by
224  // if time stepping fails ( < 1 )
225
226  MultiJac* m_jac; // Jacobian evaluator
227  MultiNewton* m_newt; // Newton iterator
228  doublereal m_rdt; // reciprocal of time step
229  bool m_jac_ok; // if true, Jacobian is current
230
231  //! number of domains
232  size_t m_nd;
233
234  size_t m_bw; // Jacobian bandwidth
235  size_t m_size; // solution vector size
236
237  std::vector<Domain1D*> m_dom, m_connect, m_bulk;
238
239  bool m_init;
240  std::vector<size_t> m_nvars;
241  std::vector<size_t> m_loc;
243  size_t m_pts;
244  doublereal m_solve_time;
245
246  // options
247  int m_ss_jac_age, m_ts_jac_age;
248
249 private:
250
251  // statistics
252  int m_nevals;
253  doublereal m_evaltime;
254  std::vector<size_t> m_gridpts;
255  vector_int m_jacEvals;
256  vector_fp m_jacElapsed;
257  vector_int m_funcEvals;
258  vector_fp m_funcElapsed;
259
260
261 };
262
263 }
264
265 #endif
266
267