Cantera  2.4.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
4 // at http://www.cantera.org/license.txt for license and copyright information.
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;
136  updateState(m_integ->solution());
137 }
138 
139 double ReactorNet::step()
140 {
141  if (!m_init) {
142  initialize();
143  } else if (!m_integrator_init) {
144  reinitialize();
145  }
146  m_time = m_integ->step(m_time + 1.0);
147  updateState(m_integ->solution());
148  return m_time;
149 }
150 
151 void ReactorNet::addReactor(Reactor& r)
152 {
153  r.setNetwork(this);
154  m_reactors.push_back(&r);
155 }
156 
157 void ReactorNet::eval(doublereal t, doublereal* y,
158  doublereal* ydot, doublereal* p)
159 {
160  updateState(y);
161  for (size_t n = 0; n < m_reactors.size(); n++) {
162  m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p);
163  }
164  checkFinite("ydot", ydot, m_nv);
165 }
166 
167 double ReactorNet::sensitivity(size_t k, size_t p)
168 {
169  if (!m_init) {
170  initialize();
171  }
172  if (p >= m_sens_params.size()) {
173  throw IndexError("ReactorNet::sensitivity",
174  "m_sens_params", p, m_sens_params.size()-1);
175  }
176  double denom = m_integ->solution(k);
177  if (denom == 0.0) {
178  denom = SmallNumber;
179  }
180  return m_integ->sensitivity(k, p) / denom;
181 }
182 
183 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
184  doublereal* ydot, doublereal* p, Array2D* j)
185 {
186  //evaluate the unperturbed ydot
187  eval(t, y, ydot, p);
188  for (size_t n = 0; n < m_nv; n++) {
189  // perturb x(n)
190  double ysave = y[n];
191  double dy = m_atol[n] + fabs(ysave)*m_rtol;
192  y[n] = ysave + dy;
193  dy = y[n] - ysave;
194 
195  // calculate perturbed residual
196  eval(t, y, m_ydot.data(), p);
197 
198  // compute nth column of Jacobian
199  for (size_t m = 0; m < m_nv; m++) {
200  j->value(m,n) = (m_ydot[m] - ydot[m])/dy;
201  }
202  y[n] = ysave;
203  }
204 }
205 
206 void ReactorNet::updateState(doublereal* y)
207 {
208  checkFinite("y", y, m_nv);
209  for (size_t n = 0; n < m_reactors.size(); n++) {
210  m_reactors[n]->updateState(y + m_start[n]);
211  }
212 }
213 
214 void ReactorNet::getState(double* y)
215 {
216  for (size_t n = 0; n < m_reactors.size(); n++) {
217  m_reactors[n]->getState(y + m_start[n]);
218  }
219 }
220 
221 size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
222 {
223  if (!m_init) {
224  initialize();
225  }
226  return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
227 }
228 
229 std::string ReactorNet::componentName(size_t i) const
230 {
231  for (auto r : m_reactors) {
232  if (i < r->neq()) {
233  return r->name() + ": " + r->componentName(i);
234  } else {
235  i -= r->neq();
236  }
237  }
238  throw CanteraError("ReactorNet::componentName", "Index out of bounds");
239 }
240 
241 size_t ReactorNet::registerSensitivityParameter(
242  const std::string& name, double value, double scale)
243 {
244  if (m_integrator_init) {
245  throw CanteraError("ReactorNet::registerSensitivityParameter",
246  "Sensitivity parameters cannot be added after the"
247  "integrator has been initialized.");
248  }
249  m_paramNames.push_back(name);
250  m_sens_params.push_back(value);
251  m_paramScales.push_back(scale);
252  return m_sens_params.size() - 1;
253 }
254 
255 }
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:108
Header file for class Wall.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
STL namespace.
virtual void initialize(doublereal t0=0.0)
Initialize the reactor.
Definition: Reactor.cpp:75
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:231
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:135
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:83
An array index is out of range.
Definition: ctexceptions.h:158
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
Class Reactor is a general-purpose class for stirred reactors.
Definition: Reactor.h:37