From b4afcd3e8eaf77ea0c57f3a199ee803a218dd4e8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 11 Nov 2014 00:11:25 +0000 Subject: [PATCH] [Kinetics] Implement addReaction using Reaction objects for GasKinetics This covers elementary, third-body, and falloff reaction types --- include/cantera/kinetics/BulkKinetics.h | 4 + include/cantera/kinetics/GasKinetics.h | 7 ++ include/cantera/kinetics/Kinetics.h | 13 +++ include/cantera/kinetics/RateCoeffMgr.h | 11 ++ include/cantera/kinetics/ReactionStoichMgr.h | 1 + src/kinetics/BulkKinetics.cpp | 30 ++++++ src/kinetics/GasKinetics.cpp | 81 +++++++++++++++ src/kinetics/Kinetics.cpp | 68 ++++++++++++ src/kinetics/ReactionStoichMgr.cpp | 1 + test/data/kineticsfromscratch.cti | 20 ++++ test/kinetics/kineticsFromScratch.cpp | 103 +++++++++++++++++++ 11 files changed, 339 insertions(+) create mode 100644 test/data/kineticsfromscratch.cti create mode 100644 test/kinetics/kineticsFromScratch.cpp diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 935900341..9ba6849f4 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -12,6 +12,8 @@ namespace Cantera { +class ElementaryReaction; + //! Partial specialization of Kinetics for chemistry in a single bulk phase class BulkKinetics : public Kinetics { @@ -33,12 +35,14 @@ public: bool doIrreversible = false); virtual void addReaction(ReactionData& r); + virtual void addReaction(shared_ptr r); virtual void init(); virtual void finalize(); virtual bool ready() const; protected: virtual void addElementaryReaction(ReactionData& r); + virtual void addElementaryReaction(ElementaryReaction& r); Rate1 m_rates; std::vector m_revindex; //!< Indices of reversible reactions diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 4acdcdc41..d7aee795c 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -12,6 +12,7 @@ #include "BulkKinetics.h" #include "ThirdBodyMgr.h" #include "FalloffMgr.h" +#include "Reaction.h" namespace Cantera { @@ -52,6 +53,7 @@ public: //! @{ virtual void init(); virtual void addReaction(ReactionData& r); + virtual void addReaction(shared_ptr r); virtual void finalize(); virtual bool ready() const; //@} @@ -105,6 +107,11 @@ protected: void addPlogReaction(ReactionData& r); void addChebyshevReaction(ReactionData& r); + void addThreeBodyReaction(ThirdBodyReaction& r); + void addFalloffReaction(FalloffReaction& r); + void addPlogReaction(PlogReaction& r); + void addChebyshevReaction(ChebyshevReaction& r); + //! Update the equilibrium constants in molar units. void updateKc(); diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 45f94a9f8..5b2b39786 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -14,12 +14,14 @@ #include "StoichManager.h" #include "cantera/thermo/mix_defs.h" #include "cantera/base/global.h" +#include "cantera/base/smart_ptr.h" namespace Cantera { // forward references class ReactionData; +class Reaction; /** * @defgroup chemkinetics Chemical Kinetics @@ -796,6 +798,14 @@ public: */ virtual void addReaction(ReactionData& r); + /** + * Add a single reaction to the mechanism. Derived classes should call the + * base class method in addition to handling their own specialized behavior. + * + * @param r Pointer to the Reaction object to be added. + */ + virtual void addReaction(shared_ptr r); + //! @deprecated To be removed after Cantera 2.2. No longer called as part //! of addReaction. virtual void installReagents(const ReactionData& r) { @@ -909,6 +919,9 @@ protected: /// progress vector. It is initialized to one. vector_fp m_perturb; + //! Vector of Reaction objects represented by this Kinetics manager + std::vector > m_reactions; + /** * This is a vector of vectors containing the reactants for * each reaction. The outer vector is over the number of diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h index 7d267607d..1dfeabfbf 100644 --- a/include/cantera/kinetics/RateCoeffMgr.h +++ b/include/cantera/kinetics/RateCoeffMgr.h @@ -45,6 +45,17 @@ public: return m_rates.size() - 1; } + /** + * Install a rate coefficient calculator. + * @param rxnNumber the reaction number + * @param rdata rate coefficient specification for the reaction + */ + void install(size_t rxnNumber, const R& rate) { + m_rxn.push_back(rxnNumber); + m_rates.push_back(rate); + } + + /** * Update the concentration-dependent parts of the rate * coefficient, if any. Used by class SurfaceArrhenius to diff --git a/include/cantera/kinetics/ReactionStoichMgr.h b/include/cantera/kinetics/ReactionStoichMgr.h index 4ad750e56..736a5f8b7 100644 --- a/include/cantera/kinetics/ReactionStoichMgr.h +++ b/include/cantera/kinetics/ReactionStoichMgr.h @@ -12,6 +12,7 @@ namespace Cantera { class ReactionData; +class Reaction; /** * Reaction mechanism stoichiometry manager. This is an internal class used diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 8cb260c98..eb138a823 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -1,4 +1,5 @@ #include "cantera/kinetics/BulkKinetics.h" +#include "cantera/kinetics/Reaction.h" namespace Cantera { @@ -127,11 +128,40 @@ void BulkKinetics::addReaction(ReactionData& r) Kinetics::addReaction(r); } +void BulkKinetics::addReaction(shared_ptr r) +{ + double dn = 0.0; + for (Composition::const_iterator iter = r->products.begin(); + iter != r->products.end(); + ++iter) { + dn += iter->second; + } + for (Composition::const_iterator iter = r->reactants.begin(); + iter != r->reactants.end(); + ++iter) { + dn -= iter->second; + } + + m_dn.push_back(dn); + + if (r->reversible) { + m_revindex.push_back(nReactions()); + } else { + m_irrev.push_back(nReactions()); + } + Kinetics::addReaction(r); +} + void BulkKinetics::addElementaryReaction(ReactionData& r) { m_rates.install(nReactions(), r); } +void BulkKinetics::addElementaryReaction(ElementaryReaction& r) +{ + m_rates.install(nReactions(), r.rate); +} + void BulkKinetics::init() { m_kk = thermo().nSpecies(); diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index accddd8e4..3145a728e 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -267,6 +267,34 @@ void GasKinetics::addReaction(ReactionData& r) BulkKinetics::addReaction(r); } +void GasKinetics::addReaction(shared_ptr r) +{ + switch (r->reaction_type) { + case ELEMENTARY_RXN: + addElementaryReaction(dynamic_cast(*r)); + break; + case THREE_BODY_RXN: + addThreeBodyReaction(dynamic_cast(*r)); + break; + case FALLOFF_RXN: + case CHEMACT_RXN: + addFalloffReaction(dynamic_cast(*r)); + break; + case PLOG_RXN: + addPlogReaction(dynamic_cast(*r)); + break; + case CHEBYSHEV_RXN: + addChebyshevReaction(dynamic_cast(*r)); + break; + default: + throw CanteraError("GasKinetics::addReaction", + "Unknown reaction type specified: " + int2str(r->reaction_type)); + } + + // operations common to all reaction types + BulkKinetics::addReaction(r); +} + void GasKinetics::addFalloffReaction(ReactionData& r) { // install high and low rate coeff calculators @@ -310,6 +338,59 @@ void GasKinetics::addChebyshevReaction(ReactionData& r) m_cheb_rates.install(nReactions(), r); } +void GasKinetics::addFalloffReaction(FalloffReaction& r) +{ + // install high and low rate coeff calculators + // and extend the high and low rate coeff value vectors + m_falloff_high_rates.install(m_nfall, r.high_rate); + m_rfn_high.push_back(0.0); + m_falloff_low_rates.install(m_nfall, r.low_rate); + m_rfn_low.push_back(0.0); + + // add this reaction number to the list of falloff reactions + m_fallindx.push_back(nReactions()); + + // install the enhanced third-body concentration calculator + map efficiencies; + for (Composition::const_iterator iter = r.third_body.efficiencies.begin(); + iter != r.third_body.efficiencies.end(); + ++iter) { + efficiencies[kineticsSpeciesIndex(iter->first)] = iter->second; + } + m_falloff_concm.install(m_nfall, efficiencies, + r.third_body.default_efficiency); + + // install the falloff function calculator for this reaction + m_falloffn.install(m_nfall, r.falloff_type, r.reaction_type, + r.falloff_parameters); + + // increment the falloff reaction counter + ++m_nfall; +} + +void GasKinetics::addThreeBodyReaction(ThirdBodyReaction& r) +{ + m_rates.install(nReactions(), r.rate); + map efficiencies; + for (Composition::const_iterator iter = r.third_body.efficiencies.begin(); + iter != r.third_body.efficiencies.end(); + ++iter) { + efficiencies[kineticsSpeciesIndex(iter->first)] = iter->second; + } + m_3b_concm.install(nReactions(), efficiencies, + r.third_body.default_efficiency); +} + +void GasKinetics::addPlogReaction(PlogReaction& r) +{ + m_plog_rates.install(nReactions(), r.rate); +} + +void GasKinetics::addChebyshevReaction(ChebyshevReaction& r) +{ + m_cheb_rates.install(nReactions(), r.rate); +} + void GasKinetics::init() { BulkKinetics::init(); diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 5cdb2cc59..e8330ade9 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -9,6 +9,7 @@ #include "cantera/kinetics/Kinetics.h" #include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/Reaction.h" #include "cantera/base/stringUtils.h" using namespace std; @@ -50,6 +51,7 @@ Kinetics& Kinetics::operator=(const Kinetics& right) m_ii = right.m_ii; m_kk = right.m_kk; m_perturb = right.m_perturb; + m_reactions = right.m_reactions; m_reactants = right.m_reactants; m_products = right.m_products; m_rrxn = right.m_rrxn; @@ -471,6 +473,72 @@ void Kinetics::addReaction(ReactionData& r) { m_ropnet.push_back(0.0); } +void Kinetics::addReaction(shared_ptr r) +{ + size_t irxn = nReactions(); + std::vector rk, pk; + vector_fp rstoich, pstoich; + for (Composition::const_iterator iter = r->reactants.begin(); + iter != r->reactants.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + rk.push_back(k); + rstoich.push_back(iter->second); + m_rrxn[k][irxn] = iter->second; + } + m_reactants.push_back(rk); + + for (Composition::const_iterator iter = r->products.begin(); + iter != r->products.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + pk.push_back(k); + pstoich.push_back(iter->second); + m_prxn[k][irxn] = iter->second; + } + m_products.push_back(pk); + + vector_fp rorder = rstoich; + for (Composition::const_iterator iter = r->orders.begin(); + iter != r->orders.end(); + ++iter) { + size_t k = kineticsSpeciesIndex(iter->first); + vector::iterator rloc = std::find(rk.begin(), rk.end(), k); + if (rloc != rk.end()) { + rorder[rloc - rk.begin()] = iter->second; + } else { + // If the reaction order involves a non-reactant species, add an + // extra term to the reactants with zero stoichiometry so that the + // stoichiometry manager can be used to compute the global forward + // reaction rate. + rk.push_back(k); + rstoich.push_back(0.0); + rorder.push_back(iter->second); + } + } + + m_reactantStoich.add(irxn, rk, rorder, rstoich); + // product orders = product stoichiometric coefficients + if (r->reversible) { + m_revProductStoich.add(irxn, pk, pstoich, pstoich); + } else { + m_irrevProductStoich.add(irxn, pk, pstoich, pstoich); + } + + incrementRxnCount(); + m_reactions.push_back(r); + m_rxneqn.push_back(r->equation()); + m_reactantStrings.push_back(r->reactantString()); + m_productStrings.push_back(r->productString()); + m_rxntype.push_back(r->reaction_type); + m_rfn.push_back(0.0); + m_rkcn.push_back(0.0); + m_ropf.push_back(0.0); + m_ropr.push_back(0.0); + m_ropnet.push_back(0.0); +} + + void Kinetics::installGroups(size_t irxn, const vector& r, const vector& p) { diff --git a/src/kinetics/ReactionStoichMgr.cpp b/src/kinetics/ReactionStoichMgr.cpp index 191c1010d..a34f734a0 100644 --- a/src/kinetics/ReactionStoichMgr.cpp +++ b/src/kinetics/ReactionStoichMgr.cpp @@ -9,6 +9,7 @@ #include "cantera/base/ctexceptions.h" #include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/Reaction.h" #include diff --git a/test/data/kineticsfromscratch.cti b/test/data/kineticsfromscratch.cti new file mode 100644 index 000000000..f2e6b5511 --- /dev/null +++ b/test/data/kineticsfromscratch.cti @@ -0,0 +1,20 @@ +units(length='m', time='s', quantity='kmol', act_energy='cal/mol') + +ideal_gas(name = "ohmech", + elements = " O H Ar ", + species = """ h2o2: H2 H O O2 OH H2O HO2 H2O2 AR""", + reactions = "all", + transport = "None", + initial_state = state(temperature = 300.0, + pressure = OneAtm) ) + +reaction('O + H2 <=> H + OH', [3.870000e+01, 2.7, 6260.0]) + +three_body_reaction('2 O + M <=> O2 + M', [1.200000e+11, -1.0, 0.0], + efficiencies='AR:0.83 H2:2.4 H2O:15.4') + +falloff_reaction('2 OH (+ M) <=> H2O2 (+ M)', + kf=[7.400000e+10, -0.37, 0.0], + kf0=[2.300000e+12, -0.9, -1700.0], + efficiencies='AR:0.7 H2:2.0 H2O:6.0', + falloff=Troe(A=0.7346, T3=94.0, T1=1756.0, T2=5182.0)) diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp new file mode 100644 index 000000000..08a46b939 --- /dev/null +++ b/test/kinetics/kineticsFromScratch.cpp @@ -0,0 +1,103 @@ +#include "gtest/gtest.h" +#include "cantera/kinetics/importKinetics.h" +#include "cantera/thermo/IdealGasPhase.h" +#include "cantera/kinetics/GasKinetics.h" + +using namespace Cantera; + +class KineticsFromScratch : public testing::Test +{ +public: + KineticsFromScratch() + : p("../data/kineticsfromscratch.cti") + , p_ref("../data/kineticsfromscratch.cti") + { + std::vector th; + th.push_back(&p_ref); + importKinetics(p_ref.xml(), th, &kin_ref); + kin.addPhase(p); + kin.init(); + } + + IdealGasPhase p; + IdealGasPhase p_ref; + GasKinetics kin; + GasKinetics kin_ref; + + //! iRef is the index of the corresponding reaction in the reference mech + void check_rates(int iRef) { + ASSERT_EQ((size_t) 1, kin.nReactions()); + + std::string X = "O:0.02 H2:0.2 O2:0.7 H:0.03 OH:0.05"; + p.setState_TPX(1200, 5*OneAtm, X); + p_ref.setState_TPX(1200, 5*OneAtm, X); + + vector_fp k(1), k_ref(kin_ref.nReactions()); + + kin.getFwdRateConstants(&k[0]); + kin_ref.getFwdRateConstants(&k_ref[0]); + EXPECT_FLOAT_EQ(k_ref[iRef], k[0]); + + kin.getRevRateConstants(&k[0]); + kin_ref.getRevRateConstants(&k_ref[0]); + EXPECT_FLOAT_EQ(k_ref[iRef], k[0]); + } +}; + +TEST_F(KineticsFromScratch, add_elementary_reaction) +{ + // reaction 0: + // reaction('O + H2 <=> H + OH', [3.870000e+01, 2.7, 6260.0]) + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + Arrhenius rate(3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + shared_ptr R(new ElementaryReaction(reac, prod, rate)); + + kin.addReaction(R); + kin.finalize(); + check_rates(0); +} + +TEST_F(KineticsFromScratch, add_three_body_reaction) +{ + // reaction 1: + // three_body_reaction('2 O + M <=> O2 + M', [1.200000e+11, -1.0, 0.0], + // efficiencies='AR:0.83 H2:2.4 H2O:15.4') + Composition reac = parseCompString("O:2"); + Composition prod = parseCompString("O2:1"); + Arrhenius rate(1.2e11, -1.0, 0.0); + ThirdBody tbody; + tbody.efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4"); + shared_ptr R(new ThirdBodyReaction(reac, prod, rate, tbody)); + + kin.addReaction(R); + kin.finalize(); + check_rates(1); +} + +TEST_F(KineticsFromScratch, add_falloff_reaction) +{ + // reaction 2: + // falloff_reaction('2 OH (+ M) <=> H2O2 (+ M)', + // kf=[7.400000e+10, -0.37, 0.0], + // kf0=[2.300000e+12, -0.9, -1700.0], + // efficiencies='AR:0.7 H2:2.0 H2O:6.0', + // falloff=Troe(A=0.7346, T3=94.0, T1=1756.0, T2=5182.0)) + Composition reac = parseCompString("OH:2"); + Composition prod = parseCompString("H2O2:1"); + Arrhenius high_rate(7.4e10, -0.37, 0.0); + Arrhenius low_rate(2.3e12, -0.9, -1700.0 / GasConst_cal_mol_K); + vector_fp falloff_params; + falloff_params.push_back(0.7346); + falloff_params.push_back(94.0); + falloff_params.push_back(1756.0); + falloff_params.push_back(5182.0); + ThirdBody tbody; + tbody.efficiencies = parseCompString("AR:0.7 H2:2.0 H2O:6.0"); + shared_ptr R(new FalloffReaction(reac, prod, low_rate, + high_rate, tbody, TROE_FALLOFF, + falloff_params)); + kin.addReaction(R); + kin.finalize(); + check_rates(2); +}