diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 61ef46370..6489c4d51 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -151,6 +151,12 @@ protected: //! Reset the reaction rate multipliers virtual void resetSensitivity(double* params); + //! Return the index in the solution vector for this reactor of the species + //! named *nm*, in either the homogeneous phase or a surface phase, relative + //! to the start of the species terms. Used to implement componentIndex for + //! specific reactor implementations. + virtual size_t speciesIndex(const std::string& nm) const; + //! Pointer to the homogeneous Kinetics object that handles the reactions Kinetics* m_kin; diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index 9b7ab83df..e0cd2ebac 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -617,16 +617,16 @@ class TestConstPressureReactor(utilities.CanteraTest): mfc2 = ct.MassFlowController(env, self.r2, mdot=0.05) if add_surf: - interface1 = ct.Interface('diamond.xml', 'diamond_100', + self.interface1 = ct.Interface('diamond.xml', 'diamond_100', (self.gas1, solid)) - interface2 = ct.Interface('diamond.xml', 'diamond_100', + self.interface2 = ct.Interface('diamond.xml', 'diamond_100', (self.gas2, solid)) - C = np.zeros(interface1.n_species) + C = np.zeros(self.interface1.n_species) C[0] = 0.3 C[4] = 0.7 - self.w1.left.kinetics = interface1 - self.w2.left.kinetics = interface2 + self.w1.left.kinetics = self.interface1 + self.w2.left.kinetics = self.interface2 self.w1.left.coverages = C self.w2.left.coverages = C @@ -636,6 +636,19 @@ class TestConstPressureReactor(utilities.CanteraTest): self.net2.set_max_time_step(0.05) self.net2.max_err_test_fails = 10 + def test_component_index(self): + self.create_reactors(add_surf=True) + for (gas,net,iface,r) in ((self.gas1, self.net1, self.interface1, self.r1), + (self.gas2, self.net2, self.interface2, self.r2)): + net.step(1.0) + + N0 = net.n_vars - gas.n_species - iface.n_species + N1 = net.n_vars - iface.n_species + for i, name in enumerate(gas.species_names): + self.assertEqual(i + N0, r.component_index(name)) + for i, name in enumerate(iface.species_names): + self.assertEqual(i + N1, r.component_index(name)) + def integrate(self, surf=False): for t in np.arange(0.5, 50, 1.0): self.net1.advance(t) diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index e1d0d5135..3f50f169c 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -223,34 +223,16 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, size_t ConstPressureReactor::componentIndex(const string& nm) const { - if (nm == "m") { - return 0; - } - if (nm == "H") { - return 1; - } - // check for a gas species name - size_t k = m_thermo->speciesIndex(nm); + size_t k = speciesIndex(nm); if (k != npos) { return k + 2; + } else if (nm == "m" || nm == "mass") { + return 0; + } else if (nm == "H" || nm == "enthalpy") { + return 1; + } else { + return npos; } - - // check for a wall species - size_t walloffset = 0, kp = 0; - thermo_t* th; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_wall[m]->kinetics(m_lr[m])) { - kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex(); - th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp); - k = th->speciesIndex(nm); - if (k != npos) { - return k + 2 + m_nsp + walloffset; - } else { - walloffset += th->nSpecies(); - } - } - } - return npos; } } diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp index 20aed33a4..297541997 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/src/zeroD/FlowReactor.cpp @@ -118,16 +118,14 @@ void FlowReactor::evalEqs(doublereal time, doublereal* y, size_t FlowReactor::componentIndex(const string& nm) const { - if (nm == "X") { - return 0; - } - if (nm == "U") { - return 1; - } // check for a gas species name size_t k = m_thermo->speciesIndex(nm); if (k != npos) { return k + 2; + } else if (nm == "X" || nm == "distance") { + return 0; + } else if (nm == "U" || nm == "velocity") { + return 1; } else { return npos; } diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index e59739643..5852c81d8 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -231,12 +231,16 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y, size_t IdealGasConstPressureReactor::componentIndex(const string& nm) const { - if (nm == "T") { + size_t k = speciesIndex(nm); + if (k != npos) { + return k + 2; + } else if (nm == "m" || nm == "mass") { + return 0; + } else if (nm == "T" || nm == "temperature") { return 1; } else { - return ConstPressureReactor::componentIndex(nm); + return npos; } - } } diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 012798609..40d1b5b87 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -251,10 +251,17 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, size_t IdealGasReactor::componentIndex(const string& nm) const { - if (nm == "T") { + size_t k = speciesIndex(nm); + if (k != npos) { + return k + 3; + } else if (nm == "m" || nm == "mass") { + return 0; + } else if (nm == "V" || nm == "volume") { + return 1; + } else if (nm == "T" || nm == "temperature") { return 2; } else { - return Reactor::componentIndex(nm); + return npos; } } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 73dbd5d97..33f24afb8 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -318,22 +318,12 @@ std::vector > Reactor::getSensitivityOrder() const return order; } -size_t Reactor::componentIndex(const string& nm) const +size_t Reactor::speciesIndex(const string& nm) const { - if (nm == "m") { - return 0; - } - if (nm == "V") { - return 1; - } - if (nm == "U") { - return 2; - } - // check for a gas species name size_t k = m_thermo->speciesIndex(nm); if (k != npos) { - return k + 3; + return k; } // check for a wall species @@ -345,7 +335,7 @@ size_t Reactor::componentIndex(const string& nm) const th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp); k = th->speciesIndex(nm); if (k != npos) { - return k + 3 + m_nsp + walloffset; + return k + m_nsp + walloffset; } else { walloffset += th->nSpecies(); } @@ -354,6 +344,22 @@ size_t Reactor::componentIndex(const string& nm) const return npos; } +size_t Reactor::componentIndex(const string& nm) const +{ + size_t k = speciesIndex(nm); + if (k != npos) { + return k + 3; + } else if (nm == "m" || nm == "mass") { + return 0; + } else if (nm == "V" || nm == "volume") { + return 1; + } else if (nm == "U" || nm == "int_energy") { + return 2; + } else { + return npos; + } +} + void Reactor::applySensitivity(double* params) { if (!params) {