From 1ba5f6b879871afe507a1391fad448587f2f2bff Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 16 Apr 2016 00:38:23 -0400 Subject: [PATCH] [Kinetics] Allow adding species after reactions for homogeneous kinetics This enables incremental mechanism construction for gas phase kinetics. For surface kinetics, adding new species changes the kinetics species index of existing species in other phases, so this feature is disabled. --- include/cantera/kinetics/BulkKinetics.h | 2 +- include/cantera/kinetics/InterfaceKinetics.h | 1 + include/cantera/kinetics/Kinetics.h | 18 ++-- src/kinetics/BulkKinetics.cpp | 4 +- src/kinetics/InterfaceKinetics.cpp | 33 +++--- src/kinetics/Kinetics.cpp | 28 +++-- test/kinetics/kineticsFromScratch.cpp | 108 ++++++++++++++++++- 7 files changed, 158 insertions(+), 36 deletions(-) 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"); +}