From be585aab3e0adaa02adef90820e2ddda9d464c52 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 11 Nov 2014 00:11:33 +0000 Subject: [PATCH] [Kinetics] Implement direct constructor for ChebyshevRate --- include/cantera/kinetics/RxnRates.h | 14 +++++++++++ src/kinetics/RxnRates.cpp | 25 +++++++++++++++++++ test/data/kineticsfromscratch.cti | 7 ++++++ test/kinetics/kineticsFromScratch.cpp | 36 ++++++++++++++++++++++++++- 4 files changed, 81 insertions(+), 1 deletion(-) diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 6c4f7ef40..1cae3d027 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -17,6 +17,7 @@ namespace Cantera { +class Array2D; class Plog; //! Arrhenius reaction rate type depends only on temperature @@ -441,6 +442,19 @@ public: //! Constructor from ReactionData. explicit ChebyshevRate(const ReactionData& rdata); + //! Constructor directly from coefficient array + /* + * @param Pmin Minimum pressure [Pa] + * @param Pmax Maximum pressure [Pa] + * @param Tmin Minimum temperature [K] + * @param Tmax Maximum temperature [K] + * @param coeffs Coefficient array dimensioned `nT` by `nP` where `nT` and + * `nP` are the number of temperatures and pressures used in the fit, + * respectively. + */ + ChebyshevRate(double Pmin, double Pmax, double Tmin, double Tmax, + const Array2D& coeffs); + //! Update concentration-dependent parts of the rate coefficient. //! @param c base-10 logarithm of the pressure in Pa void update_C(const doublereal* c) { diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp index 89fd2e036..29d0d5617 100644 --- a/src/kinetics/RxnRates.cpp +++ b/src/kinetics/RxnRates.cpp @@ -1,6 +1,7 @@ //! @file RxnRates.cpp #include "cantera/kinetics/RxnRates.h" +#include "cantera/base/Array.h" namespace Cantera { @@ -277,4 +278,28 @@ ChebyshevRate::ChebyshevRate(const ReactionData& rdata) PrDen_ = 1.0 / (logPmax - logPmin); } +ChebyshevRate::ChebyshevRate(double Pmin, double Pmax, double Tmin, double Tmax, + const Array2D& coeffs) + : nP_(coeffs.nColumns()) + , nT_(coeffs.nRows()) + , chebCoeffs_(coeffs.nColumns() * coeffs.nRows(), 0.0) + , dotProd_(coeffs.nRows()) +{ + double logPmin = std::log10(Pmin); + double logPmax = std::log10(Pmax); + double TminInv = 1.0 / Tmin; + double TmaxInv = 1.0 / Tmax; + + TrNum_ = - TminInv - TmaxInv; + TrDen_ = 1.0 / (TmaxInv - TminInv); + PrNum_ = - logPmin - logPmax; + PrDen_ = 1.0 / (logPmax - logPmin); + + for (size_t t = 0; t < nT_; t++) { + for (size_t p = 0; p < nP_; p++) { + chebCoeffs_[nP_*t + p] = coeffs(t,p); + } + } +} + } diff --git a/test/data/kineticsfromscratch.cti b/test/data/kineticsfromscratch.cti index 38add76cc..eaba17b13 100644 --- a/test/data/kineticsfromscratch.cti +++ b/test/data/kineticsfromscratch.cti @@ -24,3 +24,10 @@ pdep_arrhenius('H2 + O2 <=> 2 OH', [(1.0, 'atm'), 4.910800e+31, -4.8507, 24772.8], [(10.0, 'atm'), 1.286600e+47, -9.0246, 39796.5], [(100.0, 'atm'), 5.963200e+56, -11.529, 52599.6]) + +chebyshev_reaction('HO2 <=> OH + O', + Tmin=290.0, Tmax=3000.0, + Pmin=(0.0098692326671601278, 'atm'), Pmax=(98.692326671601279, 'atm'), + coeffs=[[ 8.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02], + [ 1.97640e+00, 1.00370e+00, 7.28650e-03, -3.04320e-02], + [ 3.17700e-01, 2.68890e-01, 9.48060e-02, -7.63850e-03]]) diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 7ba8c2335..dcd1fd792 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -2,6 +2,7 @@ #include "cantera/kinetics/importKinetics.h" #include "cantera/thermo/IdealGasPhase.h" #include "cantera/kinetics/GasKinetics.h" +#include "cantera/base/Array.h" using namespace Cantera; @@ -28,7 +29,7 @@ public: void check_rates(int iRef) { ASSERT_EQ((size_t) 1, kin.nReactions()); - std::string X = "O:0.02 H2:0.2 O2:0.7 H:0.03 OH:0.05"; + std::string X = "O:0.02 H2:0.2 O2:0.5 H:0.03 OH:0.05 H2O:0.1 HO2:0.01"; p.setState_TPX(1200, 5*OneAtm, X); p_ref.setState_TPX(1200, 5*OneAtm, X); @@ -124,3 +125,36 @@ TEST_F(KineticsFromScratch, add_plog_reaction) kin.finalize(); check_rates(3); } + +TEST_F(KineticsFromScratch, add_chebyshev_reaction) +{ + // reaction 4: + // chebyshev_reaction( + // 'HO2 <=> OH + O', + // Tmin=290.0, Tmax=3000.0, + // Pmin=(0.0098692326671601278, 'atm'), Pmax=(98.692326671601279, 'atm'), + // coeffs=[[ 8.2883e+00, -1.1397e+00, -1.2059e-01, 1.6034e-02], + // [ 1.9764e+00, 1.0037e+00, 7.2865e-03, -3.0432e-02], + // [ 3.1770e-01, 2.6889e-01, 9.4806e-02, -7.6385e-03]]) + Composition reac = parseCompString("HO2:1"); + Composition prod = parseCompString("OH:1 O:1"); + Array2D coeffs(3, 4); + coeffs(0,0) = 8.2883e+00; + coeffs(0,1) = -1.1397e+00; + coeffs(0,2) = -1.2059e-01; + coeffs(0,3) = 1.6034e-02; + coeffs(1,0) = 1.9764e+00; + coeffs(1,1) = 1.0037e+00; + coeffs(1,2) = 7.2865e-03; + coeffs(1,3) = -3.0432e-02; + coeffs(2,0) = 3.1770e-01; + coeffs(2,1) = 2.6889e-01; + coeffs(2,2) = 9.4806e-02; + coeffs(2,3) = -7.6385e-03; + ChebyshevRate rate(1000.0, 10000000.0, 290, 3000, coeffs); + + shared_ptr R(new ChebyshevReaction(reac, prod, rate)); + kin.addReaction(R); + kin.finalize(); + check_rates(4); +}