diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 3ff7590cf..ee8b48ba6 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/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 { diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 598cae870..0afbed026 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/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; diff --git a/test/kinetics/pdep.cpp b/test/kinetics/pdep.cpp index 8dee653ca..22980bf11 100644 --- a/test/kinetics/pdep.cpp +++ b/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)