Cantera  3.1.0a1
Loading...
Searching...
No Matches
ReactorNet.h
Go to the documentation of this file.
1//! @file ReactorNet.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_REACTORNET_H
7#define CT_REACTORNET_H
8
9#include "Reactor.h"
11
12
13namespace Cantera
14{
15
16class Array2D;
17class Integrator;
18class PreconditionerBase;
19
20//! A class representing a network of connected reactors.
21/*!
22 * This class is used to integrate the governing equations for a network of reactors
23 * that are time dependent (Reactor, ConstPressureReactor) connected by various
24 * means, for example Wall, MassFlowController, Valve, or PressureController; or
25 * reactors dependent on a single spatial variable (FlowReactor).
26 *
27 * @ingroup zerodGroup
28 */
29class ReactorNet : public FuncEval
30{
31public:
32 ReactorNet();
33 ~ReactorNet() override;
34 ReactorNet(const ReactorNet&) = delete;
35 ReactorNet& operator=(const ReactorNet&) = delete;
36
37 //! @name Methods to set up a simulation
38 //! @{
39
40 //! Set the type of linear solver used in the integration.
41 //! @param linSolverType type of linear solver. Default type: "DENSE"
42 //! Other options include: "DIAG", "DENSE", "GMRES", "BAND"
43 void setLinearSolverType(const string& linSolverType="DENSE");
44
45 //! Set preconditioner used by the linear solver
46 //! @param preconditioner preconditioner object used for the linear solver
47 void setPreconditioner(shared_ptr<PreconditionerBase> preconditioner);
48
49 //! Set the initial value of the independent variable (typically time).
50 //! Default = 0.0 s. Restarts integration from this value using the current mixture
51 //! state as the initial condition.
52 void setInitialTime(double time);
53
54 //! Get the initial value of the independent variable (typically time).
55 /*!
56 * @since New in %Cantera 3.0.
57 */
58 double getInitialTime() const {
59 return m_initial_time;
60 }
61
62 //! Get the maximum integrator step.
63 double maxTimeStep() const {
64 return m_maxstep;
65 }
66
67 //! Set the maximum integrator step.
68 void setMaxTimeStep(double maxstep);
69
70 //! Set the maximum number of error test failures permitted by the CVODES
71 //! integrator in a single step.
72 void setMaxErrTestFails(int nmax);
73
74 //! Set the relative and absolute tolerances for the integrator.
75 void setTolerances(double rtol, double atol);
76
77 //! Set the relative and absolute tolerances for integrating the
78 //! sensitivity equations.
79 void setSensitivityTolerances(double rtol, double atol);
80
81 //! Current value of the simulation time [s], for reactor networks that are solved
82 //! in the time domain.
83 double time();
84
85 //! Current position [m] along the length of the reactor network, for reactors that
86 //! are solved as a function of space.
87 double distance();
88
89 //! Relative tolerance.
90 double rtol() {
91 return m_rtol;
92 }
93
94 //! Absolute integration tolerance
95 double atol() {
96 return m_atols;
97 }
98
99 //! Relative sensitivity tolerance
100 double rtolSensitivity() const {
101 return m_rtolsens;
102 }
103
104 //! Absolute sensitivity tolerance
105 double atolSensitivity() const {
106 return m_atolsens;
107 }
108
109 //! Problem type of integrator
110 string linearSolverType() const;
111
112 //! Returns the maximum number of internal integration steps the
113 //! integrator will take before reaching the next output point
114 int maxSteps();
115
116 /**
117 * Advance the state of all reactors in the independent variable (time or space).
118 * Take as many internal steps as necessary to reach *t*.
119 * @param t Time/distance to advance to (s or m).
120 */
121 void advance(double t);
122
123 /**
124 * Advance the state of all reactors in the independent variable (time or space).
125 * Take as many internal steps as necessary towards *t*. If *applylimit* is true,
126 * the advance step will be automatically reduced if needed to stay within limits
127 * (set by setAdvanceLimit).
128 * Returns the time/distance at the end of integration.
129 * @param t Time/distance to advance to (s or m).
130 * @param applylimit Limit advance step (boolean).
131 */
132 double advance(double t, bool applylimit);
133
134 //! Advance the state of all reactors with respect to the independent variable
135 //! (time or space). Returns the new value of the independent variable [s or m].
136 double step();
137
138 //! Add the reactor *r* to this reactor network.
139 void addReactor(Reactor& r);
140
141 //! Return a reference to the *n*-th reactor in this network. The reactor
142 //! indices are determined by the order in which the reactors were added
143 //! to the reactor network.
144 Reactor& reactor(int n) {
145 return *m_reactors[n];
146 }
147
148 //! Returns `true` if verbose logging output is enabled.
149 bool verbose() const {
150 return m_verbose;
151 }
152
153 //! Enable or disable verbose logging while setting up and integrating the
154 //! reactor network.
155 void setVerbose(bool v = true) {
156 m_verbose = v;
157 suppressErrors(!m_verbose);
158 }
159
160 //! Return a reference to the integrator. Only valid after adding at least one
161 //! reactor to the network.
163
164 //! Update the state of all the reactors in the network to correspond to
165 //! the values in the solution vector *y*.
166 void updateState(double* y);
167
168 //! Return the sensitivity of the *k*-th solution component with respect to
169 //! the *p*-th sensitivity parameter.
170 /*!
171 * The sensitivity coefficient @f$ S_{ki} @f$ of solution variable @f$ y_k
172 * @f$ with respect to sensitivity parameter @f$ p_i @f$ is defined as:
173 *
174 * @f[ S_{ki} = \frac{1}{y_k} \frac{\partial y_k}{\partial p_i} @f]
175 *
176 * For reaction sensitivities, the parameter is a multiplier on the forward
177 * rate constant (and implicitly on the reverse rate constant for
178 * reversible reactions) which has a nominal value of 1.0, and the
179 * sensitivity is nondimensional.
180 *
181 * For species enthalpy sensitivities, the parameter is a perturbation to
182 * the molar enthalpy of formation, such that the dimensions of the
183 * sensitivity are kmol/J.
184 */
185 double sensitivity(size_t k, size_t p);
186
187 //! Return the sensitivity of the component named *component* with respect to
188 //! the *p*-th sensitivity parameter.
189 //! @copydetails ReactorNet::sensitivity(size_t, size_t)
190 double sensitivity(const string& component, size_t p, int reactor=0) {
191 size_t k = globalComponentIndex(component, reactor);
192 return sensitivity(k, p);
193 }
194
195 //! Evaluate the Jacobian matrix for the reactor network.
196 /*!
197 * @param[in] t Time/distance at which to evaluate the Jacobian
198 * @param[in] y Global state vector at *t*
199 * @param[out] ydot Derivative of the state vector evaluated at *t*, with respect
200 * to *t*.
201 * @param[in] p sensitivity parameter vector (unused?)
202 * @param[out] j Jacobian matrix, size neq() by neq().
203 */
204 void evalJacobian(double t, double* y,
205 double* ydot, double* p, Array2D* j);
206
207 // overloaded methods of class FuncEval
208 size_t neq() const override {
209 return m_nv;
210 }
211
212 size_t nReactors() const {
213 return m_reactors.size();
214 }
215
216 void eval(double t, double* y, double* ydot, double* p) override;
217
218 //! eval coupling for IDA / DAEs
219 void evalDae(double t, double* y, double* ydot, double* p,
220 double* residual) override;
221
222 void getState(double* y) override;
223 void getStateDae(double* y, double* ydot) override;
224
225 //! Return k-th derivative at the current state of the system
226 virtual void getDerivative(int k, double* dky);
227
228 void getConstraints(double* constraints) override;
229
230 size_t nparams() const override {
231 return m_sens_params.size();
232 }
233
234 //! Return the index corresponding to the component named *component* in the
235 //! reactor with index *reactor* in the global state vector for the
236 //! reactor network.
237 size_t globalComponentIndex(const string& component, size_t reactor=0);
238
239 //! Return the name of the i-th component of the global state vector. The
240 //! name returned includes both the name of the reactor and the specific
241 //! component, for example `'reactor1: CH4'`.
242 string componentName(size_t i) const;
243
244 //! Used by Reactor and Wall objects to register the addition of
245 //! sensitivity parameters so that the ReactorNet can keep track of the
246 //! order in which sensitivity parameters are added.
247 //! @param name A name describing the parameter, for example the reaction string
248 //! @param value The nominal value of the parameter
249 //! @param scale A scaling factor to be applied to the sensitivity
250 //! coefficient
251 //! @returns the index of this parameter in the vector of sensitivity
252 //! parameters (global across all reactors)
253 size_t registerSensitivityParameter(const string& name, double value, double scale);
254
255 //! The name of the p-th sensitivity parameter added to this ReactorNet.
256 const string& sensitivityParameterName(size_t p) const {
257 return m_paramNames.at(p);
258 }
259
260 //! Initialize the reactor network. Called automatically the first time
261 //! advance or step is called.
262 void initialize();
263
264 //! Reinitialize the integrator. Used to solve a new problem (different
265 //! initial conditions) but with the same configuration of the reactor
266 //! network. Can be called manually, or automatically after calling
267 //! setInitialTime or modifying a reactor's contents.
268 void reinitialize();
269
270 //! Called to trigger integrator reinitialization before further
271 //! integration.
273 m_integrator_init = false;
274 }
275
276 //! Set the maximum number of internal integration steps the
277 //! integrator will take before reaching the next output point
278 //! @param nmax The maximum number of steps, setting this value
279 //! to zero disables this option.
280 virtual void setMaxSteps(int nmax);
281
282 //! Set absolute step size limits during advance
283 void setAdvanceLimits(const double* limits);
284
285 //! Check whether ReactorNet object uses advance limits
286 bool hasAdvanceLimits() const;
287
288 //! Retrieve absolute step size limits during advance
289 bool getAdvanceLimits(double* limits) const;
290
291 void preconditionerSetup(double t, double* y, double gamma) override;
292
293 void preconditionerSolve(double* rhs, double* output) override;
294
295 //! Get solver stats from integrator
296 AnyMap solverStats() const;
297
298 //! Set derivative settings of all reactors
299 //! @param settings the settings map propagated to all reactors and kinetics objects
300 virtual void setDerivativeSettings(AnyMap& settings);
301
302protected:
303 //! Check that preconditioning is supported by all reactors in the network
304 virtual void checkPreconditionerSupported() const;
305
306 void updatePreconditioner(double gamma) override;
307
308 //! Estimate a future state based on current derivatives.
309 //! The function is intended for internal use by ReactorNet::advance
310 //! and deliberately not exposed in external interfaces.
311 virtual void getEstimate(double time, int k, double* yest);
312
313 //! Returns the order used for last solution step of the ODE integrator
314 //! The function is intended for internal use by ReactorNet::advance
315 //! and deliberately not exposed in external interfaces.
316 virtual int lastOrder() const;
317
318 vector<Reactor*> m_reactors;
319 unique_ptr<Integrator> m_integ;
320
321 //! The independent variable in the system. May be either time or space depending
322 //! on the type of reactors in the network.
323 double m_time = 0.0;
324
325 //! The initial value of the independent variable in the system.
326 double m_initial_time = 0.0;
327
328 bool m_init = false;
329 bool m_integrator_init = false; //!< True if integrator initialization is current
330 size_t m_nv = 0;
331
332 //! m_start[n] is the starting point in the state vector for reactor n
333 vector<size_t> m_start;
334
335 vector<double> m_atol;
336 double m_rtol = 1.0e-9;
337 double m_rtolsens = 1.0e-4;
338 double m_atols = 1.0e-15;
339 double m_atolsens = 1.0e-6;
340 shared_ptr<PreconditionerBase> m_precon;
341 string m_linearSolverType;
342
343 //! Maximum integrator internal timestep. Default of 0.0 means infinity.
344 double m_maxstep = 0.0;
345
346 bool m_verbose = false;
347
348 //! Indicates whether time or space is the independent variable
350
351 //! Names corresponding to each sensitivity parameter
352 vector<string> m_paramNames;
353
354 vector<double> m_ydot;
355 vector<double> m_yest;
356 vector<double> m_advancelimits;
357 //! m_LHS is a vector representing the coefficients on the
358 //! "left hand side" of each governing equation
359 vector<double> m_LHS;
360 vector<double> m_RHS;
361};
362}
363
364#endif
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
Virtual base class for ODE/DAE right-hand-side function evaluators.
Definition FuncEval.h:32
bool suppressErrors() const
Get current state of error suppression.
Definition FuncEval.h:166
vector< double > m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
Definition FuncEval.h:181
Abstract base class for ODE system integrators.
Definition Integrator.h:44
A class representing a network of connected reactors.
Definition ReactorNet.h:30
void setPreconditioner(shared_ptr< PreconditionerBase > preconditioner)
Set preconditioner used by the linear solver.
void setLinearSolverType(const string &linSolverType="DENSE")
Set the type of linear solver used in the integration.
void preconditionerSetup(double t, double *y, double gamma) override
Evaluate the setup processes for the Jacobian preconditioner.
double step()
Advance the state of all reactors with respect to the independent variable (time or space).
virtual int lastOrder() const
Returns the order used for last solution step of the ODE integrator The function is intended for inte...
const string & sensitivityParameterName(size_t p) const
The name of the p-th sensitivity parameter added to this ReactorNet.
Definition ReactorNet.h:256
void eval(double t, double *y, double *ydot, double *p) override
Evaluate the right-hand-side ODE function.
size_t nparams() const override
Number of sensitivity parameters.
Definition ReactorNet.h:230
void initialize()
Initialize the reactor network.
void advance(double t)
Advance the state of all reactors in the independent variable (time or space).
size_t neq() const override
Number of equations.
Definition ReactorNet.h:208
vector< size_t > m_start
m_start[n] is the starting point in the state vector for reactor n
Definition ReactorNet.h:333
vector< double > m_LHS
m_LHS is a vector representing the coefficients on the "left hand side" of each governing equation
Definition ReactorNet.h:359
double m_initial_time
The initial value of the independent variable in the system.
Definition ReactorNet.h:326
void evalJacobian(double t, double *y, double *ydot, double *p, Array2D *j)
Evaluate the Jacobian matrix for the reactor network.
double time()
Current value of the simulation time [s], for reactor networks that are solved in the time domain.
double getInitialTime() const
Get the initial value of the independent variable (typically time).
Definition ReactorNet.h:58
void setNeedsReinit()
Called to trigger integrator reinitialization before further integration.
Definition ReactorNet.h:272
void getConstraints(double *constraints) override
Given a vector of length neq(), mark which variables should be considered algebraic constraints.
double m_time
The independent variable in the system.
Definition ReactorNet.h:323
AnyMap solverStats() const
Get solver stats from integrator.
Reactor & reactor(int n)
Return a reference to the n-th reactor in this network.
Definition ReactorNet.h:144
virtual void setMaxSteps(int nmax)
Set the maximum number of internal integration steps the integrator will take before reaching the nex...
string componentName(size_t i) const
Return the name of the i-th component of the global state vector.
void addReactor(Reactor &r)
Add the reactor r to this reactor network.
void getStateDae(double *y, double *ydot) override
Fill in the vectors y and ydot with the current state of the system.
void setInitialTime(double time)
Set the initial value of the independent variable (typically time).
virtual void getDerivative(int k, double *dky)
Return k-th derivative at the current state of the system.
void setMaxErrTestFails(int nmax)
Set the maximum number of error test failures permitted by the CVODES integrator in a single step.
size_t registerSensitivityParameter(const string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
double maxTimeStep() const
Get the maximum integrator step.
Definition ReactorNet.h:63
double m_maxstep
Maximum integrator internal timestep. Default of 0.0 means infinity.
Definition ReactorNet.h:344
double distance()
Current position [m] along the length of the reactor network, for reactors that are solved as a funct...
double atolSensitivity() const
Absolute sensitivity tolerance.
Definition ReactorNet.h:105
void setSensitivityTolerances(double rtol, double atol)
Set the relative and absolute tolerances for integrating the sensitivity equations.
int maxSteps()
Returns the maximum number of internal integration steps the integrator will take before reaching the...
virtual void setDerivativeSettings(AnyMap &settings)
Set derivative settings of all reactors.
double sensitivity(size_t k, size_t p)
Return the sensitivity of the k-th solution component with respect to the p-th sensitivity parameter.
void updateState(double *y)
Update the state of all the reactors in the network to correspond to the values in the solution vecto...
void getState(double *y) override
Fill in the vector y with the current state of the system.
void setAdvanceLimits(const double *limits)
Set absolute step size limits during advance.
double rtol()
Relative tolerance.
Definition ReactorNet.h:90
bool verbose() const
Returns true if verbose logging output is enabled.
Definition ReactorNet.h:149
size_t globalComponentIndex(const string &component, size_t reactor=0)
Return the index corresponding to the component named component in the reactor with index reactor in ...
bool m_timeIsIndependent
Indicates whether time or space is the independent variable.
Definition ReactorNet.h:349
double atol()
Absolute integration tolerance.
Definition ReactorNet.h:95
bool hasAdvanceLimits() const
Check whether ReactorNet object uses advance limits.
double sensitivity(const string &component, size_t p, int reactor=0)
Return the sensitivity of the component named component with respect to the p-th sensitivity paramete...
Definition ReactorNet.h:190
void setMaxTimeStep(double maxstep)
Set the maximum integrator step.
double rtolSensitivity() const
Relative sensitivity tolerance.
Definition ReactorNet.h:100
void evalDae(double t, double *y, double *ydot, double *p, double *residual) override
eval coupling for IDA / DAEs
virtual void checkPreconditionerSupported() const
Check that preconditioning is supported by all reactors in the network.
bool m_integrator_init
True if integrator initialization is current.
Definition ReactorNet.h:329
void reinitialize()
Reinitialize the integrator.
Integrator & integrator()
Return a reference to the integrator.
bool getAdvanceLimits(double *limits) const
Retrieve absolute step size limits during advance.
void setVerbose(bool v=true)
Enable or disable verbose logging while setting up and integrating the reactor network.
Definition ReactorNet.h:155
string linearSolverType() const
Problem type of integrator.
void updatePreconditioner(double gamma) override
Update the preconditioner based on already computed jacobian values.
void preconditionerSolve(double *rhs, double *output) override
Evaluate the linear system Ax=b where A is the preconditioner.
vector< string > m_paramNames
Names corresponding to each sensitivity parameter.
Definition ReactorNet.h:352
void setTolerances(double rtol, double atol)
Set the relative and absolute tolerances for the integrator.
virtual void getEstimate(double time, int k, double *yest)
Estimate a future state based on current derivatives.
Class Reactor is a general-purpose class for stirred reactors.
Definition Reactor.h:47
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition utilities.h:104
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564