Cantera  2.1.2
ReactorNet.cpp
Go to the documentation of this file.
1 //! @file ReactorNet.cpp
5 #include "cantera/zeroD/Wall.h"
6 
7 #include <cstdio>
8 
9 using namespace std;
10 
11 namespace Cantera
12 {
13 
14 ReactorNet::ReactorNet() : Cantera::FuncEval(), m_nr(0), m_nreactors(0),
15  m_integ(0), m_time(0.0), m_init(false),
16  m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4),
17  m_atols(1.0e-15), m_atolsens(1.0e-4),
18  m_maxstep(-1.0), m_maxErrTestFails(0),
19  m_verbose(false), m_ntotpar(0)
20 {
21 #ifdef DEBUG_MODE
22  m_verbose = true;
23 #endif
24  m_integ = newIntegrator("CVODE");// CVodeInt;
25 
26  // use backward differencing, with a full Jacobian computed
27  // numerically, and use a Newton linear iterator
28 
29  m_integ->setMethod(BDF_Method);
30  m_integ->setProblemType(DENSE + NOJAC);
31  m_integ->setIterator(Newton_Iter);
32 }
33 
34 ReactorNet::~ReactorNet()
35 {
36  for (size_t n = 0; n < m_nr; n++) {
37  if (m_iown[n]) {
38  delete m_r[n];
39  }
40  m_r[n] = 0;
41  }
42  m_r.clear();
43  m_reactors.clear();
44  delete m_integ;
45 }
46 
48 {
49  size_t n, nv;
50  char buf[100];
51  m_nv = 0;
52  m_reactors.clear();
53  m_nreactors = 0;
54  writelog("Initializing reactor network.\n", m_verbose);
55  if (m_nr == 0)
56  throw CanteraError("ReactorNet::initialize",
57  "no reactors in network!");
58  size_t sensParamNumber = 0;
59  for (n = 0; n < m_nr; n++) {
60  if (m_r[n]->type() >= ReactorType) {
61  m_r[n]->initialize(m_time);
62  Reactor* r = (Reactor*)m_r[n];
63  m_reactors.push_back(r);
64  nv = r->neq();
65  m_size.push_back(nv);
66  m_nparams.push_back(r->nSensParams());
67  std::vector<std::pair<void*, int> > sens_objs = r->getSensitivityOrder();
68  for (size_t i = 0; i < sens_objs.size(); i++) {
69  std::map<size_t, size_t>& s = m_sensOrder[sens_objs[i]];
70  for (std::map<size_t, size_t>::iterator iter = s.begin();
71  iter != s.end();
72  ++iter) {
73  m_sensIndex.resize(std::max(iter->second + 1, m_sensIndex.size()));
74  m_sensIndex[iter->second] = sensParamNumber++;
75  }
76  }
77  m_nv += nv;
78  m_nreactors++;
79 
80  if (m_verbose) {
81  sprintf(buf,"Reactor %s: %s variables.\n",
82  int2str(n).c_str(), int2str(nv).c_str());
83  writelog(buf);
84  sprintf(buf," %s sensitivity params.\n",
85  int2str(r->nSensParams()).c_str());
86  writelog(buf);
87  }
88  if (m_r[n]->type() == FlowReactorType && m_nr > 1) {
89  throw CanteraError("ReactorNet::initialize",
90  "FlowReactors must be used alone.");
91  }
92  }
93  }
94 
95  m_connect.resize(m_nr*m_nr,0);
96  m_ydot.resize(m_nv,0.0);
97  size_t i, j, nin, nout, nw;
98  ReactorBase* r, *rj;
99  for (i = 0; i < m_nr; i++) {
100  r = m_reactors[i];
101  for (j = 0; j < m_nr; j++) {
102  if (i == j) {
103  connect(i,j);
104  } else {
105  rj = m_reactors[j];
106  nin = rj->nInlets();
107  for (n = 0; n < nin; n++) {
108  if (&rj->inlet(n).out() == r) {
109  connect(i,j);
110  }
111  }
112  nout = rj->nOutlets();
113  for (n = 0; n < nout; n++) {
114  if (&rj->outlet(n).in() == r) {
115  connect(i,j);
116  }
117  }
118  nw = rj->nWalls();
119  for (n = 0; n < nw; n++) {
120  if (&rj->wall(n).left() == rj
121  && &rj->wall(n).right() == r) {
122  connect(i,j);
123  } else if (&rj->wall(n).left() == r
124  && &rj->wall(n).right() == rj) {
125  connect(i,j);
126  }
127  }
128  }
129  }
130  }
131 
132  m_atol.resize(neq());
133  fill(m_atol.begin(), m_atol.end(), m_atols);
134  m_integ->setTolerances(m_rtol, neq(), DATA_PTR(m_atol));
135  m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
136  m_integ->setMaxStepSize(m_maxstep);
137  m_integ->setMaxErrTestFails(m_maxErrTestFails);
138  if (m_verbose) {
139  sprintf(buf, "Number of equations: %s\n", int2str(neq()).c_str());
140  writelog(buf);
141  sprintf(buf, "Maximum time step: %14.6g\n", m_maxstep);
142  writelog(buf);
143  }
144  m_integ->initialize(m_time, *this);
145  m_init = true;
146 }
147 
148 void ReactorNet::advance(doublereal time)
149 {
150  if (!m_init) {
151  if (m_maxstep < 0.0) {
152  m_maxstep = time - m_time;
153  }
154  initialize();
155  }
156  m_integ->integrate(time);
157  m_time = time;
158  updateState(m_integ->solution());
159 }
160 
161 double ReactorNet::step(doublereal time)
162 {
163  if (!m_init) {
164  if (m_maxstep < 0.0) {
165  m_maxstep = time - m_time;
166  }
167  initialize();
168  }
169  m_time = m_integ->step(time);
170  updateState(m_integ->solution());
171  return m_time;
172 }
173 
175 {
176  r->setNetwork(this);
177  if (r->type() >= ReactorType) {
178  m_r.push_back(r);
179  m_iown.push_back(iown);
180  m_nr++;
181  writelog("Adding reactor "+r->name()+"\n", m_verbose);
182  } else {
183  writelog("Not adding reactor "+r->name()+
184  ", since type = "+int2str(r->type())+"\n", m_verbose);
185  }
186 }
187 
188 void ReactorNet::eval(doublereal t, doublereal* y,
189  doublereal* ydot, doublereal* p)
190 {
191  size_t n;
192  size_t start = 0;
193  size_t pstart = 0;
194 
195  updateState(y);
196  for (n = 0; n < m_nreactors; n++) {
197  m_reactors[n]->evalEqs(t, y + start,
198  ydot + start, p + pstart);
199  start += m_size[n];
200  pstart += m_nparams[n];
201  }
202 }
203 
204 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
205  doublereal* ydot, doublereal* p, Array2D* j)
206 {
207  doublereal ysave, dy;
208  Array2D& jac = *j;
209 
210  // use a try... catch block, since exceptions are not passed
211  // through CVODE, since it is C code
212  try {
213  //evaluate the unperturbed ydot
214  eval(t, y, ydot, p);
215  for (size_t n = 0; n < m_nv; n++) {
216 
217  // perturb x(n)
218  ysave = y[n];
219  dy = m_atol[n] + fabs(ysave)*m_rtol;
220  y[n] = ysave + dy;
221  dy = y[n] - ysave;
222 
223  // calculate perturbed residual
224  eval(t, y, DATA_PTR(m_ydot), p);
225 
226  // compute nth column of Jacobian
227  for (size_t m = 0; m < m_nv; m++) {
228  jac(m,n) = (m_ydot[m] - ydot[m])/dy;
229  }
230  y[n] = ysave;
231  }
232  } catch (CanteraError& err) {
233  std::cerr << err.what() << std::endl;
234  error("Terminating execution.");
235  }
236 }
237 
238 void ReactorNet::updateState(doublereal* y)
239 {
240  size_t start = 0;
241  for (size_t n = 0; n < m_nreactors; n++) {
242  m_reactors[n]->updateState(y + start);
243  start += m_size[n];
244  }
245 }
246 
248  size_t leny, doublereal* y)
249 {
250  size_t start = 0;
251  for (size_t n = 0; n < m_nreactors; n++) {
252  m_reactors[n]->getInitialConditions(t0, m_size[n], y + start);
253  start += m_size[n];
254  }
255 }
256 
257 size_t ReactorNet::globalComponentIndex(const string& species, size_t reactor)
258 {
259  if (!m_init) {
260  initialize();
261  }
262  size_t start = 0;
263  size_t n;
264  for (n = 0; n < reactor; n++) {
265  start += m_size[n];
266  }
267  return start + m_reactors[n]->componentIndex(species);
268 }
269 
271  size_t reactionIndex, const std::string& name, int leftright)
272 {
273  std::pair<void*, int> R = std::make_pair(reactor, leftright);
274  if (m_sensOrder.count(R) &&
275  m_sensOrder[R].count(reactionIndex)) {
276  throw CanteraError("ReactorNet::registerSensitivityReaction",
277  "Attempted to register duplicate sensitivity reaction");
278  }
279  m_paramNames.push_back(name);
280  m_sensOrder[R][reactionIndex] = m_ntotpar;
281  m_ntotpar++;
282 }
283 
284 }
Backward Differentiation.
Definition: Integrator.h:33
ReactorBase & reactor(int n)
Return a reference to the n-th reactor in this network.
Definition: ReactorNet.h:122
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:245
double step(doublereal time)
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:161
std::vector< std::string > m_paramNames
Names corresponding to each sensitivity parameter.
Definition: ReactorNet.h:236
const ReactorBase & right()
Return a reference to the Reactor or Reservoir to the right of the wall.
Definition: Wall.h:138
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
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:238
virtual size_t nSensParams()
Number of sensitivity parameters associated with this reactor (including walls)
Definition: Reactor.cpp:102
Header file for class Wall.
const ReactorBase & out() const
Return a const reference to the downstream reactor.
Definition: FlowDevice.h:82
virtual size_t neq()
Number of equations.
Definition: ReactorNet.h:174
size_t globalComponentIndex(const std::string &species, size_t reactor=0)
Return the index corresponding to the species named species in the reactor with index reactor in the ...
Definition: ReactorNet.cpp:257
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
Definition: Integrator.h:176
virtual void integrate(doublereal tout)
Integrate the system of equations.
Definition: Integrator.h:122
void initialize()
Initialize the reactor network.
Definition: ReactorNet.cpp:47
void evalJacobian(doublereal t, doublereal *y, doublereal *ydot, doublereal *p, Array2D *j)
Evaluate the Jacobian matrix for the reactor network.
Definition: ReactorNet.cpp:204
virtual void setTolerances(doublereal reltol, size_t n, doublereal *abstol)
Set or reset the number of equations.
Definition: Integrator.h:73
size_t nInlets()
Return the number of inlet FlowDevice objects connected to this reactor.
Definition: ReactorBase.h:86
virtual void initialize(doublereal t0, FuncEval &func)
Initialize the integrator for a new problem.
Definition: Integrator.h:109
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:33
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:247
const char * what() const
Get a description of the error.
void addReactor(ReactorBase *r, bool iown=false)
Add the reactor r to this reactor network.
Definition: ReactorNet.cpp:174
size_t nWalls()
Return the number of Wall objects connected to this reactor.
Definition: ReactorBase.h:97
Wall & wall(size_t n)
Return a reference to the n-th Wall connected to this reactor.
Definition: ReactorBase.cpp:67
void error(const std::string &msg)
Write an error message and quit.
Definition: global.cpp:71
ReactorBase & in() const
Return a reference to the upstream reactor.
Definition: FlowDevice.h:77
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
Base class for stirred reactors.
Definition: ReactorBase.h:30
virtual void eval(doublereal t, doublereal *y, doublereal *ydot, doublereal *p)
Evaluate the right-hand-side function.
Definition: ReactorNet.cpp:188
std::string name() const
Return the name of this reactor.
Definition: ReactorBase.h:42
FlowDevice & inlet(size_t n=0)
Return a reference to the n-th inlet FlowDevice connected to this reactor.
Definition: ReactorBase.cpp:97
void setNetwork(ReactorNet *net)
Set the ReactorNet that this reactor belongs to.
Definition: ReactorBase.cpp:82
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:270
virtual int type() const
Return a constant indicating the type of this Reactor.
Definition: ReactorBase.h:37
virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol)
Set the sensitivity error tolerances.
Definition: Integrator.h:92
virtual doublereal & solution(size_t k)
The current value of the solution of equation k.
Definition: Integrator.h:137
Newton Iteration.
Definition: Integrator.h:43
virtual doublereal step(doublereal tout)
Integrate the system of equations.
Definition: Integrator.h:131
doublereal time()
Current value of the simulation time.
Definition: ReactorNet.h:79
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
std::vector< std::pair< void *, int > > getSensitivityOrder() const
Return a vector specifying the ordering of objects to use when determining sensitivity parameter indi...
Definition: Reactor.cpp:338
size_t nOutlets()
Return the number of outlet FlowDevice objects connected to this reactor.
Definition: ReactorBase.h:92
virtual size_t neq()
Number of equations (state variables) for this reactor.
Definition: Reactor.h:93
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:240
virtual void setMaxErrTestFails(int n)
Set the maximum permissible number of error test failures.
Definition: Integrator.h:186
void writelog(const std::string &msg)
Write a message to the screen.
Definition: global.cpp:43
FlowDevice & outlet(size_t n=0)
Return a reference to the n-th outlet FlowDevice connected to this reactor.
ReactorBase & left() const
Return a reference to the Reactor or Reservoir to the left of the wall.
Definition: Wall.h:132
void advance(doublereal time)
Advance the state of all reactors in time.
Definition: ReactorNet.cpp:148
Class Reactor is a general-purpose class for stirred reactors.
Definition: Reactor.h:39