[Reactor] Implement ExtensibleReactorSurface

Author: speth

doc/sphinx/python/zerodim.rst

@@ -123,6 +123,14 @@ ReactorSurface
 ^^^^^^^^^^^^^^
 .. autoclass:: ReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
 
+ExtensibleReactorSurface
+^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
+
+ExtensibleMoleReactorSurface
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleMoleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
+
 Flow Controllers
 ----------------
 

include/cantera/zeroD/Reactor.h

@@ -140,7 +140,7 @@ class Reactor : public ReactorBase
     //! Evaluate terms related to Walls. Calculates #m_vdot and #m_Qdot based on
     //! wall movement and heat transfer.
     //! @param t     the current time
-    virtual void evalWalls(double t);
+    void evalWalls(double t) override;
 
     //! Pointer to the homogeneous Kinetics object that handles the reactions
     Kinetics* m_kin = nullptr;

include/cantera/zeroD/ReactorBase.h

@@ -125,6 +125,12 @@ class ReactorBase : public std::enable_shared_from_this<ReactorBase>
             "Area is undefined for reactors of type '{}'.", type());
     }
 
+    //! Evaluate contributions from walls connected to this reactor
+    virtual void evalWalls(double t) {
+        throw NotImplementedError("ReactorBase::evalWalls",
+            "Wall evaluation is undefined for reactors of type '{}'.", type());
+    }
+
     //! Set an area associated with a reactor [m²].
     //! Examples: surface area of ReactorSurface or cross section area of FlowReactor.
     virtual void setArea(double a) {
@@ -213,6 +219,18 @@ class ReactorBase : public std::enable_shared_from_this<ReactorBase>
         return m_surfaces.size();
     }
 
+    //! Production rates on surfaces.
+    //!
+    //! For bulk reactors, this contains the total production rates [kmol/s] of bulk
+    //! phase species due to reactions on all adjacent species.
+    //!
+    //! For surfaces, this contains the production rates [kmol/m²/s] of species on the
+    //! surface and all adjacent phases, in the order defined by the InterfaceKinetics
+    //! object.
+    const vector<double>& surfaceProductionRates() const {
+        return m_sdot;
+    }
+
     /**
      * Initialize the reactor. Called automatically by ReactorNet::initialize.
      */
@@ -482,6 +500,7 @@ class ReactorBase : public std::enable_shared_from_this<ReactorBase>
 
     vector<WallBase*> m_wall;
     vector<ReactorSurface*> m_surfaces;
+    vector<double> m_sdot; //!< species production rates on surfaces
 
     //! Vector of length nWalls(), indicating whether this reactor is on the left (0)
     //! or right (1) of each wall.

include/cantera/zeroD/ReactorDelegator.h

@@ -44,6 +44,9 @@ class ReactorAccessor
     //! by Reactor::eval, this method should be called in either a "replace" or "after"
     //! delegate for Reactor::evalWalls().
     virtual void setHeatRate(double q) = 0;
+
+    //! @copydoc ReactorBase::surfaceProductionRates
+    virtual vector<double>& surfaceProductionRates() = 0;
 };
 
 //! Delegate methods of the Reactor class to external functions
@@ -52,8 +55,9 @@ template <class R>
 class ReactorDelegator : public Delegator, public R, public ReactorAccessor
 {
 public:
-    ReactorDelegator(shared_ptr<Solution> phase, bool clone, const string& name="(none)")
-        : R(phase, clone, name)
+    template <class... Args>
+    ReactorDelegator(Args&&... args)
+        : R(std::forward<Args>(args)...)
     {
         install("initialize", m_initialize, [this](double t0) { R::initialize(t0); });
         install("getState", m_getState,
@@ -67,7 +71,9 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
                 R::eval(t, LHS, RHS);
             }
         );
-        install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); });
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); });
+        }
         install("componentName", m_componentName,
             [this](size_t k) { return R::componentName(k); });
         install("componentIndex", m_componentIndex,
@@ -104,7 +110,11 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
     }
 
     void evalWalls(double t) override {
-        m_evalWalls(t);
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            m_evalWalls(t);
+        } else {
+            ReactorBase::evalWalls(t);
+        }
     }
 
     string componentName(size_t k) override {
@@ -122,19 +132,43 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
     }
 
     double expansionRate() const override {
-        return R::m_vdot;
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            return R::m_vdot;
+        } else {
+            throw NotImplementedError("ReactorDelegator::expansionRate",
+                "Expansion rate is undefined for reactors of type '{}'.", type());
+        }
     }
 
     void setExpansionRate(double v) override {
-        R::m_vdot = v;
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            R::m_vdot = v;
+        } else {
+            throw NotImplementedError("ReactorDelegator::setExpansionRate",
+                "Expansion rate is undefined for reactors of type '{}'.", type());
+        }
     }
 
     double heatRate() const override {
-        return R::m_Qdot;
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            return R::m_Qdot;
+        } else {
+            throw NotImplementedError("ReactorDelegator::heatRate",
+                "Heat rate is undefined for reactors of type '{}'.", type());
+        }
     }
 
     void setHeatRate(double q) override {
-        R::m_Qdot = q;
+        if constexpr (std::is_base_of<Reactor, R>::value) {
+            R::m_Qdot = q;
+        } else {
+            throw NotImplementedError("ReactorDelegator::setHeatRate",
+                "Heat rate is undefined for reactors of type '{}'.", type());
+        }
+    }
+
+    vector<double>& surfaceProductionRates() override {
+        return R::m_sdot;
     }
 
 private:

include/cantera/zeroD/ReactorSurface.h

@@ -91,10 +91,6 @@ class ReactorSurface : public ReactorBase
     //! number of surface species.
     void getCoverages(double* cov) const;
 
-    const vector<double>& netProductionRates() const {
-        return m_sdot;
-    }
-
     void getState(double* y) override;
     void initialize(double t0=0.0) override;
     void updateState(double* y) override;
@@ -120,7 +116,6 @@ class ReactorSurface : public ReactorBase
     shared_ptr<SurfPhase> m_surf;
     shared_ptr<InterfaceKinetics> m_kinetics;
     vector<ReactorBase*> m_reactors;
-    vector<double> m_sdot; //!< species production rates for all phases
 };
 
 //! A surface where the state variables are the total number of moles of each species.

interfaces/cython/cantera/reactor.pxd

@@ -216,6 +216,7 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera":
         void setExpansionRate(double)
         double heatRate()
         void setHeatRate(double)
+        vector[double]& surfaceProductionRates()
 
 
 ctypedef CxxReactorAccessor* CxxReactorAccessorPtr
@@ -272,6 +273,11 @@ cdef class FlowReactor(Reactor):
 cdef class ExtensibleReactor(Reactor):
     cdef public _delegates
     cdef CxxReactorAccessor* accessor
+    cdef public double[:] surface_production_rates
+    """
+    An array containing the total production rates [kmol/s] of bulk species on
+    surfaces. Updated from adjacent surfaces as part of ``Reactor.eval``.
+    """
 
 cdef class ReactorSurface(ReactorBase):
     cdef CxxReactorSurface* surface
@@ -280,6 +286,15 @@ cdef class ReactorSurface(ReactorBase):
 cdef class FlowReactorSurface(ReactorSurface):
     pass
 
+cdef class ExtensibleReactorSurface(ReactorSurface):
+    cdef public _delegates
+    cdef public double[:] surface_production_rates
+    """
+    An array containing the net production rates [kmol/m²/s] of surface and bulk species
+    on the surface, with the order corresponding to the `InterfaceKinetics` object.
+    Updated as part of ``ReactorSurface.update_state``.
+    """
+
 cdef class ConnectorNode:
     cdef shared_ptr[CxxConnectorNode] _node
     cdef CxxConnectorNode* node

interfaces/cython/cantera/reactor.pyx

@@ -567,6 +567,9 @@ cdef class ExtensibleReactor(Reactor):
 
     def __cinit__(self, *args, **kwargs):
         self.accessor = dynamic_cast[CxxReactorAccessorPtr](self.rbase)
+        cdef vector[double]* sdot = \
+            &dynamic_cast[CxxReactorAccessorPtr](self.rbase).surfaceProductionRates()
+        self.surface_production_rates = <double[:sdot.size()]> sdot.data()
 
     def __init__(self, *args, **kwargs):
         assign_delegates(self, dynamic_cast[CxxDelegatorPtr](self.rbase))
@@ -728,6 +731,9 @@ cdef class ReactorSurface(ReactorBase):
             raise TypeError("Parameter 'r' should be a ReactorBase object or a list "
                             "of ReactorBase objects.")
 
+        if kind is None and self.reactor_type != "ReactorSurface":
+            kind = self.reactor_type
+
         if kind is not None:
             self._rbase = CxxNewReactorSurface(stringify(kind), phase._base, cxx_adj,
                                                clone, stringify(name))
@@ -893,6 +899,58 @@ cdef class FlowReactorSurface(ReactorSurface):
         (<CxxFlowReactorSurface*>self.surface).setInitialMaxErrorFailures(nsteps)
 
 
+cdef class ExtensibleReactorSurface(ReactorSurface):
+    """
+    A base class for a reactor surface with delegated methods where the base
+    functionality corresponds to the `ReactorSurface` class.
+
+    Methods of the base class can be modified in the same way as in class
+    `ExtensibleReactor`. The methods that can be modified are::
+
+    - ``initialize(self, t0: double) -> None``
+    - ``get_state(self, y : double[:]) -> None``
+    - ``update_state(self, y : double[:]) -> None``
+    - ``eval(self, t : double, LHS : double[:], RHS : double[:]) -> None``
+    - ``component_name(i : int) -> string``
+    - ``component_index(name: string) -> int``
+    """
+
+    reactor_type = "ExtensibleReactorSurface"
+
+    delegatable_methods = {
+        'initialize': ('initialize', 'void(double)'),
+        'get_state': ('getState', 'void(double*)'),
+        'update_state': ('updateState', 'void(double*)'),
+        'eval': ('eval', 'void(double, double*, double*)'),
+        'component_name': ('componentName', 'string(size_t)'),
+        'component_index': ('componentIndex', 'size_t(string)'),
+    }
+
+    def __init__(self, *args, **kwargs):
+        assign_delegates(self, dynamic_cast[CxxDelegatorPtr](self.rbase))
+        cdef vector[double]* sdot = \
+            &dynamic_cast[CxxReactorAccessorPtr](self.rbase).surfaceProductionRates()
+        self.surface_production_rates = <double[:sdot.size()]> sdot.data()
+        super().__init__(*args, **kwargs)
+
+    property n_vars:
+        """
+        Get/Set the number of state variables in the reactor.
+        """
+        def __get__(self):
+            return self.reactor.neq()
+        def __set__(self, n):
+            dynamic_cast[CxxReactorAccessorPtr](self.rbase).setNEq(n)
+
+
+cdef class ExtensibleMoleReactorSurface(ExtensibleReactorSurface):
+    """
+    A variant of `ExtensibleReactorSurface` where the base behavior corresponds to the
+    :ct:`MoleReactorSurface` class.
+    """
+    reactor_type = "ExtensibleMoleReactorSurface"
+
+
 cdef class ConnectorNode:
     """
     Common base class for walls and flow devices.

src/zeroD/IdealGasConstPressureMoleReactor.cpp

@@ -120,10 +120,6 @@ void IdealGasConstPressureMoleReactor::getJacobianElements(
     // double molarVol = m_thermo->molarVolume();
     for (int k = 0; k < dnk_dnj.outerSize(); k++) {
         for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
-            // gas phase species need the addition of  V / N * omega_dot
-            // if (static_cast<size_t>(it.row()) < m_nsp) {
-            //     it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
-            // }
             trips.emplace_back(it.row() + offset, it.col() + offset, it.value());
         }
     }

src/zeroD/Reactor.cpp

@@ -293,7 +293,7 @@ void Reactor::updateSurfaceProductionRates()
 {
     m_sdot.assign(m_nsp, 0.0);
     for (auto& S : m_surfaces) {
-        const auto& sdot = S->netProductionRates();
+        const auto& sdot = S->surfaceProductionRates();
         size_t offset = S->kinetics()->kineticsSpeciesIndex(m_thermo->speciesName(0));
         for (size_t k = 0; k < m_nsp; k++) {
             m_sdot[k] += sdot[offset + k] * S->area();

src/zeroD/ReactorFactory.cpp

@@ -127,6 +127,18 @@ ReactorSurfaceFactory::ReactorSurfaceFactory()
         [](shared_ptr<Solution> surf, const vector<shared_ptr<ReactorBase>>& reactors,
            bool clone, const string& name)
         { return new FlowReactorSurface(surf, reactors, clone, name); });
+    reg("ExtensibleReactorSurface",
+        [](shared_ptr<Solution> surf, const vector<shared_ptr<ReactorBase>>& reactors,
+           bool clone, const string& name)
+        { return new ReactorDelegator<ReactorSurface>(surf, reactors, clone, name); });
+    reg("ExtensibleMoleReactorSurface",
+        [](shared_ptr<Solution> surf, const vector<shared_ptr<ReactorBase>>& reactors,
+           bool clone, const string& name)
+        { return new ReactorDelegator<MoleReactorSurface>(surf, reactors, clone, name); });
+    reg("ExtensibleFlowReactorSurface",
+        [](shared_ptr<Solution> surf, const vector<shared_ptr<ReactorBase>>& reactors,
+           bool clone, const string& name)
+        { return new ReactorDelegator<FlowReactorSurface>(surf, reactors, clone, name); });
 }
 
 // ---------- free functions ----------