19ReactorNet::ReactorNet()
24ReactorNet::~ReactorNet()
72 throw CanteraError(
"ReactorNet::time",
"Time is not the independent variable"
73 " for this reactor network.");
81 throw CanteraError(
"ReactorNet::distance",
"Distance is not the independent"
82 " variable for this reactor network.");
89 debuglog(
"Initializing reactor network.\n", m_verbose);
90 if (m_reactors.empty()) {
92 "no reactors in network!");
95 for (
size_t n = 0; n < m_reactors.size(); n++) {
103 writelog(
"Reactor {:d}: {:d} variables.\n", n, nv);
106 if (r.
type() ==
"FlowReactor" && m_reactors.size() > 1) {
108 "FlowReactors must be used alone.");
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);
123 m_integ->setPreconditioner(m_precon);
125 m_integ->initialize(
m_time, *
this);
130 if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
140 debuglog(
"Re-initializing reactor network.\n", m_verbose);
141 m_integ->reinitialize(
m_time, *
this);
142 if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
153 m_linearSolverType = linSolverType;
159 m_precon = preconditioner;
180 m_integ->integrate(
time);
213 double t =
time, delta;
214 double* y = m_integ->solution();
218 bool exceeded =
false;
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]) ) {
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]);
258 double* cvode_dky = m_integ->solution();
259 for (
size_t j = 0; j < m_nv; j++) {
260 yest[j] = cvode_dky[j];
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];
278 return m_integ->lastOrder();
286 for (
auto current : m_reactors) {
287 if (current->isOde() != r.
isOde()) {
289 "Cannot mix Reactor types using both ODEs and DAEs ({} and {})",
290 current->type(), r.
type());
294 "Cannot mix Reactor types using time and space as independent variables"
295 "\n({} and {})", current->type(), r.
type());
300 m_reactors.push_back(&r);
306 m_integ->setLinearSolverType(
"DENSE");
311 if (m_integ ==
nullptr) {
313 "Integrator has not been instantiated. Add one or more reactors first.");
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);
328 if (n == m_reactors.size() - 1) {
333 for (
size_t i =
m_start[n]; i < yEnd; i++) {
334 ydot[i] = m_RHS[i] /
m_LHS[i];
336 m_reactors[n]->resetSensitivity(p);
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);
355 for (
size_t n = 0; n < m_reactors.size(); n++) {
356 m_reactors[n]->getConstraints(constraints +
m_start[n]);
369 double denom = m_integ->solution(k);
373 return m_integ->sensitivity(k, p) / denom;
380 for (
size_t n = 0; n < m_nv; n++) {
383 double dy = m_atol[n] + fabs(ysave)*m_rtol;
388 eval(t, y, m_ydot.data(), p);
391 for (
size_t m = 0; m < m_nv; m++) {
392 j->
value(m,n) = (m_ydot[m] - ydot[m])/dy;
401 for (
size_t n = 0; n < m_reactors.size(); n++) {
402 m_reactors[n]->updateState(y +
m_start[n]);
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];
422 for (
size_t n = 0; n < m_reactors.size(); n++) {
423 m_reactors[n]->setAdvanceLimits(limits +
m_start[n]);
429 bool has_limit =
false;
430 for (
size_t n = 0; n < m_reactors.size(); n++) {
431 has_limit |= m_reactors[n]->hasAdvanceLimits();
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]);
447 for (
size_t n = 0; n < m_reactors.size(); n++) {
448 m_reactors[n]->getState(y +
m_start[n]);
454 for (
size_t n = 0; n < m_reactors.size(); n++) {
469 for (
auto r : m_reactors) {
471 return r->name() +
": " + r->componentName(i);
476 throw CanteraError(
"ReactorNet::componentName",
"Index out of bounds");
480 const string& name,
double value,
double scale)
483 throw CanteraError(
"ReactorNet::registerSensitivityParameter",
484 "Sensitivity parameters cannot be added after the"
485 "integrator has been initialized.");
496 for (
size_t i = 0; i < m_reactors.size(); i++) {
497 m_reactors[i]->setDerivativeSettings(settings);
504 return m_integ->solverStats();
513 return m_integ->linearSolverType();
523 "Must only be called after ReactorNet is initialized.");
525 m_integ->preconditionerSolve(m_nv, rhs, output);
533 auto precon = m_integ->preconditioner();
537 precon->setGamma(gamma);
539 vector<double> yCopy(m_nv);
543 precon->stateAdjustment(yCopy);
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) {
564 "Must only be called after ReactorNet is initialized.");
566 auto precon = m_integ->preconditioner();
567 precon->setGamma(gamma);
568 precon->updatePreconditioner();
573 for (
auto reactor : m_reactors) {
575 throw CanteraError(
"ReactorNet::checkPreconditionerSupported",
576 "Preconditioning is only supported for type *MoleReactor,\n"
577 "Reactor type given: '{}'.",
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.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
double & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
Base class for exceptions thrown by Cantera classes.
vector< double > m_paramScales
Scaling factors for each sensitivity parameter.
bool suppressErrors() const
Get current state of error suppression.
vector< double > m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
An array index is out of range.
Abstract base class for ODE system integrators.
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
virtual void setMaxSteps(int nmax)
Set the maximum number of time-steps the integrator can take before reaching the next output time.
virtual int maxSteps()
Returns the maximum number of time-steps the integrator can take before reaching the next output time...
virtual void setMaxErrTestFails(int n)
Set the maximum permissible number of error test failures.
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.
vector< size_t > m_start
m_start[n] is the starting point in the state vector for reactor n
vector< double > m_LHS
m_LHS is a vector representing the coefficients on the "left hand side" of each governing equation
double m_initial_time
The initial value of the independent variable in the system.
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.
AnyMap solverStats() const
Get solver stats from integrator.
Reactor & reactor(int n)
Return a reference to the n-th reactor in this network.
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.
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.
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.
double atol()
Absolute integration tolerance.
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.
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.
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.
size_t neq()
Number of equations (state variables) for this reactor.
string type() const override
String indicating the reactor model implemented.
virtual size_t nSensParams() const
Number of sensitivity parameters associated with this reactor (including walls)
void initialize(double t0=0.0) override
Initialize the reactor.
virtual bool preconditionerSupported() const
Return a false if preconditioning is not supported or true otherwise.
virtual bool isOde() const
Indicate whether the governing equations for this reactor type are a system of ODEs or DAEs.
virtual bool timeIsIndependent() const
Indicates whether the governing equations for this reactor are functions of time or a spatial variabl...
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Integrator * newIntegrator(const string &itype)
Create new Integrator object.
Namespace for the Cantera kernel.
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
@ BDF_Method
Backward Differentiation.
const double SmallNumber
smallest number to compare to zero.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...