[Kinetics] Validate rate coefficient data for Reaction objects

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;
 };
 

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

src/kinetics/Kinetics.cpp

@@ -651,6 +651,8 @@ void Kinetics::addReaction(ReactionData& r) {
 
 void Kinetics::addReaction(shared_ptr<Reaction> r)
 {
+    r->validateRateConstant();
+
     // Check for undeclared species
     for (Composition::const_iterator iter = r->reactants.begin();
          iter != r->reactants.end();

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");

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<double, Arrhenius>& 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]));
             }
         }

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<double, Arrhenius> rates;
+    typedef std::multimap<double, Arrhenius>::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<PlogReaction> 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<ElementaryReaction> 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<ElementaryReaction> 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
 {