[Input] Create InterfaceKinetics and InterfaceReaction objects from YAML

src/kinetics/KineticsFactory.cpp

@@ -43,6 +43,7 @@ KineticsFactory::KineticsFactory() {
     reg("gaskinetics", []() { return new GasKinetics(); });
     m_synonyms["gas"] = "gaskinetics";
     reg("interface", []() { return new InterfaceKinetics(); });
+    m_synonyms["surface"] = "interface";
     reg("edge", []() { return new EdgeKinetics(); });
 }
 

src/kinetics/Reaction.cpp

@@ -287,6 +287,12 @@ Arrhenius readArrhenius(const XML_Node& arrhenius_node)
 Units rateCoeffUnits(const Reaction& R, const Kinetics& kin,
                      int pressure_dependence=0)
 {
+    if (R.reaction_type == INTERFACE_RXN
+        && dynamic_cast<const InterfaceReaction&>(R).is_sticking_coefficient) {
+        // Sticking coefficients are dimensionless
+        return Units();
+    }
+
     // Determine the units of the rate coefficient
     double reaction_phase_ndim = static_cast<double>(
         kin.thermo(kin.reactionPhaseIndex()).nDim());
@@ -783,6 +789,44 @@ void setupInterfaceReaction(InterfaceReaction& R, const XML_Node& rxn_node)
     setupElementaryReaction(R, rxn_node);
 }
 
+void setupInterfaceReaction(InterfaceReaction& R, const AnyMap& node,
+                            const Kinetics& kin)
+{
+    setupReaction(R, node);
+    R.allow_negative_pre_exponential_factor = node.getBool("negative-A", false);
+
+    if (node.hasKey("rate-constant")) {
+        R.rate = readArrhenius(R, node["rate-constant"], kin, node.units());
+    } else if (node.hasKey("sticking-coefficient")) {
+        R.is_sticking_coefficient = true;
+        R.rate = readArrhenius(R, node["sticking-coefficient"], kin, node.units());
+        R.use_motz_wise_correction = node.getBool("Motz-Wise",
+            kin.thermo().input().getBool("Motz-Wise", false));
+        R.sticking_species = node.getString("sticking-species", "");
+    } else {
+        throw CanteraError("setupInterfaceReaction", "Reaction must include "
+            "either a 'rate-constant' or 'sticking-coefficient' node.");
+    }
+
+    if (node.hasKey("coverage-dependencies")) {
+        for (const auto& item : node["coverage-dependencies"].as<AnyMap>()) {
+            double a, E, m;
+            if (item.second.is<AnyMap>()) {
+                auto& cov_map = item.second.as<AnyMap>();
+                a = cov_map["a"].asDouble();
+                m = cov_map["m"].asDouble();
+                E = node.units().convertActivationEnergy(cov_map["E"], "K");
+            } else {
+                auto& cov_vec = item.second.asVector<AnyValue>(3);
+                a = cov_vec[0].asDouble();
+                m = cov_vec[1].asDouble();
+                E = node.units().convertActivationEnergy(cov_vec[2], "K");
+            }
+            R.coverage_deps[item.first] = CoverageDependency(a, E, m);
+        }
+    }
+}
+
 void setupElectrochemicalReaction(ElectrochemicalReaction& R,
                                   const XML_Node& rxn_node)
 {
@@ -935,6 +979,12 @@ unique_ptr<Reaction> newReaction(const AnyMap& node, const Kinetics& kin)
         type = node["type"].asString();
     }
 
+    if (kin.thermo().nDim() < 3) {
+        unique_ptr<InterfaceReaction> R(new InterfaceReaction());
+        setupInterfaceReaction(*R, node, kin);
+        return unique_ptr<Reaction>(move(R));
+    }
+
     if (type == "elementary") {
         unique_ptr<ElementaryReaction> R(new ElementaryReaction());
         setupElementaryReaction(*R, node, kin);

test/data/surface-phases.yaml

@@ -1,21 +1,52 @@
 phases:
+- name: gas
+  thermo: ideal-gas
+  species: [{gas-species: all}]
+  state: {T: 900 K, P: 1 atm, X: {H2: 0.4, Ar: 0.6}}
+
 - name: Pt-surf
   thermo: ideal-surface
+  kinetics: surface
   species: [{Pt-surf-species: all}]
   site-density: 2.7063e-9 mol/cm^2
   state: {T: 900 K, P: 1 atm, coverages: {Pt(s): 0.5, H(s): 0.4, O(s): 0.1}}
 
 - name: TPB
   thermo: edge
-  species: [(tpb)]
+  species: [{tpb-species: [(tpb)]}]
   site-density: 5e-17 mol/cm
   state: {T: 1073.15 K, P: 1 atm, coverages: {(tpb): 1.0}}
 
-species:
-- name: (tpb)
-  composition: {}
+
+gas-species:
+- name: H2
+  composition: {H: 2}
   thermo:
-    model: constant-cp
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [2.34433112, 0.00798052075, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
+      -917.935173, 0.683010238]
+    - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14,
+      -950.158922, -3.20502331]
+- name: H
+  composition: {H: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22,
+      25473.6599, -0.446682853]
+    - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22,
+      25473.6599, -0.446682914]
+- name: Ar
+  composition: {Ar: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 5000.0]
+    data:
+    - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366]
+
 
 Pt-surf-species:
 - name: Pt(s)
@@ -42,3 +73,21 @@ Pt-surf-species:
        3.3787992E-12, -1.3209912E+04, 3.6137905E+00]
     - [1.9454180E+00, 9.1761647E-04, -1.1226719E-07, -9.9099624E-11,
        2.4307699E-14, -1.4005187E+04, -1.1531663E+01]
+
+
+tpb-species:
+- name: (tpb)
+  composition: {}
+  thermo:
+    model: constant-cp
+
+
+reactions:
+- equation: H2 + 2 Pt(s) => 2 H(s)
+  rate-constant: {A: 4.4579e10 cm^3/mol/s, b: 0.5, Ea: 0.0}
+  orders: {H2: 1.0, Pt(s): 1.0}
+- equation: 2 H(s) => H2 + 2 Pt(s)
+  rate-constant: {A: 3.7e21 cm^2/mol/s, b: 0, Ea: 67400 J/mol}
+  coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}
+- equation: H + Pt(s) => H(s)
+  sticking-coefficient: [1.0, 0.0, 0.0]

test/kinetics/kineticsFromYaml.cpp

@@ -1,6 +1,7 @@
 #include "gtest/gtest.h"
 #include "cantera/base/Units.h"
 #include "cantera/IdealGasMix.h"
+#include "cantera/thermo/SurfPhase.h"
 #include "cantera/kinetics/KineticsFactory.h"
 #include "cantera/thermo/ThermoFactory.h"
 
@@ -197,3 +198,28 @@ TEST(Kinetics, GasKineticsFromYaml2)
     auto kin = newKinetics(phases, *phaseNode, infile);
     EXPECT_EQ(kin->nReactions(), (size_t) 3);
 }
+
+TEST(Kinetics, InterfaceKineticsFromYaml)
+{
+    shared_ptr<ThermoPhase> gas(newPhase("surface-phases.yaml", "gas"));
+    shared_ptr<ThermoPhase> surf_tp(newPhase("surface-phases.yaml", "Pt-surf"));
+    shared_ptr<SurfPhase> surf = std::dynamic_pointer_cast<SurfPhase>(surf_tp);
+    std::vector<ThermoPhase*> phases{surf_tp.get(), gas.get()};
+    auto kin = newKinetics(phases, "surface-phases.yaml", "Pt-surf");
+    EXPECT_EQ(kin->nReactions(), (size_t) 3);
+    EXPECT_EQ(kin->nTotalSpecies(), (size_t) 6);
+    auto R1 = kin->reaction(0);
+    auto IR1 = std::dynamic_pointer_cast<InterfaceReaction>(R1);
+    EXPECT_DOUBLE_EQ(R1->orders["Pt(s)"], 1.0);
+    EXPECT_DOUBLE_EQ(IR1->rate.preExponentialFactor(), 4.4579e7);
+
+    auto R2 = kin->reaction(1);
+    auto IR2 = std::dynamic_pointer_cast<InterfaceReaction>(R2);
+    EXPECT_DOUBLE_EQ(IR2->rate.preExponentialFactor(), 3.7e20);
+    EXPECT_DOUBLE_EQ(IR2->coverage_deps["H(s)"].E, -6e6 / GasConstant);
+    EXPECT_FALSE(IR2->is_sticking_coefficient);
+
+    auto R3 = kin->reaction(2);
+    auto IR3 = std::dynamic_pointer_cast<InterfaceReaction>(R3);
+    EXPECT_TRUE(IR3->is_sticking_coefficient);
+}