Cantera  3.2.0a2
Loading...
Searching...
No Matches
SteadyStateSystem.h
Go to the documentation of this file.
1//! @file SteadyStateSystem.h
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
6#ifndef CT_STEADYSTATESYSTEM_H
7#define CT_STEADYSTATESYSTEM_H
8
10#include "SystemJacobian.h"
11
12namespace Cantera
13{
14
15class Func1;
16class MultiNewton;
17
18//! Base class for representing a system of differential-algebraic equations and solving
19//! for its steady-state response.
20//!
21//! @since New in %Cantera 3.2.
22//! @ingroup numerics
24{
25public:
27 virtual ~SteadyStateSystem();
28 SteadyStateSystem(const SteadyStateSystem&) = delete;
29 SteadyStateSystem& operator=(const SteadyStateSystem&) = delete;
30
31 //! Evaluate the residual function
32 //!
33 //! @param[in] x State vector
34 //! @param[out] r On return, contains the residual vector
35 //! @param rdt Reciprocal of the time step. if omitted, then the internally stored
36 //! value (accessible using the rdt() method) is used.
37 //! @param count Set to zero to omit this call from the statistics
38 virtual void eval(double* x, double* r, double rdt=-1.0, int count=1) = 0;
39
40 //! Evaluates the Jacobian at x0 using finite differences.
41 //!
42 //! The Jacobian is computed by perturbing each component of `x0`, evaluating the
43 //! residual function, and then estimating the partial derivatives numerically using
44 //! finite differences to determine the corresponding column of the Jacobian.
45 //!
46 //! @param x0 State vector at which to evaluate the Jacobian
47 virtual void evalJacobian(double* x0) = 0;
48
49 //! Compute the weighted norm of `step`.
50 //!
51 //! The purpose is to measure the "size" of the step vector @f$ \Delta x @f$ in a
52 //! way that takes into account the relative importance or scale of different
53 //! solution components. Each component of the step vector is normalized by a weight
54 //! that depends on both the current magnitude of the solution vector and specified
55 //! tolerances. This makes the norm dimensionless and scaled appropriately, avoiding
56 //! issues where some components dominate due to differences in their scales. See
57 //! OneDim::weightedNorm() for a representative implementation.
58 virtual double weightedNorm(const double* step) const = 0;
59
60 //! Solve @f$ F(x) = 0 @f$, where @f$ F(x) @f$ is the residual function.
61 //!
62 //! @param[in] x0 Starting estimate of solution.
63 //! @param[out] x1 Final solution satisfying @f$ F(x1) = 0 @f$.
64 //! @param loglevel Controls amount of diagnostic output.
65 //!
66 //! @returns
67 //! - 1 for success
68 //! - -2 failure (maximum number of damping steps was reached)
69 //! - -3 failure (solution was up against the bounds
70 //!
71 //! @deprecated To be removed after Cantera 3.2. Use setInitialGuess(), solve(int),
72 //! and getState() instead.
73 int solve(double* x0, double* x1, int loglevel);
74
75 //! Set the initial guess. Should be called before solve().
76 void setInitialGuess(const double* x);
77
78 //! Solve the steady-state problem, taking internal timesteps as necessary until
79 //! the Newton solver can converge for the steady problem.
80 //! @param loglevel Controls amount of diagnostic output.
81 void solve(int loglevel=0);
82
83 //! Get the converged steady-state solution after calling solve().
84 void getState(double* x) const;
85
86 //! Steady-state max norm (infinity norm) of the residual evaluated using solution
87 //! x. On return, array r contains the steady-state residual values.
88 double ssnorm(double* x, double* r);
89
90 //! Total solution vector length;
91 size_t size() const {
92 return m_size;
93 }
94
95 //! Call to set the size of internal data structures after first defining the system
96 //! or if the problem size changes, for example after grid refinement.
97 //!
98 //! The #m_size variable should be updated before calling this method
99 virtual void resize();
100
101 //! Jacobian bandwidth.
102 size_t bandwidth() const {
103 return m_bw;
104 }
105
106 //! Get the name of the i-th component of the state vector
107 virtual string componentName(size_t i) const { return fmt::format("{}", i); }
108
109 //! Get header lines describing the column names included in a component label.
110 //! Headings should be aligned with items formatted in the output of
111 //! componentTableLabel(). Two lines are allowed. If only one is needed, the first
112 //! line should be blank.
113 virtual pair<string, string> componentTableHeader() const {
114 return {"", "component name"};
115 }
116
117 //! Get elements of the component name, aligned with the column headings given
118 //! by componentTableHeader().
119 virtual string componentTableLabel(size_t i) const { return componentName(i); }
120
121 //! Get the upper bound for global component *i* in the state vector
122 virtual double upperBound(size_t i) const = 0;
123
124 //! Get the lower bound for global component *i* in the state vector
125 virtual double lowerBound(size_t i) const = 0;
126
127 //! Return a reference to the Newton iterator.
129
130 //! Set the linear solver used to hold the Jacobian matrix and solve linear systems
131 //! as part of each Newton iteration.
132 void setLinearSolver(shared_ptr<SystemJacobian> solver);
133
134 //! Get the the linear solver being used to hold the Jacobian matrix and solve
135 //! linear systems as part of each Newton iteration.
136 shared_ptr<SystemJacobian> linearSolver() const { return m_jac; }
137
138 //! Reciprocal of the time step.
139 double rdt() const {
140 return m_rdt;
141 }
142
143 //! True if transient mode.
144 bool transient() const {
145 return (m_rdt != 0.0);
146 }
147
148 //! True if steady mode.
149 bool steady() const {
150 return (m_rdt == 0.0);
151 }
152
153 //! Prepare for time stepping beginning with solution *x* and timestep *dt*.
154 virtual void initTimeInteg(double dt, double* x);
155
156 //! Prepare to solve the steady-state problem. After invoking this method,
157 //! subsequent calls to solve() will solve the steady-state problem. Sets the
158 //! reciprocal of the time step to zero, and, if it was previously non-zero, signals
159 //! that a new Jacobian will be needed.
160 virtual void setSteadyMode();
161
162 //! Access the vector indicating which equations contain a transient term.
163 //! Elements are 1 for equations with a transient terms and 0 otherwise.
164 vector<int>& transientMask() {
165 return m_mask;
166 }
167
168 //! Take time steps using Backward Euler.
169 //!
170 //! @param nsteps number of steps
171 //! @param dt initial step size
172 //! @param x current solution vector
173 //! @param r solution vector after time stepping
174 //! @param loglevel controls amount of printed diagnostics
175 //! @returns size of last timestep taken
176 double timeStep(int nsteps, double dt, double* x, double* r, int loglevel);
177
178 //! Reset values such as negative species concentrations
179 virtual void resetBadValues(double* x) {}
180
181 //! @name Options
182 //! @{
183
184 //! Set the number of time steps to try when the steady Newton solver is
185 //! unsuccessful.
186 //! @param stepsize Initial time step size [s]
187 //! @param n Length of `tsteps` array
188 //! @param tsteps A sequence of time step counts to take after subsequent failures
189 //! of the steady-state solver. The last value in `tsteps` will be used again
190 //! after further unsuccessful solution attempts.
191 void setTimeStep(double stepsize, size_t n, const int* tsteps);
192
193 //! Set the minimum time step allowed during time stepping
194 void setMinTimeStep(double tmin) {
195 m_tmin = tmin;
196 }
197
198 //! Set the maximum time step allowed during time stepping
199 void setMaxTimeStep(double tmax) {
200 m_tmax = tmax;
201 }
202
203 //! Sets a factor by which the time step is reduced if the time stepping fails.
204 //! The default value is 0.5.
205 //!
206 //! @param tfactor factor time step is multiplied by if time stepping fails
207 void setTimeStepFactor(double tfactor) {
208 m_tfactor = tfactor;
209 }
210
211 //! Set the maximum number of timeteps allowed before successful steady-state solve
212 void setMaxTimeStepCount(int nmax) {
213 m_nsteps_max = nmax;
214 }
215
216 //! Get the maximum number of timeteps allowed before successful steady-state solve
217 int maxTimeStepCount() const {
218 return m_nsteps_max;
219 }
220 //! @}
221
222 //! Set the maximum number of steps that can be taken using the same Jacobian
223 //! before it must be re-evaluated.
224 //! @param ss_age Age limit during steady-state mode
225 //! @param ts_age Age limit during time stepping mode. If not specified, the
226 //! steady-state age is also used during time stepping.
227 void setJacAge(int ss_age, int ts_age=-1);
228
229 //! Set a function that will be called every time #eval is called.
230 //! Can be used to provide keyboard interrupt support in the high-level
231 //! language interfaces.
232 void setInterrupt(Func1* interrupt) {
233 m_interrupt = interrupt;
234 }
235
236 //! Set a function that will be called after each successful timestep. The
237 //! function will be called with the size of the timestep as the argument.
238 //! Intended to be used for observing solver progress for debugging
239 //! purposes.
240 void setTimeStepCallback(Func1* callback) {
241 m_time_step_callback = callback;
242 }
243
244 //! Configure perturbations used to evaluate finite difference Jacobian
245 //! @param relative Relative perturbation (multiplied by the absolute value of
246 //! each component). Default `1.0e-5`.
247 //! @param absolute Absolute perturbation (independent of component value).
248 //! Default `1.0e-10`.
249 //! @param threshold Threshold below which to exclude elements from the Jacobian
250 //! Default `0.0`.
251 void setJacobianPerturbation(double relative, double absolute, double threshold) {
252 m_jacobianRelPerturb = relative;
253 m_jacobianAbsPerturb = absolute;
254 m_jacobianThreshold = threshold;
255 }
256
257 //! Write solver debugging based on the specified log level.
258 //!
259 //! @see Sim1D::writeDebugInfo for a specific implementation of this capability.
260 virtual void writeDebugInfo(const string& header_suffix, const string& message,
261 int loglevel, int attempt_counter) {}
262
263 //! Delete debug output file that may be created by writeDebugInfo() when solving
264 //! with high `loglevel`.
265 virtual void clearDebugFile() {}
266
267protected:
268 //! Evaluate the steady-state Jacobian, accessible via linearSolver()
269 //! @param[in] x Current state vector, length size()
270 //! @param[out] rsd Storage for the residual, length size()
271 void evalSSJacobian(double* x, double* rsd);
272
273 //! Array of number of steps to take after each unsuccessful steady-state solve
274 //! before re-attempting the steady-state solution. For subsequent time stepping
275 //! calls, the final value is reused. See setTimeStep().
276 vector<int> m_steps = { 10 };
277
278 double m_tstep = 1.0e-5; //!< Initial timestep
279 double m_tmin = 1e-16; //!< Minimum timestep size
280 double m_tmax = 1e+08; //!< Maximum timestep size
281
282 //! Factor time step is multiplied by if time stepping fails ( < 1 )
283 double m_tfactor = 0.5;
284
285 shared_ptr<vector<double>> m_state; //!< Solution vector
286
287 //! Work array used to hold the residual or the new solution
288 vector<double> m_xnew;
289
290 //! State vector after the last successful set of time steps
291 vector<double> m_xlast_ts;
292
293 unique_ptr<MultiNewton> m_newt; //!< Newton iterator
294 double m_rdt = 0.0; //!< Reciprocal of time step
295
296 shared_ptr<SystemJacobian> m_jac; //!< Jacobian evaluator
297 bool m_jac_ok = false; //!< If `true`, Jacobian is current
298
299 size_t m_bw = 0; //!< Jacobian bandwidth
300 size_t m_size = 0; //!< %Solution vector size
301
302 //! Work arrays used during Jacobian evaluation
303 vector<double> m_work1, m_work2;
304
305 vector<int> m_mask; //!< Transient mask. See transientMask()
306 int m_ss_jac_age = 20; //!< Maximum age of the Jacobian in steady-state mode.
307 int m_ts_jac_age = 20; //!< Maximum age of the Jacobian in time-stepping mode.
308
309 //! Counter used to manage the number of states stored in the debug log file
310 //! generated by writeDebugInfo()
312
313 //! Constant that determines the maximum number of states stored in the debug log
314 //! file generated by writeDebugInfo()
316
317 //! Function called at the start of every call to #eval.
318 Func1* m_interrupt = nullptr;
319
320 //! User-supplied function called after each successful timestep.
322
323 int m_nsteps = 0; //!< Number of time steps taken in the current call to solve()
324 int m_nsteps_max = 500; //!< Maximum number of timesteps allowed per call to solve()
325
326 //! Threshold for ignoring small elements in Jacobian
328 //! Relative perturbation of each component in finite difference Jacobian
329 double m_jacobianRelPerturb = 1e-5;
330 //! Absolute perturbation of each component in finite difference Jacobian
331 double m_jacobianAbsPerturb = 1e-10;
332};
333
334}
335
336#endif
Declarations for class SystemJacobian.
Base class for 'functor' classes that evaluate a function of one variable.
Definition Func1.h:75
Newton iterator for multi-domain, one-dimensional problems.
Definition MultiNewton.h:25
Base class for representing a system of differential-algebraic equations and solving for its steady-s...
int solve(double *x0, double *x1, int loglevel)
Solve , where is the residual function.
int m_nsteps
Number of time steps taken in the current call to solve()
size_t m_size
Solution vector size
int m_nsteps_max
Maximum number of timesteps allowed per call to solve()
virtual void resize()
Call to set the size of internal data structures after first defining the system or if the problem si...
shared_ptr< SystemJacobian > linearSolver() const
Get the the linear solver being used to hold the Jacobian matrix and solve linear systems as part of ...
vector< double > m_xnew
Work array used to hold the residual or the new solution.
unique_ptr< MultiNewton > m_newt
Newton iterator.
virtual void resetBadValues(double *x)
Reset values such as negative species concentrations.
double m_jacobianAbsPerturb
Absolute perturbation of each component in finite difference Jacobian.
size_t size() const
Total solution vector length;.
virtual double weightedNorm(const double *step) const =0
Compute the weighted norm of step.
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
vector< int > & transientMask()
Access the vector indicating which equations contain a transient term.
virtual double upperBound(size_t i) const =0
Get the upper bound for global component i in the state vector.
double rdt() const
Reciprocal of the time step.
virtual void initTimeInteg(double dt, double *x)
Prepare for time stepping beginning with solution x and timestep dt.
virtual string componentName(size_t i) const
Get the name of the i-th component of the state vector.
virtual string componentTableLabel(size_t i) const
Get elements of the component name, aligned with the column headings given by componentTableHeader().
size_t bandwidth() const
Jacobian bandwidth.
double m_rdt
Reciprocal of time step.
double m_jacobianThreshold
Threshold for ignoring small elements in Jacobian.
void setMaxTimeStep(double tmax)
Set the maximum time step allowed during time stepping.
shared_ptr< SystemJacobian > m_jac
Jacobian evaluator.
shared_ptr< vector< double > > m_state
Solution vector.
void setMinTimeStep(double tmin)
Set the minimum time step allowed during time stepping.
vector< int > m_mask
Transient mask.
void setTimeStepCallback(Func1 *callback)
Set a function that will be called after each successful timestep.
Func1 * m_interrupt
Function called at the start of every call to eval.
virtual double lowerBound(size_t i) const =0
Get the lower bound for global component i in the state vector.
bool transient() const
True if transient mode.
void setMaxTimeStepCount(int nmax)
Set the maximum number of timeteps allowed before successful steady-state solve.
void setTimeStepFactor(double tfactor)
Sets a factor by which the time step is reduced if the time stepping fails.
vector< int > m_steps
Array of number of steps to take after each unsuccessful steady-state solve before re-attempting the ...
size_t m_bw
Jacobian bandwidth.
double m_tfactor
Factor time step is multiplied by if time stepping fails ( < 1 )
bool m_jac_ok
If true, Jacobian is current.
double m_tstep
Initial timestep.
int maxTimeStepCount() const
Get the maximum number of timeteps allowed before successful steady-state solve.
int m_ts_jac_age
Maximum age of the Jacobian in time-stepping mode.
virtual void eval(double *x, double *r, double rdt=-1.0, int count=1)=0
Evaluate the residual function.
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
void setJacAge(int ss_age, int ts_age=-1)
Set the maximum number of steps that can be taken using the same Jacobian before it must be re-evalua...
void evalSSJacobian(double *x, double *rsd)
Evaluate the steady-state Jacobian, accessible via linearSolver()
virtual void writeDebugInfo(const string &header_suffix, const string &message, int loglevel, int attempt_counter)
Write solver debugging based on the specified log level.
double m_tmin
Minimum timestep size.
void setJacobianPerturbation(double relative, double absolute, double threshold)
Configure perturbations used to evaluate finite difference Jacobian.
double m_tmax
Maximum timestep size.
int m_ss_jac_age
Maximum age of the Jacobian in steady-state mode.
virtual void setSteadyMode()
Prepare to solve the steady-state problem.
virtual void clearDebugFile()
Delete debug output file that may be created by writeDebugInfo() when solving with high loglevel.
int m_attempt_counter
Counter used to manage the number of states stored in the debug log file generated by writeDebugInfo(...
void setLinearSolver(shared_ptr< SystemJacobian > solver)
Set the linear solver used to hold the Jacobian matrix and solve linear systems as part of each Newto...
void setTimeStep(double stepsize, size_t n, const int *tsteps)
Set the number of time steps to try when the steady Newton solver is unsuccessful.
int m_max_history
Constant that determines the maximum number of states stored in the debug log file generated by write...
void setInitialGuess(const double *x)
Set the initial guess. Should be called before solve().
void getState(double *x) const
Get the converged steady-state solution after calling solve().
Func1 * m_time_step_callback
User-supplied function called after each successful timestep.
virtual pair< string, string > componentTableHeader() const
Get header lines describing the column names included in a component label.
virtual void evalJacobian(double *x0)=0
Evaluates the Jacobian at x0 using finite differences.
double m_jacobianRelPerturb
Relative perturbation of each component in finite difference Jacobian.
vector< double > m_xlast_ts
State vector after the last successful set of time steps.
bool steady() const
True if steady mode.
vector< double > m_work1
Work arrays used during Jacobian evaluation.
MultiNewton & newton()
Return a reference to the Newton iterator.
void setInterrupt(Func1 *interrupt)
Set a function that will be called every time eval is called.
This file contains definitions of constants, types and terms that are used in internal routines and a...
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595