[Reactor] Set m_nv in constructor

speth · speth · commit f7534b88903b · 2026-01-08T21:46:33.000-05:00
diff --git a/include/cantera/zeroD/ConstPressureMoleReactor.h b/include/cantera/zeroD/ConstPressureMoleReactor.h
@@ -20,16 +20,15 @@ namespace Cantera
 class ConstPressureMoleReactor : public MoleReactor
 {
 public:
-    using MoleReactor::MoleReactor; // inherit constructors
+    ConstPressureMoleReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    ConstPressureMoleReactor(shared_ptr<Solution> sol, bool clone,
+                             const string& name="(none)");
 
     string type() const override {
         return "ConstPressureMoleReactor";
     };
 
     void getState(double* y) override;
-
-    void initialize(double t0=0.0) override;
-
     void eval(double t, double* LHS, double* RHS) override;
 
     vector<size_t> steadyConstraints() const override {
diff --git a/include/cantera/zeroD/ConstPressureReactor.h b/include/cantera/zeroD/ConstPressureReactor.h
@@ -23,15 +23,15 @@ namespace Cantera
 class ConstPressureReactor : public Reactor
 {
 public:
-    using Reactor::Reactor; // inherit constructors
+    ConstPressureReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    ConstPressureReactor(shared_ptr<Solution> sol, bool clone,
+                         const string& name="(none)");
 
     string type() const override {
         return "ConstPressureReactor";
     }
 
     void getState(double* y) override;
-
-    void initialize(double t0=0.0) override;
     void eval(double t, double* LHS, double* RHS) override;
     vector<size_t> steadyConstraints() const override;
 
diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h
@@ -16,7 +16,8 @@ namespace Cantera
 class FlowReactor : public IdealGasReactor
 {
 public:
-    using IdealGasReactor::IdealGasReactor; // inherit constructors
+    FlowReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    FlowReactor(shared_ptr<Solution> sol, bool clone, const string& name="(none)");
 
     string type() const override {
         return "FlowReactor";
@@ -36,7 +37,6 @@ class FlowReactor : public IdealGasReactor
     }
 
     void getStateDae(double* y, double* ydot) override;
-    void initialize(double t0=0.0) override;
     void updateState(double* y) override;
 
     //! Not implemented; FlowReactor implements evalDae() instead.
diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h
@@ -20,20 +20,16 @@ namespace Cantera
 class MoleReactor : public Reactor
 {
 public:
-    using Reactor::Reactor; // inherit constructors
+    MoleReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    MoleReactor(shared_ptr<Solution> sol, bool clone, const string& name="(none)");
 
     string type() const override {
         return "MoleReactor";
     }
 
-    void initialize(double t0=0.0) override;
-
     void getState(double* y) override;
-
     void updateState(double* y) override;
-
     void eval(double t, double* LHS, double* RHS) override;
-
     size_t componentIndex(const string& nm) const override;
     string componentName(size_t k) override;
     double upperBound(size_t k) const override;
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
@@ -501,6 +501,10 @@ cdef class ExtensibleReactor(Reactor):
     A base class for a reactor with delegated methods where the base
     functionality corresponds to the `Reactor` class.
 
+    The ``__init__`` method of the derived class should allocate and size any
+    internal variables and set the total number of state variables associated with this
+    reactor, `n_vars` (if it is different from the base class).
+
     The following methods of the C++ :ct:`Reactor` class can be modified by a
     Python class which inherits from this class. For each method, the name below
     should be prefixed with ``before_``, ``after_``, or ``replace_``, indicating
@@ -515,9 +519,8 @@ cdef class ExtensibleReactor(Reactor):
     from the supplied method and the base class method.
 
     ``initialize(self, t0: double) -> None``
-        Responsible for allocating and setting the sizes of any internal
-        variables, initializing attached walls, and setting the total number of
-        state variables associated with this reactor, `n_vars`.
+        Responsible for initialization that can only be performed after connecting the
+        elements of the reactor network, such as initializing attached walls.
 
         Called once before the start of time integration.
 
diff --git a/samples/python/reactors/custom2.py b/samples/python/reactors/custom2.py
@@ -31,8 +31,7 @@ def __init__(self, *args, neighbor, **kwargs):
         self.k_wall = 1e-2  # proportionality constant, a_wall = k_wall * delta P
         self.neighbor = neighbor
 
-    def after_initialize(self, t0):
-        # The initialize function for the base Reactor class will have set
+        # The constructor for the base Reactor class will have set
         # n_vars to already include the volume, internal energy, mass, and mass
         # fractions of all the species. Increase this by one to account for
         # the added variable of the wall velocity.
diff --git a/samples/python/reactors/porous_media_burner.py b/samples/python/reactors/porous_media_burner.py
@@ -195,8 +195,6 @@ def __init__(self, *args, props, **kwargs):
         self.V = self.A * props.length             # reactor volume (m^3)
         self.molecular_weights = None
         self.species_offset = None
-
-    def after_initialize(self, t0):
         self.n_vars += 1                # additional equation for the solid temperature
         self.index_Ts = self.n_vars - 1
         self.species_offset = self.component_index(self.phase.species_name(0))
diff --git a/src/zeroD/ConstPressureMoleReactor.cpp b/src/zeroD/ConstPressureMoleReactor.cpp
@@ -15,6 +15,19 @@
 namespace Cantera
 {
 
+ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr<Solution> sol,
+                                                   const string& name)
+    : ConstPressureMoleReactor(sol, true, name)
+{
+}
+
+ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr<Solution> sol, bool clone,
+                                                   const string& name)
+    : MoleReactor(sol, clone, name)
+{
+    m_nv = 1 + m_nsp; // enthalpy and moles of each species
+}
+
 void ConstPressureMoleReactor::getState(double* y)
 {
     // set mass to be used in getMoles function
@@ -25,12 +38,6 @@ void ConstPressureMoleReactor::getState(double* y)
     getMoles(y + 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
diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp
@@ -14,6 +14,19 @@
 namespace Cantera
 {
 
+ConstPressureReactor::ConstPressureReactor(shared_ptr<Solution> sol,
+                                           const string& name)
+    : ConstPressureReactor(sol, true, name)
+{
+}
+
+ConstPressureReactor::ConstPressureReactor(shared_ptr<Solution> sol, bool clone,
+                                           const string& name)
+    : Reactor(sol, clone, name)
+{
+    m_nv = 2 + m_nsp; // mass, enthalpy, and mass fractions of each species
+}
+
 void ConstPressureReactor::getState(double* y)
 {
     // set the first component to the total mass
@@ -27,12 +40,6 @@ void ConstPressureReactor::getState(double* y)
 
 }
 
-void ConstPressureReactor::initialize(double t0)
-{
-    Reactor::initialize(t0);
-    m_nv -= 1; // Constant pressure reactor has one fewer state variable
-}
-
 void ConstPressureReactor::updateState(double* y)
 {
     // The components of y are [0] the total mass, [1] the total enthalpy,
diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp
@@ -17,6 +17,21 @@
 namespace Cantera
 {
 
+FlowReactor::FlowReactor(shared_ptr<Solution> sol, const string& name)
+    : FlowReactor(sol, true, name)
+{
+}
+
+FlowReactor::FlowReactor(shared_ptr<Solution> sol, bool clone, const string& name)
+    : IdealGasReactor(sol, clone, name)
+{
+    m_nv = 4 + m_nsp; // rho, u, P, T, and species mass fractions
+    m_rho = m_thermo->density();
+    // resize temporary arrays
+    m_wdot.resize(m_nsp);
+    m_hk.resize(m_nsp);
+}
+
 void FlowReactor::getStateDae(double* y, double* ydot)
 {
     m_thermo->getMassFractions(y+m_offset_Y);
@@ -132,20 +147,6 @@ void FlowReactor::getStateDae(double* y, double* ydot)
     solve(a, ydot, 1, 0);
 }
 
-void FlowReactor::initialize(double t0)
-{
-    Reactor::initialize(t0);
-    // initialize state
-    m_rho = m_thermo->density();
-    // resize temporary arrays
-    m_wdot.resize(m_nsp);
-    m_hk.resize(m_nsp);
-    // set number of variables to the number of non-species equations
-    // i.e., density, velocity, pressure and temperature
-    // plus the number of species in the gas phase
-    m_nv = m_offset_Y + m_nsp;
-}
-
 void FlowReactor::updateState(double* y)
 {
     // Set the mass fractions and density of the mixture.
diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp
@@ -105,10 +105,6 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH
 void IdealGasConstPressureMoleReactor::getJacobianElements(
     vector<Eigen::Triplet<double>>& trips)
 {
-    if (m_nv == 0) {
-        throw CanteraError("IdealGasConstPressureMoleReactor::getJacobianElements",
-                           "Reactor must be initialized first.");
-    }
     // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
     // d (dot(omega)) / d c_j, it is later transformed appropriately.
     Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp
@@ -171,10 +171,6 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS)
 
 void IdealGasMoleReactor::getJacobianElements(vector<Eigen::Triplet<double>>& trips)
 {
-    if (m_nv == 0) {
-        throw CanteraError("IdealGasMoleReactor::jacobian",
-                           "Reactor must be initialized first.");
-    }
     // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
     // d (dot(omega)) / d c_j, it is later transformed appropriately.
     Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp
@@ -18,10 +18,15 @@ namespace bmt = boost::math::tools;
 namespace Cantera
 {
 
-void MoleReactor::initialize(double t0)
+MoleReactor::MoleReactor(shared_ptr<Solution> sol, const string& name)
+    : MoleReactor(sol, true, name)
 {
-    Reactor::initialize(t0);
-    m_nv -= 1; // moles gives the state one fewer variables
+}
+
+MoleReactor::MoleReactor(shared_ptr<Solution> sol, bool clone, const string& name)
+    : Reactor(sol, clone, name)
+{
+    m_nv = 2 + m_nsp; // internal energy, volume, and moles of each species
 }
 
 void MoleReactor::getMoles(double* y)
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
@@ -34,6 +34,8 @@ Reactor::Reactor(shared_ptr<Solution> sol, bool clone, const string& name)
     setChemistryEnabled(m_kin->nReactions() > 0);
     m_vol = 1.0; // By default, the volume is set to 1.0 m^3.
     m_sdot.resize(m_nsp, 0.0);
+    m_wdot.resize(m_nsp, 0.0);
+    m_nv = 3 + m_nsp; // mass, volume, internal energy, and species mass fractions
 }
 
 void Reactor::setDerivativeSettings(AnyMap& settings)
@@ -63,15 +65,12 @@ void Reactor::initialize(double t0)
         throw CanteraError("Reactor::initialize", "Reactor contents not set"
                 " for reactor '" + m_name + "'.");
     }
-    m_wdot.resize(m_nsp, 0.0);
     updateConnected(true);
 
     for (size_t n = 0; n < m_wall.size(); n++) {
         WallBase* W = m_wall[n];
         W->initialize();
     }
-
-    m_nv = m_nsp + 3;
 }
 
 void Reactor::updateState(double* y)
@@ -214,10 +213,6 @@ vector<size_t> Reactor::steadyConstraints() const
 
 Eigen::SparseMatrix<double> Reactor::finiteDifferenceJacobian()
 {
-    if (m_nv == 0) {
-        throw CanteraError("Reactor::finiteDifferenceJacobian",
-                           "Reactor must be initialized first.");
-    }
     vector<Eigen::Triplet<double>> trips;
     Eigen::ArrayXd yCurrent(m_nv);
     getState(yCurrent.data());
@@ -429,18 +424,6 @@ bool Reactor::getAdvanceLimits(double *limits) const
 void Reactor::setAdvanceLimit(const string& nm, const double limit)
 {
     size_t k = componentIndex(nm);
-    if (m_nv == 0) {
-        if (m_net == 0) {
-            throw CanteraError("Reactor::setAdvanceLimit",
-                               "Cannot set limit on a reactor that is not "
-                               "assigned to a ReactorNet object.");
-        } else {
-            m_net->initialize();
-        }
-    } else if (k > m_nv) {
-        throw CanteraError("Reactor::setAdvanceLimit",
-                           "Index out of bounds.");
-    }
     m_advancelimits.resize(m_nv, -1.0);
     m_advancelimits[k] = limit;
 
diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py
@@ -3014,8 +3014,6 @@ def __init__(self, *args, neighbor, **kwargs):
                 self.v_wall = 0
                 self.k_wall = 1e-5
                 self.neighbor = neighbor
-
-            def after_initialize(self, t0):
                 self.n_vars += 1
                 self.i_wall = self.n_vars - 1
 

Original file line number	Diff line number	Diff line change
`@@ -16,7 +16,8 @@ namespace Cantera`
`16`	`16`	`class FlowReactor : public IdealGasReactor`
`17`	`17`	`{`
`18`	`18`	`public:`
`19`		`- using IdealGasReactor::IdealGasReactor; // inherit constructors`
	`19`	`+ FlowReactor(shared_ptr<Solution> sol, const string& name="(none)");`
	`20`	`+ FlowReactor(shared_ptr<Solution> sol, bool clone, const string& name="(none)");`
`20`	`21`
`21`	`22`	`string type() const override {`
`22`	`23`	`return "FlowReactor";`
`@@ -36,7 +37,6 @@ class FlowReactor : public IdealGasReactor`
`36`	`37`	`}`
`37`	`38`
`38`	`39`	`void getStateDae(double* y, double* ydot) override;`
`39`		`- void initialize(double t0=0.0) override;`
`40`	`40`	`void updateState(double* y) override;`
`41`	`41`
`42`	`42`	`//! Not implemented; FlowReactor implements evalDae() instead.`