Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ReactorNet.h
Go to the documentation of this file.
1 /**
2  * @file ReactorNet.h
3  */
4 // Copyright 2004 California Institute of Technology
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  */
24 {
25 public:
26  ReactorNet();
27  virtual ~ReactorNet();
28 
29  /** @name Methods to set up a simulation. */
30  //@{
31 
32  /**
33  * Set initial time. Default = 0.0 s. Restarts integration
34  * from this time using the current mixture state as the
35  * initial condition.
36  */
37  void setInitialTime(doublereal time) {
38  m_time = time;
39  m_integrator_init = false;
40  }
41 
42  //! Set the maximum time step.
43  void setMaxTimeStep(double maxstep) {
44  m_maxstep = maxstep;
45  m_init = false;
46  }
47 
48  //! Set the maximum number of error test failures permitted by the CVODES
49  //! integrator in a single time step.
50  void setMaxErrTestFails(int nmax) {
51  m_maxErrTestFails = nmax;
52  m_init = false;
53  }
54 
55  //! Set the relative and absolute tolerances for the integrator.
56  void setTolerances(doublereal rtol, doublereal atol) {
57  if (rtol >= 0.0) {
58  m_rtol = rtol;
59  }
60  if (atol >= 0.0) {
61  m_atols = atol;
62  }
63  m_init = false;
64  }
65 
66  //! Set the relative and absolute tolerances for integrating the
67  //! sensitivity equations.
68  void setSensitivityTolerances(doublereal rtol, doublereal atol) {
69  if (rtol >= 0.0) {
70  m_rtolsens = rtol;
71  }
72  if (atol >= 0.0) {
73  m_atolsens = atol;
74  }
75  m_init = false;
76  }
77 
78  //! Current value of the simulation time.
79  doublereal time() {
80  return m_time;
81  }
82 
83  //! Relative tolerance.
84  doublereal rtol() {
85  return m_rtol;
86  }
87 
88  //! Absolute integration tolerance
89  doublereal atol() {
90  return m_atols;
91  }
92 
93  //! Relative sensitivity tolerance
94  doublereal rtolSensitivity() const {
95  return m_rtolsens;
96  }
97 
98  //! Absolute sensitivity tolerance
99  doublereal atolSensitivity() const {
100  return m_atolsens;
101  }
102 
103  /**
104  * Advance the state of all reactors in time. Take as many internal
105  * timesteps as necessary to reach *time*.
106  * @param time Time to advance to (s).
107  */
108  void advance(doublereal time);
109 
110  //! Advance the state of all reactors in time. Take a single timestep
111  //! toward *time*.
112  double step(doublereal time);
113 
114  //@}
115 
116  //! Add the reactor *r* to this reactor network.
117  void addReactor(Reactor& r);
118 
119  //! Add the reactor *r* to this reactor network.
120  /**
121  * @deprecated To be removed after Cantera 2.2. Use addReactor(Reactor&)
122  * instead.
123  */
124  void addReactor(Reactor* r, bool iown = false);
125 
126  //! Return a reference to the *n*-th reactor in this network. The reactor
127  //! indices are determined by the order in which the reactors were added
128  //! to the reactor network.
129  Reactor& reactor(int n) {
130  return *m_reactors[n];
131  }
132 
133  //! Returns `true` if verbose logging output is enabled.
134  bool verbose() const {
135  return m_verbose;
136  }
137 
138  //! Enable or disable verbose logging while setting up and integrating the
139  //! reactor network.
140  void setVerbose(bool v = true) {
141  m_verbose = v;
142  }
143 
144  //! Return a reference to the integrator.
146  return *m_integ;
147  }
148 
149  //! Update the state of all the reactors in the network to correspond to
150  //! the values in the solution vector *y*.
151  void updateState(doublereal* y);
152 
153  //! Return the sensitivity of the *k*-th solution component with respect to
154  //! the *p*-th sensitivity parameter.
155  /*!
156  * The normalized sensitivity coefficient \f$ S_{ki} \f$ of solution
157  * variable \f$ y_k \f$ with respect to sensitivity parameter \f$ p_i \f$
158  * is defined as:
159  *
160  * \f[ S_{ki} = \frac{p_i}{y_k} \frac{\partial y_k}{\partial p_i} \f]
161  *
162  * For reaction sensitivities, the parameter is a multiplier on the forward
163  * rate constant (and implicitly on the reverse rate constant for
164  * reversible reactions).
165  */
166  double sensitivity(size_t k, size_t p) {
167  if (!m_init) {
168  initialize();
169  }
170  if (p >= m_sensIndex.size()) {
171  throw IndexError("ReactorNet::sensitivity",
172  "m_sensIndex", p, m_sensIndex.size()-1);
173  }
174  return m_integ->sensitivity(k, m_sensIndex[p])/m_integ->solution(k);
175  }
176 
177  //! Return the sensitivity of the component named *component* with respect to
178  //! the *p*-th sensitivity parameter.
179  //! @copydetails ReactorNet::sensitivity(size_t, size_t)
180  double sensitivity(const std::string& component, size_t p, int reactor=0) {
181  size_t k = globalComponentIndex(component, reactor);
182  return sensitivity(k, p);
183  }
184 
185  //! Evaluate the Jacobian matrix for the reactor network.
186  /*!
187  * @param[in] t Time at which to evaluate the Jacobian
188  * @param[in] y Global state vector at time *t*
189  * @param[out] ydot Time derivative of the state vector evaluated at *t*.
190  * @param[in] p sensitivity parameter vector (unused?)
191  * @param[out] j Jacobian matrix, size neq() by neq().
192  */
193  void evalJacobian(doublereal t, doublereal* y,
194  doublereal* ydot, doublereal* p, Array2D* j);
195 
196  // overloaded methods of class FuncEval
197  virtual size_t neq() {
198  return m_nv;
199  }
200  virtual void eval(doublereal t, doublereal* y,
201  doublereal* ydot, doublereal* p);
202  virtual void getInitialConditions(doublereal t0, size_t leny,
203  doublereal* y);
204  virtual size_t nparams() {
205  return m_ntotpar;
206  }
207 
208  //! Return the index corresponding to the component named *component* in the
209  //! reactor with index *reactor* in the global state vector for the
210  //! reactor network.
211  size_t globalComponentIndex(const std::string& component, size_t reactor=0);
212 
213  //! Used by Reactor and Wall objects to register the addition of
214  //! sensitivity reactions so that the ReactorNet can keep track of the
215  //! order in which sensitivity parameters are added.
216  void registerSensitivityReaction(void* reactor, size_t reactionIndex,
217  const std::string& name, int leftright=0);
218 
219  //! The name of the p-th sensitivity parameter added to this ReactorNet.
220  const std::string& sensitivityParameterName(size_t p) {
221  return m_paramNames.at(p);
222  }
223 
224  //! Reinitialize the integrator. Used to solve a new problem (different
225  //! initial conditions) but with the same configuration of the reactor
226  //! network. Can be called manually, or automatically after calling
227  //! setInitialTime or modifying a reactor's contents.
228  void reinitialize();
229 
230  //! Called to trigger integrator reinitialization before further
231  //! integration.
232  void setNeedsReinit() {
233  m_integrator_init = false;
234  }
235 
236 protected:
237  /**
238  * Initialize the reactor network. Called automatically the first time
239  * advance or step is called.
240  */
241  void initialize();
242 
243  std::vector<Reactor*> m_reactors;
244  Integrator* m_integ;
245  doublereal m_time;
246  bool m_init;
247  bool m_integrator_init; //! True if integrator initialization is current
248  size_t m_nv;
249 
250  //! m_start[n] is the starting point in the state vector for reactor n
251  std::vector<size_t> m_start;
252 
253  vector_fp m_atol;
254  doublereal m_rtol, m_rtolsens;
255  doublereal m_atols, m_atolsens;
256  doublereal m_maxstep;
257  int m_maxErrTestFails;
258  bool m_verbose;
259  size_t m_ntotpar;
260  std::vector<size_t> m_nparams;
261 
262  //! Names corresponding to each sensitivity parameter
263  std::vector<std::string> m_paramNames;
264 
265  //! Structure used to determine the order of sensitivity parameters
266  //! m_sensOrder[Reactor or Wall, leftright][reaction number] = parameter index
267  std::map<std::pair<void*, int>, std::map<size_t, size_t> > m_sensOrder;
268 
269  //! Mapping from the order in which sensitivity parameters were added to
270  //! the ReactorNet to the order in which they occur in the integrator
271  //! output.
272  std::vector<size_t> m_sensIndex;
273 
274  vector_fp m_ydot;
275 
276  std::vector<bool> m_iown;
277 };
278 }
279 
280 #endif
std::vector< size_t > m_sensIndex
Mapping from the order in which sensitivity parameters were added to the ReactorNet to the order in w...
Definition: ReactorNet.h:272
double step(doublereal time)
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:128
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:226
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:162
std::vector< std::string > m_paramNames
Names corresponding to each sensitivity parameter.
Definition: ReactorNet.h:263
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
std::vector< size_t > m_start
m_start[n] is the starting point in the state vector for reactor n
Definition: ReactorNet.h:251
void setTolerances(doublereal rtol, doublereal atol)
Set the relative and absolute tolerances for the integrator.
Definition: ReactorNet.h:56
doublereal atolSensitivity() const
Absolute sensitivity tolerance.
Definition: ReactorNet.h:99
doublereal rtolSensitivity() const
Relative sensitivity tolerance.
Definition: ReactorNet.h:94
virtual size_t neq()
Number of equations.
Definition: ReactorNet.h:197
void reinitialize()
Reinitialize the integrator.
Definition: ReactorNet.cpp:102
void initialize()
Initialize the reactor network.
Definition: ReactorNet.cpp:41
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:145
void setNeedsReinit()
Called to trigger integrator reinitialization before further integration.
Definition: ReactorNet.h:232
size_t m_nv
True if integrator initialization is current.
Definition: ReactorNet.h:248
void setMaxTimeStep(double maxstep)
Set the maximum time step.
Definition: ReactorNet.h:43
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:29
virtual void getInitialConditions(doublereal t0, size_t leny, doublereal *y)
Fill the solution vector with the initial conditions at initial time t0.
Definition: ReactorNet.cpp:217
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:180
void setVerbose(bool v=true)
Enable or disable verbose logging while setting up and integrating the reactor network.
Definition: ReactorNet.h:140
doublereal atol()
Absolute integration tolerance.
Definition: ReactorNet.h:89
void setMaxErrTestFails(int nmax)
Set the maximum number of error test failures permitted by the CVODES integrator in a single time ste...
Definition: ReactorNet.h:50
Abstract base class for ODE system integrators.
Definition: Integrator.h:52
virtual void eval(doublereal t, doublereal *y, doublereal *ydot, doublereal *p)
Evaluate the right-hand-side function.
Definition: ReactorNet.cpp:169
Reactor & reactor(int n)
Return a reference to the n-th reactor in this network.
Definition: ReactorNet.h:129
void setSensitivityTolerances(doublereal rtol, doublereal atol)
Set the relative and absolute tolerances for integrating the sensitivity equations.
Definition: ReactorNet.h:68
virtual size_t nparams()
Number of sensitivity parameters.
Definition: ReactorNet.h:204
void registerSensitivityReaction(void *reactor, size_t reactionIndex, const std::string &name, int leftright=0)
Used by Reactor and Wall objects to register the addition of sensitivity reactions so that the Reacto...
Definition: ReactorNet.cpp:234
doublereal rtol()
Relative tolerance.
Definition: ReactorNet.h:84
bool verbose() const
Returns true if verbose logging output is enabled.
Definition: ReactorNet.h:134
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
virtual doublereal & solution(size_t k)
The current value of the solution of equation k.
Definition: Integrator.h:136
doublereal time()
Current value of the simulation time.
Definition: ReactorNet.h:79
std::map< std::pair< void *, int >, std::map< size_t, size_t > > m_sensOrder
Structure used to determine the order of sensitivity parameters m_sensOrder[Reactor or Wall...
Definition: ReactorNet.h:267
void setInitialTime(doublereal time)
Set initial time.
Definition: ReactorNet.h:37
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:23
An array index is out of range.
Definition: ctexceptions.h:184
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.h:166
const std::string & sensitivityParameterName(size_t p)
The name of the p-th sensitivity parameter added to this ReactorNet.
Definition: ReactorNet.h:220
void advance(doublereal time)
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:113
Class Reactor is a general-purpose class for stirred reactors.
Definition: Reactor.h:39