diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h index 1c2c1d5c856..1332336f017 100644 --- a/include/cantera/numerics/CVodesIntegrator.h +++ b/include/cantera/numerics/CVodesIntegrator.h @@ -3,7 +3,7 @@ */ // 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 #define CT_CVODESWRAPPER_H @@ -41,6 +41,8 @@ class CVodesIntegrator : public Integrator virtual doublereal step(double tout); virtual double& solution(size_t k); virtual double* solution(); + virtual double* derivative(double tout, int n); + virtual int lastOrder() const; virtual int nEquations() const { return static_cast(m_neq); } @@ -89,6 +91,7 @@ class CVodesIntegrator : public Integrator double m_t0; double m_time; //!< The current integrator time N_Vector m_y, m_abstol; + N_Vector m_dky; int m_type; int m_itol; int m_method; diff --git a/include/cantera/numerics/Integrator.h b/include/cantera/numerics/Integrator.h index 6e0cdfc2c97..3c25847209f 100644 --- a/include/cantera/numerics/Integrator.h +++ b/include/cantera/numerics/Integrator.h @@ -7,7 +7,7 @@ */ // 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 #define CT_INTEGRATOR_H @@ -141,6 +141,18 @@ class Integrator 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("nEquations"); + return 0; + } + //! The number of equations. virtual int nEquations() const { warn("nEquations"); diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index dedd8e1d25c..d8fd6a6a552 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -1,7 +1,7 @@ //! @file ReactorNet.h // 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 #define CT_REACTORNET_H @@ -81,9 +81,21 @@ class ReactorNet : public FuncEval */ void advance(doublereal time); + /** + * Advance the state of all reactors in time. Take as many internal + * timesteps as necessary towards *time*, which will be automatically + * reduced to stay within an absolute tolerance (if specified). + * Returns the time at the end of integration. + * @param time Time to advance to (s). + */ + double advanceTowards(double time); + //! Advance the state of all reactors in time. double step(); + //! Estimate a future state based on current derivatives. + virtual void getEstimate(double time, int k, double* yest); + //@} //! Add the reactor *r* to this reactor network. @@ -164,6 +176,9 @@ class ReactorNet : public FuncEval virtual void getState(doublereal* y); + //! Return k-th derivative at the current time + virtual void getDerivative(int k, double* dky); + virtual size_t nparams() { return m_sens_params.size(); } @@ -222,6 +237,20 @@ class ReactorNet : public FuncEval return m_integ->maxSteps(); } + //! Returns the order used for last solution step of the ODE integrator + virtual int lastOrder() { + return m_integ->lastOrder(); + } + + //! Set absolute tolerances during advanceTowards + virtual void setAtolAdvance(const double* atols); + + //! Retrieve absolute tolerances during advanceTowards + virtual void getAtolAdvance(double* atols); + + //! Set an individual absolute tolerance during advanceTowards + virtual void setAtolAdvance(size_t k, double atol); + protected: //! Initialize the reactor network. Called automatically the first time //! advance or step is called. @@ -246,11 +275,14 @@ class ReactorNet : public FuncEval int m_maxErrTestFails; bool m_verbose; + bool m_limitadvance; //! Names corresponding to each sensitivity parameter std::vector m_paramNames; vector_fp m_ydot; + vector_fp m_yest; + vector_fp m_atoladvance; }; } diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 1c304ee3928..119f716e19f 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -581,7 +581,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": CxxReactorNet() void addReactor(CxxReactor&) void advance(double) except +translate_exception + double advanceTowards(double) except +translate_exception double step() except +translate_exception + void getEstimate(double, int, double *) except +translate_exception void reinitialize() except +translate_exception double time() void setInitialTime(double) @@ -596,6 +598,11 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void setVerbose(cbool) size_t neq() void getState(double*) + void getDerivative(int, double *) except +translate_exception + int lastOrder() except +translate_exception + void setAtolAdvance(double*) + void getAtolAdvance(double*) + void setAtolAdvance(int, double) except +translate_exception string componentName(size_t) except +translate_exception void setSensitivityTolerances(double, double) diff --git a/interfaces/cython/cantera/examples/reactors/reactor1.py b/interfaces/cython/cantera/examples/reactors/reactor1.py index a7347a3dc2c..a2a55dd6a51 100644 --- a/interfaces/cython/cantera/examples/reactors/reactor1.py +++ b/interfaces/cython/cantera/examples/reactors/reactor1.py @@ -1,5 +1,7 @@ """ Constant-pressure, adiabatic kinetics simulation. +Simulation is equivalent to reactor1.py, except for an automatically +reduced step size when large changes of temperature are anticipated. """ import sys @@ -15,13 +17,22 @@ time = 0.0 states = ct.SolutionArray(gas, extra=['t']) -print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]')) +# limit advance when temperature difference is exceeded +ix = sim.component_index(r.name, 'temperature') +delta_T_max = 20. +sim.set_atol_advance(ix, delta_T_max) +sim.verbose = True + +print('{:10s} {:10s} {:10s} {:14s}'.format('t [s]','T [K]','P [Pa]','u [J/kg]')) for n in range(100): time += 1.e-5 - sim.advance(time) + while sim.advance_towards(time) < time: + states.append(r.thermo.state, t=sim.time*1e3) + print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(sim.time, r.T, + r.thermo.P, r.thermo.u)) states.append(r.thermo.state, t=time*1e3) - print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T, - r.thermo.P, r.thermo.u)) + print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(sim.time, r.T, + r.thermo.P, r.thermo.u)) # Plot the results if matplotlib is installed. # See http://matplotlib.org/ to get it. diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 92302359bd3..87e8f533f7f 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -866,6 +866,17 @@ cdef class ReactorNet: """ self.net.advance(t) + def advance_towards(self, double t): + """ + Advance the state of the reactor network in time from the current + time towards time *t* [s]. If an absolute advance tolerance is + specified, the reactor state at the end of the timestep is estimated + prior to advancing. If the difference exceed limits, the step size is + reduced by half until the projected end state remains within the + limits. Returns the time reached at the end of integration. + """ + return self.net.advanceTowards(t) + def step(self): """ Take a single internal time step. The time after taking the step is @@ -873,6 +884,21 @@ cdef class ReactorNet: """ return self.net.step() + def get_estimate(self, double t, int k=-1): + """ + Estimate the state of the reactor network at a future time *t* [s] + based on derivatives at the current time. The estimate uses a Taylor + expansion with order *k*. If k<1, the last order of the integration + step is used. + """ + if not self.n_vars: + raise CanteraError('ReactorNet empty or not initialized.') + if k < 1: + k = self.net.lastOrder() + cdef np.ndarray[np.double_t, ndim = 1] yest = np.zeros(self.n_vars) + self.net.getEstimate(t, k, & yest[0]) + return yest + def reinitialize(self): """ Reinitialize the integrator after making changing to the state of the @@ -966,6 +992,26 @@ cdef class ReactorNet: def __set__(self, pybool v): self.net.setVerbose(v) + def component_index(self, reactor_name, name): + """ + Returns the index of a component named *name* of a reactor named + *reactor_name* within the global state vector. I.e. this determines the + (absolute) index of the component. *reactor_name* is the name of the + reactor that holds the component, e.g. 'reactor1'. *name* is either a + species name or the name of a reactor state variable, e.g. 'int_energy', + 'temperature', depending on the reactor's equations. + """ + k = 0 + for R in self._reactors: + if R.name == reactor_name: + k += R.component_index(name) + break + else: + k += R.n_vars + if k == self.n_vars: + raise IndexError('No such component: {!r} in {!r}'.format(name, reactor_name)) + return k + def component_name(self, int i): """ Return the name of the i-th component of the global state vector. The @@ -1071,6 +1117,48 @@ cdef class ReactorNet: self.net.getState(&y[0]) 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 last_order: + """ + get the integration order of the last step evaluation + """ + def __get__(self): + return self.net.lastOrder() + + property atol_advance: + """ + Get or set the absolute tolerance a state during `advance_towards` + (positive values are considered; negative values are ignored). + """ + def __get__(self): + cdef np.ndarray[np.double_t, ndim=1] atols = np.empty(self.n_vars) + self.net.getAtolAdvance(&atols[0]) + return atols + + def __set__(self, atols): + if len(atols) != self.n_vars: + raise ValueError('array must be of length n_vars') + + cdef np.ndarray[np.double_t, ndim=1] data = \ + np.ascontiguousarray(atols, dtype=np.double) + self.net.setAtolAdvance(&data[0]) + + def set_atol_advance(self, int index, double atol): + """ + Limit absolute tolerance of a state component during `advance_towards`. + (positive *atol* are considered; negative *atol* are ignored). + """ + self.net.setAtolAdvance(index, atol) + def advance_to_steady_state(self, int max_steps=10000, double residual_threshold=0., double atol=0., pybool return_residuals=False): diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index a7f7f201e6d..b7e509d3c4b 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -121,6 +121,45 @@ def test_disjoint2(self): self.assertNear(P1, self.r1.thermo.P) self.assertNear(P2, self.r2.thermo.P) + def test_integration_order(self): + T1, P1 = 300, 101325 + + self.make_reactors(n_reactors=1, T1=T1, P1=P1) + self.net.advance(1.0) + + # CVode integrator order is non-zero + self.assertTrue(self.net.last_order>0) + + 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_estimate(self): + T1, P1 = 300, 101325 + + self.make_reactors(n_reactors=1, T1=T1, P1=P1) + self.net.advance(1.0) + + # estimate should be sufficiently accurate + t_next = 1. + 1.e-4 + yest = self.net.get_estimate(t_next) + self.net.advance(t_next) + y = self.net.get_state() + for i in range(self.net.n_vars): + self.assertNear(yest[i], y[i]) + def test_timestepping(self): self.make_reactors() @@ -205,6 +244,40 @@ def integrate(atol, rtol): self.assertTrue(n_baseline > n_rtol) self.assertTrue(n_baseline > n_atol) + def test_advance_towards(self): + def integrate(atol_H2 = None): + 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) + self.net.step() + comp = {self.net.component_name(i).split()[-1]: i + for i in range(self.net.n_vars)} + ix = comp['H2'] + if atol_H2 is not None: + self.net.set_atol_advance(ix, atol_H2) + self.assertEqual(self.net.atol_advance[ix], atol_H2) + + tEnd = 1.0 + tStep = 1.e-3 + nSteps = 0 + + t = tStep + while t < tEnd: + t_curr = self.net.advance_towards(t) + nSteps += 1 + if t_curr == t: + t += tStep + + return nSteps + + n_baseline = integrate() + n_advance_coarse = integrate(.01) + n_advance_fine = integrate(.001) + + self.assertTrue(n_advance_coarse > n_baseline) + self.assertTrue(n_advance_fine > n_advance_coarse) + def test_heat_transfer1(self): # Connected reactors reach thermal equilibrium after some time self.make_reactors(T1=300, T2=1000) diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index d1930aae66e..806c67dd45e 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -1,7 +1,7 @@ //! @file CVodesIntegrator.cpp // 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/base/stringUtils.h" @@ -85,6 +85,7 @@ CVodesIntegrator::CVodesIntegrator() : m_t0(0.0), m_y(0), m_abstol(0), + m_dky(0), m_type(DENSE+NOJAC), m_itol(CV_SS), m_method(CV_BDF), @@ -126,6 +127,9 @@ CVodesIntegrator::~CVodesIntegrator() if (m_abstol) { N_VDestroy_Serial(m_abstol); } + if (m_dky) { + N_VDestroy_Serial(m_dky); + } if (m_yS) { N_VDestroyVectorArray_Serial(m_yS, static_cast(m_np)); } @@ -274,6 +278,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) } m_y = N_VNew_Serial(static_cast(m_neq)); // allocate solution vector 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(m_neq)); // allocate derivative vector + N_VConst(0.0, m_dky); // check abs tolerance array size if (m_itol == CV_SV && m_nabs < m_neq) { throw CanteraError("CVodesIntegrator::initialize", @@ -484,6 +493,29 @@ double CVodesIntegrator::step(double tout) 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 { long int ne; diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 590d51acbc6..8d742c559a8 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -1,7 +1,7 @@ //! @file ReactorNet.cpp // 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/FlowDevice.h" @@ -20,7 +20,7 @@ ReactorNet::ReactorNet() : m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4), m_atols(1.0e-15), m_atolsens(1.0e-6), m_maxstep(0.0), m_maxErrTestFails(0), - m_verbose(false) + m_verbose(false), m_limitadvance(false) { suppressErrors(true); @@ -98,6 +98,8 @@ void ReactorNet::initialize() } m_ydot.resize(m_nv,0.0); + m_yest.resize(m_nv,0.0); + m_atoladvance.resize(m_nv,-1.0); m_atol.resize(neq()); fill(m_atol.begin(), m_atol.end(), m_atols); m_integ->setTolerances(m_rtol, neq(), m_atol.data()); @@ -136,6 +138,49 @@ void ReactorNet::advance(doublereal time) updateState(m_integ->solution()); } +double ReactorNet::advanceTowards(double time) +{ + if (!m_limitadvance) { + // take full step + advance(time); + return time; + } else if (!m_init || !m_integrator_init) { + step(); + } + + // ensure 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_atoladvance[j] > 0.) && ( delta > m_atoladvance[j]) ) { + exceeded = true; + if (m_verbose) { + writelog(" Limiting index {:d} (dt = {:9.4g}):" + "{:11.6g} > {:9.4g}\n", + j, t - m_time, delta, m_atoladvance[j]); + } + } + } + if (!exceeded) { + break; + } + t = .5 * (m_time + t); + } + advance(t); + return t; +} + double ReactorNet::step() { if (!m_init) { @@ -148,6 +193,26 @@ double ReactorNet::step() 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) { r.setNetwork(this); @@ -218,6 +283,53 @@ 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::setAtolAdvance(const double *atols) +{ + if (!m_init) { + initialize(); + } + m_limitadvance = false; + for (size_t j = 0; j < m_nv; j++) { + if (atols[j] > 0.) { + m_limitadvance = true; + } + m_atoladvance[j] = atols[j]; + } +} + +void ReactorNet::setAtolAdvance(size_t k, double atol) +{ + if (!m_init) { + initialize(); + } + if (k > m_nv) { + throw CanteraError("ReactorNet::setAtolAdvance", + "Index out of bounds."); + } + m_limitadvance = false; + m_atoladvance[k] = atol; + for (size_t j = 0; j < m_nv; j++) { + if (m_atoladvance[j] > 0.) { + m_limitadvance = true; + } + } +} + +void ReactorNet::getAtolAdvance(double *atols) +{ + for (size_t j = 0; j < m_nv; j++) { + atols[j] = m_atoladvance[j]; + } +} + size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor) { if (!m_init) {