From 8b502d4065dd8f1f40b27d06e92ac4b7b03fdd76 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Wed, 5 Dec 2018 18:43:00 -0500
Subject: [PATCH] [Input] Parse YAML entries for three body reactions

---
 src/kinetics/Reaction.cpp          | 31 +++++++++++++++++++++++++++++-
 test/kinetics/kineticsFromYaml.cpp | 31 ++++++++++++++++++++++++++++++
 2 files changed, 61 insertions(+), 1 deletion(-)

diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp
index dfcb384e9..acb386a3c 100644
--- a/src/kinetics/Reaction.cpp
+++ b/src/kinetics/Reaction.cpp
@@ -293,7 +293,10 @@ Arrhenius readArrhenius(const Reaction& R, const AnyValue& rate_node,
     for (const auto& stoich : R.reactants) {
         // Order for each reactant is the reactant stoichiometric coefficient,
         // unless already overridden by user-specified orders
-        if (R.orders.find(stoich.first) == R.orders.end()) {
+        if (stoich.first == "M") {
+            len_dim += kin.thermo(kin.reactionPhaseIndex()).nDim();
+            quantity_dim -= 1.0;
+        } else if (R.orders.find(stoich.first) == R.orders.end()) {
             len_dim += stoich.second * kin.speciesPhase(stoich.first).nDim();
             quantity_dim -= stoich.second;
         }
@@ -360,6 +363,14 @@ void readEfficiencies(ThirdBody& tbody, const XML_Node& rc_node)
     tbody.efficiencies = parseCompString(eff_node.value());
 }
 
+void readEfficiencies(ThirdBody& tbody, const AnyMap& node)
+{
+    tbody.default_efficiency = node.getDouble("default-efficiency", 1.0);
+    if (node.hasKey("efficiencies")) {
+        tbody.efficiencies = node.at("efficiencies").asMap<double>();
+    }
+}
+
 void setupReaction(Reaction& R, const XML_Node& rxn_node)
 {
     // Reactant and product stoichiometries
@@ -473,6 +484,20 @@ void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node)
     setupElementaryReaction(R, rxn_node);
 }
 
+void setupThreeBodyReaction(ThreeBodyReaction& R, const AnyMap& node,
+                            const Kinetics& kin, const UnitSystem& units)
+{
+    setupElementaryReaction(R, node, kin, units);
+    if (R.reactants.count("M") != 1 || R.products.count("M") != 1) {
+        throw CanteraError("setupThreeBodyReaction",
+            "Reaction equation '{}' does not contain third body 'M'",
+            node.at("equation").asString());
+    }
+    R.reactants.erase("M");
+    R.products.erase("M");
+    readEfficiencies(R.third_body, node);
+}
+
 void setupFalloffReaction(FalloffReaction& R, const XML_Node& rxn_node)
 {
     XML_Node& rc_node = rxn_node.child("rateCoeff");
@@ -761,6 +786,10 @@ unique_ptr<Reaction> newReaction(const AnyMap& node, const Kinetics& kin,
         unique_ptr<ElementaryReaction> R(new ElementaryReaction());
         setupElementaryReaction(*R, node, kin, units);
         return unique_ptr<Reaction>(move(R));
+    } else if (type == "three-body") {
+        unique_ptr<ThreeBodyReaction> R(new ThreeBodyReaction());
+        setupThreeBodyReaction(*R, node, kin, units);
+        return unique_ptr<Reaction>(move(R));
     } else {
         throw CanteraError("newReaction", "Unknown reaction type '{}'", type);
     }
diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp
index 0c73225aa..eacc7d59f 100644
--- a/test/kinetics/kineticsFromYaml.cpp
+++ b/test/kinetics/kineticsFromYaml.cpp
@@ -26,3 +26,34 @@ TEST(Reaction, ElementaryFromYaml)
     EXPECT_TRUE(ER.allow_negative_pre_exponential_factor);
     EXPECT_FALSE(ER.allow_negative_orders);
 }
+
+TEST(Reaction, ThreeBodyFromYaml1)
+{
+    IdealGasMix gas("gri30.xml");
+    AnyMap rxn = AnyMap::fromYamlString(
+        "{equation: 2 O + M <=> O2 + M,"
+        " type: three-body,"
+        " rate-constant: [1.20000E+17 cm^6/mol^2/s, -1, 0],"
+        " efficiencies: {AR: 0.83, H2O: 5}}");
+
+    UnitSystem U;
+    auto R = newReaction(rxn, gas, U);
+    EXPECT_EQ(R->reactants.count("M"), (size_t) 0);
+
+    auto TBR = dynamic_cast<ThreeBodyReaction&>(*R);
+    EXPECT_DOUBLE_EQ(TBR.rate.preExponentialFactor(), 1.2e11);
+    EXPECT_DOUBLE_EQ(TBR.third_body.efficiencies["H2O"], 5.0);
+    EXPECT_DOUBLE_EQ(TBR.third_body.default_efficiency, 1.0);
+}
+
+TEST(Reaction, ThreeBodyFromYaml2)
+{
+    IdealGasMix gas("gri30.xml");
+    AnyMap rxn = AnyMap::fromYamlString(
+        "{equation: 2 O <=> O2," // Missing "M" on each side of the equation
+        " type: three-body,"
+        " rate-constant: [1.20000E+17, -1, 0]}");
+
+    UnitSystem U;
+    EXPECT_THROW(newReaction(rxn, gas, U), CanteraError);
+}