Add a default evaluation function for reactors

It can sometimes be useful to add on to the default evaluation in such a
way that the standard before or after delegations cannot handle. A default
evaluation function allows for this.

include/cantera/zeroD/ReactorDelegator.h

@@ -52,6 +52,10 @@ class ReactorAccessor
     //! Set the state of the thermo object for surface *n* to correspond to the
     //! state of that surface
     virtual void restoreSurfaceState(size_t n) = 0;
+
+    //! Public access to the default evaluation function so it can be used in
+    //! replace functions
+    virtual void defaultEval(double t, double* LHS, double* RHS) = 0;
 };
 
 //! Delegate methods of the Reactor class to external functions
@@ -161,6 +165,10 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         m_build_jacobian(jacVector);
     }
 
+    void defaultEval(double t, double* LHS, double* RHS) override {
+        R::eval(t, LHS, RHS);
+    }
+
     // Public access to protected Reactor variables needed by derived classes
 
     void setNEq(size_t n) override {

interfaces/cython/cantera/reactor.pxd

@@ -215,6 +215,7 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera":
         void setHeatRate(double)
         void restoreThermoState() except +translate_exception
         void restoreSurfaceState(size_t) except +translate_exception
+        void defaultEval(double time, double* LHS, double* RHS)
 
 
 ctypedef CxxReactorAccessor* CxxReactorAccessorPtr

interfaces/cython/cantera/reactor.pyx

@@ -2,6 +2,7 @@
 # at https://cantera.org/license.txt for license and copyright information.
 
 import warnings
+import numpy as np
 from collections import defaultdict as _defaultdict
 import numbers as _numbers
 from cython.operator cimport dereference as deref
@@ -687,6 +688,16 @@ cdef class ExtensibleReactor(Reactor):
         """
         self.accessor.restoreSurfaceState(n)
 
+    def default_eval(self, time, LHS, RHS):
+        """
+        Evaluation of the base reactors `eval` function to be used in `replace`
+        functions and maintain original functionality.
+        """
+        assert len(LHS) == self.n_vars and len(RHS) == self.n_vars
+        cdef np.ndarray[np.double_t, ndim=1, mode="c"] rhs = np.frombuffer(RHS)
+        cdef np.ndarray[np.double_t, ndim=1, mode="c"] lhs = np.frombuffer(LHS)
+        self.accessor.defaultEval(time, &lhs[0], &rhs[0])
+
 
 cdef class ExtensibleIdealGasReactor(ExtensibleReactor):
     """

test/python/test_reactor.py

@@ -3047,3 +3047,26 @@ def replace_build_jacobian(self, jac_vector):
         # test that jacobian wall element is hard coded value
         jac = r.jacobian
         assert np.sum(jac) == 0
+
+    def test_replace_with_default_eval(self):
+        class ReplaceEvalReactor(ct.ExtensibleIdealGasConstPressureMoleReactor):
+
+            def replace_eval(self, t, LHS, RHS):
+                self.default_eval(t, LHS, RHS)
+
+        # setup thermo object
+        gas = ct.Solution("h2o2.yaml", "ohmech")
+        gas.set_equivalence_ratio(0.5, "H2:1.0", "O2:1.0")
+        gas.equilibrate("HP")
+        # replacement reactor
+        r = ReplaceEvalReactor(gas)
+        r.volume = 1.0
+        # default reactor
+        rstd = ct.IdealGasConstPressureMoleReactor(gas)
+        rstd.volume = r.volume
+        # network of both reactors
+        net = ct.ReactorNet([r, rstd])
+        net.preconditioner = ct.AdaptivePreconditioner()
+        net.advance_to_steady_state()
+        # reactors should have the same solution because the default is used
+        self.assertArrayNear(r.get_state(), rstd.get_state())