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