Cantera  2.3.0
ReactorNet.cpp
Go to the documentation of this file.
1 //! @file ReactorNet.cpp
2
3 // This file is part of Cantera. See License.txt in the top-level directory or
5
8 #include "cantera/zeroD/Wall.h"
9
10 #include <cstdio>
11
12 using namespace std;
13
14 namespace Cantera
15 {
16
17 ReactorNet::ReactorNet() :
18  m_integ(newIntegrator("CVODE")),
19  m_time(0.0), m_init(false), m_integrator_init(false),
20  m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4),
21  m_atols(1.0e-15), m_atolsens(1.0e-6),
22  m_maxstep(0.0), m_maxErrTestFails(0),
23  m_verbose(false)
24 {
25  suppressErrors(true);
26
27  // use backward differencing, with a full Jacobian computed
28  // numerically, and use a Newton linear iterator
29  m_integ->setMethod(BDF_Method);
30  m_integ->setProblemType(DENSE + NOJAC);
31  m_integ->setIterator(Newton_Iter);
32 }
33
34 void ReactorNet::setInitialTime(double time)
35 {
36  m_time = time;
37  m_integrator_init = false;
38 }
39
40 void ReactorNet::setMaxTimeStep(double maxstep)
41 {
42  m_maxstep = maxstep;
43  m_init = false;
44 }
45
46 void ReactorNet::setMaxErrTestFails(int nmax)
47 {
48  m_maxErrTestFails = nmax;
49  m_init = false;
50 }
51
52 void ReactorNet::setTolerances(double rtol, double atol)
53 {
54  if (rtol >= 0.0) {
55  m_rtol = rtol;
56  }
57  if (atol >= 0.0) {
58  m_atols = atol;
59  }
60  m_init = false;
61 }
62
63 void ReactorNet::setSensitivityTolerances(double rtol, double atol)
64 {
65  if (rtol >= 0.0) {
66  m_rtolsens = rtol;
67  }
68  if (atol >= 0.0) {
69  m_atolsens = atol;
70  }
71  m_init = false;
72 }
73
74 void ReactorNet::initialize()
75 {
76  m_nv = 0;
77  debuglog("Initializing reactor network.\n", m_verbose);
78  if (m_reactors.empty()) {
79  throw CanteraError("ReactorNet::initialize",
80  "no reactors in network!");
81  }
82  m_start.assign(1, 0);
83  for (size_t n = 0; n < m_reactors.size(); n++) {
84  Reactor& r = *m_reactors[n];
85  r.initialize(m_time);
86  size_t nv = r.neq();
87  m_nv += nv;
88  m_start.push_back(m_nv);
89
90  if (m_verbose) {
91  writelog("Reactor {:d}: {:d} variables.\n", n, nv);
92  writelog(" {:d} sensitivity params.\n", r.nSensParams());
93  }
94  if (r.type() == FlowReactorType && m_reactors.size() > 1) {
95  throw CanteraError("ReactorNet::initialize",
96  "FlowReactors must be used alone.");
97  }
98  }
99
100  m_ydot.resize(m_nv,0.0);
101  m_atol.resize(neq());
102  fill(m_atol.begin(), m_atol.end(), m_atols);
103  m_integ->setTolerances(m_rtol, neq(), m_atol.data());
104  m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
105  m_integ->setMaxStepSize(m_maxstep);
106  m_integ->setMaxErrTestFails(m_maxErrTestFails);
107  if (m_verbose) {
108  writelog("Number of equations: {:d}\n", neq());
109  writelog("Maximum time step: {:14.6g}\n", m_maxstep);
110  }
111  m_integ->initialize(m_time, *this);
112  m_integrator_init = true;
113  m_init = true;
114 }
115
116 void ReactorNet::reinitialize()
117 {
118  if (m_init) {
119  debuglog("Re-initializing reactor network.\n", m_verbose);
120  m_integ->reinitialize(m_time, *this);
121  m_integrator_init = true;
122  } else {
123  initialize();
124  }
125 }
126
127 void ReactorNet::advance(doublereal time)
128 {
129  if (!m_init) {
130  initialize();
131  } else if (!m_integrator_init) {
132  reinitialize();
133  }
134  m_integ->integrate(time);
135  m_time = time;
137 }
138
139 double ReactorNet::step(doublereal time)
140 {
141  if (time != -999) {
142  warn_deprecated("ReactorNet::step(t)", "The argument to this function"
143  " is deprecated and will be removed after Cantera 2.3.");
144  }
145  if (!m_init) {
146  initialize();
147  } else if (!m_integrator_init) {
148  reinitialize();
149  }
150  m_time = m_integ->step(m_time + 1.0);
152  return m_time;
153 }
154
155 void ReactorNet::addReactor(Reactor& r)
156 {
157  r.setNetwork(this);
158  m_reactors.push_back(&r);
159 }
160
161 void ReactorNet::eval(doublereal t, doublereal* y,
162  doublereal* ydot, doublereal* p)
163 {
165  for (size_t n = 0; n < m_reactors.size(); n++) {
166  m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p);
167  }
168  checkFinite("ydot", ydot, m_nv);
169 }
170
171 double ReactorNet::sensitivity(size_t k, size_t p)
172 {
173  if (!m_init) {
174  initialize();
175  }
176  if (p >= m_sens_params.size()) {
177  throw IndexError("ReactorNet::sensitivity",
178  "m_sens_params", p, m_sens_params.size()-1);
179  }
180  double denom = m_integ->solution(k);
181  if (denom == 0.0) {
182  denom = SmallNumber;
183  }
184  return m_integ->sensitivity(k, p) / denom;
185 }
186
187 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
188  doublereal* ydot, doublereal* p, Array2D* j)
189 {
190  //evaluate the unperturbed ydot
191  eval(t, y, ydot, p);
192  for (size_t n = 0; n < m_nv; n++) {
193  // perturb x(n)
194  double ysave = y[n];
195  double dy = m_atol[n] + fabs(ysave)*m_rtol;
196  y[n] = ysave + dy;
197  dy = y[n] - ysave;
198
199  // calculate perturbed residual
200  eval(t, y, m_ydot.data(), p);
201
202  // compute nth column of Jacobian
203  for (size_t m = 0; m < m_nv; m++) {
204  j->value(m,n) = (m_ydot[m] - ydot[m])/dy;
205  }
206  y[n] = ysave;
207  }
208 }
209
210 void ReactorNet::updateState(doublereal* y)
211 {
212  checkFinite("y", y, m_nv);
213  for (size_t n = 0; n < m_reactors.size(); n++) {
214  m_reactors[n]->updateState(y + m_start[n]);
215  }
216 }
217
218 void ReactorNet::getInitialConditions(double t0, size_t leny, double* y)
219 {
220  warn_deprecated("ReactorNet::getInitialConditions",
221  "Use getState instead. To be removed after Cantera 2.3.");
222  getState(y);
223 }
224
225 void ReactorNet::getState(double* y)
226 {
227  for (size_t n = 0; n < m_reactors.size(); n++) {
228  m_reactors[n]->getState(y + m_start[n]);
229  }
230 }
231
232 size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
233 {
234  if (!m_init) {
235  initialize();
236  }
237  return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
238 }
239
240 std::string ReactorNet::componentName(size_t i) const
241 {
242  for (auto r : m_reactors) {
243  if (i < r->neq()) {
244  return r->name() + ": " + r->componentName(i);
245  } else {
246  i -= r->neq();
247  }
248  }
249  throw CanteraError("ReactorNet::componentName", "Index out of bounds");
250 }
251
252 size_t ReactorNet::registerSensitivityParameter(
253  const std::string& name, double value, double scale)
254 {
255  if (m_integrator_init) {
256  throw CanteraError("ReactorNet::registerSensitivityParameter",
257  "Sensitivity parameters cannot be added after the"
258  "integrator has been initialized.");
259  }
260  m_paramNames.push_back(name);
261  m_sens_params.push_back(value);
262  m_paramScales.push_back(scale);
263  return m_sens_params.size() - 1;
264 }
265
266 }
Backward Differentiation.
Definition: Integrator.h:33
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
Definition: checkFinite.cpp:15
virtual size_t nSensParams()
Number of sensitivity parameters associated with this reactor (including walls)
Definition: Reactor.cpp:118
Header file for class Wall.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
STL namespace.
virtual void initialize(doublereal t0=0.0)
Initialize the reactor.
Definition: Reactor.cpp:82
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
doublereal & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition: Array.h:253
void setNetwork(ReactorNet *net)
Set the ReactorNet that this reactor belongs to.
Definition: ReactorBase.cpp:95
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:161
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
Newton Iteration.
Definition: Integrator.h:43
virtual int type() const
Return a constant indicating the type of this Reactor.
Definition: Reactor.h:42
virtual size_t neq()
Number of equations (state variables) for this reactor.
Definition: Reactor.h:99
An array index is out of range.
Definition: ctexceptions.h:164
Namespace for the Cantera kernel.
Definition: application.cpp:29
Class Reactor is a general-purpose class for stirred reactors.
Definition: Reactor.h:37