implemented automatic reduction of advance via ReactorNet::advanceTow…

…ards

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

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

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

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)

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.

interfaces/cython/cantera/reactor.pyx

@@ -866,13 +866,39 @@ 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 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):

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)