Cantera  2.0
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 class ReactorNet : public Cantera::FuncEval
18 {
19 
20 public:
21 
22  //! Constructor
23  ReactorNet();
24 
25  //! Destructor
26  virtual ~ReactorNet();
27 
28  //-----------------------------------------------------
29 
30  /** @name Methods to set up a simulation. */
31  //@{
32 
33  /**
34  * Set initial time. Default = 0.0 s. Restarts integration
35  * from this time using the current mixture state as the
36  * initial condition.
37  */
38  void setInitialTime(doublereal time) {
39  m_time = time;
40  m_init = false;
41  }
42 
43  /// Set the maximum time step.
44  void setMaxTimeStep(double maxstep) {
45  m_maxstep = maxstep;
46  m_init = false;
47  }
48 
49  void setTolerances(doublereal rtol, doublereal atol) {
50  if (rtol >= 0.0) {
51  m_rtol = rtol;
52  }
53  if (atol >= 0.0) {
54  m_atols = atol;
55  }
56  m_init = false;
57  }
58 
59  void setSensitivityTolerances(doublereal rtol, doublereal atol) {
60  if (rtol >= 0.0) {
61  m_rtolsens = rtol;
62  }
63  if (atol >= 0.0) {
64  m_atolsens = atol;
65  }
66  m_init = false;
67  }
68 
69  /// Current value of the simulation time.
70  doublereal time() {
71  return m_time;
72  }
73 
74  /// Relative tolerance.
75  doublereal rtol() {
76  return m_rtol;
77  }
78  doublereal atol() {
79  return m_atols;
80  }
81 
82  /**
83  * Initialize the reactor network.
84  */
85  void initialize(doublereal t0 = 0.0);
86 
87  /**
88  * Advance the state of all reactors in time.
89  * @param time Time to advance to (s).
90  */
91  void advance(doublereal time);
92 
93  double step(doublereal time);
94 
95  //@}
96 
97  void addReactor(ReactorBase* r, bool iown = false);
98 
99  ReactorBase& reactor(int n) {
100  return *m_r[n];
101  }
102 
103  bool verbose() const {
104  return m_verbose;
105  }
106  void setVerbose(bool v = true) {
107  m_verbose = v;
108  }
109 
110  /// Return a reference to the integrator.
111  Integrator& integrator() {
112  return *m_integ;
113  }
114 
115  void updateState(doublereal* y);
116 
117  double sensitivity(size_t k, size_t p) {
118  return m_integ->sensitivity(k, p)/m_integ->solution(k);
119  }
120 
121  double sensitivity(std::string species, size_t p, int reactor=0) {
122  size_t k = globalComponentIndex(species, reactor);
123  return sensitivity(k, p);
124  }
125 
126  void evalJacobian(doublereal t, doublereal* y,
127  doublereal* ydot, doublereal* p, Array2D* j);
128 
129  //-----------------------------------------------------
130 
131  // overloaded methods of class FuncEval
132  virtual size_t neq() {
133  return m_nv;
134  }
135  virtual void eval(doublereal t, doublereal* y,
136  doublereal* ydot, doublereal* p);
137  virtual void getInitialConditions(doublereal t0, size_t leny,
138  doublereal* y);
139  virtual size_t nparams() {
140  return m_ntotpar;
141  }
142 
143  size_t globalComponentIndex(std::string species, size_t reactor=0);
144 
145  void connect(size_t i, size_t j) {
146  m_connect[j*m_nr + i] = 1;
147  m_connect[i*m_nr + j] = 1;
148  }
149 
150  bool connected(size_t i, size_t j) {
151  return (m_connect[m_nr*i + j] == 1);
152  }
153 
154 protected:
155 
156  std::vector<ReactorBase*> m_r;
157  std::vector<Reactor*> m_reactors;
158  size_t m_nr;
159  size_t m_nreactors;
160  Integrator* m_integ;
161  doublereal m_time;
162  bool m_init;
163  size_t m_nv;
164  std::vector<size_t> m_size;
165  vector_fp m_atol;
166  doublereal m_rtol, m_rtolsens;
167  doublereal m_atols, m_atolsens;
168  doublereal m_maxstep;
169  bool m_verbose;
170  size_t m_ntotpar;
171  std::vector<size_t> m_nparams;
172  vector_int m_connect;
173  vector_fp m_ydot;
174 
175  std::vector<bool> m_iown;
176 
177 private:
178 
179 };
180 }
181 
182 #endif
183