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