From 82fbff7a059aad182f10fe263702aa42305c1056 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 15 Nov 2014 00:47:09 +0000 Subject: [PATCH] [Kinetics] Validate rate coefficient data for Reaction objects --- include/cantera/kinetics/Reaction.h | 10 +++++--- include/cantera/kinetics/RxnRates.h | 2 +- src/kinetics/Kinetics.cpp | 2 ++ src/kinetics/Reaction.cpp | 30 ++++++++++++++++++++-- src/kinetics/RxnRates.cpp | 9 +++---- test/kinetics/kineticsFromScratch.cpp | 37 +++++++++++++++++++++++++++ 6 files changed, 79 insertions(+), 11 deletions(-) diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index d8ac3364f..6d588c15c 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -28,6 +28,9 @@ public: virtual std::string productString() const; std::string equation() const; + //! Ensure that the rate constant for this reaction is valid. + virtual void validateRateConstant() {} + //! Type of the reaction. The valid types are listed in the file, //! reaction_defs.h, with constants ending in `RXN`. int reaction_type; @@ -50,8 +53,6 @@ public: //! True if the current reaction is reversible. False otherwise bool reversible; - bool validate; //!< Perform validation of the rate coefficient data - //! True if the current reaction is marked as duplicate bool duplicate; }; @@ -65,8 +66,10 @@ public: ElementaryReaction(); ElementaryReaction(const Composition& reactants, const Composition products, const Arrhenius& rate); + virtual void validateRateConstant(); Arrhenius rate; + bool allow_negative_pre_exponential_factor; }; //! A class for managing third-body efficiencies, including default values @@ -112,6 +115,7 @@ public: const vector_fp& falloff_params); virtual std::string reactantString() const; virtual std::string productString() const; + virtual void validateRateConstant(); Arrhenius low_rate; Arrhenius high_rate; @@ -144,7 +148,7 @@ public: PlogReaction(); PlogReaction(const Composition& reactants, const Composition& products, const Plog& rate); - + virtual void validateRateConstant(); Plog rate; }; diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index ed717c44c..5a2abacb9 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -409,7 +409,7 @@ public: //! temperatures at each interpolation pressure. This is potentially an //! issue when one of the Arrhenius expressions at a particular pressure //! has a negative pre-exponential factor. - void validate(const ReactionData& rdata); + void validate(const std::string& equation); protected: //! log(p) to (index range) in A_, n, Ea vectors diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 8ba3d715f..bff58de0e 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -651,6 +651,8 @@ void Kinetics::addReaction(ReactionData& r) { void Kinetics::addReaction(shared_ptr r) { + r->validateRateConstant(); + // Check for undeclared species for (Composition::const_iterator iter = r->reactants.begin(); iter != r->reactants.end(); diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index b5aba0764..fd1acc62f 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -15,7 +15,6 @@ namespace Cantera Reaction::Reaction(int type) : reaction_type(type) , reversible(true) - , validate(true) , duplicate(false) { } @@ -26,7 +25,6 @@ Reaction::Reaction(int type, const Composition& reactants_, , reactants(reactants_) , products(products_) , reversible(true) - , validate(true) , duplicate(false) { } @@ -79,14 +77,26 @@ ElementaryReaction::ElementaryReaction(const Composition& reactants_, const Arrhenius& rate_) : Reaction(ELEMENTARY_RXN, reactants_, products_) , rate(rate_) + , allow_negative_pre_exponential_factor(false) { } ElementaryReaction::ElementaryReaction() : Reaction(ELEMENTARY_RXN) + , allow_negative_pre_exponential_factor(false) { } +void ElementaryReaction::validateRateConstant() +{ + if (!allow_negative_pre_exponential_factor && + rate.preExponentialFactor() < 0) { + throw CanteraError("ElementaryReaction::validateRateConstant", + "Undeclared negative pre-exponential factor found in reaction '" + + equation() + "'"); + } +} + ThirdBody::ThirdBody(double default_eff) : default_efficiency(default_eff) { @@ -155,6 +165,14 @@ std::string FalloffReaction::productString() const { } } +void FalloffReaction::validateRateConstant() { + if (low_rate.preExponentialFactor() < 0 || + high_rate.preExponentialFactor() < 0) { + throw CanteraError("FalloffReaction::validateRateConstant", "Negative " + "pre-exponential factor found for reaction '" + equation() + "'"); + } +} + ChemicallyActivatedReaction::ChemicallyActivatedReaction() { reaction_type = CHEMACT_RXN; @@ -321,6 +339,9 @@ void setupElementaryReaction(ElementaryReaction& R, const XML_Node& rxn_node) } else { throw CanteraError("setupElementaryReaction", "Couldn't find Arrhenius node"); } + if (rxn_node["negative_A"] == "yes") { + R.allow_negative_pre_exponential_factor = true; + } setupReaction(R, rxn_node); } @@ -400,6 +421,11 @@ void setupPlogReaction(PlogReaction& R, const XML_Node& rxn_node) setupReaction(R, rxn_node); } +void PlogReaction::validateRateConstant() +{ + rate.validate(equation()); +} + void setupChebyshevReaction(ChebyshevReaction& R, const XML_Node& rxn_node) { XML_Node& rc = rxn_node.child("rateCoeff"); diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 4ca7d7bf8..c707673d3 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -191,7 +191,7 @@ Plog::Plog(const ReactionData& rdata) Ea2_.resize(maxRates_); if (rdata.validate) { - validate(rdata); + validate(rdata.equation); } } @@ -251,7 +251,7 @@ Plog::Plog(const std::multimap& rates) Ea2_.resize(maxRates_); } -void Plog::validate(const ReactionData& rdata) +void Plog::validate(const std::string& equation) { double T[] = {200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0}; for (pressureIter iter = pressures_.begin(); @@ -265,9 +265,8 @@ void Plog::validate(const ReactionData& rdata) // message will correctly indicate that the problematic rate // expression is at the higher of the adjacent pressures. throw CanteraError("Plog::validate", - "Invalid rate coefficient for reaction #" + - int2str(rdata.number) + ":\n" + rdata.equation + "\n" + - "at P = " + fp2str(std::exp((++iter)->first)) + + "Invalid rate coefficient for reaction '" + equation + + "'\nat P = " + fp2str(std::exp((++iter)->first)) + ", T = " + fp2str(T[i])); } } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 1e3885646..9f9cc285d 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -155,6 +155,21 @@ TEST_F(KineticsFromScratch, add_plog_reaction) check_rates(3); } +TEST_F(KineticsFromScratch, plog_invalid_rate) +{ + Composition reac = parseCompString("H2:1, O2:1"); + Composition prod = parseCompString("OH:2"); + std::multimap rates; + typedef std::multimap::value_type item; + rates.insert(item(0.01*101325, Arrhenius(1.2124e+16, -0.5779, 10872.7 / GasConst_cal_mol_K))); + rates.insert(item(10.0*101325, Arrhenius(1e15, -1, 10000 / GasConst_cal_mol_K))); + rates.insert(item(10.0*101325, Arrhenius(-2e20, -2.0, 20000 / GasConst_cal_mol_K))); + rates.insert(item(100.0*101325, Arrhenius(5.9632e+56, -11.529, 52599.6 / GasConst_cal_mol_K))); + + shared_ptr R(new PlogReaction(reac, prod, Plog(rates))); + ASSERT_THROW(kin.addReaction(R), CanteraError); +} + TEST_F(KineticsFromScratch, add_chebyshev_reaction) { // reaction 4: @@ -211,6 +226,28 @@ TEST_F(KineticsFromScratch, skip_undeclared_species) ASSERT_EQ(0, kin.nReactions()); } +TEST_F(KineticsFromScratch, negative_A_error) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + Arrhenius rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + shared_ptr R(new ElementaryReaction(reac, prod, rate)); + + ASSERT_THROW(kin.addReaction(R), CanteraError); + ASSERT_EQ(0, kin.nReactions()); +} + +TEST_F(KineticsFromScratch, allow_negative_A) +{ + Composition reac = parseCompString("O:1 H2:1"); + Composition prod = parseCompString("H:1 OH:1"); + Arrhenius rate(-3.87e1, 2.7, 6260.0 / GasConst_cal_mol_K); + shared_ptr R(new ElementaryReaction(reac, prod, rate)); + R->allow_negative_pre_exponential_factor = true; + + kin.addReaction(R); + ASSERT_EQ((size_t) 1, kin.nReactions()); +} class InterfaceKineticsFromScratch : public testing::Test {