[Kinetics] Implement addReaction from Reaction for InterfaceKinetics

include/cantera/kinetics/InterfaceKinetics.h

@@ -11,7 +11,7 @@
 #include "cantera/thermo/mix_defs.h"
 #include "Kinetics.h"
 #include "cantera/kinetics/RxnMolChange.h"
-
+#include "Reaction.h"
 #include "cantera/base/utilities.h"
 #include "RateCoeffMgr.h"
 
@@ -210,6 +210,7 @@ public:
 
     virtual void init();
     virtual void addReaction(ReactionData& r);
+    virtual void addReaction(shared_ptr<Reaction> r);
     virtual void finalize();
     virtual bool ready() const;
     //! @}
@@ -389,8 +390,12 @@ public:
     int phaseStability(const size_t iphase) const;
 
     virtual void determineFwdOrdersBV(ReactionData& rdata, vector_fp& fwdFullorders);
+    virtual void determineFwdOrdersBV(ElectrochemicalReaction& r, vector_fp& fwdFullorders);
 
 protected:
+    void addElementaryReaction(InterfaceReaction& rdata);
+    void addGlobalReaction(InterfaceReaction& r);
+
     //! Temporary work vector of length m_kk
     vector_fp m_grt;
 

include/cantera/kinetics/Reaction.h

@@ -209,6 +209,8 @@ public:
 
     //! Forward value of the apparent Electrochemical transfer coefficient
     doublereal beta;
+
+    bool exchange_current_density_formulation;
 };
 
 }

src/kinetics/InterfaceKinetics.cpp

@@ -896,6 +896,122 @@ void InterfaceKinetics::addReaction(ReactionData& r)
     }
 }
 
+void InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
+{
+    InterfaceReaction& r = dynamic_cast<InterfaceReaction&>(*r_base);
+    // Create a SurfaceArrhenius rate calculator and set the coverage dependencies
+    SurfaceArrhenius rate(r.rate);
+
+    // Turn on the global flag indicating surface coverage dependence
+    if (!r.coverage_deps.empty()) {
+        m_has_coverage_dependence = true;
+    }
+
+    for (map<string, CoverageDependency>::const_iterator iter = r.coverage_deps.begin();
+         iter != r.coverage_deps.end();
+         ++iter) {
+        size_t k = kineticsSpeciesIndex(iter->first);
+        rate.addCoverageDependence(k, iter->second.a, iter->second.m, iter->second.E);
+    }
+
+    m_rates.install(m_ii, rate);
+
+    // Store activation energy
+    m_E.push_back(rate.activationEnergy_R());
+
+    ElectrochemicalReaction* re = dynamic_cast<ElectrochemicalReaction*>(&r);
+    if (re) {
+        m_has_electrochem_rxns = true;
+        m_beta.push_back(re->beta);
+        m_ctrxn.push_back(m_ii);
+        if (re->exchange_current_density_formulation) {
+            m_has_exchange_current_density_formulation = true;
+            m_ctrxn_ecdf.push_back(1);
+        } else {
+            m_ctrxn_ecdf.push_back(0);
+        }
+        m_ctrxn_resistivity_.push_back(re->film_resistivity);
+
+        if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN ||
+            r.reaction_type == BUTLERVOLMER_RXN ||
+            r.reaction_type == SURFACEAFFINITY_RXN ||
+            r.reaction_type == GLOBAL_RXN) {
+            //   Specify alternative forms of the electrochemical reaction
+            if (r.reaction_type == BUTLERVOLMER_RXN) {
+                m_ctrxn_BVform.push_back(1);
+            } else if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) {
+                m_ctrxn_BVform.push_back(2);
+            } else {
+                // set the default to be the normal forward / reverse calculation method
+                m_ctrxn_BVform.push_back(0);
+            }
+            if (!r.orders.empty()) {
+                vector_fp orders(nTotalSpecies(), 0.0);
+                for (Composition::const_iterator iter = r.orders.begin();
+                     iter != r.orders.end();
+                     ++iter) {
+                    orders[kineticsSpeciesIndex(iter->first)] = iter->second;
+                }
+                RxnOrders* ro = new RxnOrders();
+                ro->fill(orders);
+                m_ctrxn_ROPOrdersList_.push_back(ro);
+                m_ctrxn_FwdOrdersList_.push_back(0);
+
+                // Fill in the Fwd Orders dependence here for B-V reactions
+                if (r.reaction_type == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN ||
+                    r.reaction_type == BUTLERVOLMER_RXN) {
+                    vector_fp fwdFullorders(m_kk, 0.0);
+                    determineFwdOrdersBV(*re, fwdFullorders);
+                    RxnOrders* ro = new RxnOrders();
+                    ro->fill(fwdFullorders);
+                    m_ctrxn_FwdOrdersList_[m_ii] = ro;
+                }
+            } else {
+                m_ctrxn_ROPOrdersList_.push_back(0);
+                m_ctrxn_FwdOrdersList_.push_back(0);
+            }
+
+        } else {
+            m_ctrxn_BVform.push_back(0);
+            m_ctrxn_ROPOrdersList_.push_back(0);
+            m_ctrxn_FwdOrdersList_.push_back(0);
+            if (re->film_resistivity > 0.0) {
+                throw CanteraError("InterfaceKinetics::addReaction()",
+                                   "film resistivity set for elementary reaction");
+            }
+        }
+    }
+
+    if (r.reversible) {
+        m_revindex.push_back(nReactions());
+        m_nrev++;
+    } else {
+        m_irrev.push_back(nReactions());
+        m_nirrev++;
+    }
+    Kinetics::addReaction(r_base);
+
+    m_rxnPhaseIsReactant.push_back(std::vector<bool>(nPhases(), false));
+    m_rxnPhaseIsProduct.push_back(std::vector<bool>(nPhases(), false));
+
+    size_t i = m_ii - 1;
+    for (Composition::const_iterator iter = r.reactants.begin();
+         iter != r.reactants.end();
+         ++iter) {
+        size_t k = kineticsSpeciesIndex(iter->first);
+        size_t p = speciesPhaseIndex(k);
+        m_rxnPhaseIsReactant[i][p] = true;
+    }
+    for (Composition::const_iterator iter = r.products.begin();
+         iter != r.products.end();
+         ++iter) {
+        size_t k = kineticsSpeciesIndex(iter->first);
+        size_t p = speciesPhaseIndex(k);
+        m_rxnPhaseIsProduct[i][p] = true;
+    }
+}
+
+
 void InterfaceKinetics::setIOFlag(int ioFlag)
 {
     m_ioFlag = ioFlag;
@@ -1102,6 +1218,49 @@ void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector<do
     }
 }
 
+void InterfaceKinetics::determineFwdOrdersBV(ElectrochemicalReaction& r, std::vector<doublereal>& fwdFullOrders)
+{
+    // Start out with the full ROP orders vector.
+    // This vector will have the BV exchange current density orders in it.
+    fwdFullOrders.assign(nTotalSpecies(), 0.0);
+    for (Composition::const_iterator iter = r.orders.begin();
+         iter != r.orders.end();
+         ++iter) {
+        fwdFullOrders[kineticsSpeciesIndex(iter->first)] = iter->second;
+    }
+
+    //   forward and reverse beta values
+    double betaf = r.beta;
+    double betar = 1.0 - betaf;
+
+    // Loop over the reactants doing away with the BV terms.
+    // This should leave the reactant terms only, even if they are non-mass action.
+    for (Composition::const_iterator iter = r.reactants.begin();
+         iter != r.reactants.end();
+         ++iter) {
+        size_t k = kineticsSpeciesIndex(iter->first);
+        fwdFullOrders[k] += betaf * iter->second;
+        // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
+        if (abs(fwdFullOrders[k]) < 0.00001) {
+            fwdFullOrders[k] = 0.0;
+        }
+    }
+
+    // Loop over the products doing away with the BV terms.
+    // This should leave the reactant terms only, even if they are non-mass action.
+    for (Composition::const_iterator iter = r.products.begin();
+         iter != r.products.end();
+         ++iter) {
+        size_t k = kineticsSpeciesIndex(iter->first);
+        fwdFullOrders[k] -= betaf * iter->second;
+        // just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
+        if (abs(fwdFullOrders[k]) < 0.00001) {
+            fwdFullOrders[k] = 0.0;
+        }
+    }
+}
+
+
 void EdgeKinetics::finalize()
 {
     //  Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below.

src/kinetics/Reaction.cpp

@@ -109,6 +109,7 @@ ElectrochemicalReaction::ElectrochemicalReaction(const Composition& reactants_,
     , equilibrium_constant_power(1.0)
     , affinity_power(1.0)
     , beta(0.0)
+    , exchange_current_density_formulation(false)
 {
 }
 

test/kinetics/kineticsFromScratch.cpp

@@ -1,7 +1,9 @@
 #include "gtest/gtest.h"
 #include "cantera/kinetics/importKinetics.h"
 #include "cantera/thermo/IdealGasPhase.h"
+#include "cantera/thermo/SurfPhase.h"
 #include "cantera/kinetics/GasKinetics.h"
+#include "cantera/kinetics/InterfaceKinetics.h"
 #include "cantera/base/Array.h"
 
 using namespace Cantera;
@@ -158,3 +160,68 @@ TEST_F(KineticsFromScratch, add_chebyshev_reaction)
     kin.finalize();
     check_rates(4);
 }
+
+class InterfaceKineticsFromScratch : public testing::Test
+{
+public:
+    InterfaceKineticsFromScratch()
+        : gas("../data/sofc-test.xml", "gas")
+        , gas_ref("../data/sofc-test.xml", "gas")
+        , surf("../data/sofc-test.xml", "metal_surface")
+        , surf_ref("../data/sofc-test.xml", "metal_surface")
+    {
+        std::vector<ThermoPhase*> th;
+        th.push_back(&surf_ref);
+        th.push_back(&gas_ref);
+        importKinetics(surf_ref.xml(), th, &kin_ref);
+        kin.addPhase(surf);
+        kin.addPhase(gas);
+        kin.init();
+    }
+
+    IdealGasPhase gas;
+    IdealGasPhase gas_ref;
+    SurfPhase surf;
+    SurfPhase surf_ref;
+    InterfaceKinetics kin;
+    InterfaceKinetics 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 = "H2:0.2 O2:0.5 H2O:0.1 N2:0.2";
+        std::string Xs = "H(m):0.1 O(m):0.2 OH(m):0.3 (m):0.4";
+        gas.setState_TPX(1200, 5*OneAtm, X);
+        gas_ref.setState_TPX(1200, 5*OneAtm, X);
+        surf.setState_TP(1200, 5*OneAtm);
+        surf_ref.setState_TP(1200, 5*OneAtm);
+        surf.setCoveragesByName(Xs);
+        surf_ref.setCoveragesByName(Xs);
+
+        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(InterfaceKineticsFromScratch, add_surface_reaction)
+{
+    // Reaction 3 on the metal surface
+    // surface_reaction( "H(m) + O(m) <=> OH(m) + (m)",
+    //                   [5.00000E+22, 0, 100.0], id = 'metal-rxn4')
+    Composition reac = parseCompString("H(m):1 O(m):1");
+    Composition prod = parseCompString("OH(m):1 (m):1");
+    Arrhenius rate(5e21, 0, 100.0e6 / GasConstant); // kJ/mol -> J/kmol
+
+    shared_ptr<InterfaceReaction>R(new InterfaceReaction(reac, prod, rate));
+    kin.addReaction(R);
+    kin.finalize();
+    check_rates(3);
+}