Implemented Chebyshev rate expressions

include/cantera/kinetics/RxnRates.h

@@ -660,27 +660,61 @@ public:
     ChebyshevRate() {}
 
     //! Constructor from ReactionData.
-    explicit ChebyshevRate(const ReactionData& rdata)
+    explicit ChebyshevRate(const ReactionData& rdata) :
+        nP_(rdata.chebDegreeP),
+        nT_(rdata.chebDegreeT),
+        chebCoeffs_(rdata.chebCoeffs),
+        dotProd_(rdata.chebDegreeT)
     {
+        double logPmin = log10(rdata.chebPmin);
+        double logPmax = log10(rdata.chebPmax);
+        double TminInv = 1.0 / rdata.chebTmin;
+        double TmaxInv = 1.0 / rdata.chebTmax;
+
+        TrNum_ = - TminInv - TmaxInv;
+        TrDen_ = 1.0 / (TmaxInv - TminInv);
+        PrNum_ = - logPmin - logPmax;
+        PrDen_ = 1.0 / (logPmax - logPmin);
     }
 
     //! Update concentration-dependent parts of the rate coefficient.
-    //! @param c natural log of the pressure in Pa
+    //! @param c base-10 logarithm of the pressure in Pa
     void update_C(const doublereal* c)
     {
+        double Pr = (2 * c[0] + PrNum_) * PrDen_;
+        double Cnm1 = 1;
+        double Cn = Pr;
+        double Cnp1;
+        for (size_t j = 0; j < nT_; j++) {
+            dotProd_[j] = chebCoeffs_[nP_*j] + Pr * chebCoeffs_[nP_*j+1];
+        }
+        for (size_t i = 2; i < nP_; i++) {
+            Cnp1 = 2 * Pr * Cn - Cnm1;
+            for (size_t j = 0; j < nT_; j++) {
+                dotProd_[j] += Cnp1 * chebCoeffs_[nP_*j + i];
+            }
+            Cnm1 = Cn;
+            Cn = Cnp1;
+        }
     }
 
     /**
-     * Update the value of the logarithm of the rate constant.
-     *
-     * Note, this function should never be called for negative A values.
-     * If it does then it will produce a negative overflow result, and
-     * a zero net forwards reaction rate, instead of a negative reaction
-     * rate constant that is the expected result.
+     * Update the value of the base-10 logarithm of the rate constant.
      */
     doublereal update(doublereal logT, doublereal recipT) const
     {
-        return 0.0;
+        double Tr = (2 * recipT + TrNum_) * TrDen_;
+        double Cnm1 = 1;
+        double Cn = Tr;
+        double Cnp1;
+        double logk = dotProd_[0] + Tr * dotProd_[1];
+        for (size_t i = 2; i < nT_; i++) {
+            Cnp1 = 2 * Tr * Cn - Cnm1;
+            logk += Cnp1 * dotProd_[i];
+            Cnm1 = Cn;
+            Cn = Cnp1;
+        }
+        return logk;
     }
 
     /**
@@ -689,7 +723,7 @@ public:
      * This function returns the actual value of the rate constant.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return exp(update(logT, recipT));
+        return pow(10, update(logT, recipT));
     }
 
     doublereal activationEnergy_R() const {
@@ -701,6 +735,13 @@ public:
     }
 
 protected:
+    double TrNum_, TrDen_; //!< terms appearing in the reduced temperature
+    double PrNum_, PrDen_; //!< terms appearing in the reduced pressure
+
+    size_t nP_; //!< number of points in the pressure direction
+    size_t nT_; //!< number of points in the temperature direction
+    vector_fp chebCoeffs_; //!< Chebyshev coefficients, length nP * nT
+    vector_fp dotProd_; //!< dot product of chebCoeffs with the reduced pressure polynomial
 };
 
 //     class LandauTeller {

src/kinetics/GasKinetics.cpp

@@ -206,16 +206,16 @@ _update_rates_C()
         m_falloff_concm.update(m_conc, ctot, &concm_falloff_values[0]);
     }
 
-    double logP = log(thermo().pressure());
-
     // P-log reactions
     if (m_plog_rates.nReactions()) {
+        double logP = log(thermo().pressure());
         m_plog_rates.update_C(&logP);
     }
 
     // Chebyshev reactions
     if (m_cheb_rates.nReactions()) {
-        m_cheb_rates.update_C(&logP);
+        double log10P = log10(thermo().pressure());
+        m_cheb_rates.update_C(&log10P);
     }
 
     m_ROP_ok = false;

test/kinetics/pdep.cpp

@@ -156,6 +156,62 @@ TEST_F(PdepTest, PlogIntermediatePressure3) {
     EXPECT_NEAR(2.224601e+07, ropf[2], 1e+3);
     EXPECT_NEAR(1.007440e+07, ropf[3], 1e+3);
 }
+
+TEST_F(PdepTest, ChebyshevIntermediate1) {
+    // Test Chebyshev rates in the normal interpolation region
+    vector_fp kf(6);
+
+    set_TP(1100.0, 20 * 101325);
+    kin_->getFwdRateConstants(&kf[0]);
+    // Expected rates computed using RMG-py
+    EXPECT_NEAR(3.130698657e+06, kf[4], 1e-1);
+    EXPECT_NEAR(1.187949573e+00, kf[5], 1e-7);
+}
+
+TEST_F(PdepTest, ChebyshevIntermediate2) {
+    // Test Chebyshev rates in the normal interpolation region
+    vector_fp kf(6);
+
+    set_TP(400.0, 0.1 * 101325);
+    kin_->getFwdRateConstants(&kf[0]);
+    // Expected rates computed using RMG-py
+    EXPECT_NEAR(1.713599902e+05, kf[4], 1e-3);
+    EXPECT_NEAR(9.581780687e-24, kf[5], 1e-31);
+}
+
+TEST_F(PdepTest, ChebyshevIntermediateROP) {
+    set_TP(1100.0, 30 * 101325);
+    vector_fp ropf(6);
+    // Expected rates computed using Chemkin
+    kin_->getFwdRatesOfProgress(&ropf[0]);
+    EXPECT_NEAR(4.552930e+03, ropf[4], 1e-1);
+    EXPECT_NEAR(4.877390e-02, ropf[5], 1e-5);
+}
+
+TEST_F(PdepTest, ChebyshevEdgeCases) {
+    vector_fp kf(6);
+
+    // Minimum P
+    set_TP(500.0, 1000.0);
+    kin_->getFwdRateConstants(&kf[0]);
+    EXPECT_NEAR(1.225785655e+06, kf[4], 1e-2);
+
+    // Maximum P
+    set_TP(500.0, 1.0e7);
+    kin_->getFwdRateConstants(&kf[0]);
+    EXPECT_NEAR(1.580981157e+03, kf[4], 1e-5);
+
+    // Minimum T
+    set_TP(300.0, 101325);
+    kin_->getFwdRateConstants(&kf[0]);
+    EXPECT_NEAR(5.405987017e+03, kf[4], 1e-5);
+
+    // Maximum T
+    set_TP(2000.0, 101325);
+    kin_->getFwdRateConstants(&kf[0]);
+    EXPECT_NEAR(3.354054351e+07, kf[4], 1e-1);
+}
+
 } // namespace Cantera
 
 int main(int argc, char** argv)