From ecce98c1dc4fe381e71fc60be47c1035c3390412 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 13 Sep 2019 14:24:07 -0400 Subject: [PATCH] [Kinetics] Fix Chebyshev rate evaluation with only 1 point in T or P --- include/cantera/kinetics/RxnRates.h | 16 ++++---- .../cython/cantera/test/test_kinetics.py | 41 +++++++++++++++++++ 2 files changed, 49 insertions(+), 8 deletions(-) diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 1a9db4492..ddd18d9d3 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -364,13 +364,13 @@ public: //! @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 Cnm1 = Pr; + double Cn = 1; double Cnp1; for (size_t j = 0; j < nT_; j++) { - dotProd_[j] = chebCoeffs_[nP_*j] + Pr * chebCoeffs_[nP_*j+1]; + dotProd_[j] = chebCoeffs_[nP_*j]; } - for (size_t i = 2; i < nP_; i++) { + for (size_t i = 1; i < nP_; i++) { Cnp1 = 2 * Pr * Cn - Cnm1; for (size_t j = 0; j < nT_; j++) { dotProd_[j] += Cnp1 * chebCoeffs_[nP_*j + i]; @@ -387,11 +387,11 @@ public: */ doublereal updateRC(doublereal logT, doublereal recipT) const { double Tr = (2 * recipT + TrNum_) * TrDen_; - double Cnm1 = 1; - double Cn = Tr; + double Cnm1 = Tr; + double Cn = 1; double Cnp1; - double logk = dotProd_[0] + Tr * dotProd_[1]; - for (size_t i = 2; i < nT_; i++) { + double logk = dotProd_[0]; + for (size_t i = 1; i < nT_; i++) { Cnp1 = 2 * Tr * Cn - Cnm1; logk += Cnp1 * dotProd_[i]; Cnm1 = Cn; diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 65b42b2d6..c9b220cc1 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -943,6 +943,47 @@ class TestReaction(utilities.CanteraTest): self.assertNear(gas2.net_rates_of_progress[0], gas1.net_rates_of_progress[4]) + def test_chebyshev_single_P(self): + species = ct.Species.listFromFile('pdep-test.cti') + r = ct.ChebyshevReaction() + r.reactants = 'R5:1, H:1' + r.products = 'P5A:1, P5B:1' + r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + coeffs=[[ 5.28830e+00], + [ 1.97640e+00], + [ 3.17700e-01], + [-3.12850e-02]]) + + gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=species, reactions=[r]) + + # rate constant should be pressure independent + for T in [300, 500, 1500]: + gas.TP = T, 1e4 + k1 = gas.forward_rate_constants[0] + gas.TP = T, 1e6 + k2 = gas.forward_rate_constants[0] + self.assertNear(k1, k2) + + def test_chebyshev_single_T(self): + species = ct.Species.listFromFile('pdep-test.cti') + r = ct.ChebyshevReaction() + r.reactants = 'R5:1, H:1' + r.products = 'P5A:1, P5B:1' + r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + coeffs=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02]]) + + gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=species, reactions=[r]) + + # rate constant should be temperature independent + for P in [1e4, 2e5, 8e6]: + gas.TP = 400, P + k1 = gas.forward_rate_constants[0] + gas.TP = 1700, P + k2 = gas.forward_rate_constants[0] + self.assertNear(k1, k2) + def test_chebyshev_rate(self): gas1 = ct.Solution('pdep-test.cti') gas1.TP = 800, 2*ct.one_atm