Fix issues for undeclared third-bodies and implicit third-body reactions

include/cantera/kinetics/Kinetics.h

@@ -1625,6 +1625,9 @@ class Kinetics
     //! See skipUndeclaredThirdBodies()
     bool m_skipUndeclaredThirdBodies = false;
 
+    //! Flag indicating whether reactions include undeclared third bodies
+    bool m_hasUndeclaredThirdBodies = false;
+
     //! reference to Solution
     std::weak_ptr<Solution> m_root;
 };

interfaces/cython/cantera/kinetics.pxd

@@ -35,8 +35,8 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera":
         void addThermo(shared_ptr[CxxThermoPhase]) except +translate_exception
         void init() except +translate_exception
         void skipUndeclaredThirdBodies(cbool)
-        void addReaction(shared_ptr[CxxReaction]) except +translate_exception
-        void addReaction(shared_ptr[CxxReaction], cbool) except +translate_exception
+        cbool addReaction(shared_ptr[CxxReaction]) except +translate_exception
+        cbool addReaction(shared_ptr[CxxReaction], cbool) except +translate_exception
         void modifyReaction(int, shared_ptr[CxxReaction]) except +translate_exception
         void invalidateCache() except +translate_exception
         void resizeReactions()

interfaces/cython/cantera/reaction.pyx

@@ -1502,10 +1502,19 @@
         species names and the values, are the stoichiometric coefficients, for example
         ``{'CH4':1, 'OH':1}``, or as a composition string, for example
         ``'CH4:1, OH:1'``.
+
+        .. deprecated:: 3.0
+
+            Setter for reactants is deprecated and will be removed after Cantera 3.0.
+            Use constructor instead.
         """
         def __get__(self):
             return comp_map_to_dict(self.reaction.reactants)
         def __set__(self, reactants):
+            warnings.warn(
+                "'Reaction.reactants' setter is deprecated and will be removed after "
+                "Cantera 3.0.\nInstantiate using constructor instead.",
+                DeprecationWarning)
             self.reaction.reactants = comp_map(reactants)
 
     property products:
@@ -1514,10 +1523,19 @@
         species names and the values, are the stoichiometric coefficients, for example
         ``{'CH3':1, 'H2O':1}``, or as a composition string, for example
         ``'CH3:1, H2O:1'``.
+
+        .. deprecated:: 3.0
+
+            Setter for products is deprecated and will be removed after Cantera 3.0.
+            Use constructor instead.
         """
         def __get__(self):
             return comp_map_to_dict(self.reaction.products)
         def __set__(self, products):
+            warnings.warn(
+                "'Reaction.products' setter is deprecated and will be removed after "
+                "Cantera 3.0.\nInstantiate using constructor instead.",
+                DeprecationWarning)
             self.reaction.products = comp_map(products)
 
     def __contains__(self, species):
@@ -1655,6 +1673,18 @@
             if self.reaction.usesThirdBody():
                 return ThirdBody.wrap(self.reaction.thirdBody())
 
+    @property
+    def third_body_name(self):
+        """
+        Returns name of `ThirdBody` collider if `Reaction` uses a third body collider,
+        and ``None`` otherwise.
+
+        .. versionadded:: 3.0
+        """
+        if self.reaction.usesThirdBody():
+            return ThirdBody.wrap(self.reaction.thirdBody()).name
+        return None
+
     property efficiencies:
         """
         Get/Set a `dict` defining non-default third-body efficiencies for this reaction,

src/kinetics/BulkKinetics.cpp

@@ -88,6 +88,8 @@ void BulkKinetics::addThirdBody(shared_ptr<Reaction> r)
             throw CanteraError("BulkKinetics::addThirdBody", "Found third-body"
                 " efficiency for undefined species '{}' while adding reaction '{}'",
                 name, r->equation());
+        } else {
+            m_hasUndeclaredThirdBodies = true;
         }
     }
     m_multi_concm.install(nReactions() - 1, efficiencies,

src/kinetics/Kinetics.cpp

@@ -643,6 +643,9 @@ AnyMap Kinetics::parameters()
         if (nReactions() == 0) {
             out["reactions"] = "none";
         }
+        if (m_hasUndeclaredThirdBodies) {
+            out["skip-undeclared-third-bodies"] = true;
+        }
     }
     return out;
 }

src/kinetics/Reaction.cpp

@@ -35,21 +35,35 @@ Reaction::Reaction(const Composition& reactants_,
     , m_from_composition(true)
     , m_third_body(tbody_)
 {
+    if (reactants.count("M") || products.count("M")) {
+        throw CanteraError("Reaction::Reaction",
+            "Third body 'M' must not be included in either reactant or product maps.");
+    }
     setRate(rate_);
 
     // set flags ensuring correct serialization output
-    bool count = 0;
+    Composition third;
     for (const auto& [name, stoich] : reactants) {
         if (products.count(name)) {
-            count = true;
-            break;
+            third[name] = products.at(name) - stoich;
         }
     }
-    if (count) {
-        if (tbody_ && tbody_->name() != "M") {
+    if (tbody_) {
+        string name = tbody_->name();
+        if (reactants.count(name) && products.count(name)) {
+            throw CanteraError("Reaction::Reaction",
+                "'{}' not acting as third body collider must not be included in both "
+                "reactant and product maps.", name);
+        }
+        if (name != "M") {
+            m_third_body->explicit_3rd = true;
+        }
+    } else if (!tbody_ && third.size() == 1) {
+        // implicit third body
+        string name = third.begin()->first;
+        m_third_body = make_shared<ThirdBody>(name);
+        if (name != "M") {
             m_third_body->explicit_3rd = true;
-        } else if (!tbody_) {
-            m_explicit_type = true;
         }
     }
 }
@@ -61,7 +75,7 @@ Reaction::Reaction(const string& equation,
 {
     setEquation(equation);
     setRate(rate_);
-    if (tbody_ && tbody_->name() != "M") {
+    if (m_third_body && m_third_body->name() != "M") {
         m_third_body->explicit_3rd = true;
     }
 }
@@ -822,7 +836,9 @@ bool ThirdBody::checkSpecies(const Reaction& rxn, const Kinetics& kin) const
             throw InputFileError("ThirdBody::checkSpecies", rxn.input, "Reaction '{}'\n"
                 "is a three-body reaction with undeclared species: '{}'",
                 rxn.equation(), ba::join(undeclared, "', '"));
-        } else if (kin.skipUndeclaredSpecies() && m_name != "M") {
+        } else if (kin.skipUndeclaredThirdBodies() && m_name != "M") {
+            // Prevent addition of reaction silently as "skip-undeclared-third-bodies"
+            // is set to true
             return false;
         }
     }

test/kinetics/kineticsFromScratch.cpp

@@ -102,6 +102,11 @@ TEST_F(KineticsFromScratch, add_three_body_reaction1)
 
     kin.addReaction(R);
     check_rates(1);
+
+    reac = parseCompString("O:2, M:1");
+    prod = parseCompString("O2:1, M:1");
+    ASSERT_THROW(make_shared<Reaction>(reac, prod, rate, tbody), CanteraError);
+
 }
 
 TEST_F(KineticsFromScratch, add_three_body_reaction2)
@@ -116,6 +121,8 @@ TEST_F(KineticsFromScratch, add_three_body_reaction2)
     auto tbody = make_shared<ThirdBody>();
     tbody->efficiencies = parseCompString("AR:0.83 H2:2.4 H2O:15.4");
     auto R = make_shared<Reaction>(equation, rate, tbody);
+    auto reac = R->reactants;
+    EXPECT_EQ(reac.count("M"), (size_t) 0);
 
     kin.addReaction(R);
     check_rates(1);
@@ -160,6 +167,9 @@ TEST_F(KineticsFromScratch, multiple_third_bodies4)
     auto R = make_shared<Reaction>(equation, rate, tbody);
     EXPECT_EQ(R->thirdBody()->name(), "O2");
     EXPECT_EQ(R->thirdBody()->default_efficiency, 0.);
+    EXPECT_EQ(R->reactants.count("H2"), (size_t) 1);
+    EXPECT_EQ(R->reactants.count("O2"), (size_t) 0);
+    EXPECT_EQ(R->reactants.count("M"), (size_t) 0);
 
     AnyMap input = R->parameters(false);
     EXPECT_FALSE(input.hasKey("type"));
@@ -187,6 +197,9 @@ TEST_F(KineticsFromScratch, multiple_third_bodies6)
     auto R = make_shared<Reaction>(reac, prod, rate, tbody);
     EXPECT_EQ(R->thirdBody()->name(), "O2");
     EXPECT_EQ(R->thirdBody()->default_efficiency, 0.);
+    EXPECT_EQ(R->reactants.count("H2"), (size_t) 1);
+    EXPECT_EQ(R->reactants.count("O2"), (size_t) 0);
+    EXPECT_EQ(R->reactants.count("M"), (size_t) 0);
 
     AnyMap input = R->parameters(false);
     EXPECT_FALSE(input.hasKey("type"));
@@ -235,6 +248,24 @@ TEST_F(KineticsFromScratch, multiple_third_bodies8)
     EXPECT_EQ(input["default-efficiency"].asDouble(), 0.);
 }
 
+TEST_F(KineticsFromScratch, multiple_third_bodies9)
+{
+    Composition reac = parseCompString("H2:1, O2:1");
+    Composition prod = parseCompString("H2:1, O:2");
+    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");
+    EXPECT_EQ(R->thirdBody()->default_efficiency, 0.);
+    EXPECT_EQ(R->reactants.count("H2"), (size_t) 1);
+    EXPECT_EQ(R->reactants.count("O2"), (size_t) 1);
+    EXPECT_EQ(R->reactants.count("M"), (size_t) 0);
+
+    reac = parseCompString("H2:1, O2:1");
+    prod = parseCompString("H2:1, O2:1");
+    ASSERT_THROW(make_shared<Reaction>(reac, prod, rate, tbody), CanteraError);
+}
+
 TEST_F(KineticsFromScratch, add_two_temperature_plasma)
 {
     string equation = "O + H => O + H";
@@ -269,6 +300,35 @@ TEST_F(KineticsFromScratch, skip_undefined_third_body)
     ASSERT_EQ((size_t) 1, kin.nReactions());
 }
 
+TEST_F(KineticsFromScratch, skip_explicit_third_body)
+{
+    string equation = "2 O + CO2 <=> O2 + CO2";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto R = make_shared<Reaction>(equation, rate);
+    EXPECT_EQ(R->thirdBody()->name(), "CO2");
+
+    ASSERT_THROW(kin.addReaction(R), CanteraError);
+    kin.skipUndeclaredThirdBodies(true);
+    kin.addReaction(R);
+    ASSERT_EQ((size_t) 0, kin.nReactions());
+}
+
+TEST_F(KineticsFromScratch, third_body_composition)
+{
+    string equation = "2 O + H2O <=> O2 + H2O";
+    auto rate = make_shared<ArrheniusRate>(1.2e11, -1.0, 0.0);
+    auto R = make_shared<Reaction>(equation, rate);
+    EXPECT_EQ(R->thirdBody()->name(), "H2O");
+    EXPECT_TRUE(R->thirdBody()->explicit_3rd);
+
+    Composition reac = R->reactants;
+    EXPECT_EQ(reac.count("H2O"), (size_t) 0);
+    EXPECT_EQ(reac.count("M"), (size_t) 0);
+    Composition prod = R->products;
+    EXPECT_EQ(prod.count("H2O"), (size_t) 0);
+    EXPECT_EQ(reac.count("M"), (size_t) 0);
+}
+
 TEST_F(KineticsFromScratch, add_falloff_reaction1)
 {
     // reaction 2:

test/python/test_kinetics.py

@@ -304,7 +304,7 @@ def test_coverage_dependence_flags(self):
         surf.net_rates_of_progress_ddCi
 
     def test_electrochemistry_flags(self):
-             # Phases
+        # Phases
         mech = "lithium_ion_battery.yaml"
         anode, cathode, metal, electrolyte = ct.import_phases(
             mech, ["anode", "cathode", "electron", "electrolyte"])
@@ -315,6 +315,41 @@ def test_electrochemistry_flags(self):
         anode_int.derivative_settings = {"skip-electrochemistry": True}
         anode_int.net_rates_of_progress_ddCi
 
+    def test_submechanism(self):
+        # Simplified samples/python/kinetics/extract_submechanism.py
+        gri30 = ct.Solution('gri30.yaml', transport_model=None)
+        h2o2 = ct.Solution('h2o2.yaml', transport_model=None)
+
+        reactions_plus = []
+        reactions = []
+        dest_species = set(h2o2.species_names)
+        colliders = dest_species.union([None, "M"])
+        for R in gri30.reactions():
+            if not all(S in dest_species for S in R.reactants):
+                continue
+            if not all(S in dest_species for S in R.products):
+                continue
+            reactions_plus.append(R)
+            if R.third_body_name not in colliders:
+                continue
+            reactions.append(R)
+
+        gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
+                          species=h2o2.species(), reactions=reactions_plus)
+        # there is one third-body reaction with an undeclared third body species
+        gas.n_reactions < len(reactions_plus)
+        assert gas.n_species == len(h2o2.species_names)
+
+        gas = ct.Solution(thermo='ideal-gas', kinetics='gas', species=h2o2.species(),
+                          reactions=reactions)
+        assert gas.n_reactions == len(reactions)
+        assert gas.n_species == len(h2o2.species_names)
+
+        gas.write_yaml("reduced.yaml")
+        restored = ct.Solution("reduced.yaml")
+        assert gas.species_names == restored.species_names
+        assert gas.reaction_equations() == restored.reaction_equations()
+
 
 class KineticsRepeatability(utilities.CanteraTest):
     """