implemented limited advance step

This commit is contained in:
Ingmar Schoegl 2019-04-22 15:37:54 -05:00 committed by Ray Speth
parent bbab606a20
commit ae792dde00
11 changed files with 421 additions and 23 deletions

View file

@ -3,7 +3,7 @@
*/ */
// This file is part of Cantera. See License.txt in the top-level directory or // This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information. // at https://cantera.org/license.txt for license and copyright information.
#ifndef CT_CVODESWRAPPER_H #ifndef CT_CVODESWRAPPER_H
#define CT_CVODESWRAPPER_H #define CT_CVODESWRAPPER_H
@ -41,6 +41,8 @@ public:
virtual doublereal step(double tout); virtual doublereal step(double tout);
virtual double& solution(size_t k); virtual double& solution(size_t k);
virtual double* solution(); virtual double* solution();
virtual double* derivative(double tout, int n);
virtual int lastOrder() const;
virtual int nEquations() const { virtual int nEquations() const {
return static_cast<int>(m_neq); return static_cast<int>(m_neq);
} }
@ -89,6 +91,7 @@ private:
double m_t0; double m_t0;
double m_time; //!< The current integrator time double m_time; //!< The current integrator time
N_Vector m_y, m_abstol; N_Vector m_y, m_abstol;
N_Vector m_dky;
int m_type; int m_type;
int m_itol; int m_itol;
int m_method; int m_method;

View file

@ -7,7 +7,7 @@
*/ */
// This file is part of Cantera. See License.txt in the top-level directory or // This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information. // at https://cantera.org/license.txt for license and copyright information.
#ifndef CT_INTEGRATOR_H #ifndef CT_INTEGRATOR_H
#define CT_INTEGRATOR_H #define CT_INTEGRATOR_H
@ -141,6 +141,18 @@ public:
return 0; return 0;
} }
//! n-th derivative of the output function at time tout.
virtual double* derivative(double tout, int n) {
warn("derivative");
return 0;
}
//! Order used during the last solution step
virtual int lastOrder() const {
warn("lastOrder");
return 0;
}
//! The number of equations. //! The number of equations.
virtual int nEquations() const { virtual int nEquations() const {
warn("nEquations"); warn("nEquations");

View file

@ -132,6 +132,20 @@ public:
//! @see componentIndex() //! @see componentIndex()
virtual std::string componentName(size_t k); virtual std::string componentName(size_t k);
//! Set absolute step size limits during advance
//! @param limits array of step size limits with length neq
virtual void setAdvanceLimits(const double* limits);
//! Retrieve absolute step size limits during advance
//! @param[out] limits array of step size limits with length neq
//! @returns True if at least one limit is set, False otherwise
virtual bool getAdvanceLimits(double* limits);
//! Set individual step size limit for compoment name *nm*
//! @param nm component name
//! @param limit value for step size limit
virtual void setAdvanceLimit(const std::string& nm, const double limit);
protected: protected:
//! Set reaction rate multipliers based on the sensitivity variables in //! Set reaction rate multipliers based on the sensitivity variables in
//! *params*. //! *params*.
@ -180,6 +194,8 @@ protected:
bool m_energy; bool m_energy;
size_t m_nv; size_t m_nv;
vector_fp m_advancelimits; //!< Advance step limit
// Data associated each sensitivity parameter // Data associated each sensitivity parameter
std::vector<SensitivityParameter> m_sensParams; std::vector<SensitivityParameter> m_sensParams;
}; };

View file

@ -1,7 +1,7 @@
//! @file ReactorNet.h //! @file ReactorNet.h
// This file is part of Cantera. See License.txt in the top-level directory or // This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information. // at https://cantera.org/license.txt for license and copyright information.
#ifndef CT_REACTORNET_H #ifndef CT_REACTORNET_H
#define CT_REACTORNET_H #define CT_REACTORNET_H
@ -81,6 +81,17 @@ public:
*/ */
void advance(doublereal time); void advance(doublereal time);
/**
* Advance the state of all reactors in time. Take as many internal
* timesteps as necessary towards *time*. If *applylimit* is true,
* the advance step will be automatically reduced if needed to
* stay within limits (set by setAdvanceLimit).
* Returns the time at the end of integration.
* @param time Time to advance to (s).
* @param applylimit Limit advance step (boolean).
*/
double advance(double time, bool applylimit);
//! Advance the state of all reactors in time. //! Advance the state of all reactors in time.
double step(); double step();
@ -164,6 +175,9 @@ public:
virtual void getState(doublereal* y); virtual void getState(doublereal* y);
//! Return k-th derivative at the current time
virtual void getDerivative(int k, double* dky);
virtual size_t nparams() { virtual size_t nparams() {
return m_sens_params.size(); return m_sens_params.size();
} }
@ -195,6 +209,10 @@ public:
return m_paramNames.at(p); return m_paramNames.at(p);
} }
//! Initialize the reactor network. Called automatically the first time
//! advance or step is called.
void initialize();
//! Reinitialize the integrator. Used to solve a new problem (different //! Reinitialize the integrator. Used to solve a new problem (different
//! initial conditions) but with the same configuration of the reactor //! initial conditions) but with the same configuration of the reactor
//! network. Can be called manually, or automatically after calling //! network. Can be called manually, or automatically after calling
@ -222,10 +240,25 @@ public:
return m_integ->maxSteps(); return m_integ->maxSteps();
} }
//! Set absolute step size limits during advance
virtual void setAdvanceLimits(const double* limits);
//! Retrieve absolute step size limits during advance
virtual bool getAdvanceLimits(double* limits);
protected: protected:
//! Initialize the reactor network. Called automatically the first time
//! advance or step is called. //! Estimate a future state based on current derivatives.
void initialize(); //! The function is intended for internal use by ReactorNet::advance
//! and deliberately not exposed in external interfaces.
virtual void getEstimate(double time, int k, double* yest);
//! Returns the order used for last solution step of the ODE integrator
//! The function is intended for internal use by ReactorNet::advance
//! and deliberately not exposed in external interfaces.
virtual int lastOrder() {
return m_integ->lastOrder();
}
std::vector<Reactor*> m_reactors; std::vector<Reactor*> m_reactors;
std::unique_ptr<Integrator> m_integ; std::unique_ptr<Integrator> m_integ;
@ -251,6 +284,8 @@ protected:
std::vector<std::string> m_paramNames; std::vector<std::string> m_paramNames;
vector_fp m_ydot; vector_fp m_ydot;
vector_fp m_yest;
vector_fp m_advancelimits;
}; };
} }

View file

@ -524,7 +524,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
size_t neq() size_t neq()
void getState(double*) void getState(double*)
void addSurface(CxxReactorSurface*) void addSurface(CxxReactorSurface*)
void setAdvanceLimit(string&, double) except +translate_exception
void addSensitivityReaction(size_t) except +translate_exception void addSensitivityReaction(size_t) except +translate_exception
void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception
size_t nSensParams() size_t nSensParams()
@ -602,8 +602,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
cdef cppclass CxxReactorNet "Cantera::ReactorNet": cdef cppclass CxxReactorNet "Cantera::ReactorNet":
CxxReactorNet() CxxReactorNet()
void addReactor(CxxReactor&) void addReactor(CxxReactor&)
void advance(double) except +translate_exception double advance(double, cbool) except +translate_exception
double step() except +translate_exception double step() except +translate_exception
void initialize() except +translate_exception
void reinitialize() except +translate_exception void reinitialize() except +translate_exception
double time() double time()
void setInitialTime(double) void setInitialTime(double)
@ -618,8 +619,11 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
void setVerbose(cbool) void setVerbose(cbool)
size_t neq() size_t neq()
void getState(double*) void getState(double*)
void getDerivative(int, double *) except +translate_exception
void setAdvanceLimits(double*)
cbool getAdvanceLimits(double*)
string componentName(size_t) except +translate_exception string componentName(size_t) except +translate_exception
size_t globalComponentIndex(string&, int) except +translate_exception
void setSensitivityTolerances(double, double) void setSensitivityTolerances(double, double)
double rtolSensitivity() double rtolSensitivity()
double atolSensitivity() double atolSensitivity()

View file

@ -12,16 +12,22 @@ gas.TPX = 1001.0, ct.one_atm, 'H2:2,O2:1,N2:4'
r = ct.IdealGasConstPressureReactor(gas) r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r]) sim = ct.ReactorNet([r])
time = 0.0 sim.verbose = True
# limit advance when temperature difference is exceeded
delta_T_max = 20.
r.set_advance_limit('temperature', delta_T_max)
dt_max = 1.e-5
t_end = 100 * dt_max
states = ct.SolutionArray(gas, extra=['t']) states = ct.SolutionArray(gas, extra=['t'])
print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]')) print('{:10s} {:10s} {:10s} {:14s}'.format('t [s]','T [K]','P [Pa]','u [J/kg]'))
for n in range(100): while sim.time < t_end:
time += 1.e-5 sim.advance(sim.time + dt_max)
sim.advance(time) states.append(r.thermo.state, t=sim.time*1e3)
states.append(r.thermo.state, t=time*1e3) print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(sim.time, r.T,
print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T, r.thermo.P, r.thermo.u))
r.thermo.P, r.thermo.u))
# Plot the results if matplotlib is installed. # Plot the results if matplotlib is installed.
# See http://matplotlib.org/ to get it. # See http://matplotlib.org/ to get it.

View file

@ -315,6 +315,17 @@ cdef class Reactor(ReactorBase):
self.reactor.getState(&y[0]) self.reactor.getState(&y[0])
return y return y
def set_advance_limit(self, name, limit):
"""
Limit absolute change of component *name* during `ReactorNet.advance`.
(positive *limit* values are considered; negative values disable a
previously set advance limit for a solution component). Note that
limits are disabled by default (with individual values set to -1.).
"""
if limit is None:
limit = -1.
self.reactor.setAdvanceLimit(stringify(name), limit)
cdef class Reservoir(ReactorBase): cdef class Reservoir(ReactorBase):
""" """
@ -859,12 +870,17 @@ cdef class ReactorNet:
self._reactors.append(r) self._reactors.append(r)
self.net.addReactor(deref(r.reactor)) self.net.addReactor(deref(r.reactor))
def advance(self, double t): def advance(self, double t, pybool apply_limit=True):
""" """
Advance the state of the reactor network in time from the current Advance the state of the reactor network in time from the current time
time to time *t* [s], taking as many integrator timesteps as necessary. towards time *t* [s], taking as many integrator timesteps as necessary.
If *apply_limit* is true and an advance limit is specified, the reactor
state at the end of the timestep is estimated prior to advancing. If
the difference exceed limits, the end time is reduced by half until
the projected end state remains within specified limits.
Returns the time reached at the end of integration.
""" """
self.net.advance(t) return self.net.advance(t, apply_limit)
def step(self): def step(self):
""" """
@ -873,6 +889,12 @@ cdef class ReactorNet:
""" """
return self.net.step() return self.net.step()
def initialize(self):
"""
Force initialization of the integrator after initial setup.
"""
self.net.initialize()
def reinitialize(self): def reinitialize(self):
""" """
Reinitialize the integrator after making changing to the state of the Reinitialize the integrator after making changing to the state of the
@ -966,6 +988,17 @@ cdef class ReactorNet:
def __set__(self, pybool v): def __set__(self, pybool v):
self.net.setVerbose(v) self.net.setVerbose(v)
def global_component_index(self, name, int reactor):
"""
Returns the index of a component named *name* of a reactor with index
*reactor* within the global state vector. I.e. this determines the
(absolute) index of the component, where *reactor* is the index of the
reactor that holds the component. *name* is either a species name or the
name of a reactor state variable, e.g. 'int_energy', 'temperature', etc.
depending on the reactor's equations.
"""
return self.net.globalComponentIndex(stringify(name), reactor)
def component_name(self, int i): def component_name(self, int i):
""" """
Return the name of the i-th component of the global state vector. The Return the name of the i-th component of the global state vector. The
@ -1071,6 +1104,38 @@ cdef class ReactorNet:
self.net.getState(&y[0]) self.net.getState(&y[0])
return y return y
def get_derivative(self, k):
"""
Get the k-th time derivative of the state vector of the reactor network.
"""
if not self.n_vars:
raise CanteraError('ReactorNet empty or not initialized.')
cdef np.ndarray[np.double_t, ndim = 1] dky = np.zeros(self.n_vars)
self.net.getDerivative(k, & dky[0])
return dky
property advance_limits:
"""
Get or set absolute limits for state changes during `ReactorNet.advance`
(positive values are considered; negative values disable a previously
set advance limit for a solution component). Note that limits are
disabled by default (with individual values set to -1.).
"""
def __get__(self):
cdef np.ndarray[np.double_t, ndim=1] limits = np.empty(self.n_vars)
self.net.getAdvanceLimits(&limits[0])
return limits
def __set__(self, limits):
if limits is None:
limits = -1. * np.ones([self.n_vars])
elif len(limits) != self.n_vars:
raise ValueError('array must be of length n_vars')
cdef np.ndarray[np.double_t, ndim=1] data = \
np.ascontiguousarray(limits, dtype=np.double)
self.net.setAdvanceLimits(&data[0])
def advance_to_steady_state(self, int max_steps=10000, def advance_to_steady_state(self, int max_steps=10000,
double residual_threshold=0., double atol=0., double residual_threshold=0., double atol=0.,
pybool return_residuals=False): pybool return_residuals=False):

View file

@ -121,6 +121,22 @@ class TestReactor(utilities.CanteraTest):
self.assertNear(P1, self.r1.thermo.P) self.assertNear(P1, self.r1.thermo.P)
self.assertNear(P2, self.r2.thermo.P) self.assertNear(P2, self.r2.thermo.P)
def test_derivative(self):
T1, P1 = 300, 101325
self.make_reactors(n_reactors=1, T1=T1, P1=P1)
self.net.advance(1.0)
# compare cvode derivative to numerical derivative
dydt = self.net.get_derivative(1)
dt = -self.net.time
dy = -self.net.get_state()
self.net.step()
dt += self.net.time
dy += self.net.get_state()
for i in range(self.net.n_vars):
self.assertNear(dydt[i], dy[i]/dt)
def test_timestepping(self): def test_timestepping(self):
self.make_reactors() self.make_reactors()
@ -205,6 +221,63 @@ class TestReactor(utilities.CanteraTest):
self.assertTrue(n_baseline > n_rtol) self.assertTrue(n_baseline > n_rtol)
self.assertTrue(n_baseline > n_atol) self.assertTrue(n_baseline > n_atol)
def test_advance_limits(self):
P0 = 10 * ct.one_atm
T0 = 1100
X0 = 'H2:1.0, O2:0.5, AR:8.0'
self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0)
limit_H2 = .01
ix = self.net.global_component_index('H2', 0)
self.r1.set_advance_limit('H2', limit_H2)
self.assertEqual(self.net.advance_limits[ix], limit_H2)
self.r1.set_advance_limit('H2', None)
self.assertEqual(self.net.advance_limits[ix], -1.)
self.r1.set_advance_limit('H2', limit_H2)
self.net.advance_limits = None
self.assertEqual(self.net.advance_limits[ix], -1.)
self.r1.set_advance_limit('H2', limit_H2)
self.net.advance_limits = 0 * self.net.advance_limits - 1.
self.assertEqual(self.net.advance_limits[ix], -1.)
def test_advance_with_limits(self):
def integrate(limit_H2 = None, apply=True):
P0 = 10 * ct.one_atm
T0 = 1100
X0 = 'H2:1.0, O2:0.5, AR:8.0'
self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0)
if limit_H2 is not None:
self.r1.set_advance_limit('H2', limit_H2)
ix = self.net.global_component_index('H2', 0)
self.assertEqual(self.net.advance_limits[ix], limit_H2)
tEnd = 1.0
tStep = 1.e-3
nSteps = 0
t = tStep
while t < tEnd:
t_curr = self.net.advance(t, apply_limit=apply)
nSteps += 1
if t_curr == t:
t += tStep
return nSteps
n_baseline = integrate()
n_advance_coarse = integrate(.01)
n_advance_fine = integrate(.001)
n_advance_negative = integrate(-.001)
n_advance_override = integrate(.001, False)
self.assertTrue(n_advance_coarse > n_baseline)
self.assertTrue(n_advance_fine > n_advance_coarse)
self.assertTrue(n_advance_negative == n_baseline)
self.assertTrue(n_advance_override == n_baseline)
def test_heat_transfer1(self): def test_heat_transfer1(self):
# Connected reactors reach thermal equilibrium after some time # Connected reactors reach thermal equilibrium after some time
self.make_reactors(T1=300, T2=1000) self.make_reactors(T1=300, T2=1000)

View file

@ -1,7 +1,7 @@
//! @file CVodesIntegrator.cpp //! @file CVodesIntegrator.cpp
// This file is part of Cantera. See License.txt in the top-level directory or // This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information. // at https://cantera.org/license.txt for license and copyright information.
#include "cantera/numerics/CVodesIntegrator.h" #include "cantera/numerics/CVodesIntegrator.h"
#include "cantera/base/stringUtils.h" #include "cantera/base/stringUtils.h"
@ -85,6 +85,7 @@ CVodesIntegrator::CVodesIntegrator() :
m_t0(0.0), m_t0(0.0),
m_y(0), m_y(0),
m_abstol(0), m_abstol(0),
m_dky(0),
m_type(DENSE+NOJAC), m_type(DENSE+NOJAC),
m_itol(CV_SS), m_itol(CV_SS),
m_method(CV_BDF), m_method(CV_BDF),
@ -126,6 +127,9 @@ CVodesIntegrator::~CVodesIntegrator()
if (m_abstol) { if (m_abstol) {
N_VDestroy_Serial(m_abstol); N_VDestroy_Serial(m_abstol);
} }
if (m_dky) {
N_VDestroy_Serial(m_dky);
}
if (m_yS) { if (m_yS) {
N_VDestroyVectorArray_Serial(m_yS, static_cast<sd_size_t>(m_np)); N_VDestroyVectorArray_Serial(m_yS, static_cast<sd_size_t>(m_np));
} }
@ -274,6 +278,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
} }
m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate solution vector m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate solution vector
N_VConst(0.0, m_y); N_VConst(0.0, m_y);
if (m_dky) {
N_VDestroy_Serial(m_dky); // free derivative vector if already allocated
}
m_dky = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate derivative vector
N_VConst(0.0, m_dky);
// check abs tolerance array size // check abs tolerance array size
if (m_itol == CV_SV && m_nabs < m_neq) { if (m_itol == CV_SV && m_nabs < m_neq) {
throw CanteraError("CVodesIntegrator::initialize", throw CanteraError("CVodesIntegrator::initialize",
@ -484,6 +493,29 @@ double CVodesIntegrator::step(double tout)
return m_time; return m_time;
} }
double* CVodesIntegrator::derivative(double tout, int n)
{
int flag = CVodeGetDky(m_cvode_mem, tout, n, m_dky);
if (flag != CV_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
f_errs = "Exceptions caught evaluating derivative:\n" + f_errs;
}
throw CanteraError("CVodesIntegrator::derivative",
"CVodes error encountered. Error code: {}\n{}\n"
"{}",
flag, m_error_message, f_errs);
}
return NV_DATA_S(m_dky);
}
int CVodesIntegrator::lastOrder() const
{
int ord;
CVodeGetLastOrder(m_cvode_mem, &ord);
return ord;
}
int CVodesIntegrator::nEvals() const int CVodesIntegrator::nEvals() const
{ {
long int ne; long int ne;

View file

@ -103,6 +103,7 @@ void Reactor::initialize(doublereal t0)
} }
} }
m_work.resize(maxnt); m_work.resize(maxnt);
m_advancelimits.resize(m_nv, -1.0);
} }
size_t Reactor::nSensParams() size_t Reactor::nSensParams()
@ -429,4 +430,50 @@ void Reactor::resetSensitivity(double* params)
} }
} }
void Reactor::setAdvanceLimits(const double *limits)
{
if (m_thermo == 0) {
throw CanteraError("getState",
"Error: reactor is empty.");
}
for (size_t j = 0; j < m_nv; j++) {
m_advancelimits[j] = limits[j];
}
}
bool Reactor::getAdvanceLimits(double *limits)
{
bool has_limit = false;
for (size_t j = 0; j < m_nv; j++) {
limits[j] = m_advancelimits[j];
has_limit |= (limits[j] > 0.);
}
return has_limit;
}
void Reactor::setAdvanceLimit(const string& nm, const double limit)
{
size_t k = componentIndex(nm);
if (m_thermo == 0) {
throw CanteraError("getState",
"Error: reactor is empty.");
}
if (m_nv == 0) {
if (m_net == 0) {
throw CanteraError("Reactor::setAdvanceLimit",
"Cannot set limit on a reactor that is not "
"assigned to a ReactorNet object.");
} else {
m_net->initialize();
}
} else if (k > m_nv) {
throw CanteraError("Reactor::setAdvanceLimit",
"Index out of bounds.");
}
m_advancelimits[k] = limit;
}
} }

View file

@ -1,7 +1,7 @@
//! @file ReactorNet.cpp //! @file ReactorNet.cpp
// This file is part of Cantera. See License.txt in the top-level directory or // This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information. // at https://cantera.org/license.txt for license and copyright information.
#include "cantera/zeroD/ReactorNet.h" #include "cantera/zeroD/ReactorNet.h"
#include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/FlowDevice.h"
@ -98,6 +98,8 @@ void ReactorNet::initialize()
} }
m_ydot.resize(m_nv,0.0); m_ydot.resize(m_nv,0.0);
m_yest.resize(m_nv,0.0);
m_advancelimits.resize(m_nv,-1.0);
m_atol.resize(neq()); m_atol.resize(neq());
fill(m_atol.begin(), m_atol.end(), m_atols); fill(m_atol.begin(), m_atol.end(), m_atols);
m_integ->setTolerances(m_rtol, neq(), m_atol.data()); m_integ->setTolerances(m_rtol, neq(), m_atol.data());
@ -136,6 +138,60 @@ void ReactorNet::advance(doublereal time)
updateState(m_integ->solution()); updateState(m_integ->solution());
} }
double ReactorNet::advance(double time, bool applylimit)
{
if (!m_init) {
initialize();
} else if (!m_integrator_init) {
reinitialize();
}
if (!applylimit) {
// take full step
advance(time);
return time;
}
bool limitadvance = getAdvanceLimits(m_advancelimits.data());
if (!limitadvance) {
// take full step
advance(time);
return time;
}
// ensure that gradient is available
while (lastOrder() < 1) {
step();
}
int k = lastOrder();
double t = time, delta;
double* y = m_integ->solution();
// reduce time step if limits are exceeded
while (true) {
bool exceeded = false;
getEstimate(t, k, &m_yest[0]);
for (size_t j = 0; j < m_nv; j++) {
delta = abs(m_yest[j] - y[j]);
if ( (m_advancelimits[j] > 0.) && ( delta > m_advancelimits[j]) ) {
exceeded = true;
if (m_verbose) {
writelog(" Limiting global state vector component {:d} (dt = {:9.4g}):"
"{:11.6g} > {:9.4g}\n",
j, t - m_time, delta, m_advancelimits[j]);
}
}
}
if (!exceeded) {
break;
}
t = .5 * (m_time + t);
}
advance(t);
return t;
}
double ReactorNet::step() double ReactorNet::step()
{ {
if (!m_init) { if (!m_init) {
@ -148,6 +204,26 @@ double ReactorNet::step()
return m_time; return m_time;
} }
void ReactorNet::getEstimate(double time, int k, double* yest)
{
// initialize
double* cvode_dky = m_integ->solution();
for (size_t j = 0; j < m_nv; j++) {
yest[j] = cvode_dky[j];
}
// Taylor expansion
double factor = 1.;
double deltat = time - m_time;
for (int n = 1; n <= k; n++) {
factor *= deltat / n;
cvode_dky = m_integ->derivative(m_time, n);
for (size_t j = 0; j < m_nv; j++) {
yest[j] += factor * cvode_dky[j];
}
}
}
void ReactorNet::addReactor(Reactor& r) void ReactorNet::addReactor(Reactor& r)
{ {
r.setNetwork(this); r.setNetwork(this);
@ -218,6 +294,35 @@ void ReactorNet::getState(double* y)
} }
} }
void ReactorNet::getDerivative(int k, double* dky)
{
double* cvode_dky = m_integ->derivative(m_time, k);
for (size_t j = 0; j < m_nv; j++) {
dky[j] = cvode_dky[j];
}
}
void ReactorNet::setAdvanceLimits(const double *limits)
{
if (!m_init) {
initialize();
}
for (size_t n = 0; n < m_reactors.size(); n++) {
m_reactors[n]->setAdvanceLimits(limits + m_start[n]);
}
}
bool ReactorNet::getAdvanceLimits(double *limits)
{
bool has_limit = false;
for (size_t n = 0; n < m_reactors.size(); n++) {
has_limit |= m_reactors[n]->getAdvanceLimits(limits + m_start[n]);
}
return has_limit;
}
size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor) size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
{ {
if (!m_init) { if (!m_init) {