diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 620cb89c792..a03fb3d2291 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/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 diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 410839569fc..6922d4930b9 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/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, 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, 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, bool> undeclaredThirdBodies( + const Kinetics& kin) const; }; //! A reaction where the rate decreases as pressure increases due to collisional diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 88968872871..55000f3c039 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/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') diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index d5fd3024786..753ce16a3e5 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/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) diff --git a/src/base/ctexceptions.cpp b/src/base/ctexceptions.cpp index 307ced49f23..6461af24f6c 100644 --- a/src/base/ctexceptions.cpp +++ b/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) diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index c0b94b1fa71..b89505d2b5b 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/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, diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 97b9486f52d..4874d7d70ec 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -17,6 +17,7 @@ #include "cantera/base/utilities.h" #include "cantera/base/global.h" #include +#include using namespace std; @@ -231,44 +232,9 @@ double Kinetics::checkDuplicateStoich(std::map& 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 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 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 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 05933324981..7e783631e2c 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -86,20 +86,31 @@ void Reaction::validate() if (!allow_nonreactant_orders) { for (const auto& order : orders) { if (reactants.find(order.first) == reactants.end()) { - throw CanteraError("Reaction::validate", "Reaction order " - "specified for non-reactant species '" + order.first + "'"); - } + throw InputFileError("Reaction::validate", input, + "Reaction order specified for non-reactant species '{}'", + order.first); + } } } if (!allow_negative_orders) { for (const auto& order : orders) { if (order.second < 0.0) { - throw CanteraError("Reaction::validate", "Negative reaction " - "order specified for species '" + order.first + "'"); + throw InputFileError("Reaction::validate", input, + "Negative reaction order specified for species '{}'", + order.first); } } } + + // 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 (reversible && !orders.empty()) { + throw InputFileError("Reaction::validate", input, + "Reaction orders may only be given for irreversible reactions"); + } } AnyMap Reaction::parameters(bool withInput) const @@ -211,6 +222,127 @@ void Reaction::calculateRateCoeffUnits(const Kinetics& kin) } } +void updateUndeclared(std::vector& undeclared, + const Composition& comp, const Kinetics& kin) +{ + for (const auto& sp: comp) { + if (kin.kineticsSpeciesIndex(sp.first) == npos) { + undeclared.emplace_back(sp.first); + } + } +} + +std::pair, bool> Reaction::undeclaredThirdBodies( + const Kinetics& kin) const +{ + std::vector undeclared; + return std::make_pair(undeclared, false); +} + +void Reaction::checkBalance(const Kinetics& kin) const +{ + Composition balr, balp; + + // iterate over products and reactants + for (const auto& sp : products) { + const ThermoPhase& ph = kin.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 : reactants) { + const ThermoPhase& ph = kin.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); + } + } + + std::string msg; + bool ok = true; + for (const auto& el : balr) { + const std::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("Reaction::checkBalance", input, + "The following reaction is unbalanced: {}\n" + " Element Reactants Products\n{}", + equation(), msg); + } +} + +bool Reaction::checkSpecies(const Kinetics& kin) const +{ + // Check for undeclared species + std::vector undeclared; + updateUndeclared(undeclared, reactants, kin); + updateUndeclared(undeclared, products, kin); + if (!undeclared.empty()) { + if (kin.skipUndeclaredSpecies()) { + return false; + } else { + throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" + "contains undeclared species: '{}'", + equation(), boost::algorithm::join(undeclared, "', '")); + } + } + + undeclared.clear(); + updateUndeclared(undeclared, orders, kin); + if (!undeclared.empty()) { + if (kin.skipUndeclaredSpecies()) { + return false; + } else { + if (input.hasKey("orders")) { + throw InputFileError("Reaction::checkSpecies", input["orders"], + "Reaction '{}'\n" + "defines reaction orders for undeclared species: '{}'", + equation(), boost::algorithm::join(undeclared, "', '")); + } + // Error for empty input AnyMap (e.g. XML) + throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" + "defines reaction orders for undeclared species: '{}'", + equation(), boost::algorithm::join(undeclared, "', '")); + } + } + + // Use helper function while there is no uniform handling of third bodies + auto third = undeclaredThirdBodies(kin); + undeclared = third.first; + bool specified_collision_partner_ = third.second; + if (!undeclared.empty()) { + if (!kin.skipUndeclaredThirdBodies()) { + if (input.hasKey("efficiencies")) { + throw InputFileError("Reaction::checkSpecies", input["efficiencies"], + "Reaction '{}'\n" + "defines third-body efficiencies for undeclared species: '{}'", + equation(), boost::algorithm::join(undeclared, "', '")); + } + // Error for specified ThirdBody or empty input AnyMap + throw InputFileError("Reaction::checkSpecies", input, "Reaction '{}'\n" + "is a three-body reaction with undeclared species: '{}'", + equation(), boost::algorithm::join(undeclared, "', '")); + } else if (kin.skipUndeclaredSpecies() && specified_collision_partner_) { + return false; + } + } + + checkBalance(kin); + + return true; +} + ElementaryReaction::ElementaryReaction(const Composition& reactants_, const Composition products_, const Arrhenius& rate_) @@ -233,7 +365,7 @@ void ElementaryReaction::validate() Reaction::validate(); if (!allow_negative_pre_exponential_factor && rate.preExponentialFactor() < 0) { - throw CanteraError("ElementaryReaction::validate", + throw InputFileError("ElementaryReaction::validate", input, "Undeclared negative pre-exponential factor found in reaction '" + equation() + "'"); } @@ -275,7 +407,8 @@ ThreeBodyReaction::ThreeBodyReaction(const Composition& reactants_, reaction_type = THREE_BODY_RXN; } -std::string ThreeBodyReaction::reactantString() const { +std::string ThreeBodyReaction::reactantString() const +{ if (specified_collision_partner) { return ElementaryReaction::reactantString() + " + " + third_body.efficiencies.begin()->first; @@ -284,7 +417,8 @@ std::string ThreeBodyReaction::reactantString() const { } } -std::string ThreeBodyReaction::productString() const { +std::string ThreeBodyReaction::productString() const +{ if (specified_collision_partner) { return ElementaryReaction::productString() + " + " + third_body.efficiencies.begin()->first; @@ -327,6 +461,14 @@ void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const } } +std::pair, bool> ThreeBodyReaction::undeclaredThirdBodies( + const Kinetics& kin) const +{ + std::vector undeclared; + updateUndeclared(undeclared, third_body.efficiencies, kin); + return std::make_pair(undeclared, specified_collision_partner); +} + FalloffReaction::FalloffReaction() : Reaction() , falloff(new Falloff()) @@ -351,7 +493,8 @@ FalloffReaction::FalloffReaction( reaction_type = FALLOFF_RXN; } -std::string FalloffReaction::reactantString() const { +std::string FalloffReaction::reactantString() const +{ if (third_body.default_efficiency == 0 && third_body.efficiencies.size() == 1) { return Reaction::reactantString() + " (+" + @@ -361,7 +504,8 @@ std::string FalloffReaction::reactantString() const { } } -std::string FalloffReaction::productString() const { +std::string FalloffReaction::productString() const +{ if (third_body.default_efficiency == 0 && third_body.efficiencies.size() == 1) { return Reaction::productString() + " (+" + @@ -371,16 +515,17 @@ std::string FalloffReaction::productString() const { } } -void FalloffReaction::validate() { +void FalloffReaction::validate() +{ Reaction::validate(); if (!allow_negative_pre_exponential_factor && (low_rate.preExponentialFactor() < 0 || high_rate.preExponentialFactor() < 0)) { - throw CanteraError("FalloffReaction::validate", "Negative " + throw InputFileError("FalloffReaction::validate", input, "Negative " "pre-exponential factor found for reaction '" + equation() + "'"); } if (low_rate.preExponentialFactor() * high_rate.preExponentialFactor() < 0) { - throw CanteraError("FalloffReaction::validate", "High and " + throw InputFileError("FalloffReaction::validate", input, "High and " "low rate pre-exponential factors must have the same sign." "Reaction: '{}'", equation()); } @@ -413,6 +558,14 @@ void FalloffReaction::getParameters(AnyMap& reactionNode) const } } +std::pair, bool> FalloffReaction::undeclaredThirdBodies( + const Kinetics& kin) const +{ + std::vector undeclared; + updateUndeclared(undeclared, third_body.efficiencies, kin); + return std::make_pair(undeclared, false); +} + ChemicallyActivatedReaction::ChemicallyActivatedReaction() { reaction_type = CHEMACT_RXN; @@ -657,7 +810,7 @@ void BlowersMaselReaction::validate() Reaction::validate(); if (!allow_negative_pre_exponential_factor && rate.preExponentialFactor() < 0) { - throw CanteraError("BlowersMaselReaction::validate", + throw InputFileError("BlowersMaselReaction::validate", input, "Undeclared negative pre-exponential factor found in reaction '" + equation() + "'"); } @@ -1447,13 +1600,11 @@ std::vector> getReactions(const AnyValue& items, std::vector> all_reactions; for (const auto& node : items.asVector()) { shared_ptr R(newReaction(node, kinetics)); - if (R->valid()) { + R->validate(); + if (R->valid() && R->checkSpecies(kinetics)) { all_reactions.emplace_back(R); - } else if (!kinetics.skipUndeclaredSpecies()) { - throw InputFileError("getReactions", node, - "Reaction '{}' contains undeclared species.", R->equation()); } - }; + } return all_reactions; }