Cantera  3.0.0
Loading...
Searching...
No Matches
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 https://cantera.org/license.txt for license and copyright information.
5
10#include "cantera/base/Array.h"
13
14#include <cstdio>
15
16namespace Cantera
17{
18
19ReactorNet::ReactorNet()
20{
21 suppressErrors(true);
22}
23
24ReactorNet::~ReactorNet()
25{
26}
27
29{
30 m_time = time;
32 m_integrator_init = false;
33}
34
35void ReactorNet::setMaxTimeStep(double maxstep)
36{
37 m_maxstep = maxstep;
39}
40
42{
44}
45
46void ReactorNet::setTolerances(double rtol, double atol)
47{
48 if (rtol >= 0.0) {
49 m_rtol = rtol;
50 }
51 if (atol >= 0.0) {
52 m_atols = atol;
53 }
54 m_init = false;
55}
56
57void ReactorNet::setSensitivityTolerances(double rtol, double atol)
58{
59 if (rtol >= 0.0) {
60 m_rtolsens = rtol;
61 }
62 if (atol >= 0.0) {
63 m_atolsens = atol;
64 }
65 m_init = false;
66}
67
70 return m_time;
71 } else {
72 throw CanteraError("ReactorNet::time", "Time is not the independent variable"
73 " for this reactor network.");
74 }
75}
76
79 return m_time;
80 } else {
81 throw CanteraError("ReactorNet::distance", "Distance is not the independent"
82 " variable for this reactor network.");
83 }
84}
85
87{
88 m_nv = 0;
89 debuglog("Initializing reactor network.\n", m_verbose);
90 if (m_reactors.empty()) {
91 throw CanteraError("ReactorNet::initialize",
92 "no reactors in network!");
93 }
94 m_start.assign(1, 0);
95 for (size_t n = 0; n < m_reactors.size(); n++) {
96 Reactor& r = *m_reactors[n];
98 size_t nv = r.neq();
99 m_nv += nv;
100 m_start.push_back(m_nv);
101
102 if (m_verbose) {
103 writelog("Reactor {:d}: {:d} variables.\n", n, nv);
104 writelog(" {:d} sensitivity params.\n", r.nSensParams());
105 }
106 if (r.type() == "FlowReactor" && m_reactors.size() > 1) {
107 throw CanteraError("ReactorNet::initialize",
108 "FlowReactors must be used alone.");
109 }
110 }
111
112 m_ydot.resize(m_nv,0.0);
113 m_yest.resize(m_nv,0.0);
114 m_advancelimits.resize(m_nv,-1.0);
115 m_atol.resize(neq());
116 fill(m_atol.begin(), m_atol.end(), m_atols);
117 m_integ->setTolerances(m_rtol, neq(), m_atol.data());
118 m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
119 if (!m_linearSolverType.empty()) {
120 m_integ->setLinearSolverType(m_linearSolverType);
121 }
122 if (m_precon) {
123 m_integ->setPreconditioner(m_precon);
124 }
125 m_integ->initialize(m_time, *this);
126 if (m_verbose) {
127 writelog("Number of equations: {:d}\n", neq());
128 writelog("Maximum time step: {:14.6g}\n", m_maxstep);
129 }
130 if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
132 }
133 m_integrator_init = true;
134 m_init = true;
135}
136
138{
139 if (m_init) {
140 debuglog("Re-initializing reactor network.\n", m_verbose);
141 m_integ->reinitialize(m_time, *this);
142 if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
144 }
145 m_integrator_init = true;
146 } else {
147 initialize();
148 }
149}
150
151void ReactorNet::setLinearSolverType(const string& linSolverType)
152{
153 m_linearSolverType = linSolverType;
154 m_integrator_init = false;
155}
156
157void ReactorNet::setPreconditioner(shared_ptr<PreconditionerBase> preconditioner)
158{
159 m_precon = preconditioner;
160 m_integrator_init = false;
161}
162
164{
165 integrator().setMaxSteps(nmax);
166}
167
169{
170 return integrator().maxSteps();
171}
172
173void ReactorNet::advance(double time)
174{
175 if (!m_init) {
176 initialize();
177 } else if (!m_integrator_init) {
178 reinitialize();
179 }
180 m_integ->integrate(time);
181 m_time = time;
182 updateState(m_integ->solution());
183}
184
185double ReactorNet::advance(double time, bool applylimit)
186{
187 if (!m_init) {
188 initialize();
189 } else if (!m_integrator_init) {
190 reinitialize();
191 }
192
193 if (!applylimit) {
194 // take full step
195 advance(time);
196 return time;
197 }
198
199 if (!hasAdvanceLimits()) {
200 // take full step
201 advance(time);
202 return time;
203 }
204
205 getAdvanceLimits(m_advancelimits.data());
206
207 // ensure that gradient is available
208 while (lastOrder() < 1) {
209 step();
210 }
211
212 int k = lastOrder();
213 double t = time, delta;
214 double* y = m_integ->solution();
215
216 // reduce time step if limits are exceeded
217 while (true) {
218 bool exceeded = false;
219 getEstimate(t, k, &m_yest[0]);
220 for (size_t j = 0; j < m_nv; j++) {
221 delta = abs(m_yest[j] - y[j]);
222 if ( (m_advancelimits[j] > 0.) && ( delta > m_advancelimits[j]) ) {
223 exceeded = true;
224 if (m_verbose) {
225 writelog(" Limiting global state vector component {:d} (dt = {:9.4g}):"
226 "{:11.6g} > {:9.4g}\n",
227 j, t - m_time, delta, m_advancelimits[j]);
228 }
229 }
230 }
231 if (!exceeded) {
232 break;
233 }
234 t = .5 * (m_time + t);
235 }
236 advance(t);
237 return t;
238}
239
241{
242 if (!m_init) {
243 initialize();
244 } else if (!m_integrator_init) {
245 reinitialize();
246 }
247 m_time = m_integ->step(m_time + 1.0);
248 updateState(m_integ->solution());
249 return m_time;
250}
251
252void ReactorNet::getEstimate(double time, int k, double* yest)
253{
254 if (!m_init) {
255 initialize();
256 }
257 // initialize
258 double* cvode_dky = m_integ->solution();
259 for (size_t j = 0; j < m_nv; j++) {
260 yest[j] = cvode_dky[j];
261 }
262
263 // Taylor expansion
264 double factor = 1.;
265 double deltat = time - m_time;
266 for (int n = 1; n <= k; n++) {
267 factor *= deltat / n;
268 cvode_dky = m_integ->derivative(m_time, n);
269 for (size_t j = 0; j < m_nv; j++) {
270 yest[j] += factor * cvode_dky[j];
271 }
272 }
273}
274
276{
277 if (m_integ) {
278 return m_integ->lastOrder();
279 } else {
280 return 0;
281 }
282}
283
285{
286 for (auto current : m_reactors) {
287 if (current->isOde() != r.isOde()) {
288 throw CanteraError("ReactorNet::addReactor",
289 "Cannot mix Reactor types using both ODEs and DAEs ({} and {})",
290 current->type(), r.type());
291 }
292 if (current->timeIsIndependent() != r.timeIsIndependent()) {
293 throw CanteraError("ReactorNet::addReactor",
294 "Cannot mix Reactor types using time and space as independent variables"
295 "\n({} and {})", current->type(), r.type());
296 }
297 }
299 r.setNetwork(this);
300 m_reactors.push_back(&r);
301 if (!m_integ) {
302 m_integ.reset(newIntegrator(r.isOde() ? "CVODE" : "IDA"));
303 // use backward differencing, with a full Jacobian computed
304 // numerically, and use a Newton linear iterator
305 m_integ->setMethod(BDF_Method);
306 m_integ->setLinearSolverType("DENSE");
307 }
308}
309
311 if (m_integ == nullptr) {
312 throw CanteraError("ReactorNet::integrator",
313 "Integrator has not been instantiated. Add one or more reactors first.");
314 }
315 return *m_integ;
316}
317
318void ReactorNet::eval(double t, double* y, double* ydot, double* p)
319{
320 m_time = t;
321 updateState(y);
322 m_LHS.assign(m_nv, 1);
323 m_RHS.assign(m_nv, 0);
324 for (size_t n = 0; n < m_reactors.size(); n++) {
325 m_reactors[n]->applySensitivity(p);
326 m_reactors[n]->eval(t, m_LHS.data() + m_start[n], m_RHS.data() + m_start[n]);
327 size_t yEnd = 0;
328 if (n == m_reactors.size() - 1) {
329 yEnd = m_RHS.size();
330 } else {
331 yEnd = m_start[n + 1];
332 }
333 for (size_t i = m_start[n]; i < yEnd; i++) {
334 ydot[i] = m_RHS[i] / m_LHS[i];
335 }
336 m_reactors[n]->resetSensitivity(p);
337 }
338 checkFinite("ydot", ydot, m_nv);
339}
340
341void ReactorNet::evalDae(double t, double* y, double* ydot, double* p, double* residual)
342{
343 m_time = t;
344 updateState(y);
345 for (size_t n = 0; n < m_reactors.size(); n++) {
346 m_reactors[n]->applySensitivity(p);
347 m_reactors[n]->evalDae(t, y, ydot, residual);
348 m_reactors[n]->resetSensitivity(p);
349 }
350 checkFinite("ydot", ydot, m_nv);
351}
352
353void ReactorNet::getConstraints(double* constraints)
354{
355 for (size_t n = 0; n < m_reactors.size(); n++) {
356 m_reactors[n]->getConstraints(constraints + m_start[n]);
357 }
358}
359
360double ReactorNet::sensitivity(size_t k, size_t p)
361{
362 if (!m_init) {
363 initialize();
364 }
365 if (p >= m_sens_params.size()) {
366 throw IndexError("ReactorNet::sensitivity",
367 "m_sens_params", p, m_sens_params.size()-1);
368 }
369 double denom = m_integ->solution(k);
370 if (denom == 0.0) {
371 denom = SmallNumber;
372 }
373 return m_integ->sensitivity(k, p) / denom;
374}
375
376void ReactorNet::evalJacobian(double t, double* y, double* ydot, double* p, Array2D* j)
377{
378 //evaluate the unperturbed ydot
379 eval(t, y, ydot, p);
380 for (size_t n = 0; n < m_nv; n++) {
381 // perturb x(n)
382 double ysave = y[n];
383 double dy = m_atol[n] + fabs(ysave)*m_rtol;
384 y[n] = ysave + dy;
385 dy = y[n] - ysave;
386
387 // calculate perturbed residual
388 eval(t, y, m_ydot.data(), p);
389
390 // compute nth column of Jacobian
391 for (size_t m = 0; m < m_nv; m++) {
392 j->value(m,n) = (m_ydot[m] - ydot[m])/dy;
393 }
394 y[n] = ysave;
395 }
396}
397
399{
400 checkFinite("y", y, m_nv);
401 for (size_t n = 0; n < m_reactors.size(); n++) {
402 m_reactors[n]->updateState(y + m_start[n]);
403 }
404}
405
406void ReactorNet::getDerivative(int k, double* dky)
407{
408 if (!m_init) {
409 initialize();
410 }
411 double* cvode_dky = m_integ->derivative(m_time, k);
412 for (size_t j = 0; j < m_nv; j++) {
413 dky[j] = cvode_dky[j];
414 }
415}
416
417void ReactorNet::setAdvanceLimits(const double *limits)
418{
419 if (!m_init) {
420 initialize();
421 }
422 for (size_t n = 0; n < m_reactors.size(); n++) {
423 m_reactors[n]->setAdvanceLimits(limits + m_start[n]);
424 }
425}
426
428{
429 bool has_limit = false;
430 for (size_t n = 0; n < m_reactors.size(); n++) {
431 has_limit |= m_reactors[n]->hasAdvanceLimits();
432 }
433 return has_limit;
434}
435
436bool ReactorNet::getAdvanceLimits(double *limits) const
437{
438 bool has_limit = false;
439 for (size_t n = 0; n < m_reactors.size(); n++) {
440 has_limit |= m_reactors[n]->getAdvanceLimits(limits + m_start[n]);
441 }
442 return has_limit;
443}
444
445void ReactorNet::getState(double* y)
446{
447 for (size_t n = 0; n < m_reactors.size(); n++) {
448 m_reactors[n]->getState(y + m_start[n]);
449 }
450}
451
452void ReactorNet::getStateDae(double* y, double* ydot)
453{
454 for (size_t n = 0; n < m_reactors.size(); n++) {
455 m_reactors[n]->getStateDae(y + m_start[n], ydot + m_start[n]);
456 }
457}
458
459size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
460{
461 if (!m_init) {
462 initialize();
463 }
464 return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
465}
466
467string ReactorNet::componentName(size_t i) const
468{
469 for (auto r : m_reactors) {
470 if (i < r->neq()) {
471 return r->name() + ": " + r->componentName(i);
472 } else {
473 i -= r->neq();
474 }
475 }
476 throw CanteraError("ReactorNet::componentName", "Index out of bounds");
477}
478
480 const string& name, double value, double scale)
481{
482 if (m_integrator_init) {
483 throw CanteraError("ReactorNet::registerSensitivityParameter",
484 "Sensitivity parameters cannot be added after the"
485 "integrator has been initialized.");
486 }
487 m_paramNames.push_back(name);
488 m_sens_params.push_back(value);
489 m_paramScales.push_back(scale);
490 return m_sens_params.size() - 1;
491}
492
494{
495 // Apply given settings to all reactors
496 for (size_t i = 0; i < m_reactors.size(); i++) {
497 m_reactors[i]->setDerivativeSettings(settings);
498 }
499}
500
502{
503 if (m_integ) {
504 return m_integ->solverStats();
505 } else {
506 return AnyMap();
507 }
508}
509
511{
512 if (m_integ) {
513 return m_integ->linearSolverType();
514 } else {
515 return "";
516 }
517}
518
519void ReactorNet::preconditionerSolve(double* rhs, double* output)
520{
521 if (!m_integ) {
522 throw CanteraError("ReactorNet::preconditionerSolve",
523 "Must only be called after ReactorNet is initialized.");
524 }
525 m_integ->preconditionerSolve(m_nv, rhs, output);
526}
527
528void ReactorNet::preconditionerSetup(double t, double* y, double gamma)
529{
530 // ensure state is up to date.
531 updateState(y);
532 // get the preconditioner
533 auto precon = m_integ->preconditioner();
534 // Reset preconditioner
535 precon->reset();
536 // Set gamma value for M =I - gamma*J
537 precon->setGamma(gamma);
538 // Make a copy of state to adjust it for preconditioner
539 vector<double> yCopy(m_nv);
540 // Get state of reactor
541 getState(yCopy.data());
542 // transform state based on preconditioner rules
543 precon->stateAdjustment(yCopy);
544 // update network with adjusted state
545 updateState(yCopy.data());
546 // Get jacobians and give elements to preconditioners
547 for (size_t i = 0; i < m_reactors.size(); i++) {
548 Eigen::SparseMatrix<double> rJac = m_reactors[i]->jacobian();
549 for (int k=0; k<rJac.outerSize(); ++k) {
550 for (Eigen::SparseMatrix<double>::InnerIterator it(rJac, k); it; ++it) {
551 precon->setValue(it.row() + m_start[i], it.col() + m_start[i],
552 it.value());
553 }
554 }
555 }
556 // post reactor setup operations
557 precon->setup();
558}
559
561{
562 if (!m_integ) {
563 throw CanteraError("ReactorNet::updatePreconditioner",
564 "Must only be called after ReactorNet is initialized.");
565 }
566 auto precon = m_integ->preconditioner();
567 precon->setGamma(gamma);
568 precon->updatePreconditioner();
569}
570
572 // check for non-mole-based reactors and throw an error otherwise
573 for (auto reactor : m_reactors) {
575 throw CanteraError("ReactorNet::checkPreconditionerSupported",
576 "Preconditioning is only supported for type *MoleReactor,\n"
577 "Reactor type given: '{}'.",
578 reactor->type());
579 }
580 }
581}
582
583}
Header file for class Cantera::Array2D.
Header file for base class WallBase.
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Definition Array.h:176
Base class for exceptions thrown by Cantera classes.
vector< double > m_paramScales
Scaling factors for each sensitivity parameter.
Definition FuncEval.h:184
bool suppressErrors() const
Get current state of error suppression.
Definition FuncEval.h:166
vector< double > m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
Definition FuncEval.h:181
An array index is out of range.
Abstract base class for ODE system integrators.
Definition Integrator.h:51
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
Definition Integrator.h:242
virtual void setMaxSteps(int nmax)
Set the maximum number of time-steps the integrator can take before reaching the next output time.
Definition Integrator.h:260
virtual int maxSteps()
Returns the maximum number of time-steps the integrator can take before reaching the next output time...
Definition Integrator.h:266
virtual void setMaxErrTestFails(int n)
Set the maximum permissible number of error test failures.
Definition Integrator.h:252
void setNetwork(ReactorNet *net)
Set the ReactorNet that this reactor belongs to.
void setPreconditioner(shared_ptr< PreconditionerBase > preconditioner)
Set preconditioner used by the linear solver.
void setLinearSolverType(const string &linSolverType="DENSE")
Set the type of linear solver used in the integration.
void preconditionerSetup(double t, double *y, double gamma) override
Evaluate the setup processes for the Jacobian preconditioner.
double step()
Advance the state of all reactors with respect to the independent variable (time or space).
virtual int lastOrder() const
Returns the order used for last solution step of the ODE integrator The function is intended for inte...
void eval(double t, double *y, double *ydot, double *p) override
Evaluate the right-hand-side ODE function.
void initialize()
Initialize the reactor network.
void advance(double t)
Advance the state of all reactors in the independent variable (time or space).
size_t neq() const override
Number of equations.
Definition ReactorNet.h:208
vector< size_t > m_start
m_start[n] is the starting point in the state vector for reactor n
Definition ReactorNet.h:333
vector< double > m_LHS
m_LHS is a vector representing the coefficients on the "left hand side" of each governing equation
Definition ReactorNet.h:359
double m_initial_time
The initial value of the independent variable in the system.
Definition ReactorNet.h:326
void evalJacobian(double t, double *y, double *ydot, double *p, Array2D *j)
Evaluate the Jacobian matrix for the reactor network.
double time()
Current value of the simulation time [s], for reactor networks that are solved in the time domain.
void getConstraints(double *constraints) override
Given a vector of length neq(), mark which variables should be considered algebraic constraints.
double m_time
The independent variable in the system.
Definition ReactorNet.h:323
AnyMap solverStats() const
Get solver stats from integrator.
Reactor & reactor(int n)
Return a reference to the n-th reactor in this network.
Definition ReactorNet.h:144
virtual void setMaxSteps(int nmax)
Set the maximum number of internal integration steps the integrator will take before reaching the nex...
string componentName(size_t i) const
Return the name of the i-th component of the global state vector.
void addReactor(Reactor &r)
Add the reactor r to this reactor network.
void getStateDae(double *y, double *ydot) override
Fill in the vectors y and ydot with the current state of the system.
void setInitialTime(double time)
Set the initial value of the independent variable (typically time).
virtual void getDerivative(int k, double *dky)
Return k-th derivative at the current state of the system.
void setMaxErrTestFails(int nmax)
Set the maximum number of error test failures permitted by the CVODES integrator in a single step.
size_t registerSensitivityParameter(const string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
double m_maxstep
Maximum integrator internal timestep. Default of 0.0 means infinity.
Definition ReactorNet.h:344
double distance()
Current position [m] along the length of the reactor network, for reactors that are solved as a funct...
void setSensitivityTolerances(double rtol, double atol)
Set the relative and absolute tolerances for integrating the sensitivity equations.
int maxSteps()
Returns the maximum number of internal integration steps the integrator will take before reaching the...
virtual void setDerivativeSettings(AnyMap &settings)
Set derivative settings of all reactors.
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.
void updateState(double *y)
Update the state of all the reactors in the network to correspond to the values in the solution vecto...
void getState(double *y) override
Fill in the vector y with the current state of the system.
void setAdvanceLimits(const double *limits)
Set absolute step size limits during advance.
double rtol()
Relative tolerance.
Definition ReactorNet.h:90
size_t globalComponentIndex(const string &component, size_t reactor=0)
Return the index corresponding to the component named component in the reactor with index reactor in ...
bool m_timeIsIndependent
Indicates whether time or space is the independent variable.
Definition ReactorNet.h:349
double atol()
Absolute integration tolerance.
Definition ReactorNet.h:95
bool hasAdvanceLimits() const
Check whether ReactorNet object uses advance limits.
void setMaxTimeStep(double maxstep)
Set the maximum integrator step.
void evalDae(double t, double *y, double *ydot, double *p, double *residual) override
eval coupling for IDA / DAEs
virtual void checkPreconditionerSupported() const
Check that preconditioning is supported by all reactors in the network.
bool m_integrator_init
True if integrator initialization is current.
Definition ReactorNet.h:329
void reinitialize()
Reinitialize the integrator.
Integrator & integrator()
Return a reference to the integrator.
bool getAdvanceLimits(double *limits) const
Retrieve absolute step size limits during advance.
string linearSolverType() const
Problem type of integrator.
void updatePreconditioner(double gamma) override
Update the preconditioner based on already computed jacobian values.
void preconditionerSolve(double *rhs, double *output) override
Evaluate the linear system Ax=b where A is the preconditioner.
vector< string > m_paramNames
Names corresponding to each sensitivity parameter.
Definition ReactorNet.h:352
void setTolerances(double rtol, double atol)
Set the relative and absolute tolerances for the integrator.
virtual void getEstimate(double time, int k, double *yest)
Estimate a future state based on current derivatives.
Class Reactor is a general-purpose class for stirred reactors.
Definition Reactor.h:44
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:103
string type() const override
String indicating the reactor model implemented.
Definition Reactor.h:48
virtual size_t nSensParams() const
Number of sensitivity parameters associated with this reactor (including walls)
Definition Reactor.cpp:112
void initialize(double t0=0.0) override
Initialize the reactor.
Definition Reactor.cpp:84
virtual bool preconditionerSupported() const
Return a false if preconditioning is not supported or true otherwise.
Definition Reactor.h:238
virtual bool isOde() const
Indicate whether the governing equations for this reactor type are a system of ODEs or DAEs.
Definition Reactor.h:55
virtual bool timeIsIndependent() const
Indicates whether the governing equations for this reactor are functions of time or a spatial variabl...
Definition Reactor.h:61
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:158
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition utilities.h:104
Integrator * newIntegrator(const string &itype)
Create new Integrator object.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
@ BDF_Method
Backward Differentiation.
Definition Integrator.h:31
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...