Update three body logic

include/cantera/kinetics/Reaction.h

@@ -195,7 +195,10 @@ class Reaction
     bool m_valid = true;
 
     //! Flag indicating that serialization uses explicit type
-    bool m_explicit_rate = false;
+    bool m_explicit_type = false;
+
+    //! Flag indicating that third body collider is ambiguous
+    bool m_explicit_3rd = false;
 
     //! Flag indicating that object was instantiated from reactant/product compositions
     bool m_from_composition = false;
@@ -239,8 +242,10 @@ class ThirdBody
     void setParameters(const AnyMap& node);
 
     //! Get third-body efficiencies from AnyMap *node*
+    //! @param node  AnyMap receiving serialized parameters
+    //! @param explicit_3rd  Flag triggering explicit third body efficiency output
     //! @since  New in Cantera 3.0
-    void getParameters(AnyMap& node) const;
+    void getParameters(AnyMap& node, bool explicit_3rd=false) const;
 
     //! Get the third-body efficiency for species *k*
     double efficiency(const std::string& k) const;

interfaces/cython/cantera/reaction.pyx

@@ -1146,7 +1146,7 @@ cdef class Reaction:
 
     def __cinit__(self, reactants=None, products=None, rate=None, *,
                   equation=None, init=True, efficiencies=None,
-                  Kinetics kinetics=None, third_body=None, **kwargs):
+                  Kinetics kinetics=None, third_body=None):
         if kinetics:
             warnings.warn(
                 "Parameter 'kinetics' is no longer used and will be removed after "

src/kinetics/Reaction.cpp

@@ -37,6 +37,22 @@ Reaction::Reaction(const Composition& reactants_,
     , m_third_body(tbody_)
 {
     setRate(rate_);
+
+    // ensure safe serialization
+    bool count = 0;
+    for (const auto& reac : reactants) {
+        if (products.count(reac.first)) {
+            count = true;
+            break;
+        }
+    }
+    if (count) {
+        if (tbody_ && tbody_->name() != "M") {
+            m_explicit_3rd = true;
+        } else if (!tbody_) {
+            m_explicit_type = true;
+        }
+    }
 }
 
 Reaction::Reaction(const std::string& equation,
@@ -46,6 +62,9 @@ Reaction::Reaction(const std::string& equation,
 {
     setEquation(equation);
     setRate(rate_);
+    if (tbody_ && tbody_->name() != "M") {
+        m_explicit_3rd = true;
+    }
 }
 
 Reaction::Reaction(const AnyMap& node, const Kinetics& kin)
@@ -59,7 +78,13 @@ Reaction::Reaction(const AnyMap& node, const Kinetics& kin)
     setParameters(node, kin);
     size_t nDim = kin.thermo(kin.reactionPhaseIndex()).nDim();
     if (nDim == 3) {
-        setRate(newReactionRate(node, calculateRateCoeffUnits(kin)));
+        if (ba::starts_with(rate_type, "three-body-")) {
+            AnyMap rateNode = node;
+            rateNode["type"] = rate_type.substr(11, rate_type.size() - 11);
+            setRate(newReactionRate(rateNode, calculateRateCoeffUnits(kin)));
+        } else {
+            setRate(newReactionRate(node, calculateRateCoeffUnits(kin)));
+        }
     } else {
         AnyMap rateNode = node;
         if (rateNode.hasKey("rate-constant")) {
@@ -162,24 +187,26 @@ void Reaction::getParameters(AnyMap& reactionNode) const
     reactionNode.update(m_rate->parameters());
 
     // strip information not needed for reconstruction
-    std::string type = reactionNode["type"].asString();
-    if (type == "pressure-dependent-Arrhenius") {
+    std::string rtype = reactionNode["type"].asString();
+    if (rtype == "pressure-dependent-Arrhenius") {
         // skip
-    } else if (m_explicit_rate && ba::ends_with(type, "Arrhenius")) {
+    } else if (m_explicit_type && ba::ends_with(rtype, "Arrhenius")) {
         // retain type information
         if (m_third_body) {
             reactionNode["type"] = "three-body";
         } else {
             reactionNode["type"] = "elementary";
         }
-    } else if (ba::ends_with(type, "Arrhenius")) {
+    } else if (ba::ends_with(rtype, "Arrhenius")) {
         reactionNode.erase("type");
-    } else if (ba::ends_with(type, "Blowers-Masel")) {
+    } else if (m_explicit_type) {
+        reactionNode["type"] = type();
+    } else if (ba::ends_with(rtype, "Blowers-Masel")) {
         reactionNode["type"] = "Blowers-Masel";
     }
 
     if (m_third_body) {
-        m_third_body->getParameters(reactionNode);
+        m_third_body->getParameters(reactionNode, m_explicit_3rd);
     }
 }
 
@@ -210,7 +237,14 @@ void Reaction::setParameters(const AnyMap& node, const Kinetics& kin)
     allow_nonreactant_orders = node.getBool("nonreactant-orders", false);
 
     if (m_third_body) {
-        m_third_body->setParameters(node);
+        try {
+            m_third_body->setParameters(node);
+            if (m_third_body->name() == "M" && m_third_body->efficiencies.size() == 1) {
+                m_explicit_3rd = true;
+            }
+        } catch (CanteraError& err) {
+            throw InputFileError("Reaction::setParameters", input, err.getMessage());
+        }
     } else if (node.hasKey("default-efficiency") || node.hasKey("efficiencies")) {
         throw InputFileError("Reaction::setParameters", input,
             "Reaction '{}' specifies efficiency parameters\n"
@@ -310,12 +344,12 @@ void Reaction::setEquation(const std::string& equation, const Kinetics* kin)
     parseReactionEquation(*this, equation, input, kin);
 
     std::string rate_type = input.getString("type", "");
-    if (rate_type == "three-body") {
+    if (ba::starts_with(rate_type, "three-body")) {
         // state type when serializing
-        m_explicit_rate = true;
+        m_explicit_type = true;
     } else if (rate_type == "elementary") {
         // user override
-        m_explicit_rate = true;
+        m_explicit_type = true;
         return;
     } else if (kin && kin->thermo(kin->reactionPhaseIndex()).nDim() != 3) {
         // interface reactions
@@ -341,7 +375,7 @@ void Reaction::setEquation(const std::string& equation, const Kinetics* kin)
     }
 
     if (count == 0) {
-        if (rate_type == "three-body") {
+        if (ba::starts_with(rate_type, "three-body")) {
             throw InputFileError("Reaction::setEquation", input,
                 "Reactants for reaction '{}'\n"
                 "do not contain a valid third body collider", equation);
@@ -355,12 +389,48 @@ void Reaction::setEquation(const std::string& equation, const Kinetics* kin)
     } else if (count > 1) {
         // equations with more than one explicit third-body collider are handled as a
         // regular elementary reaction unless the equation contains a generic third body
-        if (!countM) {
+        if (countM) {
+            // generic collider 'M' is selected as third body
+            third_body = "M";
+        } else if (m_third_body) {
+            // third body is defined as explicit object
+            auto& effs = m_third_body->efficiencies;
+            if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
+                throw InputFileError("Reaction::setEquation", input,
+                    "Detected ambiguous third body colliders in reaction '{}'\n"
+                    "ThirdBody object needs to specify a single species", equation);
+            }
+            third_body = effs.begin()->first;
+            m_explicit_3rd = true;
+        } else if (input.hasKey("efficiencies") && input.hasKey("default-efficiency")) {
+            // third body is implicitly defined by efficiency
+            auto effs = input["efficiencies"].asMap<double>();
+            if (effs.size() != 1 || !reactants.count(effs.begin()->first)) {
+                throw InputFileError("Reaction::setEquation", input,
+                    "Detected ambiguous third body colliders in reaction '{}'\n"
+                    "Collision efficiencies need to specify single species", equation);
+            }
+            third_body = effs.begin()->first;
+            m_explicit_3rd = true;
+        } else if (input.hasKey("efficiencies") || input.hasKey("default-efficiency")) {
+            // insufficient disambiguation of third bodies
+            throw InputFileError("Reaction::setEquation", input,
+                "Detected ambiguous third body colliders in reaction '{}'\n"
+                "Third-body definition requires specification of efficiencies "
+                "as well as default efficiency", equation);
+        } else if (ba::starts_with(rate_type, "three-body")) {
+            // no disambiguation of third bodies
+            throw InputFileError("Reaction::setEquation", input,
+                "Detected ambiguous third body colliders in reaction '{}'\n"
+                "A valid ThirdBody or collision efficiency definition is required",
+                equation);
+        } else {
             return;
         }
-        third_body = "M";
 
-    } else if (!ba::starts_with(third_body, "(+")) {
+    } else if (third_body != "M" && !ba::starts_with(rate_type, "three-body")
+            && !ba::starts_with(third_body, "(+"))
+    {
         // check for conditions of three-body reactions:
         // - integer stoichiometric conditions
         // - either reactant or product side involves exactly three species
@@ -390,6 +460,12 @@ void Reaction::setEquation(const std::string& equation, const Kinetics* kin)
     }
 
     if (m_third_body) {
+        std::string tName = m_third_body->name();
+        if (tName != third_body && third_body != "M" && tName != "M") {
+            throw InputFileError("Reaction::setEquation", input,
+                "Detected incompatible third body colliders in reaction '{}'\n"
+                "ThirdBody definition does not match equation", equation);
+        }
         m_third_body->setName(third_body);
     } else {
         m_third_body.reset(new ThirdBody(third_body));
@@ -656,9 +732,10 @@ void ThirdBody::setName(const std::string& third_body)
     if (name == m_name) {
         return;
     }
-    if (name == "M") {
-        throw CanteraError("ThirdBody::setName",
-            "Unable to revert explicit third body '{}' to 'M'", m_name);
+    if (name == "M" && efficiencies.size() == 1) {
+        // revert from explicit name to generic collider
+        m_name = name;
+        return;
     }
     if (efficiencies.size()) {
         throw CanteraError("ThirdBody::setName",
@@ -689,11 +766,17 @@ void ThirdBody::setParameters(const AnyMap& node)
     if (node.hasKey("efficiencies")) {
         efficiencies = node["efficiencies"].asMap<double>();
     }
+    if (m_name != "M"
+        && (efficiencies.size() != 1 || efficiencies.begin()->first != m_name))
+    {
+        throw CanteraError("Reaction::setEquation",
+            "Detected incompatible third body colliders definitions");
+    }
 }
 
-void ThirdBody::getParameters(AnyMap& node) const
+void ThirdBody::getParameters(AnyMap& node, bool explicit_3rd) const
 {
-    if (m_name == "M") {
+    if (m_name == "M" || explicit_3rd) {
         if (efficiencies.size()) {
             node["efficiencies"] = efficiencies;
             node["efficiencies"].setFlowStyle();

test/kinetics/kineticsFromScratch.cpp

@@ -149,6 +149,89 @@ TEST_F(KineticsFromScratch, multiple_third_bodies3)
     EXPECT_TRUE(R->usesThirdBody());
 }
 
+TEST_F(KineticsFromScratch, multiple_third_bodies4)
+{
+    std::string equation = "H2 + O2 => H2 + O2";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto tbody = make_shared<ThirdBody>("O2");
+    auto R = make_shared<Reaction>(equation, rate, tbody);
+    EXPECT_EQ(R->thirdBody()->name(), "O2");
+
+    AnyMap input = R->parameters();
+    EXPECT_FALSE(input.hasKey("type"));
+    EXPECT_TRUE(input.hasKey("efficiencies"));
+    auto efficiencies = input["efficiencies"].asMap<double>();
+    EXPECT_EQ(efficiencies.size(), 1);
+    EXPECT_EQ(efficiencies.begin()->first, "O2");
+    EXPECT_TRUE(input.hasKey("default-efficiency"));
+    EXPECT_EQ(input["default-efficiency"].asDouble(), 0.);
+}
+
+TEST_F(KineticsFromScratch, multiple_third_bodies5)
+{
+    std::string equation = "H2 + O2 => H2 + O2";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto tbody = make_shared<ThirdBody>("AR"); // incompatible third body
+    ASSERT_THROW(Reaction(equation, rate, tbody), CanteraError);
+}
+
+TEST_F(KineticsFromScratch, multiple_third_bodies6)
+{
+    Composition reac = parseCompString("H2:1");
+    Composition prod = parseCompString("H2:1");
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto tbody = make_shared<ThirdBody>("O2");
+    auto R = make_shared<Reaction>(reac, prod, rate, tbody);
+    EXPECT_EQ(R->thirdBody()->name(), "O2");
+
+    AnyMap input = R->parameters();
+    EXPECT_FALSE(input.hasKey("type"));
+    EXPECT_TRUE(input.hasKey("efficiencies"));
+    auto efficiencies = input["efficiencies"].asMap<double>();
+    EXPECT_EQ(efficiencies.size(), 1);
+    EXPECT_EQ(efficiencies.begin()->first, "O2");
+    EXPECT_TRUE(input.hasKey("default-efficiency"));
+    EXPECT_EQ(input["default-efficiency"].asDouble(), 0.);
+}
+
+TEST_F(KineticsFromScratch, multiple_third_bodies7)
+{
+    std::string equation = "CH2OCH + M <=> CH2CHO + M";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto tbody = make_shared<ThirdBody>("O2");
+    auto R = make_shared<Reaction>(equation, rate, tbody);
+    EXPECT_EQ(R->thirdBody()->name(), "M");
+
+    AnyMap input = R->parameters();
+    EXPECT_FALSE(input.hasKey("type"));
+    EXPECT_TRUE(input.hasKey("efficiencies"));
+    auto efficiencies = input["efficiencies"].asMap<double>();
+    EXPECT_EQ(efficiencies.size(), 1);
+    EXPECT_EQ(efficiencies.begin()->first, "O2");
+    EXPECT_TRUE(input.hasKey("default-efficiency"));
+    EXPECT_EQ(input["default-efficiency"].asDouble(), 0.);
+}
+
+TEST_F(KineticsFromScratch, multiple_third_bodies8)
+{
+    std::string equation = "CH2OCH + M <=> CH2CHO + M";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto tbody = make_shared<ThirdBody>();
+    tbody->efficiencies = parseCompString("O2:1");
+    tbody->default_efficiency = 0.;
+    auto R = make_shared<Reaction>(equation, rate, tbody);
+    EXPECT_EQ(R->thirdBody()->name(), "M");
+
+    AnyMap input = R->parameters();
+    EXPECT_FALSE(input.hasKey("type"));
+    EXPECT_TRUE(input.hasKey("efficiencies"));
+    auto efficiencies = input["efficiencies"].asMap<double>();
+    EXPECT_EQ(efficiencies.size(), 1);
+    EXPECT_EQ(efficiencies.begin()->first, "O2");
+    EXPECT_TRUE(input.hasKey("default-efficiency"));
+    EXPECT_EQ(input["default-efficiency"].asDouble(), 0.);
+}
+
 TEST_F(KineticsFromScratch, add_two_temperature_plasma)
 {
     std::string equation = "O + H => O + H";