From 67395c84383c1b96b6d1743fadf3d4ab129e1598 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Tue, 24 Jul 2012 23:03:28 +0000
Subject: [PATCH] Added validation of Plog rate expressions

---
 include/cantera/kinetics/ReactionData.h |  2 ++
 include/cantera/kinetics/RxnRates.h     | 34 +++++++++++++++++++++++--
 src/kinetics/importKinetics.cpp         |  1 +
 3 files changed, 35 insertions(+), 2 deletions(-)

diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h
index 0bab6ce51..ddcf48dc4 100644
--- a/include/cantera/kinetics/ReactionData.h
+++ b/include/cantera/kinetics/ReactionData.h
@@ -19,6 +19,7 @@ public:
     //! Default constructor
     ReactionData() {
         reactionType = ELEMENTARY_RXN;
+        validate = false;
         number = 0;
         rxn_number = 0;
         reversible = true;
@@ -41,6 +42,7 @@ public:
      */
     int reactionType;
 
+    bool validate;
     int number;
     int rxn_number;
     std::vector<size_t> reactants;
diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h
index a6265454e..617ac9680 100644
--- a/include/cantera/kinetics/RxnRates.h
+++ b/include/cantera/kinetics/RxnRates.h
@@ -535,6 +535,10 @@ public:
         n2_.resize(maxRates_);
         Ea1_.resize(maxRates_);
         Ea2_.resize(maxRates_);
+
+        if (rdata.validate) {
+            validate();
+        }
     }
 
     //! Update concentration-dependent parts of the rate coefficient.
@@ -584,7 +588,7 @@ public:
         if (m1_ == 1) {
             log_k1 = A1_[0] + n1_[0] * logT - Ea1_[0] * recipT;
         } else {
-            double k = 0.0;
+            double k = 1e-300; // non-zero to make log(k) finite
             for (size_t m = 0; m < m1_; m++) {
                 k += A1_[m] * exp(n1_[m] * logT - Ea1_[m] * recipT);
             }
@@ -594,7 +598,7 @@ public:
         if (m2_ == 1) {
             log_k2 = A2_[0] + n2_[0] * logT - Ea2_[0] * recipT;
         } else {
-            double k = 0.0;
+            double k = 1e-300; // non-zero to make log(k) finite
             for (size_t m = 0; m < m2_; m++) {
                 k += A2_[m] * exp(n2_[m] * logT - Ea2_[m] * recipT);
             }
@@ -621,6 +625,32 @@ public:
         return false;
     }
 
+    //! Check to make sure that the rate expression is finite over a range of
+    //! 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() {
+        double T[] = {1.0, 10.0, 100.0, 1000.0, 10000.0};
+        for (pressureIter iter = pressures_.begin();
+             iter->first < 1000;
+             iter++)
+        {
+            update_C(&iter->first);
+            for (size_t i=0; i < 5; i++) {
+                double k = updateRC(log(T[i]), 1.0/T[i]);
+                if (!(k >= 0)) {
+                    // k is NaN. Increment the iterator so that the error
+                    // message will correctly indicate that the problematic rate
+                    // expression is at the higher of the adjacent pressures.
+                    throw CanteraError("Plog::validate",
+                                       "Invalid rate coefficient at P = " +
+                                       fp2str(exp((++iter)->first)) +
+                                       ", T = " + fp2str(T[i]));
+                }
+            }
+        }
+    }
+
 protected:
     //! log(p) to (index range) in A_, n, Ea vectors
     std::map<double, std::pair<size_t, size_t> > pressures_;
diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp
index ed877ab3d..cf750c223 100644
--- a/src/kinetics/importKinetics.cpp
+++ b/src/kinetics/importKinetics.cpp
@@ -689,6 +689,7 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& r, Kinetics& kin,
     // xml data. Then, when we have collected everything we add the reaction to
     // the kinetics object, kin, at the end of the routine.
     ReactionData rdata;
+    rdata.validate = validate_rxn;
 
     // Check to see if the reaction is specified to be a duplicate of another
     // reaction. It's an error if the reaction is a duplicate and this is not