From c18ba5324d67415717f71d43191c3fe91dd12a83 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 22 Apr 2019 15:37:54 -0500 Subject: [PATCH] implemented automatic reduction of advance via ReactorNet::advanceTowards --- include/cantera/numerics/CVodesIntegrator.h | 3 + include/cantera/numerics/Integrator.h | 12 ++ include/cantera/zeroD/ReactorNet.h | 32 ++++++ interfaces/cython/cantera/_cantera.pxd | 7 ++ .../cantera/examples/reactors/reactor3.py | 61 ++++++++++ interfaces/cython/cantera/reactor.pyx | 65 +++++++++++ src/numerics/CVodesIntegrator.cpp | 32 ++++++ src/zeroD/ReactorNet.cpp | 106 +++++++++++++++++- 8 files changed, 317 insertions(+), 1 deletion(-) create mode 100644 interfaces/cython/cantera/examples/reactors/reactor3.py diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h index 1c2c1d5c856..dea7e6890e9 100644 --- a/include/cantera/numerics/CVodesIntegrator.h +++ b/include/cantera/numerics/CVodesIntegrator.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..7cc2299088b 100644 --- a/include/cantera/numerics/Integrator.h +++ b/include/cantera/numerics/Integrator.h @@ -141,6 +141,18 @@ class Integrator return 0; } + //! n-th derivative of the output function at time tout. + virtual doublereal* 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..892299a74f8 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/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 estimate(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 during the last solution step + 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 03992e9580f..4852e5a0c21 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -573,7 +573,9 @@ cdef extern from "cantera/zeroD/ReactorNet.h": CxxReactorNet() void addReactor(CxxReactor&) void advance(double) except +translate_exception + double advanceTowards(double) except +translate_exception double step() except +translate_exception + void estimate(double, int, double *) except +translate_exception void reinitialize() except +translate_exception double time() void setInitialTime(double) @@ -588,6 +590,11 @@ cdef extern from "cantera/zeroD/ReactorNet.h": 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/reactor3.py b/interfaces/cython/cantera/examples/reactors/reactor3.py new file mode 100644 index 00000000000..46c3a7415ce --- /dev/null +++ b/interfaces/cython/cantera/examples/reactors/reactor3.py @@ -0,0 +1,61 @@ +""" +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 +import numpy as np + +import cantera as ct + +gas = ct.Solution('gri30.xml') +gas.TPX = 1001.0, ct.one_atm, 'H2:2,O2:1,N2:4' +r = ct.IdealGasConstPressureReactor(gas) + +sim = ct.ReactorNet([r]) +time = 0.0 +states = ct.SolutionArray(gas, extra=['t']) + +# limit advance when temperature difference is exceeded +ix = r.component_index('temperature') +delta_T_max = 20. +sim.set_atol_advance(ix, delta_T_max) +sim.verbose = True + +print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]')) +for n in range(100): + time += 1.e-5 + while sim.advance_towards(time) < time: + states.append(r.thermo.state, t=sim.time*1e3) + print('%10.3e %10.3f %10.3f %14.6e' % (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)) + +# Plot the results if matplotlib is installed. +# See http://matplotlib.org/ to get it. +if '--plot' in sys.argv[1:]: + import matplotlib.pyplot as plt + plt.clf() + plt.subplot(2, 2, 1) + plt.plot(states.t, states.T) + plt.xlabel('Time (ms)') + plt.ylabel('Temperature (K)') + plt.subplot(2, 2, 2) + plt.plot(states.t, states.X[:,gas.species_index('OH')]) + plt.xlabel('Time (ms)') + plt.ylabel('OH Mole Fraction') + plt.subplot(2, 2, 3) + plt.plot(states.t, states.X[:,gas.species_index('H')]) + plt.xlabel('Time (ms)') + plt.ylabel('H Mole Fraction') + plt.subplot(2, 2, 4) + plt.plot(states.t, states.X[:,gas.species_index('H2')]) + plt.xlabel('Time (ms)') + plt.ylabel('H2 Mole Fraction') + plt.tight_layout() + plt.show() +else: + print("To view a plot of these results, run this script with the option --plot") diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index db7730f1f97..3e4823f0f0d 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -854,6 +854,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 @@ -861,6 +872,18 @@ cdef class ReactorNet: """ return self.net.step() + def estimate(self, double t, int k): + """ + 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 not self.n_vars: + raise CanteraError('ReactorNet empty or not initialized.') + cdef np.ndarray[np.double_t, ndim = 1] yest = np.zeros(self.n_vars) + self.net.estimate(t, k, & yest[0]) + return yest + def reinitialize(self): """ Reinitialize the integrator after making changing to the state of the @@ -1059,6 +1082,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/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index d1930aae66e..95bbe7666d9 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -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..60741e2f651 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -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,47 @@ 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; + estimate(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 +191,28 @@ double ReactorNet::step() return m_time; } +void ReactorNet::estimate(double time, int k, double* yest) +{ + double* cvode_dky; + double deltat = time - m_time; + double factor = 1.; + + // initialize + cvode_dky = m_integ->solution(); + for (size_t j = 0; j < m_nv; j++) { + yest[j] = cvode_dky[j]; + } + + // taylor expansion + for (int n = 1; n <= k; n++) { + factor *= deltat / ((double)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,45 @@ 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) {