Cantera  2.4.0
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 http://www.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"
12 #include "cantera/base/Array.h"
13 
14 namespace Cantera
15 {
16 
17 //! A class representing a network of connected reactors.
18 /*!
19  * This class is used to integrate the time-dependent governing equations for
20  * a network of reactors (Reactor, ConstPressureReactor) connected by various
21  * means, e.g. Wall, MassFlowController, Valve, PressureController.
22  */
23 class ReactorNet : public FuncEval
24 {
25 public:
26  ReactorNet();
27  virtual ~ReactorNet() {};
28  ReactorNet(const ReactorNet&) = delete;
29  ReactorNet& operator=(const ReactorNet&) = delete;
30 
31  //! @name Methods to set up a simulation.
32  //@{
33 
34  //! Set initial time. Default = 0.0 s. Restarts integration from this time
35  //! using the current mixture state as the initial condition.
36  void setInitialTime(double time);
37 
38  //! Set the maximum time step.
39  void setMaxTimeStep(double maxstep);
40 
41  //! Set the maximum number of error test failures permitted by the CVODES
42  //! integrator in a single time step.
43  void setMaxErrTestFails(int nmax);
44 
45  //! Set the relative and absolute tolerances for the integrator.
46  void setTolerances(double rtol, double atol);
47 
48  //! Set the relative and absolute tolerances for integrating the
49  //! sensitivity equations.
50  void setSensitivityTolerances(double rtol, double atol);
51 
52  //! Current value of the simulation time.
53  doublereal time() {
54  return m_time;
55  }
56 
57  //! Relative tolerance.
58  doublereal rtol() {
59  return m_rtol;
60  }
61 
62  //! Absolute integration tolerance
63  doublereal atol() {
64  return m_atols;
65  }
66 
67  //! Relative sensitivity tolerance
68  doublereal rtolSensitivity() const {
69  return m_rtolsens;
70  }
71 
72  //! Absolute sensitivity tolerance
73  doublereal atolSensitivity() const {
74  return m_atolsens;
75  }
76 
77  /**
78  * Advance the state of all reactors in time. Take as many internal
79  * timesteps as necessary to reach *time*.
80  * @param time Time to advance to (s).
81  */
82  void advance(doublereal time);
83 
84  //! Advance the state of all reactors in time.
85  double step();
86 
87  //@}
88 
89  //! Add the reactor *r* to this reactor network.
90  void addReactor(Reactor& r);
91 
92  //! Return a reference to the *n*-th reactor in this network. The reactor
93  //! indices are determined by the order in which the reactors were added
94  //! to the reactor network.
95  Reactor& reactor(int n) {
96  return *m_reactors[n];
97  }
98 
99  //! Returns `true` if verbose logging output is enabled.
100  bool verbose() const {
101  return m_verbose;
102  }
103 
104  //! Enable or disable verbose logging while setting up and integrating the
105  //! reactor network.
106  void setVerbose(bool v = true) {
107  m_verbose = v;
108  suppressErrors(!m_verbose);
109  }
110 
111  //! Return a reference to the integrator.
113  return *m_integ;
114  }
115 
116  //! Update the state of all the reactors in the network to correspond to
117  //! the values in the solution vector *y*.
118  void updateState(doublereal* y);
119 
120  //! Return the sensitivity of the *k*-th solution component with respect to
121  //! the *p*-th sensitivity parameter.
122  /*!
123  * The sensitivity coefficient \f$ S_{ki} \f$ of solution variable \f$ y_k
124  * \f$ with respect to sensitivity parameter \f$ p_i \f$ is defined as:
125  *
126  * \f[ S_{ki} = \frac{1}{y_k} \frac{\partial y_k}{\partial p_i} \f]
127  *
128  * For reaction sensitivities, the parameter is a multiplier on the forward
129  * rate constant (and implicitly on the reverse rate constant for
130  * reversible reactions) which has a nominal value of 1.0, and the
131  * sensitivity is nondimensional.
132  *
133  * For species enthalpy sensitivities, the parameter is a perturbation to
134  * the molar enthalpy of formation, such that the dimensions of the
135  * sensitivity are kmol/J.
136  */
137  double sensitivity(size_t k, size_t p);
138 
139  //! Return the sensitivity of the component named *component* with respect to
140  //! the *p*-th sensitivity parameter.
141  //! @copydetails ReactorNet::sensitivity(size_t, size_t)
142  double sensitivity(const std::string& component, size_t p, int reactor=0) {
143  size_t k = globalComponentIndex(component, reactor);
144  return sensitivity(k, p);
145  }
146 
147  //! Evaluate the Jacobian matrix for the reactor network.
148  /*!
149  * @param[in] t Time at which to evaluate the Jacobian
150  * @param[in] y Global state vector at time *t*
151  * @param[out] ydot Time derivative of the state vector evaluated at *t*.
152  * @param[in] p sensitivity parameter vector (unused?)
153  * @param[out] j Jacobian matrix, size neq() by neq().
154  */
155  void evalJacobian(doublereal t, doublereal* y,
156  doublereal* ydot, doublereal* p, Array2D* j);
157 
158  // overloaded methods of class FuncEval
159  virtual size_t neq() {
160  return m_nv;
161  }
162  virtual void eval(doublereal t, doublereal* y,
163  doublereal* ydot, doublereal* p);
164 
165  virtual void getState(doublereal* y);
166 
167  virtual size_t nparams() {
168  return m_sens_params.size();
169  }
170 
171  //! Return the index corresponding to the component named *component* in the
172  //! reactor with index *reactor* in the global state vector for the
173  //! reactor network.
174  size_t globalComponentIndex(const std::string& component, size_t reactor=0);
175 
176  //! Return the name of the i-th component of the global state vector. The
177  //! name returned includes both the name of the reactor and the specific
178  //! component, e.g. `'reactor1: CH4'`.
179  std::string componentName(size_t i) const;
180 
181  //! Used by Reactor and Wall objects to register the addition of
182  //! sensitivity parameters so that the ReactorNet can keep track of the
183  //! order in which sensitivity parameters are added.
184  //! @param name A name describing the parameter, e.g. the reaction string
185  //! @param value The nominal value of the parameter
186  //! @param scale A scaling factor to be applied to the sensitivity
187  //! coefficient
188  //! @returns the index of this parameter in the vector of sensitivity
189  //! parameters (global across all reactors)
190  size_t registerSensitivityParameter(const std::string& name, double value,
191  double scale);
192 
193  //! The name of the p-th sensitivity parameter added to this ReactorNet.
194  const std::string& sensitivityParameterName(size_t p) {
195  return m_paramNames.at(p);
196  }
197 
198  //! Reinitialize the integrator. Used to solve a new problem (different
199  //! initial conditions) but with the same configuration of the reactor
200  //! network. Can be called manually, or automatically after calling
201  //! setInitialTime or modifying a reactor's contents.
202  void reinitialize();
203 
204  //! Called to trigger integrator reinitialization before further
205  //! integration.
206  void setNeedsReinit() {
207  m_integrator_init = false;
208  }
209 
210 protected:
211  //! Initialize the reactor network. Called automatically the first time
212  //! advance or step is called.
213  void initialize();
214 
215  std::vector<Reactor*> m_reactors;
216  std::unique_ptr<Integrator> m_integ;
217  doublereal m_time;
218  bool m_init;
219  bool m_integrator_init; //!< True if integrator initialization is current
220  size_t m_nv;
221 
222  //! m_start[n] is the starting point in the state vector for reactor n
223  std::vector<size_t> m_start;
224 
225  vector_fp m_atol;
226  doublereal m_rtol, m_rtolsens;
227  doublereal m_atols, m_atolsens;
228 
229  //! Maximum integrator internal timestep. Default of 0.0 means infinity.
230  doublereal m_maxstep;
231 
232  int m_maxErrTestFails;
233  bool m_verbose;
234 
235  //! Names corresponding to each sensitivity parameter
236  std::vector<std::string> m_paramNames;
237 
238  vector_fp m_ydot;
239 };
240 }
241 
242 #endif
size_t globalComponentIndex(const std::string &component, size_t reactor=0)
Return the index corresponding to the component named component in the reactor with index reactor in ...
Definition: ReactorNet.cpp:221
A class representing a network of connected reactors.
Definition: ReactorNet.h:23
void addReactor(Reactor &r)
Add the reactor r to this reactor network.
Definition: ReactorNet.cpp:151
std::vector< std::string > m_paramNames
Names corresponding to each sensitivity parameter.
Definition: ReactorNet.h:236
void setSensitivityTolerances(double rtol, double atol)
Set the relative and absolute tolerances for integrating the sensitivity equations.
Definition: ReactorNet.cpp:63
void updateState(doublereal *y)
Update the state of all the reactors in the network to correspond to the values in the solution vecto...
Definition: ReactorNet.cpp:206
std::vector< size_t > m_start
m_start[n] is the starting point in the state vector for reactor n
Definition: ReactorNet.h:223
doublereal rtolSensitivity() const
Relative sensitivity tolerance.
Definition: ReactorNet.h:68
bool verbose() const
Returns true if verbose logging output is enabled.
Definition: ReactorNet.h:100
virtual size_t neq()
Number of equations.
Definition: ReactorNet.h:159
void reinitialize()
Reinitialize the integrator.
Definition: ReactorNet.cpp:116
void initialize()
Initialize the reactor network.
Definition: ReactorNet.cpp:74
void evalJacobian(doublereal t, doublereal *y, doublereal *ydot, doublereal *p, Array2D *j)
Evaluate the Jacobian matrix for the reactor network.
Definition: ReactorNet.cpp:183
Integrator & integrator()
Return a reference to the integrator.
Definition: ReactorNet.h:112
void setNeedsReinit()
Called to trigger integrator reinitialization before further integration.
Definition: ReactorNet.h:206
double step()
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:139
void setMaxTimeStep(double maxstep)
Set the maximum time step.
Definition: ReactorNet.cpp:40
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
size_t registerSensitivityParameter(const std::string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
Definition: ReactorNet.cpp:241
Header file for class Cantera::Array2D.
double sensitivity(const std::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:142
void setVerbose(bool v=true)
Enable or disable verbose logging while setting up and integrating the reactor network.
Definition: ReactorNet.h:106
doublereal atol()
Absolute integration tolerance.
Definition: ReactorNet.h:63
std::string componentName(size_t i) const
Return the name of the i-th component of the global state vector.
Definition: ReactorNet.cpp:229
void setMaxErrTestFails(int nmax)
Set the maximum number of error test failures permitted by the CVODES integrator in a single time ste...
Definition: ReactorNet.cpp:46
Abstract base class for ODE system integrators.
Definition: Integrator.h:52
bool m_integrator_init
True if integrator initialization is current.
Definition: ReactorNet.h:219
virtual void eval(doublereal t, doublereal *y, doublereal *ydot, doublereal *p)
Evaluate the right-hand-side function.
Definition: ReactorNet.cpp:157
vector_fp m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
Definition: FuncEval.h:81
Reactor & reactor(int n)
Return a reference to the n-th reactor in this network.
Definition: ReactorNet.h:95
virtual void getState(doublereal *y)
Fill in the vector y with the current state of the system.
Definition: ReactorNet.cpp:214
doublereal m_maxstep
Maximum integrator internal timestep. Default of 0.0 means infinity.
Definition: ReactorNet.h:230
bool suppressErrors() const
Get current state of error suppression.
Definition: FuncEval.h:71
virtual size_t nparams()
Number of sensitivity parameters.
Definition: ReactorNet.h:167
doublereal rtol()
Relative tolerance.
Definition: ReactorNet.h:58
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 scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
doublereal time()
Current value of the simulation time.
Definition: ReactorNet.h:53
void setTolerances(double rtol, double atol)
Set the relative and absolute tolerances for the integrator.
Definition: ReactorNet.cpp:52
void setInitialTime(double time)
Set initial time.
Definition: ReactorNet.cpp:34
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:26
doublereal atolSensitivity() const
Absolute sensitivity tolerance.
Definition: ReactorNet.h:73
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...
Definition: ReactorNet.cpp:167
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
const std::string & sensitivityParameterName(size_t p)
The name of the p-th sensitivity parameter added to this ReactorNet.
Definition: ReactorNet.h:194
void advance(doublereal time)
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:127
Class Reactor is a general-purpose class for stirred reactors.
Definition: Reactor.h:37