Detect unspecified thirdbody collision partners

include/cantera/kinetics/Reaction.h

@@ -212,6 +212,8 @@ class ThreeBodyReaction : public ElementaryReaction
     //! Relative efficiencies of third-body species in enhancing the reaction
     //! rate.
     ThirdBody third_body;
+
+    bool specified = false; //!< Explicitly specified collision partner
 };
 
 //! A reaction that is first-order in [M] at low pressure, like a third-body
@@ -505,6 +507,8 @@ void parseReactionEquation(Reaction& R, const AnyValue& equation,
                            const Kinetics& kin);
 
 // declarations of setup functions
+void setupReaction(Reaction& R, const XML_Node& rxn_node);
+
 void setupElementaryReaction(ElementaryReaction&, const XML_Node&);
 //! @internal May be changed without notice in future versions
 void setupElementaryReaction(ElementaryReaction&, const AnyMap&,

interfaces/cython/cantera/test/test_kinetics.py

@@ -1221,9 +1221,11 @@ def test_modify_invalid(self):
             self.gas.modify_reaction(0, R2)
 
         # different reactants
-        R = self.gas.reaction(7)
+        R = self.gas.reaction(7) # implicitly defined three-body reaction
+        R2 = ct.ElementaryReaction(R.reactants, R.products)
+        R2.rate = R.rate
         with self.assertRaisesRegex(ct.CanteraError, 'Reactants are different'):
-            self.gas.modify_reaction(23, R)
+            self.gas.modify_reaction(23, R2)
 
         # different products
         R = self.gas.reaction(14)

interfaces/cython/cantera/test/test_reaction.py

@@ -4,6 +4,58 @@
 from . import utilities
 
 
+class TestImplicitThirdBody(utilities.CanteraTest):
+
+    @classmethod
+    def setUpClass(cls):
+        utilities.CanteraTest.setUpClass()
+        cls.gas = ct.Solution('gri30.yaml')
+
+    def test_implicit_three_body(self):
+        yaml1 = """
+            equation: H + 2 O2 <=> HO2 + O2
+            rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
+            """
+        rxn1 = ct.Reaction.fromYaml(yaml1, self.gas)
+        self.assertEqual(rxn1.reaction_type, 'three-body')
+
+        yaml2 = """
+            equation: H + O2 + M <=> HO2 + M
+            rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
+            type: three-body
+            default-efficiency: 0
+            efficiencies: {O2: 1.0}
+            """
+        rxn2 = ct.Reaction.fromYaml(yaml2, self.gas)
+        self.assertEqual(rxn1.efficiencies, rxn2.efficiencies)
+        self.assertEqual(rxn1.default_efficiency, rxn2.default_efficiency)
+
+    def test_non_integer_stoich(self):
+        yaml = """
+            equation: H + 1.5 O2 <=> HO2 + O2
+            rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
+            """
+        rxn = ct.Reaction.fromYaml(yaml, self.gas)
+        self.assertEqual(rxn.reaction_type, 'elementary')
+
+    def test_not_three_body(self):
+        yaml = """
+            equation: HCNO + H <=> H + HNCO  # Reaction 270
+            rate-constant: {A: 2.1e+15, b: -0.69, Ea: 2850.0}
+            """
+        rxn = ct.Reaction.fromYaml(yaml, self.gas)
+        self.assertEqual(rxn.reaction_type, 'elementary')
+
+    def test_user_override(self):
+        yaml = """
+            equation: H + 2 O2 <=> HO2 + O2
+            rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
+            type: elementary
+            """
+        rxn = ct.Reaction.fromYaml(yaml, self.gas)
+        self.assertEqual(rxn.reaction_type, 'elementary')
+
+
 class TestElementary(utilities.CanteraTest):
 
     _cls = ct.ElementaryReaction

src/kinetics/Kinetics.cpp

@@ -143,8 +143,13 @@ std::pair<size_t, size_t> Kinetics::checkDuplicates(bool throw_err) const
                     continue; // No overlap in third body efficiencies
                 }
             } else if (R.type() == "three-body") {
-                ThirdBody& tb1 = dynamic_cast<ThreeBodyReaction&>(R).third_body;
-                ThirdBody& tb2 = dynamic_cast<ThreeBodyReaction&>(other).third_body;
+                auto& R1 = dynamic_cast<ThreeBodyReaction&>(R);
+                auto& R2 = dynamic_cast<ThreeBodyReaction&>(other);
+                if (R1.specified != R2.specified) {
+                    continue; // No overlap as one third body is explicitly specified
+                }
+                ThirdBody& tb1 = R1.third_body;
+                ThirdBody& tb2 = R2.third_body;
                 bool thirdBodyOk = true;
                 for (size_t k = 0; k < nTotalSpecies(); k++) {
                     string s = kineticsSpeciesName(k);

src/kinetics/Reaction.cpp

@@ -276,18 +276,37 @@ ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants_,
 }
 
 std::string ThreeBodyReaction::reactantString() const {
-    return ElementaryReaction::reactantString() + " + M";
+    if (specified) {
+        return ElementaryReaction::reactantString() + " + "
+            + third_body.efficiencies.begin()->first;
+    } else {
+        return ElementaryReaction::reactantString() + " + M";
+    }
 }
 
 std::string ThreeBodyReaction::productString() const {
-    return ElementaryReaction::productString() + " + M";
+    if (specified) {
+        return ElementaryReaction::productString() + " + "
+            + third_body.efficiencies.begin()->first;
+    } else {
+        return ElementaryReaction::productString() + " + M";
+    }
 }
 
 void ThreeBodyReaction::calculateRateCoeffUnits(const Kinetics& kin)
 {
     ElementaryReaction::calculateRateCoeffUnits(kin);
-    const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
-    rate_units *= rxn_phase.standardConcentrationUnits().pow(-1);
+    bool specified = false;
+    for (const auto& reac : reactants) {
+        if (reac.first != "M" && products.count(reac.first)) {
+            // detected specified third-body collision partner
+            specified = true;
+        }
+    }
+    if (!specified) {
+        const ThermoPhase& rxn_phase = kin.thermo(kin.reactionPhaseIndex());
+        rate_units *= rxn_phase.standardConcentrationUnits().pow(-1);
+    }
 }
 
 void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const
@@ -858,6 +877,46 @@ BlowersMasel readBlowersMasel(const Reaction& R, const AnyValue& rate,
     return BlowersMasel(A, b, Ta0, w);
 }
 
+bool detectEfficiencies(ThreeBodyReaction& R)
+{
+    for (const auto& reac : R.reactants) {
+        // detect explicitly specified collision partner
+        if (R.products.count(reac.first)) {
+            R.third_body.efficiencies[reac.first] = 1.;
+        }
+    }
+
+    if (R.third_body.efficiencies.size() == 0) {
+        return false;
+    } else if (R.third_body.efficiencies.size() > 1) {
+        throw CanteraError("detectEfficiencies",
+            "Found more than one explicitly specified collision partner\n"
+            "in reaction '{}'.", R.equation());
+    }
+
+    R.third_body.default_efficiency = 0.;
+    R.specified = true;
+    auto sp = R.third_body.efficiencies.begin();
+
+    // adjust reactant coefficients
+    auto reac = R.reactants.find(sp->first);
+    if (abs(reac->second - 1) > SmallNumber) {
+        reac->second -= 1.;
+    } else {
+        R.reactants.erase(reac);
+    }
+
+    // adjust product coefficients
+    auto prod = R.products.find(sp->first);
+    if (abs(prod->second - 1) > SmallNumber) {
+        prod->second -= 1.;
+    } else {
+        R.products.erase(prod);
+    }
+
+    return true;
+}
+
 void setupReaction(Reaction& R, const XML_Node& rxn_node)
 {
     // Reactant and product stoichiometries
@@ -1000,20 +1059,26 @@ void setupThreeBodyReaction(ThreeBodyReaction& R, const XML_Node& rxn_node)
 {
     readEfficiencies(R.third_body, rxn_node.child("rateCoeff"));
     setupElementaryReaction(R, rxn_node);
+    if (R.third_body.efficiencies.size() == 0) {
+        detectEfficiencies(R);
+    }
 }
 
 void setupThreeBodyReaction(ThreeBodyReaction& R, const AnyMap& node,
                             const Kinetics& kin)
 {
     setupElementaryReaction(R, node, kin);
     if (R.reactants.count("M") != 1 || R.products.count("M") != 1) {
-        throw InputFileError("setupThreeBodyReaction", node["equation"],
-            "Reaction equation '{}' does not contain third body 'M'",
-            node["equation"].asString());
+        if (!detectEfficiencies(R)) {
+            throw InputFileError("setupThreeBodyReaction", node["equation"],
+                "Reaction equation '{}' does not contain third body 'M'",
+                node["equation"].asString());
+        }
+    } else {
+        R.reactants.erase("M");
+        R.products.erase("M");
+        readEfficiencies(R.third_body, node);
     }
-    R.reactants.erase("M");
-    R.products.erase("M");
-    readEfficiencies(R.third_body, node);
 }
 
 void setupFalloffReaction(FalloffReaction& R, const XML_Node& rxn_node)

src/kinetics/ReactionFactory.cpp

@@ -159,6 +159,44 @@ ReactionFactory::ReactionFactory()
                });
 }
 
+bool isThreeBody(Reaction& R)
+{
+    // detect explicitly specified collision partner
+    size_t found = 0;
+    for (const auto& reac : R.reactants) {
+        auto prod = R.products.find(reac.first);
+        if (prod != R.products.end() &&
+            trunc(reac.second) == reac.second && trunc(prod->second) == prod->second) {
+            // candidate species with integer stoichiometric coefficients on both sides
+            found += 1;
+        }
+    }
+    if (found != 1) {
+        return false;
+    }
+
+    // ensure that all reactants have integer stoichiometric coefficients
+    size_t nreac = 0;
+    for (const auto& reac : R.reactants) {
+       if (trunc(reac.second) != reac.second) {
+           return false;
+       }
+       nreac += reac.second;
+    }
+
+    // ensure that all products have integer stoichiometric coefficients
+    size_t nprod = 0;
+    for (const auto& prod : R.products) {
+       if (trunc(prod.second) != prod.second) {
+           return false;
+       }
+       nprod += prod.second;
+    }
+
+    // either reactant or product side involves exactly three species
+    return (nreac == 3) || (nprod == 3);
+}
+
 bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin)
 {
     vector_fp e_counter(kin.nPhases(), 0.0);
@@ -212,6 +250,16 @@ unique_ptr<Reaction> newReaction(const XML_Node& rxn_node)
         throw CanteraError("newReaction",
             "Unknown reaction type '" + rxn_node["type"] + "'");
     }
+    if (type == "") {
+        // Reaction type is not specified
+        // See if this is a three-body reaction with a specified collision partner
+        ElementaryReaction testReaction;
+        setupReaction(testReaction, rxn_node);
+        if (isThreeBody(testReaction)) {
+            type = "three-body";
+            R = ReactionFactory::factory()->create(type);
+        }
+    }
     if (type != "electrochemical") {
         type = R->type();
     }
@@ -226,6 +274,14 @@ unique_ptr<Reaction> newReaction(const AnyMap& rxn_node,
     std::string type = "elementary";
     if (rxn_node.hasKey("type")) {
         type = rxn_node["type"].asString();
+    } else if (kin.thermo().nDim() == 3) {
+        // Reaction type is not specified
+        // See if this is a three-body reaction with a specified collision partner
+        ElementaryReaction testReaction;
+        parseReactionEquation(testReaction, rxn_node["equation"], kin);
+        if (isThreeBody(testReaction)) {
+            type = "three-body";
+        }
     }
 
     if (kin.thermo().nDim() < 3 && type == "elementary") {

test_problems/clib_test/output_blessed.txt

@@ -107,15 +107,15 @@
                 CO + O2 <=> CO2 + O   1.185678e-05
             CH2O + O2 <=> HCO + HO2   2.924612e-12
              H + O2 + M <=> HO2 + M   4.077955e-05
-              H + 2 O2 <=> HO2 + O2   9.267647e-07
-         H + H2O + O2 <=> H2O + HO2   4.308072e-04
-           H + N2 + O2 <=> HO2 + N2   1.003166e-04
-           AR + H + O2 <=> AR + HO2   0.000000e+00
+           H + O2 + O2 <=> HO2 + O2   9.267647e-07
+         H + O2 + H2O <=> HO2 + H2O   4.308072e-04
+           H + O2 + N2 <=> HO2 + N2   1.003166e-04
+           H + O2 + AR <=> HO2 + AR   0.000000e+00
                   H + O2 <=> O + OH   9.566542e-02
                  2 H + M <=> H2 + M   3.752098e-06
-                  2 H + H2 <=> 2 H2   6.185839e-08
+               2 H + H2 <=> H2 + H2   6.185839e-08
              2 H + H2O <=> H2 + H2O   8.033966e-06
-             CO2 + 2 H <=> CO2 + H2   1.048274e-07
+             2 H + CO2 <=> H2 + CO2   1.048274e-07
              H + OH + M <=> H2O + M   3.496064e-04
                 H + HO2 <=> H2O + O   1.093006e-05
                 H + HO2 <=> H2 + O2   1.120905e-04
@@ -239,7 +239,7 @@
          CH3 + CH3OH <=> CH3O + CH4   7.138914e-32
           C2H4 + CH3 <=> C2H3 + CH4   3.146094e-40
           C2H6 + CH3 <=> C2H5 + CH4   3.552777e-46
-         H2O + HCO <=> CO + H + H2O   1.106338e-05
+         HCO + H2O <=> CO + H + H2O   1.106338e-05
              HCO + M <=> CO + H + M   7.088738e-06
               HCO + O2 <=> CO + HO2   5.164787e-07
           CH2OH + O2 <=> CH2O + HO2   1.869986e-14
@@ -376,7 +376,7 @@
         CH3CHO + H => CH3 + CO + H2   1.329705e-24
       CH3CHO + OH => CH3 + CO + H2O   5.882935e-24
     CH3CHO + HO2 => CH3 + CO + H2O2   1.642887e-29
-     CH3 + CH3CHO => CH3 + CH4 + CO   6.720371e-39
+     CH3CHO + CH3 => CH4 + CO + CH3   6.720371e-39
      CH2CO + H (+M) <=> CH2CHO (+M)   1.464625e-20
         CH2CHO + O => CH2 + CO2 + H   1.178961e-22
       CH2CHO + O2 => CH2O + CO + OH   1.668470e-25