From a399f034ce55cfeeeda5b0ab659a67335ef4d611 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 29 Apr 2021 14:50:53 -0500 Subject: [PATCH 1/9] [Kinetics] update detection of undeclared species --- include/cantera/kinetics/Reaction.h | 16 +++++++ src/kinetics/GasKinetics.cpp | 8 ---- src/kinetics/Kinetics.cpp | 67 +++++++++++++++++----------- src/kinetics/Reaction.cpp | 69 +++++++++++++++++++++++++---- 4 files changed, 118 insertions(+), 42 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index 410839569fc..f6af74609e8 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -75,6 +75,18 @@ class Reaction m_valid = valid; } + //! Return vector containing list of undeclared species + //! @param kin Kinetics object + std::vector undeclaredSpecies(const Kinetics& kin) const; + + //! Return vector containing list of undeclared third body species + //! @param kin Kinetics object + virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; + + //! Return vector containing list of undeclared reaction order species + //! @param kin Kinetics object + std::vector undeclaredOrders(const Kinetics& kin) const; + //! Type of the reaction. The valid types are listed in the file, //! reaction_defs.h, with constants ending in `RXN`. /*! @@ -209,6 +221,8 @@ class ThreeBodyReaction : public ElementaryReaction virtual void calculateRateCoeffUnits(const Kinetics& kin); virtual void getParameters(AnyMap& reactionNode) const; + virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; + //! Relative efficiencies of third-body species in enhancing the reaction //! rate. ThirdBody third_body; @@ -236,6 +250,8 @@ class FalloffReaction : public Reaction virtual void calculateRateCoeffUnits(const Kinetics& kin); virtual void getParameters(AnyMap& reactionNode) const; + virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; + //! The rate constant in the low-pressure limit Arrhenius low_rate; 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..b14727a40c4 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; @@ -537,38 +538,52 @@ bool Kinetics::addReaction(shared_ptr r) } // 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); - } + std::vector undeclared; + + undeclared = r->undeclaredSpecies(*this); + if (!undeclared.empty()) { + if (m_skipUndeclaredSpecies) { + return false; + } else { + throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" + "contains undeclared species: '{}'", + r->equation(), boost::algorithm::join(undeclared, "', '")); } } - 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); + + undeclared = r->undeclaredOrders(*this); + if (!undeclared.empty()) { + if (m_skipUndeclaredSpecies) { + return false; + } else { + if (r->input.hasKey("orders")) { + throw InputFileError("Kinetics::addReaction", r->input["orders"], + "Reaction '{}'\n" + "defines reaction orders for undeclared species: '{}'", + r->equation(), boost::algorithm::join(undeclared, "', '")); } + // Error for empty r->input AnyMap (e.g. XML) + throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" + "defines reaction orders for undeclared species: '{}'", + r->equation(), boost::algorithm::join(undeclared, "', '")); } } - 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); + + undeclared = r->undeclaredThirdBodies(*this); + if (!undeclared.empty()) { + if (!m_skipUndeclaredThirdBodies) { + if (r->input.hasKey("efficiencies")) { + throw InputFileError("Kinetics::addReaction", r->input["efficiencies"], + "Reaction '{}'\n" + "defines third-body efficiencies for undeclared species: '{}'", + r->equation(), boost::algorithm::join(undeclared, "', '")); } + // Error for specified ThirdBody or empty r->input AnyMap + throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" + "is a three-body reaction with undeclared species: '{}'", + r->equation(), boost::algorithm::join(undeclared, "', '")); + } else if (m_skipUndeclaredSpecies) { + return false; } } diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 05933324981..93e9eb8e0f0 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -86,17 +86,19 @@ 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); } } } @@ -211,6 +213,33 @@ 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::vector Reaction::undeclaredSpecies(const Kinetics& kin) const { + std::vector undeclared; + updateUndeclared(undeclared, reactants, kin); + updateUndeclared(undeclared, products, kin); + return undeclared; +} + +std::vector Reaction::undeclaredThirdBodies(const Kinetics& kin) const { + std::vector undeclared; + return undeclared; +} + +std::vector Reaction::undeclaredOrders(const Kinetics& kin) const { + std::vector undeclared; + updateUndeclared(undeclared, orders, kin); + return undeclared; +} + ElementaryReaction::ElementaryReaction(const Composition& reactants_, const Composition products_, const Arrhenius& rate_) @@ -327,6 +356,13 @@ void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const } } +std::vector ThreeBodyReaction::undeclaredThirdBodies( + const Kinetics& kin) const { + std::vector undeclared; + updateUndeclared(undeclared, third_body.efficiencies, kin); + return undeclared; +} + FalloffReaction::FalloffReaction() : Reaction() , falloff(new Falloff()) @@ -413,6 +449,13 @@ void FalloffReaction::getParameters(AnyMap& reactionNode) const } } +std::vector FalloffReaction::undeclaredThirdBodies( + const Kinetics& kin) const { + std::vector undeclared; + updateUndeclared(undeclared, third_body.efficiencies, kin); + return undeclared; +} + ChemicallyActivatedReaction::ChemicallyActivatedReaction() { reaction_type = CHEMACT_RXN; @@ -1447,13 +1490,23 @@ 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()) { + std::vector undeclared, more; + undeclared = R->undeclaredSpecies(kinetics); + more = R->undeclaredOrders(kinetics); + undeclared.insert(undeclared.end(), more.begin(), more.end()); + more = R->undeclaredThirdBodies(kinetics); + + if (R->valid() && undeclared.empty() && more.empty()) { all_reactions.emplace_back(R); - } else if (!kinetics.skipUndeclaredSpecies()) { + } else if (!kinetics.skipUndeclaredSpecies() && !undeclared.empty()) { throw InputFileError("getReactions", node, "Reaction '{}' contains undeclared species.", R->equation()); + } else if (!kinetics.skipUndeclaredThirdBodies() && !more.empty()) { + throw InputFileError("getReactions", node, + "Reaction '{}' contains undeclared third body species.", + R->equation()); } - }; + } return all_reactions; } From 8583b1d61b1148c85434d318c55df58a2cb859d1 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 29 Apr 2021 14:52:40 -0500 Subject: [PATCH 2/9] [Test] add unit tests for detection of undeclared species --- .../cython/cantera/test/test_kinetics.py | 87 +++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 88968872871..844dbcceaaf 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -423,6 +423,93 @@ 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_bodies(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, 2) + + 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') From 15607e7df1c905febf13ee523ac68e2a073c45c1 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 29 Apr 2021 16:25:38 -0500 Subject: [PATCH 3/9] Make divider lines for error messages 80 characters --- src/base/ctexceptions.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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) From 9a777192817d198430c57908629855ab23e1f583 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 8 May 2021 04:55:31 -0500 Subject: [PATCH 4/9] [Kinetics] Move checks of reactions from Kinetics to Reaction Also, fix logic to retain reactions with undeclared third body species. --- include/cantera/kinetics/Reaction.h | 31 +++++----- src/kinetics/Kinetics.cpp | 61 ++----------------- src/kinetics/Reaction.cpp | 92 ++++++++++++++++++++--------- 3 files changed, 86 insertions(+), 98 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index f6af74609e8..d8ed505262c 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -75,17 +75,12 @@ class Reaction m_valid = valid; } - //! Return vector containing list of undeclared species + //! 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 - std::vector undeclaredSpecies(const Kinetics& kin) const; - - //! Return vector containing list of undeclared third body species - //! @param kin Kinetics object - virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; - - //! Return vector containing list of undeclared reaction order species - //! @param kin Kinetics object - std::vector undeclaredOrders(const Kinetics& kin) const; + 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`. @@ -139,6 +134,12 @@ class Reaction //! Flag indicating whether reaction is set up correctly bool m_valid; + + //! @internal Helper function returning vector of undeclared third body species. + //! 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::vector undeclaredThirdBodies(const Kinetics& kin) const; }; @@ -221,13 +222,14 @@ class ThreeBodyReaction : public ElementaryReaction virtual void calculateRateCoeffUnits(const Kinetics& kin); virtual void getParameters(AnyMap& reactionNode) const; - virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; - //! Relative efficiencies of third-body species in enhancing the reaction //! rate. ThirdBody third_body; bool specified_collision_partner = false; //!< Input specifies collision partner + +protected: + virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; }; //! A reaction that is first-order in [M] at low pressure, like a third-body @@ -250,8 +252,6 @@ class FalloffReaction : public Reaction virtual void calculateRateCoeffUnits(const Kinetics& kin); virtual void getParameters(AnyMap& reactionNode) const; - virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; - //! The rate constant in the low-pressure limit Arrhenius low_rate; @@ -270,6 +270,9 @@ 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::vector undeclaredThirdBodies(const Kinetics& kin) const; }; //! A reaction where the rate decreases as pressure increases due to collisional diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index b14727a40c4..0d9050da32d 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -528,63 +528,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 - std::vector undeclared; - - undeclared = r->undeclaredSpecies(*this); - if (!undeclared.empty()) { - if (m_skipUndeclaredSpecies) { - return false; - } else { - throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" - "contains undeclared species: '{}'", - r->equation(), boost::algorithm::join(undeclared, "', '")); - } - } - - undeclared = r->undeclaredOrders(*this); - if (!undeclared.empty()) { - if (m_skipUndeclaredSpecies) { - return false; - } else { - if (r->input.hasKey("orders")) { - throw InputFileError("Kinetics::addReaction", r->input["orders"], - "Reaction '{}'\n" - "defines reaction orders for undeclared species: '{}'", - r->equation(), boost::algorithm::join(undeclared, "', '")); - } - // Error for empty r->input AnyMap (e.g. XML) - throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" - "defines reaction orders for undeclared species: '{}'", - r->equation(), boost::algorithm::join(undeclared, "', '")); - } - } - - undeclared = r->undeclaredThirdBodies(*this); - if (!undeclared.empty()) { - if (!m_skipUndeclaredThirdBodies) { - if (r->input.hasKey("efficiencies")) { - throw InputFileError("Kinetics::addReaction", r->input["efficiencies"], - "Reaction '{}'\n" - "defines third-body efficiencies for undeclared species: '{}'", - r->equation(), boost::algorithm::join(undeclared, "', '")); - } - // Error for specified ThirdBody or empty r->input AnyMap - throw InputFileError("Kinetics::addReaction", r->input, "Reaction '{}'\n" - "is a three-body reaction with undeclared species: '{}'", - r->equation(), boost::algorithm::join(undeclared, "', '")); - } else if (m_skipUndeclaredSpecies) { - return false; - } + // 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 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 93e9eb8e0f0..2cf611a3654 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -102,6 +102,15 @@ void Reaction::validate() } } } + + // 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 @@ -222,22 +231,64 @@ void updateUndeclared(std::vector& undeclared, } } -std::vector Reaction::undeclaredSpecies(const Kinetics& kin) const { - std::vector undeclared; - updateUndeclared(undeclared, reactants, kin); - updateUndeclared(undeclared, products, kin); - return undeclared; -} - std::vector Reaction::undeclaredThirdBodies(const Kinetics& kin) const { std::vector undeclared; return undeclared; } -std::vector Reaction::undeclaredOrders(const Kinetics& kin) const { +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); - return undeclared; + 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 + undeclared = undeclaredThirdBodies(kin); + 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, "', '")); + } + } + + return true; } ElementaryReaction::ElementaryReaction(const Composition& reactants_, @@ -262,7 +313,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() + "'"); } @@ -412,11 +463,11 @@ void FalloffReaction::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()); } @@ -700,7 +751,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() + "'"); } @@ -1490,21 +1541,8 @@ std::vector> getReactions(const AnyValue& items, std::vector> all_reactions; for (const auto& node : items.asVector()) { shared_ptr R(newReaction(node, kinetics)); - std::vector undeclared, more; - undeclared = R->undeclaredSpecies(kinetics); - more = R->undeclaredOrders(kinetics); - undeclared.insert(undeclared.end(), more.begin(), more.end()); - more = R->undeclaredThirdBodies(kinetics); - - if (R->valid() && undeclared.empty() && more.empty()) { + if (R->valid() && R->checkSpecies(kinetics)) { all_reactions.emplace_back(R); - } else if (!kinetics.skipUndeclaredSpecies() && !undeclared.empty()) { - throw InputFileError("getReactions", node, - "Reaction '{}' contains undeclared species.", R->equation()); - } else if (!kinetics.skipUndeclaredThirdBodies() && !more.empty()) { - throw InputFileError("getReactions", node, - "Reaction '{}' contains undeclared third body species.", - R->equation()); } } return all_reactions; From f0ca92c9aab6f97cbadf7869fea289191a2a0bff Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 8 May 2021 05:21:56 -0500 Subject: [PATCH 5/9] [Kinetics] Ensure that reaction balance is always checked Calling the check from within 'Reaction::checkSpecies' ensures that the reaction balance is always checked. Previously, this check was skipped for 'getReactions'. For consistency, the check itself is moved from 'Kinetics::checkReactionBalance' to 'Reaction::checkBalance'. --- include/cantera/kinetics/Kinetics.h | 4 ++ include/cantera/kinetics/Reaction.h | 6 +++ .../cython/cantera/test/test_kinetics.py | 2 +- src/kinetics/Kinetics.cpp | 42 ++--------------- src/kinetics/Reaction.cpp | 47 +++++++++++++++++++ 5 files changed, 61 insertions(+), 40 deletions(-) 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 d8ed505262c..a597ad87f19 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -75,6 +75,12 @@ 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, diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 844dbcceaaf..8d8f4d43c19 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') diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 0d9050da32d..4874d7d70ec 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -232,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, @@ -540,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 2cf611a3654..568c75a2c38 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -236,6 +236,49 @@ std::vector Reaction::undeclaredThirdBodies(const Kinetics& kin) co return undeclared; } +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 @@ -288,6 +331,9 @@ bool Reaction::checkSpecies(const Kinetics& kin) const { } } + // + checkBalance(kin); + return true; } @@ -1541,6 +1587,7 @@ std::vector> getReactions(const AnyValue& items, std::vector> all_reactions; for (const auto& node : items.asVector()) { shared_ptr R(newReaction(node, kinetics)); + R->validate(); if (R->valid() && R->checkSpecies(kinetics)) { all_reactions.emplace_back(R); } From bb781882a7253f7d53e30010ab3f6615831dd96c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 9 May 2021 03:29:49 -0500 Subject: [PATCH 6/9] [Kinetics] Refine logic for skipped reactions --- include/cantera/kinetics/Reaction.h | 12 ++++++++---- src/kinetics/Reaction.cpp | 19 ++++++++++++------- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index a597ad87f19..6922d4930b9 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -141,11 +141,13 @@ class Reaction //! Flag indicating whether reaction is set up correctly bool m_valid; - //! @internal Helper function returning vector of undeclared third body species. + //! @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::vector undeclaredThirdBodies(const Kinetics& kin) const; + virtual std::pair, bool> + undeclaredThirdBodies(const Kinetics& kin) const; }; @@ -235,7 +237,8 @@ class ThreeBodyReaction : public ElementaryReaction bool specified_collision_partner = false; //!< Input specifies collision partner protected: - virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; + virtual std::pair, bool> + undeclaredThirdBodies(const Kinetics& kin) const; }; //! A reaction that is first-order in [M] at low pressure, like a third-body @@ -278,7 +281,8 @@ class FalloffReaction : public Reaction Units low_rate_units; protected: - virtual std::vector undeclaredThirdBodies(const Kinetics& kin) const; + virtual std::pair, bool> undeclaredThirdBodies( + const Kinetics& kin) const; }; //! A reaction where the rate decreases as pressure increases due to collisional diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 568c75a2c38..4254f8e8c51 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -231,9 +231,10 @@ void updateUndeclared(std::vector& undeclared, } } -std::vector Reaction::undeclaredThirdBodies(const Kinetics& kin) const { +std::pair, bool> Reaction::undeclaredThirdBodies( + const Kinetics& kin) const { std::vector undeclared; - return undeclared; + return std::make_pair(undeclared, false); } void Reaction::checkBalance(const Kinetics& kin) const { @@ -315,7 +316,9 @@ bool Reaction::checkSpecies(const Kinetics& kin) const { } // Use helper function while there is no uniform handling of third bodies - undeclared = undeclaredThirdBodies(kin); + auto third = undeclaredThirdBodies(kin); + undeclared = third.first; + bool specified_collision_partner_ = third.second; if (!undeclared.empty()) { if (!kin.skipUndeclaredThirdBodies()) { if (input.hasKey("efficiencies")) { @@ -328,6 +331,8 @@ bool Reaction::checkSpecies(const Kinetics& kin) const { 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; } } @@ -453,11 +458,11 @@ void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const } } -std::vector ThreeBodyReaction::undeclaredThirdBodies( +std::pair, bool> ThreeBodyReaction::undeclaredThirdBodies( const Kinetics& kin) const { std::vector undeclared; updateUndeclared(undeclared, third_body.efficiencies, kin); - return undeclared; + return std::make_pair(undeclared, specified_collision_partner); } FalloffReaction::FalloffReaction() @@ -546,11 +551,11 @@ void FalloffReaction::getParameters(AnyMap& reactionNode) const } } -std::vector FalloffReaction::undeclaredThirdBodies( +std::pair, bool> FalloffReaction::undeclaredThirdBodies( const Kinetics& kin) const { std::vector undeclared; updateUndeclared(undeclared, third_body.efficiencies, kin); - return undeclared; + return std::make_pair(undeclared, false); } ChemicallyActivatedReaction::ChemicallyActivatedReaction() From e9fe939a8a29e5f60b90f419ad11da9e30ce1126 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 8 May 2021 06:11:05 -0500 Subject: [PATCH 7/9] [Test] Augment test for undeclared third bodies --- .../cython/cantera/test/test_kinetics.py | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 8d8f4d43c19..55000f3c039 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -457,7 +457,7 @@ def test_raise_undeclared_third_bodies(self): with self.assertRaisesRegex(ct.CanteraError, "three-body reaction with undeclared"): ct.Solution(yaml=gas_def) - def test_skip_undeclared_third_bodies(self): + def test_skip_undeclared_third_bodies1(self): gas_def = self._gas_def + """ reactions: @@ -466,7 +466,29 @@ def test_skip_undeclared_third_bodies(self): """ gas = ct.Solution(yaml=gas_def) - self.assertEqual(gas.n_reactions, 2) + 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): From 66b71386a7a10e39dd332f8308d9a59dfc2f8f5d Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 9 May 2021 03:48:41 -0500 Subject: [PATCH 8/9] [Test] Switch ion flame tests to ch4_ion.yaml --- interfaces/cython/cantera/test/test_onedim.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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) From eee271c843e117c49e74f0c09a04e8226f1bb33c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 10 May 2021 01:15:43 -0500 Subject: [PATCH 9/9] Update formatting --- src/kinetics/Reaction.cpp | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 4254f8e8c51..7e783631e2c 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -223,7 +223,8 @@ void Reaction::calculateRateCoeffUnits(const Kinetics& kin) } void updateUndeclared(std::vector& undeclared, - const Composition& comp, const Kinetics& kin) { + const Composition& comp, const Kinetics& kin) +{ for (const auto& sp: comp) { if (kin.kineticsSpeciesIndex(sp.first) == npos) { undeclared.emplace_back(sp.first); @@ -232,13 +233,14 @@ void updateUndeclared(std::vector& undeclared, } std::pair, bool> Reaction::undeclaredThirdBodies( - const Kinetics& kin) const { + const Kinetics& kin) const +{ std::vector undeclared; return std::make_pair(undeclared, false); } -void Reaction::checkBalance(const Kinetics& kin) const { - +void Reaction::checkBalance(const Kinetics& kin) const +{ Composition balr, balp; // iterate over products and reactants @@ -280,8 +282,8 @@ void Reaction::checkBalance(const Kinetics& kin) const { } } -bool Reaction::checkSpecies(const Kinetics& kin) const { - +bool Reaction::checkSpecies(const Kinetics& kin) const +{ // Check for undeclared species std::vector undeclared; updateUndeclared(undeclared, reactants, kin); @@ -336,7 +338,6 @@ bool Reaction::checkSpecies(const Kinetics& kin) const { } } - // checkBalance(kin); return true; @@ -406,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; @@ -415,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; @@ -459,7 +462,8 @@ void ThreeBodyReaction::getParameters(AnyMap& reactionNode) const } std::pair, bool> ThreeBodyReaction::undeclaredThirdBodies( - const Kinetics& kin) const { + const Kinetics& kin) const +{ std::vector undeclared; updateUndeclared(undeclared, third_body.efficiencies, kin); return std::make_pair(undeclared, specified_collision_partner); @@ -489,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() + " (+" + @@ -499,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() + " (+" + @@ -509,7 +515,8 @@ std::string FalloffReaction::productString() const { } } -void FalloffReaction::validate() { +void FalloffReaction::validate() +{ Reaction::validate(); if (!allow_negative_pre_exponential_factor && (low_rate.preExponentialFactor() < 0 || @@ -552,7 +559,8 @@ void FalloffReaction::getParameters(AnyMap& reactionNode) const } std::pair, bool> FalloffReaction::undeclaredThirdBodies( - const Kinetics& kin) const { + const Kinetics& kin) const +{ std::vector undeclared; updateUndeclared(undeclared, third_body.efficiencies, kin); return std::make_pair(undeclared, false);