[Kinetics] Implement addReaction using Reaction objects for GasKinetics

This covers elementary, third-body, and falloff reaction types

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<Reaction> r);
     virtual void init();
     virtual void finalize();
     virtual bool ready() const;
 
 protected:
     virtual void addElementaryReaction(ReactionData& r);
+    virtual void addElementaryReaction(ElementaryReaction& r);
 
     Rate1<Arrhenius> m_rates;
     std::vector<size_t> m_revindex; //!< Indices of reversible reactions

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<Reaction> 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();
 

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<Reaction> 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<shared_ptr<Reaction> > m_reactions;
+
     /**
      * This is a vector of vectors containing the reactants for
      * each reaction. The outer vector is over the number of

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

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

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<Reaction> 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();

src/kinetics/GasKinetics.cpp

@@ -267,6 +267,34 @@ void GasKinetics::addReaction(ReactionData& r)
     BulkKinetics::addReaction(r);
 }
 
+void GasKinetics::addReaction(shared_ptr<Reaction> r)
+{
+    switch (r->reaction_type) {
+    case ELEMENTARY_RXN:
+        addElementaryReaction(dynamic_cast<ElementaryReaction&>(*r));
+        break;
+    case THREE_BODY_RXN:
+        addThreeBodyReaction(dynamic_cast<ThirdBodyReaction&>(*r));
+        break;
+    case FALLOFF_RXN:
+    case CHEMACT_RXN:
+        addFalloffReaction(dynamic_cast<FalloffReaction&>(*r));
+        break;
+    case PLOG_RXN:
+        addPlogReaction(dynamic_cast<PlogReaction&>(*r));
+        break;
+    case CHEBYSHEV_RXN:
+        addChebyshevReaction(dynamic_cast<ChebyshevReaction&>(*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<size_t, double> 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<size_t, double> 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();

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<Reaction> r)
+{
+    size_t irxn = nReactions();
+    std::vector<size_t> 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<size_t>::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<grouplist_t>& r,
                              const vector<grouplist_t>& p)
 {

src/kinetics/ReactionStoichMgr.cpp

@@ -9,6 +9,7 @@
 
 #include "cantera/base/ctexceptions.h"
 #include "cantera/kinetics/ReactionData.h"
+#include "cantera/kinetics/Reaction.h"
 
 #include <fstream>
 

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))

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<ThermoPhase*> 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<ElementaryReaction> 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<ThirdBodyReaction> 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<FalloffReaction> R(new FalloffReaction(reac, prod, low_rate,
+                                                      high_rate, tbody, TROE_FALLOFF,
+                                                      falloff_params));
+    kin.addReaction(R);
+    kin.finalize();
+    check_rates(2);
+}