Update detection of skipped/undeclared species

include/cantera/kinetics/Kinetics.h

@@ -857,6 +857,10 @@ class Kinetics
     //! Check that the specified reaction is balanced (same number of atoms for
     //! each element in the reactants and products). Raises an exception if the
     //! reaction is not balanced.
+    /*!
+     * @deprecated To be removed in Cantera 2.6.
+     *             Replaceable by Reaction::checkBalance.
+     */
     void checkReactionBalance(const Reaction& R);
 
     //! @name Stoichiometry management

include/cantera/kinetics/Reaction.h

@@ -75,6 +75,19 @@ class Reaction
         m_valid = valid;
     }
 
+    //! Check that the specified reaction is balanced (same number of atoms for
+    //! each element in the reactants and products). Raises an exception if the
+    //! reaction is not balanced. Used by checkSpecies.
+    //! @param kin  Kinetics object
+    void checkBalance(const Kinetics& kin) const;
+
+    //! Verify that all species involved in the reaction are defined in the Kinetics
+    //! object. The function returns true if all species are found, and raises an
+    //! exception unless the kinetics object is configured to skip undeclared species,
+    //! in which case false is returned.
+    //! @param kin  Kinetics object
+    bool checkSpecies(const Kinetics& kin) const;
+
     //! Type of the reaction. The valid types are listed in the file,
     //! reaction_defs.h, with constants ending in `RXN`.
     /*!
@@ -127,6 +140,14 @@ class Reaction
 
     //! Flag indicating whether reaction is set up correctly
     bool m_valid;
+
+    //! @internal  Helper function returning vector of undeclared third body species
+    //! and a boolean expression indicating whether the third body is specified.
+    //! The function is used by the checkSpecies method and only needed as long as
+    //! there is no unified approach to handle third body collision partners.
+    //! @param kin  Kinetics object
+    virtual std::pair<std::vector<std::string>, bool>
+        undeclaredThirdBodies(const Kinetics& kin) const;
 };
 
 
@@ -214,6 +235,10 @@ class ThreeBodyReaction : public ElementaryReaction
     ThirdBody third_body;
 
     bool specified_collision_partner = false; //!< Input specifies collision partner
+
+protected:
+    virtual std::pair<std::vector<std::string>, bool>
+        undeclaredThirdBodies(const Kinetics& kin) const;
 };
 
 //! A reaction that is first-order in [M] at low pressure, like a third-body
@@ -254,6 +279,10 @@ class FalloffReaction : public Reaction
     //! The units of the low-pressure rate constant. The units of the
     //! high-pressure rate constant are stored in #rate_units.
     Units low_rate_units;
+
+protected:
+    virtual std::pair<std::vector<std::string>, bool> undeclaredThirdBodies(
+        const Kinetics& kin) const;
 };
 
 //! A reaction where the rate decreases as pressure increases due to collisional

interfaces/cython/cantera/test/test_kinetics.py

@@ -413,7 +413,7 @@ def test_pdep_err(self):
                    "at P = 1.0132e+05, T = 1000.0",
                    "Invalid rate coefficient for reaction 'CH2CHOO <=> CH3 + CO2'",
                    "at P = 1.0132e+05, T = 500.0",
-                   "InputFileError thrown by Kinetics::checkReactionBalance:",
+                   "InputFileError thrown by Reaction::checkBalance:",
                    "The following reaction is unbalanced: H2O2 + OH <=> 2 H2O + HO2")
         try:
             ct.Solution('addReactions_err_test.yaml')
@@ -423,6 +423,115 @@ def test_pdep_err(self):
             for msg in err_msg:
                 self.assertIn(msg, err_msg_list)
 
+
+class TestUndeclared(utilities.CanteraTest):
+
+    _gas_def = """
+            phases:
+            - name: gas
+              species:
+              - h2o2.yaml/species: [H, O2, H2O, HO2]
+              thermo: ideal-gas
+              kinetics: gas
+            """
+
+    def test_raise_undeclared_species(self):
+
+        gas_def = self._gas_def + """
+            reactions:
+            - equation: O + H2 <=> H + OH  # Reaction 3
+              rate-constant: {A: 3.87e+04, b: 2.7, Ea: 6260.0}
+            """
+
+        with self.assertRaisesRegex(ct.CanteraError, "contains undeclared species"):
+            ct.Solution(yaml=gas_def)
+
+    def test_raise_undeclared_third_bodies(self):
+
+        gas_def = self._gas_def + """
+            reactions:
+            - equation: H + O2 + AR <=> HO2 + AR  # Reaction 10
+              rate-constant: {A: 7.0e+17, b: -0.8, Ea: 0.0}
+            """
+
+        with self.assertRaisesRegex(ct.CanteraError, "three-body reaction with undeclared"):
+            ct.Solution(yaml=gas_def)
+
+    def test_skip_undeclared_third_bodies1(self):
+
+        gas_def = self._gas_def + """
+              reactions:
+              - h2o2.yaml/reactions: declared-species
+              skip-undeclared-third-bodies: true
+            """
+
+        gas = ct.Solution(yaml=gas_def)
+        self.assertEqual(gas.n_reactions, 3)
+
+    def test_skip_undeclared_third_bodies2(self):
+
+        gas_def = """
+            phases:
+            - name: gas
+              species:
+              - h2o2.yaml/species: [H, O2, HO2]
+              thermo: ideal-gas
+              skip-undeclared-third-bodies: true
+              kinetics: gas
+              reactions:
+              - h2o2.yaml/reactions: declared-species
+            """
+
+        gas = ct.Solution(yaml=gas_def)
+        found = False
+        for rxn in gas.reactions():
+            if rxn.equation == "H + O2 + M <=> HO2 + M":
+                found = True
+                break
+        self.assertTrue(found)
+
+    def test_skip_undeclared_orders(self):
+
+        gas_def = self._gas_def + """
+            reactions:
+            - equation: H + O2 => HO2
+              rate-constant: {A: 1.255943e+13, b: -2.0, Ea: 5000.0 cal/mol}
+              orders:
+                H2O: 0.2
+              nonreactant-orders: true
+            """
+
+        gas = ct.Solution(yaml=gas_def)
+        self.assertEqual(gas.n_reactions, 1)
+
+    def test_raise_nonreactant_orders(self):
+
+        gas_def = self._gas_def + """
+            reactions:
+            - equation: H + O2 => HO2
+              rate-constant: {A: 1.255943e+13, b: -2.0, Ea: 5000.0 cal/mol}
+              orders:
+                H2O: 0.2
+            """
+
+        with self.assertRaisesRegex(ct.CanteraError, "Reaction order specified"):
+            ct.Solution(yaml=gas_def)
+
+    def test_raise_undeclared_orders(self):
+
+        gas_def = self._gas_def + """
+            reactions:
+            - equation: H + O2 => HO2
+              rate-constant: {A: 1.255943e+13, b: -2.0, Ea: 5000.0 cal/mol}
+              orders:
+                N2: 0.2
+              nonreactant-orders: true
+            """
+
+        with self.assertRaisesRegex(ct.CanteraError, "reaction orders for undeclared"):
+            ct.Solution(yaml=gas_def)
+
+
 class TestEmptyKinetics(utilities.CanteraTest):
     def test_empty(self):
         gas = ct.Solution('air-no-reactions.xml')

interfaces/cython/cantera/test/test_onedim.py

@@ -1253,7 +1253,7 @@ def test_ion_profile(self):
         width = 0.03
 
         # Solution object used to compute mixture properties
-        self.gas = ct.Solution('ch4_ion.cti')
+        self.gas = ct.Solution('ch4_ion.yaml')
         self.gas.TPX = Tin, p, reactants
         self.sim = ct.IonFreeFlame(self.gas, width=width)
         self.sim.set_refine_criteria(ratio=4, slope=0.8, curve=1.0)
@@ -1279,7 +1279,7 @@ def test_ion_profile(self):
         width = 0.01
 
         # Solution object used to compute mixture properties
-        self.gas = ct.Solution('ch4_ion.cti')
+        self.gas = ct.Solution('ch4_ion.yaml')
         self.gas.TPX = Tburner, p, reactants
         self.sim = ct.IonBurnerFlame(self.gas, width=width)
         self.sim.set_refine_criteria(ratio=4, slope=0.1, curve=0.1)

src/base/ctexceptions.cpp

@@ -14,7 +14,8 @@ namespace Cantera
 
 // *** Exceptions ***
 
-static const char* stars = "***********************************************************************\n";
+static const char* stars = ("*****************************************"
+                            "**************************************\n");
 
 CanteraError::CanteraError(const std::string& procedure) :
     procedure_(procedure)

src/kinetics/GasKinetics.cpp

@@ -288,10 +288,6 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r)
         size_t k = kineticsSpeciesIndex(eff.first);
         if (k != npos) {
             efficiencies[k] = eff.second;
-        } else if (!m_skipUndeclaredThirdBodies) {
-            throw CanteraError("GasKinetics::addFalloffReaction", "Found "
-                "third-body efficiency for undefined species '" + eff.first +
-                "' while adding reaction '" + r.equation() + "'");
         }
     }
     m_falloff_concm.install(nfall, efficiencies,
@@ -311,10 +307,6 @@ void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r)
         size_t k = kineticsSpeciesIndex(eff.first);
         if (k != npos) {
             efficiencies[k] = eff.second;
-        } else if (!m_skipUndeclaredThirdBodies) {
-            throw CanteraError("GasKinetics::addThreeBodyReaction", "Found "
-                "third-body efficiency for undefined species '" + eff.first +
-                "' while adding reaction '" + r.equation() + "'");
         }
     }
     m_3b_concm.install(nReactions()-1, efficiencies,

src/kinetics/Kinetics.cpp

@@ -17,6 +17,7 @@
 #include "cantera/base/utilities.h"
 #include "cantera/base/global.h"
 #include <unordered_set>
+#include <boost/algorithm/string/join.hpp>
 
 using namespace std;
 
@@ -231,44 +232,9 @@ double Kinetics::checkDuplicateStoich(std::map<int, double>& r1,
 
 void Kinetics::checkReactionBalance(const Reaction& R)
 {
-    Composition balr, balp;
-    // iterate over the products
-    for (const auto& sp : R.products) {
-        const ThermoPhase& ph = speciesPhase(sp.first);
-        size_t k = ph.speciesIndex(sp.first);
-        double stoich = sp.second;
-        for (size_t m = 0; m < ph.nElements(); m++) {
-            balr[ph.elementName(m)] = 0.0; // so that balr contains all species
-            balp[ph.elementName(m)] += stoich*ph.nAtoms(k,m);
-        }
-    }
-    for (const auto& sp : R.reactants) {
-        const ThermoPhase& ph = speciesPhase(sp.first);
-        size_t k = ph.speciesIndex(sp.first);
-        double stoich = sp.second;
-        for (size_t m = 0; m < ph.nElements(); m++) {
-            balr[ph.elementName(m)] += stoich*ph.nAtoms(k,m);
-        }
-    }
-
-    string msg;
-    bool ok = true;
-    for (const auto& el : balr) {
-        const string& elem = el.first;
-        double elemsum = balr[elem] + balp[elem];
-        double elemdiff = fabs(balp[elem] - balr[elem]);
-        if (elemsum > 0.0 && elemdiff/elemsum > 1e-4) {
-            ok = false;
-            msg += fmt::format("  {}           {}           {}\n",
-                               elem, balr[elem], balp[elem]);
-        }
-    }
-    if (!ok) {
-        throw InputFileError("Kinetics::checkReactionBalance", R.input,
-            "The following reaction is unbalanced: {}\n"
-            "  Element    Reactants    Products\n{}",
-            R.equation(), msg);
-    }
+    R.checkBalance(*this);
+    warn_deprecated("Kinetics::checkReactionBalance",
+        "To be removed after Cantera 2.6. Replacable by Reaction::checkBalance.");
 }
 
 void Kinetics::selectPhase(const double* data, const ThermoPhase* phase,
@@ -527,49 +493,10 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r)
     }
     resizeSpecies();
 
-    // If reaction orders are specified, then this reaction does not follow
-    // mass-action kinetics, and is not an elementary reaction. So check that it
-    // is not reversible, since computing the reverse rate from thermochemistry
-    // only works for elementary reactions.
-    if (r->reversible && !r->orders.empty()) {
-        throw InputFileError("Kinetics::addReaction", r->input,
-            "Reaction orders may only be given for irreversible reactions");
-    }
-
-    // Check for undeclared species
-    for (const auto& sp : r->reactants) {
-        if (kineticsSpeciesIndex(sp.first) == npos) {
-            if (m_skipUndeclaredSpecies) {
-                return false;
-            } else {
-                throw InputFileError("Kinetics::addReaction", r->input,
-                    "Reaction '{}' contains the undeclared species '{}'",
-                    r->equation(), sp.first);
-            }
-        }
-    }
-    for (const auto& sp : r->products) {
-        if (kineticsSpeciesIndex(sp.first) == npos) {
-            if (m_skipUndeclaredSpecies) {
-                return false;
-            } else {
-                throw InputFileError("Kinetics::addReaction", r->input,
-                    "Reaction '{}' contains the undeclared species '{}'",
-                    r->equation(), sp.first);
-            }
-        }
-    }
-    for (const auto& sp : r->orders) {
-        if (kineticsSpeciesIndex(sp.first) == npos) {
-            if (m_skipUndeclaredSpecies) {
-                return false;
-            } else {
-                throw InputFileError("Kinetics::addReaction", r->input,
-                    "Reaction '{}' has a reaction order specified for the "
-                    "undeclared species '{}'",
-                    r->equation(), sp.first);
-            }
-        }
+    // Check validity of reaction within the context of the Kinetics object
+    if (!r->checkSpecies(*this)) {
+        // Do not add reaction
+        return false;
     }
 
     // For reactions created outside the context of a Kinetics object, the units
@@ -578,7 +505,6 @@ bool Kinetics::addReaction(shared_ptr<Reaction> r)
         r->calculateRateCoeffUnits(*this);
     }
 
-    checkReactionBalance(*r);
     size_t irxn = nReactions(); // index of the new reaction
 
     // indices of reactant and product species within this Kinetics object