From 9a2d1b2d80e8798009cf4ff424d7c7af400105d6 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Mon, 1 Jan 2024 17:32:43 -0500 Subject: [PATCH] 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 | 8 ++++++++ interfaces/cython/cantera/reactor.pxd | 1 + interfaces/cython/cantera/reactor.pyx | 11 +++++++++++ test/python/test_reactor.py | 23 +++++++++++++++++++++++ 4 files changed, 43 insertions(+) diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index a733e5e6d33..52a56388cc8 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/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 { diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 61ede6a7b57..f571897d80c 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/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 diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index d287eaec4ed..31750a19510 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/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 @@ -738,6 +739,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): """ diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 1629469c4a2..fbe97d602da 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3250,3 +3250,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())