diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 43d6e476c..0420d7b44 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -35,7 +35,7 @@ public: bool doIrreversible = false); virtual bool addReaction(shared_ptr r); - virtual void init(); + virtual void resizeSpecies(); virtual void setMultiplier(size_t i, double f); virtual void invalidateCache(); diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index d50eda2a5..bfecbaf10 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -198,6 +198,7 @@ public: virtual void addPhase(thermo_t& thermo); virtual void init(); + virtual void resizeSpecies(); virtual bool addReaction(shared_ptr r); virtual void modifyReaction(size_t i, shared_ptr rNew); //! @} diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 419b9bf4a..5a7ef34ee 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -714,15 +714,21 @@ public: virtual void addPhase(thermo_t& thermo); /** - * Prepare the class for the addition of reactions. This method is called - * by importKinetics() after all phases have been added but before any - * reactions have been. The base class method does nothing, but derived - * classes may use this to perform any initialization (allocating arrays, - * etc.) that requires knowing the phases and species, but before any - * reactions are added. + * Prepare the class for the addition of reactions, after all phases have + * been added. This method is called automatically when the first reaction + * is added. It needs to be called directly only in the degenerate case + * where there are no reactions. The base class method does nothing, but + * derived classes may use this to perform any initialization (allocating + * arrays, etc.) that requires knowing the phases. */ virtual void init() {} + /** + * Resize arrays with sizes that depend on the total number of species. + * Automatically called before adding each Reaction and Phase. + */ + virtual void resizeSpecies(); + /** * Finish adding reactions and prepare for use. This method is called by * importKinetics() after all reactions have been entered into the diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index fdc4c9ce7..21e6f785b 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -137,9 +137,9 @@ void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew) m_rates.replace(i, rNew.rate); } -void BulkKinetics::init() +void BulkKinetics::resizeSpecies() { - m_kk = thermo().nSpecies(); + Kinetics::resizeSpecies(); m_conc.resize(m_kk); m_grt.resize(m_kk); } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index d2654d4a6..58cd21d47 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -828,20 +828,6 @@ void InterfaceKinetics::addPhase(thermo_t& thermo) void InterfaceKinetics::init() { - m_kk = 0; - for (size_t n = 0; n < nPhases(); n++) { - m_kk += thermo(n).nSpecies(); - } - m_actConc.resize(m_kk); - m_conc.resize(m_kk); - m_StandardConc.resize(m_kk, 0.0); - m_mu0.resize(m_kk); - m_mu.resize(m_kk); - m_mu0_Kc.resize(m_kk); - m_grt.resize(m_kk); - m_pot.resize(m_kk, 0.0); - m_phi.resize(nPhases(), 0.0); - size_t ks = reactionPhaseIndex(); if (ks == npos) throw CanteraError("InterfaceKinetics::finalize", "no surface phase is present."); @@ -855,6 +841,25 @@ void InterfaceKinetics::init() } } +void InterfaceKinetics::resizeSpecies() +{ + size_t kOld = m_kk; + Kinetics::resizeSpecies(); + if (m_kk != kOld && nReactions()) { + throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add" + " species to InterfaceKinetics after reactions have been added."); + } + m_actConc.resize(m_kk); + m_conc.resize(m_kk); + m_StandardConc.resize(m_kk, 0.0); + m_mu0.resize(m_kk); + m_mu.resize(m_kk); + m_mu0_Kc.resize(m_kk); + m_grt.resize(m_kk); + m_pot.resize(m_kk, 0.0); + m_phi.resize(nPhases(), 0.0); +} + doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const { for (size_t i = 0; i < m_ctrxn.size(); i++) { diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 18a8075c0..fae7da91b 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -483,16 +483,6 @@ void Kinetics::getNetProductionRates(doublereal* net) void Kinetics::addPhase(thermo_t& thermo) { - // if not the first thermo object, set the start position - // to that of the last object added + the number of its species - if (m_thermo.size() > 0) { - m_start.push_back(m_start.back() - + m_thermo.back()->nSpecies()); - } else { - // otherwise start at 0 - m_start.push_back(0); - } - // the phase with lowest dimensionality is assumed to be the // phase/interface at which reactions take place if (thermo.nDim() <= m_mindim) { @@ -513,7 +503,19 @@ void Kinetics::addPhase(thermo_t& thermo) } m_thermo.push_back(&thermo); m_phaseindex[m_thermo.back()->id()] = nPhases(); - m_kk += thermo.nSpecies(); + resizeSpecies(); +} + +void Kinetics::resizeSpecies() +{ + m_kk = 0; + m_start.resize(nPhases()); + + for (size_t i = 0; i < m_thermo.size(); i++) { + m_start[i] = m_kk; // global index of first species of phase i + m_kk += m_thermo[i]->nSpecies(); + } + invalidateCache(); } void Kinetics::finalize() @@ -525,6 +527,10 @@ void Kinetics::finalize() bool Kinetics::addReaction(shared_ptr r) { r->validate(); + if (m_kk == 0) { + init(); + } + resizeSpecies(); // If reaction orders are specified, then this reaction does not follow // mass-action kinetics, and is not an elementary reaction. So check that it diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 83e19fc43..b553c7ed0 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -19,7 +19,6 @@ public: th.push_back(&p_ref); importKinetics(p_ref.xml(), th, &kin_ref); kin.addPhase(p); - kin.init(); } IdealGasPhase p; @@ -318,7 +317,6 @@ public: importKinetics(surf_ref.xml(), th, &kin_ref); kin.addPhase(surf); kin.addPhase(gas); - kin.init(); } IdealGasPhase gas; @@ -380,3 +378,109 @@ TEST_F(InterfaceKineticsFromScratch, add_sticking_reaction) kin.addReaction(R); check_rates(0); } + +class KineticsAddSpecies : public testing::Test +{ +public: + KineticsAddSpecies() + : p_ref("../data/kineticsfromscratch.cti") + { + std::vector th; + th.push_back(&p_ref); + importKinetics(p_ref.xml(), th, &kin_ref); + + p.addUndefinedElements(); + kin.addPhase(p); + + std::vector> S = getSpecies(*get_XML_File("h2o2.cti")); + for (auto sp : S) { + species[sp->name] = sp; + } + reactions = getReactions(*get_XML_File("../data/kineticsfromscratch.cti")); + } + + IdealGasPhase p; + IdealGasPhase p_ref; + GasKinetics kin; + GasKinetics kin_ref; + std::vector> reactions; + std::map> species; + + void check_rates(size_t N, const std::string& X) { + for (size_t i = 0; i < kin_ref.nReactions(); i++) { + if (i >= N) { + kin_ref.setMultiplier(i, 0); + } else { + kin_ref.setMultiplier(i, 1); + } + } + p.setState_TPX(1200, 5*OneAtm, X); + p_ref.setState_TPX(1200, 5*OneAtm, X); + + vector_fp k(kin.nReactions()), k_ref(kin_ref.nReactions()); + vector_fp w(kin.nTotalSpecies()), w_ref(kin_ref.nTotalSpecies()); + + kin.getCreationRates(w.data()); + kin_ref.getCreationRates(w_ref.data()); + for (size_t i = 0; i < kin.nTotalSpecies(); i++) { + size_t iref = p_ref.speciesIndex(p.speciesName(i)); + EXPECT_DOUBLE_EQ(w_ref[iref], w[i]) << "sp = " << p.speciesName(i) << "; N = " << N; + } + + kin.getFwdRateConstants(k.data()); + kin_ref.getFwdRateConstants(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + + kin.getRevRateConstants(k.data()); + kin_ref.getRevRateConstants(k_ref.data()); + for (size_t i = 0; i < kin.nReactions(); i++) { + EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; + } + } +}; + +TEST_F(KineticsAddSpecies, add_species_sequential) +{ + ASSERT_EQ((size_t) 0, kin.nReactions()); + + for (auto s : {"AR", "O", "H2", "H", "OH"}) { + p.addSpecies(species[s]); + } + kin.addReaction(reactions[0]); + ASSERT_EQ(5, (int) kin.nTotalSpecies()); + check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88"); + + p.addSpecies(species["O2"]); + p.addSpecies(species["H2O"]); + kin.addReaction(reactions[1]); + ASSERT_EQ(7, (int) kin.nTotalSpecies()); + ASSERT_EQ(2, (int) kin.nReactions()); + check_rates(2, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38"); + + p.addSpecies(species["H2O2"]); + kin.addReaction(reactions[2]); + kin.addReaction(reactions[3]); + check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38"); // no change + check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.35, H2O2:0.03"); + + p.addSpecies(species["HO2"]); + kin.addReaction(reactions[4]); + check_rates(5, "O:0.01, H2:0.1, H:0.02, OH:0.03, O2:0.4, AR:0.3, H2O2:0.03, HO2:0.01"); +} + +TEST_F(KineticsAddSpecies, add_species_err_first) +{ + for (auto s : {"AR", "O", "H2", "H"}) { + p.addSpecies(species[s]); + } + ASSERT_THROW(kin.addReaction(reactions[0]), CanteraError); + ASSERT_EQ((size_t) 0, kin.nReactions()); + + p.addSpecies(species["OH"]); + kin.addReaction(reactions[0]); + ASSERT_EQ(5, (int) kin.nTotalSpecies()); + ASSERT_EQ((size_t) 1, kin.nReactions()); + check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88"); +}