implemented automatic reduction of advance via ReactorNet::advanceTow…

…ards

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<int>(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;

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");

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<std::string> m_paramNames;
 
     vector_fp m_ydot;
+    vector_fp m_yest;
+    vector_fp m_atoladvance;
 };
 }
 

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)

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")

interfaces/cython/cantera/reactor.pyx

@@ -854,13 +854,36 @@ 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
         returned.
         """
         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):

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<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
     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
     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;