From 2907cdefee751d56168144fa60d3dc0dc3013fbc Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 5 Jan 2026 15:42:08 -0500 Subject: [PATCH] [Reactor] Set m_nv in constructor --- .../cantera/zeroD/ConstPressureMoleReactor.h | 7 ++--- include/cantera/zeroD/ConstPressureReactor.h | 6 ++-- include/cantera/zeroD/FlowReactor.h | 4 +-- include/cantera/zeroD/MoleReactor.h | 8 ++--- interfaces/cython/cantera/reactor.pyx | 9 ++++-- samples/python/reactors/custom2.py | 3 +- .../python/reactors/porous_media_burner.py | 2 -- src/zeroD/ConstPressureMoleReactor.cpp | 19 ++++++++---- src/zeroD/ConstPressureReactor.cpp | 19 ++++++++---- src/zeroD/FlowReactor.cpp | 29 ++++++++++--------- .../IdealGasConstPressureMoleReactor.cpp | 4 --- src/zeroD/IdealGasMoleReactor.cpp | 4 --- src/zeroD/MoleReactor.cpp | 11 +++++-- src/zeroD/Reactor.cpp | 21 ++------------ test/python/test_reactor.py | 2 -- 15 files changed, 68 insertions(+), 80 deletions(-) diff --git a/include/cantera/zeroD/ConstPressureMoleReactor.h b/include/cantera/zeroD/ConstPressureMoleReactor.h index 5031e4889..0ea9c0ce2 100644 --- 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 sol, const string& name="(none)"); + ConstPressureMoleReactor(shared_ptr 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 steadyConstraints() const override { diff --git a/include/cantera/zeroD/ConstPressureReactor.h b/include/cantera/zeroD/ConstPressureReactor.h index a443cc6a4..9566289b7 100644 --- 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 sol, const string& name="(none)"); + ConstPressureReactor(shared_ptr 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 steadyConstraints() const override; diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h index ae6a610c0..ea2f4a481 100644 --- 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 sol, const string& name="(none)"); + FlowReactor(shared_ptr sol, bool clone, const string& name="(none)"); string type() const override { return "FlowReactor"; @@ -36,7 +37,6 @@ public: } 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 index 6db8f33bb..4a9719751 100644 --- 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 sol, const string& name="(none)"); + MoleReactor(shared_ptr 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 index 14a8624c2..45e7b5e63 100644 --- 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 index 1215139b8..7094ffe36 100644 --- a/samples/python/reactors/custom2.py +++ b/samples/python/reactors/custom2.py @@ -31,8 +31,7 @@ class InertialWallReactor(ct.ExtensibleIdealGasReactor): 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 index 35a9bdb05..13a7aec72 100644 --- a/samples/python/reactors/porous_media_burner.py +++ b/samples/python/reactors/porous_media_burner.py @@ -195,8 +195,6 @@ class PMReactor(ct.ExtensibleIdealGasConstPressureReactor): 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 index 1b5ca149f..134eb9263 100644 --- a/src/zeroD/ConstPressureMoleReactor.cpp +++ b/src/zeroD/ConstPressureMoleReactor.cpp @@ -15,6 +15,19 @@ namespace Cantera { +ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr sol, + const string& name) + : ConstPressureMoleReactor(sol, true, name) +{ +} + +ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr 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 index d76bf667f..0d2eedcf1 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -14,6 +14,19 @@ namespace Cantera { +ConstPressureReactor::ConstPressureReactor(shared_ptr sol, + const string& name) + : ConstPressureReactor(sol, true, name) +{ +} + +ConstPressureReactor::ConstPressureReactor(shared_ptr 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 index f3bb257f0..7667f9c8a 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/src/zeroD/FlowReactor.cpp @@ -17,6 +17,21 @@ namespace Cantera { +FlowReactor::FlowReactor(shared_ptr sol, const string& name) + : FlowReactor(sol, true, name) +{ +} + +FlowReactor::FlowReactor(shared_ptr 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 index b382d26a5..564ad8ce9 100644 --- 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>& 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 dnk_dnj = m_kin->netProductionRates_ddCi(); diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 2997c2025..9c832210c 100644 --- 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>& 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 dnk_dnj = m_kin->netProductionRates_ddCi(); diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp index f2b627706..32c717832 100644 --- 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 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 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 index 30a5074ea..54e6bf2e5 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -34,6 +34,8 @@ Reactor::Reactor(shared_ptr 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 Reactor::steadyConstraints() const Eigen::SparseMatrix Reactor::finiteDifferenceJacobian() { - if (m_nv == 0) { - throw CanteraError("Reactor::finiteDifferenceJacobian", - "Reactor must be initialized first."); - } vector> 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 index 62bd8c796..c1d67cf39 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -3021,8 +3021,6 @@ class TestExtensibleReactor: 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