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),
30 m_integ->setProblemType(DENSE + NOJAC);
34 void ReactorNet::setInitialTime(
double time)
37 m_integrator_init =
false;
40 void ReactorNet::setMaxTimeStep(
double maxstep)
46 void ReactorNet::setMaxErrTestFails(
int nmax)
48 m_maxErrTestFails = nmax;
52 void ReactorNet::setTolerances(
double rtol,
double atol)
63 void ReactorNet::setSensitivityTolerances(
double rtol,
double atol)
74 void ReactorNet::initialize()
77 debuglog(
"Initializing reactor network.\n", m_verbose);
78 if (m_reactors.empty()) {
80 "no reactors in network!");
83 for (
size_t n = 0; n < m_reactors.size(); n++) {
88 m_start.push_back(m_nv);
91 writelog(
"Reactor {:d}: {:d} variables.\n", n, nv);
94 if (r.
type() == FlowReactorType && m_reactors.size() > 1) {
96 "FlowReactors must be used alone.");
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);
108 writelog(
"Number of equations: {:d}\n", neq());
109 writelog(
"Maximum time step: {:14.6g}\n", m_maxstep);
111 m_integ->initialize(m_time, *
this);
112 m_integrator_init =
true;
116 void ReactorNet::reinitialize()
119 debuglog(
"Re-initializing reactor network.\n", m_verbose);
120 m_integ->reinitialize(m_time, *
this);
121 m_integrator_init =
true;
127 void ReactorNet::advance(doublereal time)
131 }
else if (!m_integrator_init) {
134 m_integ->integrate(time);
136 updateState(m_integ->solution());
139 double ReactorNet::step()
143 }
else if (!m_integrator_init) {
146 m_time = m_integ->step(m_time + 1.0);
147 updateState(m_integ->solution());
154 m_reactors.push_back(&r);
157 void ReactorNet::eval(doublereal t, doublereal* y,
158 doublereal* ydot, doublereal* p)
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);
167 double ReactorNet::sensitivity(
size_t k,
size_t p)
172 if (p >= m_sens_params.size()) {
174 "m_sens_params", p, m_sens_params.size()-1);
176 double denom = m_integ->solution(k);
180 return m_integ->sensitivity(k, p) / denom;
183 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
184 doublereal* ydot, doublereal* p,
Array2D* j)
188 for (
size_t n = 0; n < m_nv; n++) {
191 double dy = m_atol[n] + fabs(ysave)*m_rtol;
196 eval(t, y, m_ydot.data(), p);
199 for (
size_t m = 0; m < m_nv; m++) {
200 j->
value(m,n) = (m_ydot[m] - ydot[m])/dy;
206 void ReactorNet::updateState(doublereal* y)
209 for (
size_t n = 0; n < m_reactors.size(); n++) {
210 m_reactors[n]->updateState(y + m_start[n]);
214 void ReactorNet::getState(
double* y)
216 for (
size_t n = 0; n < m_reactors.size(); n++) {
217 m_reactors[n]->getState(y + m_start[n]);
221 size_t ReactorNet::globalComponentIndex(
const string& component,
size_t reactor)
226 return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
229 std::string ReactorNet::componentName(
size_t i)
const 231 for (
auto r : m_reactors) {
233 return r->name() +
": " + r->componentName(i);
238 throw CanteraError(
"ReactorNet::componentName",
"Index out of bounds");
241 size_t ReactorNet::registerSensitivityParameter(
242 const std::string& name,
double value,
double scale)
244 if (m_integrator_init) {
245 throw CanteraError(
"ReactorNet::registerSensitivityParameter",
246 "Sensitivity parameters cannot be added after the" 247 "integrator has been initialized.");
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;
Backward Differentiation.
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
virtual size_t nSensParams()
Number of sensitivity parameters associated with this reactor (including walls)
Header file for class Wall.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
virtual void initialize(doublereal t0=0.0)
Initialize the reactor.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Base class for exceptions thrown by Cantera classes.
doublereal & value(size_t i, size_t j)
Returns a changeable reference to position in the matrix.
void setNetwork(ReactorNet *net)
Set the ReactorNet that this reactor belongs to.
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
const doublereal SmallNumber
smallest number to compare to zero.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
virtual int type() const
Return a constant indicating the type of this Reactor.
virtual size_t neq()
Number of equations (state variables) for this reactor.
An array index is out of range.
Namespace for the Cantera kernel.
Class Reactor is a general-purpose class for stirred reactors.