Implements MoleReactor and ConstPressureMoleReactor, tests, and docs.

This commit implements MoleReactor and ConstPressureMoleReactor classes, it adds the appropriate python interfaces for them, it also adds a test comparing surface chemistry of mole reactors to traditional reactors. test_component_names had a bug because the network was not initialized so self.net.n_vars was 0 and the loop was never entered. Adding a line to initialize the network reveal that the test failed for IdealGasMoleReactor so the appropriate lines were added to fix it. Add wrappers for Extensible...MoleReactor, update docs, comment fixes.
2025-02-25 18:55:29 -06:00 · 2022-09-01 13:26:39 -04:00
parent e3080edd57
commit a999ead68f
12 changed files with 562 additions and 70 deletions
--- a/doc/sphinx/cython/zerodim.rst
+++ b/doc/sphinx/cython/zerodim.rst
@@ -42,6 +42,10 @@ Reactor
 ^^^^^^^
 .. autoclass:: Reactor

+MoleReactor
+^^^^^^^^^^^
+.. autoclass:: MoleReactor(contents=None, *, name=None, energy='on')
+
 IdealGasReactor
 ^^^^^^^^^^^^^^^
 .. autoclass:: IdealGasReactor(contents=None, *, name=None, energy='on')
@@ -54,6 +58,10 @@ ConstPressureReactor
 ^^^^^^^^^^^^^^^^^^^^
 .. autoclass:: ConstPressureReactor(contents=None, *, name=None, energy='on')

+ConstPressureMoleReactor
+^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ConstPressureMoleReactor(contents=None, *, name=None, energy='on')
+
 IdealGasConstPressureReactor
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autoclass:: IdealGasConstPressureReactor(contents=None, *, name=None, energy='on')
@@ -82,6 +90,22 @@ ExtensibleIdealGasConstPressureReactor
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autoclass:: ExtensibleIdealGasConstPressureReactor(contents=None, *, name=None, energy='on')

+ExtensibleMoleReactor
+^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleMoleReactor(contents=None, *, name=None, energy='on')
+
+ExtensibleIdealGasMoleReactor
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleIdealGasMoleReactor(contents=None, *, name=None, energy='on')
+
+ExtensibleConstPressureMoleReactor
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleConstPressureMoleReactor(contents=None, *, name=None, energy='on')
+
+ExtensibleIdealGasConstPressureMoleReactor
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleIdealGasConstPressureMoleReactor(contents=None, *, name=None, energy='on')
+
 Walls
 -----

--- a/include/cantera/zeroD/ConstPressureMoleReactor.h
+++ b/include/cantera/zeroD/ConstPressureMoleReactor.h
@@ -0,0 +1,45 @@
+//! @file ConstPressureMoleReactor.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_CONSTPRESSMOLE_REACTOR_H
+#define CT_CONSTPRESSMOLE_REACTOR_H
+
+#include "cantera/zeroD/MoleReactor.h"
+
+namespace Cantera
+{
+
+/*!
+ * ConstPressureMoleReactor is a class for constant-pressure reactors
+ * which use a state of moles.
+ */
+class ConstPressureMoleReactor : public MoleReactor
+{
+public:
+    ConstPressureMoleReactor() {}
+
+    virtual std::string type() const {
+        return "ConstPressureMoleReactor";
+    };
+
+    virtual size_t componentIndex(const std::string& nm) const;
+
+    virtual std::string componentName(size_t k);
+
+    virtual void getState(double* y);
+
+    virtual void initialize(double t0 = 0.0);
+
+    virtual void eval(double t, double* LHS, double* RHS);
+
+    virtual void updateState(double* y);
+
+protected:
+    const size_t m_sidx = 1;
+};
+
+}
+
+#endif
--- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
+++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
@@ -6,7 +6,7 @@
 #ifndef CT_IDEALGASCONSTPRESSMOLE_REACTOR_H
 #define CT_IDEALGASCONSTPRESSMOLE_REACTOR_H

-#include "cantera/zeroD/MoleReactor.h"
+#include "cantera/zeroD/ConstPressureMoleReactor.h"

 namespace Cantera
 {
@@ -15,7 +15,7 @@ namespace Cantera
 * IdealGasConstPressureMoleReactor is a class for ideal gas constant-pressure reactors
 * which use a state of moles.
 */
-class IdealGasConstPressureMoleReactor : public MoleReactor
+class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
 {
 public:
    IdealGasConstPressureMoleReactor() {}
@@ -47,8 +47,6 @@ public:

 protected:
    vector_fp m_hk; //!< Species molar enthalpies
-
-    const size_t m_sidx = 1;
 };

 }
--- a/include/cantera/zeroD/MoleReactor.h
+++ b/include/cantera/zeroD/MoleReactor.h
@@ -13,8 +13,7 @@ namespace Cantera

 /*!
 * MoleReactor is meant to serve the same purpose as the reactor class but with a state
- * vector composed of moles. It is currently not functional and only serves as a base
- * class for other reactors.
+ * vector composed of moles. It also serves as the base class for other mole reactors.
 */
 class MoleReactor : public Reactor
 {
@@ -27,19 +26,25 @@ public:

    virtual void initialize(double t0 = 0.0);

-    virtual void eval(double t, double* LHS, double* RHS) {
-        throw NotImplementedError("MoleReactor::eval()");
-    }
+    virtual void getState(double* y);

-    size_t componentIndex(const std::string& nm) const {
-        throw NotImplementedError("MoleReactor::componentIndex");
-    }
+    virtual void updateState(double* y);

-    std::string componentName(size_t k) {
-        throw NotImplementedError("MoleReactor::componentName");
-    }
+    virtual void eval(double t, double* LHS, double* RHS);
+
+    size_t componentIndex(const std::string& nm) const;
+
+    std::string componentName(size_t k);

 protected:
+    //! Get moles of the system from mass fractions stored by thermo object
+    //! @param y vector for moles to be put into
+    virtual void getMoles(double* y);
+
+    //! Set internal mass variable based on moles given
+    //! @param y vector of moles of the system
+    virtual void setMassFromMoles(double* y);
+
    virtual void evalSurfaces(double* LHS, double* RHS, double* sdot);

    virtual void updateSurfaceState(double* y);
--- a/interfaces/cython/cantera/reactor.pxd
+++ b/interfaces/cython/cantera/reactor.pxd
@@ -193,12 +193,18 @@ cdef class Reactor(ReactorBase):
    cdef CxxReactor* reactor
    cdef object _kinetics

+cdef class MoleReactor(Reactor):
+    pass
+
 cdef class Reservoir(ReactorBase):
    pass

 cdef class ConstPressureReactor(Reactor):
    pass

+cdef class ConstPressureMoleReactor(Reactor):
+    pass
+
 cdef class IdealGasReactor(Reactor):
    pass

--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -369,6 +369,14 @@ cdef class Reactor(ReactorBase):
            limit = -1.
        self.reactor.setAdvanceLimit(stringify(name), limit)

+cdef class MoleReactor(Reactor):
+    """
+    A homogeneous zero-dimensional reactor with a mole based state vector. By default,
+    they are closed (no inlets or outlets), have fixed volume, and have adiabatic,
+    chemically-inert walls. These properties may all be changed by adding
+    appropriate components such as `Wall`, `MassFlowController` and `Valve`.
+    """
+    reactor_type = "MoleReactor"

 cdef class Reservoir(ReactorBase):
    """
@@ -386,6 +394,13 @@ cdef class ConstPressureReactor(Reactor):
    """
    reactor_type = "ConstPressureReactor"

+cdef class ConstPressureMoleReactor(Reactor):
+    """A homogeneous, constant pressure, zero-dimensional reactor with a mole based
+    state vector. The volume of the reactor changes as a function of time in order to
+    keep the pressure constant.
+    """
+    reactor_type = "ConstPressureMoleReactor"
+

 cdef class IdealGasReactor(Reactor):
    """ A constant volume, zero-dimensional reactor for ideal gas mixtures. """
@@ -613,6 +628,38 @@ cdef class ExtensibleIdealGasConstPressureReactor(ExtensibleReactor):
    reactor_type = "ExtensibleIdealGasConstPressureReactor"


+cdef class ExtensibleMoleReactor(ExtensibleReactor):
+    """
+    A variant of `ExtensibleReactor` where the base behavior corresponds to the
+    `MoleReactor` class.
+    """
+    reactor_type = "ExtensibleMoleReactor"
+
+
+cdef class ExtensibleIdealGasMoleReactor(ExtensibleReactor):
+    """
+    A variant of `ExtensibleReactor` where the base behavior corresponds to the
+    `IdealGasMoleReactor` class.
+    """
+    reactor_type = "ExtensibleIdealGasMoleReactor"
+
+
+cdef class ExtensibleConstPressureMoleReactor(ExtensibleReactor):
+    """
+    A variant of `ExtensibleReactor` where the base behavior corresponds to the
+    `ConstPressureMoleReactor` class.
+    """
+    reactor_type = "ExtensibleConstPressureMoleReactor"
+
+
+cdef class ExtensibleIdealGasConstPressureMoleReactor(ExtensibleReactor):
+    """
+    A variant of `ExtensibleReactor` where the base behavior corresponds to the
+    `IdealGasConstPressureMoleReactor` class.
+    """
+    reactor_type = "ExtensibleIdealGasConstPressureMoleReactor"
+
+
 cdef class ReactorSurface:
    """
    Represents a surface in contact with the contents of a reactor.
--- a/src/zeroD/ConstPressureMoleReactor.cpp
+++ b/src/zeroD/ConstPressureMoleReactor.cpp
@@ -0,0 +1,146 @@
+//! @file ConstPressureMoleReactor.cpp A constant pressure
+//! zero-dimensional reactor with moles as the state
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/zeroD/Wall.h"
+#include "cantera/zeroD/FlowDevice.h"
+#include "cantera/zeroD/ConstPressureMoleReactor.h"
+#include "cantera/base/utilities.h"
+#include "cantera/thermo/SurfPhase.h"
+#include "cantera/kinetics/Kinetics.h"
+
+using namespace std;
+
+namespace Cantera
+{
+
+void ConstPressureMoleReactor::getState(double* y)
+{
+    if (m_thermo == 0) {
+        throw CanteraError("ConstPressureMoleReactor::getState",
+                           "Error: reactor is empty.");
+    }
+    m_thermo->restoreState(m_state);
+    // set mass to be used in getMoles function
+    m_mass = m_thermo->density() * m_vol;
+    // set the first array element to enthalpy
+    y[0] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol;
+    // get moles of species in remaining state
+    getMoles(y + m_sidx);
+    // set the remaining components to the surface species moles on the walls
+    getSurfaceInitialConditions(y+m_nsp+m_sidx);
+}
+
+void ConstPressureMoleReactor::initialize(double t0)
+{
+    MoleReactor::initialize(t0);
+    m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor
+}
+
+void ConstPressureMoleReactor::updateState(double* y)
+{
+    // the components of y are: [0] the enthalpy, [1...K+1) are the
+    // moles of each species, and [K+1...] are the moles of surface
+    // species on each wall.
+    setMassFromMoles(y + m_sidx);
+    m_thermo->setMolesNoTruncate(y + m_sidx);
+    if (m_energy) {
+        m_thermo->setState_HP(y[0] / m_mass, m_pressure);
+    } else {
+        m_thermo->setPressure(m_pressure);
+    }
+    m_vol = m_mass / m_thermo->density();
+    updateConnected(false);
+    updateSurfaceState(y + m_nsp + m_sidx);
+}
+
+void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS)
+{
+    double* dndt = RHS + m_sidx; // kmol per s
+
+    evalWalls(time);
+
+    m_thermo->restoreState(m_state);
+
+    const vector_fp& imw = m_thermo->inverseMolecularWeights();
+
+    if (m_chem) {
+        m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
+    }
+
+    // evaluate reactor surfaces
+    evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
+
+    // external heat transfer
+    double dHdt = m_Qdot;
+
+    for (size_t n = 0; n < m_nsp; n++) {
+        // production in gas phase and from surfaces
+        dndt[n] = m_wdot[n] * m_vol + m_sdot[n];
+    }
+
+    // add terms for outlets
+    for (auto outlet : m_outlet) {
+        // determine enthalpy contribution
+        dHdt -= outlet->massFlowRate() * m_enthalpy;
+        // flow of species into system and dilution by other species
+        for (size_t n = 0; n < m_nsp; n++) {
+            dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
+        }
+    }
+
+    // add terms for inlets
+    for (auto inlet : m_inlet) {
+        // enthalpy contribution from inlets
+        dHdt += inlet->enthalpy_mass() * inlet->massFlowRate();
+        // flow of species into system and dilution by other species
+        for (size_t n = 0; n < m_nsp; n++) {
+            dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
+        }
+    }
+
+    if (m_energy) {
+        RHS[0] = dHdt;
+    } else {
+        RHS[0] = 0.0;
+    }
+}
+
+size_t ConstPressureMoleReactor::componentIndex(const string& nm) const
+{
+    size_t k = speciesIndex(nm);
+    if (k != npos) {
+        return k + m_sidx;
+    } else if (nm == "enthalpy") {
+        return 0;
+    } else {
+        return npos;
+    }
+}
+
+std::string ConstPressureMoleReactor::componentName(size_t k) {
+    if (k == 0) {
+        return "enthalpy";
+    } else if (k >= m_sidx && k < neq()) {
+        k -= m_sidx;
+        if (k < m_thermo->nSpecies()) {
+            return m_thermo->speciesName(k);
+        } else {
+            k -= m_thermo->nSpecies();
+        }
+        for (auto& S : m_surfaces) {
+            ThermoPhase* th = S->thermo();
+            if (k < th->nSpecies()) {
+                return th->speciesName(k);
+            } else {
+                k -= th->nSpecies();
+            }
+        }
+    }
+    throw CanteraError("ConstPressureMoleReactor::componentName",
+                       "Index is out of bounds.");
+}
+
+}
--- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp
+++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp
@@ -25,7 +25,7 @@ void IdealGasConstPressureMoleReactor::setThermoMgr(ThermoPhase& thermo)
        throw CanteraError("IdealGasConstPressureMoleReactor::setThermoMgr",
                           "Incompatible phase type provided");
    }
-    MoleReactor::setThermoMgr(thermo);
+    ConstPressureMoleReactor::setThermoMgr(thermo);
 }

 void IdealGasConstPressureMoleReactor::getState(double* y)
@@ -39,21 +39,15 @@ void IdealGasConstPressureMoleReactor::getState(double* y)
    m_mass = m_thermo->density() * m_vol;
    // set the first component to the temperature
    y[0] = m_thermo->temperature();
-    // Use inverse molecular weights
-    const double* Y = m_thermo->massFractions();
-    const vector_fp& imw = m_thermo->inverseMolecularWeights();
-    double *ys = y + m_sidx;
-    for (size_t i = 0; i < m_nsp; i++) {
-        ys[i] = m_mass * imw[i] * Y[i];
-    }
+    // get moles of species in remaining state
+    getMoles(y + m_sidx);
    // set the remaining components to the surface species moles on the walls
    getSurfaceInitialConditions(y + m_nsp + m_sidx);
 }

 void IdealGasConstPressureMoleReactor::initialize(double t0)
 {
-    MoleReactor::initialize(t0);
-    m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor
+    ConstPressureMoleReactor::initialize(t0);
    m_hk.resize(m_nsp, 0.0);
 }

@@ -62,14 +56,9 @@ void IdealGasConstPressureMoleReactor::updateState(double* y)
    // the components of y are: [0] the temperature, [1...K+1) are the
    // moles of each species, and [K+1...] are the moles of surface
    // species on each wall.
+    setMassFromMoles(y + m_sidx);
    m_thermo->setMolesNoTruncate(y + m_sidx);
    m_thermo->setState_TP(y[0], m_pressure);
-    // calculate mass from moles
-    const vector_fp& mw = m_thermo->molecularWeights();
-    m_mass = 0;
-    for (size_t i = 0; i < m_nsp; i++) {
-        m_mass += y[i + m_sidx] * mw[i];
-    }
    m_vol = m_mass / m_thermo->density();
    updateConnected(false);
    updateSurfaceState(y + m_nsp + m_sidx);
--- a/src/zeroD/IdealGasMoleReactor.cpp
+++ b/src/zeroD/IdealGasMoleReactor.cpp
@@ -45,13 +45,8 @@ void IdealGasMoleReactor::getState(double* y)
    // set the second component to the volume
    y[1] = m_vol;

-    // use inverse molecular weights
-    const double* Y = m_thermo->massFractions();
-    const vector_fp& imw = m_thermo->inverseMolecularWeights();
-    double* ys = y + m_sidx;
-    for (size_t i = 0; i < m_nsp; i++) {
-        ys[i] = m_mass * imw[i] * Y[i];
-    }
+    // get moles of species in remaining state
+    getMoles(y + m_sidx);
    // set the remaining components to the surface species moles on
    // the walls
    getSurfaceInitialConditions(y + m_nsp + m_sidx);
@@ -75,25 +70,9 @@ std::string IdealGasMoleReactor::componentName(size_t k)
 {
    if (k == 0) {
        return "temperature";
-    } else if (k == 1) {
-        return "volume";
-    } else if (k >= m_sidx && k < neq()) {
-        k -= m_sidx;
-        if (k < m_thermo->nSpecies()) {
-            return m_thermo->speciesName(k);
-        } else {
-            k -= m_thermo->nSpecies();
-        }
-        for (auto& S : m_surfaces) {
-            ThermoPhase* th = S->thermo();
-            if (k < th->nSpecies()) {
-                return th->speciesName(k);
-            } else {
-                k -= th->nSpecies();
-            }
-        }
+    } else {
+        return MoleReactor::componentName(k);
    }
-    throw IndexError("IdealGasMoleReactor::componentName", "components", k, neq() - 1);
 }

 void IdealGasMoleReactor::initialize(double t0)
@@ -107,12 +86,7 @@ void IdealGasMoleReactor::updateState(double* y)
    // the components of y are: [0] the temperature, [1] the volume, [2...K+1) are the
    // moles of each species, and [K+1...] are the moles of surface
    // species on each wall.
-    const vector_fp& mw = m_thermo->molecularWeights();
-    // calculate mass from moles
-    m_mass = 0;
-    for (size_t i = 0; i < m_nv - m_sidx; i++) {
-        m_mass += y[i + m_sidx] * mw[i];
-    }
+    setMassFromMoles(y + m_sidx);
    m_vol = y[1];
    // set state
    m_thermo->setMolesNoTruncate(y + m_sidx);
--- a/src/zeroD/MoleReactor.cpp
+++ b/src/zeroD/MoleReactor.cpp
@@ -5,11 +5,15 @@
 // at https://cantera.org/license.txt for license and copyright information.

 #include "cantera/zeroD/MoleReactor.h"
-#include "cantera/thermo/SurfPhase.h"
+#include "cantera/zeroD/FlowDevice.h"
 #include "cantera/zeroD/ReactorSurface.h"
+#include "cantera/thermo/SurfPhase.h"
 #include "cantera/kinetics/Kinetics.h"
+#include "cantera/base/utilities.h"
+#include <boost/math/tools/roots.hpp>

 using namespace std;
+namespace bmt = boost::math::tools;

 namespace Cantera
 {
@@ -80,4 +84,186 @@ void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot)
    }
 }

+void MoleReactor::getMoles(double* y)
+{
+    // Use inverse molecular weights to convert to moles
+    const double* Y = m_thermo->massFractions();
+    const vector_fp& imw = m_thermo->inverseMolecularWeights();
+    for (size_t i = 0; i < m_nsp; i++) {
+        y[i] = m_mass * imw[i] * Y[i];
+    }
+}
+
+void MoleReactor::setMassFromMoles(double* y)
+{
+    const vector_fp& mw = m_thermo->molecularWeights();
+    // calculate mass from moles
+    m_mass = 0;
+    for (size_t i = 0; i < m_nsp; i++) {
+        m_mass += y[i] * mw[i];
+    }
+}
+
+void MoleReactor::getState(double* y)
+{
+    if (m_thermo == 0) {
+        throw CanteraError("MoleReactor::getState",
+                           "Error: reactor is empty.");
+    }
+    m_thermo->restoreState(m_state);
+    // set the first component to the internal energy
+    m_mass = m_thermo->density() * m_vol;
+    y[0] = m_thermo->intEnergy_mass() * m_mass;
+    // set the second component to the total volume
+    y[1] = m_vol;
+    // set components y+2 ... y+K+2 to the moles of each species
+    getMoles(y + m_sidx);
+    // set the remaining components to the surface species
+    // moles on walls
+    getSurfaceInitialConditions(y+m_nsp+m_sidx);
+}
+
+
+void MoleReactor::updateState(double* y)
+{
+    // The components of y are [0] total internal energy, [1] the total volume, and
+    // [2...K+3] are the moles of each species, and [K+3...] are the moles
+    // of surface species on each wall.
+    setMassFromMoles(y + m_sidx);
+    m_vol = y[1];
+    m_thermo->setMolesNoTruncate(y + m_sidx);
+    if (m_energy) {
+        double U = y[0];
+        // Residual function: error in internal energy as a function of T
+        auto u_err = [this, U](double T) {
+            m_thermo->setState_TR(T, m_mass / m_vol);
+            return m_thermo->intEnergy_mass() * m_mass - U;
+        };
+        double T = m_thermo->temperature();
+        boost::uintmax_t maxiter = 100;
+        std::pair<double, double> TT;
+        try {
+            TT = bmt::bracket_and_solve_root(
+                u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
+        } catch (std::exception&) {
+            // Try full-range bisection if bracketing fails (for example, near
+            // temperature limits for the phase's equation of state)
+            try {
+                TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
+                    bmt::eps_tolerance<double>(48), maxiter);
+            } catch (std::exception& err2) {
+                // Set m_thermo back to a reasonable state if root finding fails
+                m_thermo->setState_TR(T, m_mass / m_vol);
+                throw CanteraError("MoleReactor::updateState",
+                    "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
+            }
+        }
+        if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
+            throw CanteraError("MoleReactor::updateState", "root finding failed");
+        }
+        m_thermo->setState_TR(TT.second, m_mass / m_vol);
+    } else {
+        m_thermo->setDensity(m_mass / m_vol);
+    }
+    updateConnected(true);
+    updateSurfaceState(y + m_nsp + m_sidx);
+}
+
+void MoleReactor::eval(double time, double* LHS, double* RHS)
+{
+    double* dndt = RHS + m_sidx; // moles per time
+
+    evalWalls(time);
+    m_thermo->restoreState(m_state);
+
+    evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
+    // inverse molecular weights for conversion
+    const vector_fp& imw = m_thermo->inverseMolecularWeights();
+    // volume equation
+    RHS[1] = m_vdot;
+
+    if (m_chem) {
+        m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
+    }
+
+    // Energy equation.
+    // \f[
+    //     \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
+    // \f]
+    if (m_energy) {
+        RHS[0] = - m_thermo->pressure() * m_vdot + m_Qdot;
+    } else {
+        RHS[0] = 0.0;
+    }
+
+    for (size_t k = 0; k < m_nsp; k++) {
+        // production in gas phase and from surfaces
+        dndt[k] = m_wdot[k] * m_vol + m_sdot[k];
+    }
+
+    // add terms for outlets
+    for (auto outlet : m_outlet) {
+        // flow of species into system and dilution by other species
+        for (size_t n = 0; n < m_nsp; n++) {
+            dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
+        }
+        // energy update based on mass flow
+        double mdot = outlet->massFlowRate();
+        if (m_energy) {
+            RHS[0] -= mdot * m_enthalpy;
+        }
+    }
+
+    // add terms for inlets
+    for (auto inlet : m_inlet) {
+        double mdot = inlet->massFlowRate();
+        for (size_t n = 0; n < m_nsp; n++) {
+            // flow of species into system and dilution by other species
+            dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
+        }
+        if (m_energy) {
+            RHS[0] += mdot * inlet->enthalpy_mass();
+        }
+    }
+}
+
+
+size_t MoleReactor::componentIndex(const string& nm) const
+{
+    size_t k = speciesIndex(nm);
+    if (k != npos) {
+        return k + m_sidx;
+    } else if (nm == "int_energy") {
+        return 0;
+    } else if (nm == "volume") {
+        return 1;
+    } else {
+        return npos;
+    }
+}
+
+std::string MoleReactor::componentName(size_t k) {
+    if (k == 0) {
+        return "int_energy";
+    } else if (k == 1) {
+        return "volume";
+    } else if (k >= m_sidx && k < neq()) {
+        k -= m_sidx;
+        if (k < m_thermo->nSpecies()) {
+            return m_thermo->speciesName(k);
+        } else {
+            k -= m_thermo->nSpecies();
+        }
+        for (auto& S : m_surfaces) {
+            ThermoPhase* th = S->thermo();
+            if (k < th->nSpecies()) {
+                return th->speciesName(k);
+            } else {
+                k -= th->nSpecies();
+            }
+        }
+    }
+    throw CanteraError("MoleReactor::componentName", "Index is out of bounds.");
+}
+
 }
--- a/src/zeroD/ReactorFactory.cpp
+++ b/src/zeroD/ReactorFactory.cpp
@@ -6,8 +6,10 @@
 #include "cantera/zeroD/ReactorFactory.h"
 #include "cantera/zeroD/Reservoir.h"
 #include "cantera/zeroD/Reactor.h"
+#include "cantera/zeroD/MoleReactor.h"
 #include "cantera/zeroD/FlowReactor.h"
 #include "cantera/zeroD/ConstPressureReactor.h"
+#include "cantera/zeroD/ConstPressureMoleReactor.h"
 #include "cantera/zeroD/IdealGasReactor.h"
 #include "cantera/zeroD/IdealGasMoleReactor.h"
 #include "cantera/zeroD/IdealGasConstPressureReactor.h"
@@ -38,8 +40,19 @@ ReactorFactory::ReactorFactory()
        []() { return new ReactorDelegator<ConstPressureReactor>(); });
    reg("ExtensibleIdealGasConstPressureReactor",
        []() { return new ReactorDelegator<IdealGasConstPressureReactor>(); });
-    reg("IdealGasConstPressureMoleReactor", []() { return new IdealGasConstPressureMoleReactor(); });
+    reg("ExtensibleMoleReactor",
+        []() { return new ReactorDelegator<MoleReactor>(); });
+    reg("ExtensibleConstPressureMoleReactor",
+        []() { return new ReactorDelegator<ConstPressureMoleReactor>(); });
+    reg("ExtensibleIdealGasMoleReactor",
+        []() { return new ReactorDelegator<IdealGasMoleReactor>(); });
+    reg("ExtensibleIdealGasConstPressureMoleReactor",
+        []() { return new ReactorDelegator<IdealGasConstPressureMoleReactor>(); });
+    reg("IdealGasConstPressureMoleReactor", []() { return new
+        IdealGasConstPressureMoleReactor(); });
    reg("IdealGasMoleReactor", []() { return new IdealGasMoleReactor(); });
+    reg("ConstPressureMoleReactor", []() { return new ConstPressureMoleReactor(); });
+    reg("MoleReactor", []() { return new MoleReactor(); });
 }

 ReactorBase* ReactorFactory::newReactor(const std::string& reactorType)
--- a/test/python/test_reactor.py
+++ b/test/python/test_reactor.py
@@ -89,6 +89,7 @@ class TestReactor(utilities.CanteraTest):

    def test_component_names(self):
        self.make_reactors(n_reactors=2)
+        self.net.initialize()
        N = self.net.n_vars // 2
        for i in range(N):
            self.assertEqual(self.r1.component_index(self.r1.component_name(i)), i)
@@ -258,10 +259,9 @@ class TestReactor(utilities.CanteraTest):

        n_baseline = integrate(1e-10, 1e-20)
        n_rtol = integrate(5e-7, 1e-20)
-        n_atol = integrate(1e-10, 1e-6)
-
-        self.assertTrue(n_baseline > n_rtol)
-        self.assertTrue(n_baseline > n_atol)
+        n_atol = integrate(1e-10, 1e-5)
+        assert n_baseline > n_rtol
+        assert n_baseline > n_atol

    def test_advance_limits(self):
        P0 = 10 * ct.one_atm
@@ -823,6 +823,56 @@ class TestReactor(utilities.CanteraTest):
        with self.assertRaises(TypeError):
            r1 = self.reactorClass(foobar=3.14)

+class TestMoleReactor(TestReactor):
+    reactorClass = ct.MoleReactor
+
+    def test_mole_reactor_surface_chem(self):
+        model = "ptcombust.yaml"
+        # initial conditions
+        T0 = 1500
+        P0 = ct.one_atm
+        equiv_ratio = 1
+        surf_area = 0.527
+        fuel = "CH4"
+        air = "O2:1.0, N2:3.773"
+        # reactor 1
+        gas1 = ct.Solution(model, transport_model=None)
+        gas1.TP = T0, P0
+        gas1.set_equivalence_ratio(equiv_ratio, fuel, air)
+        r1 = self.reactorClass(gas1)
+        # comparison reactor
+        gas2 = ct.Solution(model, transport_model=None)
+        gas2.TP = T0, P0
+        gas2.set_equivalence_ratio(equiv_ratio, fuel, air)
+        if "ConstPressure" in r1.type:
+            r2 = ct.ConstPressureReactor(gas2)
+        else:
+            r2 = ct.Reactor(gas2)
+        # surf 1
+        surf1 = ct.Interface(model, "Pt_surf", [gas1])
+        surf1.TP = T0, P0
+        surf1.coverages = {"PT(S)":1}
+        rsurf1 = ct.ReactorSurface(surf1, r1, A=surf_area)
+        # surf 2
+        surf2 = ct.Interface(model, "Pt_surf", [gas2])
+        surf2.TP = T0, P0
+        surf2.coverages = {"PT(S)":1}
+        rsurf2 = ct.ReactorSurface(surf2, r2, A=surf_area)
+        # reactor network setup
+        net1 = ct.ReactorNet([r1,])
+        net2 = ct.ReactorNet([r2,])
+        net1.rtol = net2.rtol = 1e-9
+        net1.atol = net2.atol = 1e-18
+        # steady state occurs at ~0.002 seconds
+        for i in np.linspace(0, 0.0025, 50)[1:]:
+            net1.advance(i)
+            net2.advance(i)
+            self.assertArrayNear(r1.thermo.Y, r2.thermo.Y, rtol=5e-4, atol=1e-6)
+            self.assertNear(r1.T, r2.T, rtol=5e-5)
+            self.assertNear(r1.thermo.P, r2.thermo.P, rtol=1e-6)
+            self.assertArrayNear(rsurf1.coverages, rsurf2.coverages, rtol=1e-4,
+                atol=1e-8)
+

 class TestIdealGasReactor(TestReactor):
    reactorClass = ct.IdealGasReactor
@@ -1070,12 +1120,21 @@ class TestConstPressureReactor(utilities.CanteraTest):
        self.net1.rtol = self.net2.rtol = 1e-9
        self.integrate(surf=True)

+class TestConstPressureMoleReactor(TestConstPressureReactor):
+    """
+    The constant pressure reactor should give the same results as
+    as a regular "Reactor" with a wall with a very high expansion rate
+    coefficient.
+    """
+    reactorClass = ct.ConstPressureMoleReactor
+    test_mole_reactor_surface_chem = TestMoleReactor.test_mole_reactor_surface_chem
+

 class TestIdealGasConstPressureReactor(TestConstPressureReactor):
    reactorClass = ct.IdealGasConstPressureReactor


-class TestIdealGasConstPressureMoleReactor(TestIdealGasConstPressureReactor):
+class TestIdealGasConstPressureMoleReactor(TestConstPressureMoleReactor):
    reactorClass = ct.IdealGasConstPressureMoleReactor

    def create_reactors(self, **kwargs):
@@ -1097,7 +1156,7 @@ class TestIdealGasConstPressureMoleReactor(TestIdealGasConstPressureReactor):
        with pytest.raises(ct.CanteraError):
            super().test_component_index()

-class TestIdealGasMoleReactor(TestReactor):
+class TestIdealGasMoleReactor(TestMoleReactor):
    reactorClass = ct.IdealGasMoleReactor

    def test_adaptive_precon_integration(self):