From 5ad656e34286efe5e9700a049a6b5a1b1773d06f Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 25 Feb 2019 00:08:24 -0500 Subject: [PATCH] [ck2cti] Fix Chebyshev rate constant when using local quantity units Since the leading Chebyshev coefficient has effective units like log(cm^3/kmol), it needs to be converted directly to the default units of the CTI file. Analogous to the fix for PLOG reactions in #435. --- interfaces/cython/cantera/ck2cti.py | 10 ++- .../cython/cantera/test/test_convert.py | 2 +- test/data/pdep-test.cti | 47 +++++++++++- test/data/pdep-test.inp | 23 ++++++ test/data/pdep-test.xml | 75 ++++++++++++++++++- test/kinetics/pdep.cpp | 30 +++++--- 6 files changed, 164 insertions(+), 23 deletions(-) diff --git a/interfaces/cython/cantera/ck2cti.py b/interfaces/cython/cantera/ck2cti.py index 2f9dc64a2..8f25033f4 100644 --- a/interfaces/cython/cantera/ck2cti.py +++ b/interfaces/cython/cantera/ck2cti.py @@ -65,6 +65,8 @@ ENERGY_UNITS = {'CAL/': 'cal/mol', 'KJOULES/MOL': 'kJ/mol', 'KJOULES/MOLE': 'kJ/mol'} +Avogadro = 6.02214129e23 # in molec/mol; value consistent with ct_defs.h. + _open = open if sys.version_info[0] == 2: string_types = (str, unicode) @@ -580,7 +582,7 @@ class Chebyshev(KineticsModel): """ - def __init__(self, coeffs=None, kunits='', **kwargs): + def __init__(self, coeffs=None, **kwargs): KineticsModel.__init__(self, **kwargs) if coeffs is not None: self.coeffs = np.array(coeffs, np.float64) @@ -590,7 +592,6 @@ class Chebyshev(KineticsModel): self.coeffs = None self.degreeT = 0 self.degreeP = 0 - self.kunits = kunits def isPressureDependent(self): """ @@ -1499,6 +1500,11 @@ class Parser(object): for p in range(chebyshev.degreeP): chebyshev.coeffs[t,p] = chebyshevCoeffs[index] index += 1 + if quantity_units == 'mol' and self.quantity_units == 'molec': + chebyshev.coeffs[0, 0] -= np.log10(Avogadro)*quantity_dim + elif quantity_units == 'molec' and self.quantity_units == 'mol': + chebyshev.coeffs[0, 0] += np.log10(Avogadro)*quantity_dim + reaction.kinetics = chebyshev elif pdepArrhenius is not None: reaction.kinetics = PDepArrhenius( diff --git a/interfaces/cython/cantera/test/test_convert.py b/interfaces/cython/cantera/test/test_convert.py index eefcd47c3..59f917b26 100644 --- a/interfaces/cython/cantera/test/test_convert.py +++ b/interfaces/cython/cantera/test/test_convert.py @@ -506,7 +506,7 @@ class CtmlConverterTest(utilities.CanteraTest): def test_pdep(self): gas = ct.Solution('pdep-test.cti') - self.assertEqual(gas.n_reactions, 6) + self.assertEqual(gas.n_reactions, 7) def test_invalid(self): try: diff --git a/test/data/pdep-test.cti b/test/data/pdep-test.cti index 0233a9b20..58b9327f0 100644 --- a/test/data/pdep-test.cti +++ b/test/data/pdep-test.cti @@ -2,10 +2,8 @@ units(length='cm', time='s', quantity='mol', act_energy='cal/mol') ideal_gas(name='gas', elements="H C", - species="""H R1A R1B P1 R2 - P2A P2B R3 P3A P3B - R4 P4 R5 P5A P5B - R6 P6A P6B""", + species="""H R1A R1B P1 R2 P2A P2B R3 P3A P3B R4 + P4 R5 P5A P5B R6 P6A P6B R7 P7A P7B""", reactions='all', initial_state=state(temperature=300.0, pressure=OneAtm)) @@ -211,6 +209,39 @@ species(name='P6B', 1.22292535E-09, -1.01815230E-13, -9.46834459E+03, 1.84373180E+01]))) +species(name='R7', + atoms='C:2 H:7', + thermo=(NASA([200.00, 1000.00], + [ 5.14987613E+00, -1.36709788E-02, 4.91800599E-05, + -4.84743026E-08, 1.66693956E-11, -1.02466476E+04, + -4.64130376E+00]), + NASA([1000.00, 3500.00], + [ 7.48514950E-02, 1.33909467E-02, -5.73285809E-06, + 1.22292535E-09, -1.01815230E-13, -9.46834459E+03, + 1.84373180E+01]))) + +species(name='P7A', + atoms='C:1 H:4', + thermo=(NASA([200.00, 1000.00], + [ 5.14987613E+00, -1.36709788E-02, 4.91800599E-05, + -4.84743026E-08, 1.66693956E-11, -1.02466476E+04, + -4.64130376E+00]), + NASA([1000.00, 3500.00], + [ 7.48514950E-02, 1.33909467E-02, -5.73285809E-06, + 1.22292535E-09, -1.01815230E-13, -9.46834459E+03, + 1.84373180E+01]))) + +species(name='P7B', + atoms='C:1 H:4', + thermo=(NASA([200.00, 1000.00], + [ 5.14987613E+00, -1.36709788E-02, 4.91800599E-05, + -4.84743026E-08, 1.66693956E-11, -1.02466476E+04, + -4.64130376E+00]), + NASA([1000.00, 3500.00], + [ 7.48514950E-02, 1.33909467E-02, -5.73285809E-06, + 1.22292535E-09, -1.01815230E-13, -9.46834459E+03, + 1.84373180E+01]))) + #------------------------------------------------------------------------------- # Reaction data #------------------------------------------------------------------------------- @@ -265,3 +296,11 @@ chebyshev_reaction('R6 <=> P6A + P6B', [-2.93660e-01, 2.85680e-01, -9.33730e-03, -4.01020e-03], [-2.26210e-01, 1.69190e-01, 4.85810e-03, -2.38030e-03], [-1.43220e-01, 7.71110e-02, 1.27080e-02, -6.41540e-04]]) +# Reaction 7 +chebyshev_reaction('R7 + H <=> P7A + P7B', + Tmin=300.0, Tmax=2000.0, + Pmin=(0.009869232667160128, 'atm'), Pmax=(98.69232667160128, 'atm'), + coeffs=[[ 3.20681e+01, -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], + [-3.12850e-02, -3.94120e-02, 4.43750e-02, 1.44580e-02]]) diff --git a/test/data/pdep-test.inp b/test/data/pdep-test.inp index bced1c162..48d986504 100644 --- a/test/data/pdep-test.inp +++ b/test/data/pdep-test.inp @@ -12,6 +12,7 @@ R3 P3A P3B R4 P4 R5 P5A P5B R6 P6A P6B +R7 P7A P7B END THERMO ALL @@ -95,6 +96,19 @@ P6B C 1H 4 G 200.000 3500.000 1000.000 1 7.48514950E-02 1.33909467E-02-5.73285809E-06 1.22292535E-09-1.01815230E-13 2 -9.46834459E+03 1.84373180E+01 5.14987613E+00-1.36709788E-02 4.91800599E-05 3 -4.84743026E-08 1.66693956E-11-1.02466476E+04-4.64130376E+00 4 + +R7 C 2H 7 G 200.000 3500.000 1000.000 1 + 7.48514950E-02 1.33909467E-02-5.73285809E-06 1.22292535E-09-1.01815230E-13 2 +-9.46834459E+03 1.84373180E+01 5.14987613E+00-1.36709788E-02 4.91800599E-05 3 +-4.84743026E-08 1.66693956E-11-1.02466476E+04-4.64130376E+00 4 +P7A C 1H 4 G 200.000 3500.000 1000.000 1 + 7.48514950E-02 1.33909467E-02-5.73285809E-06 1.22292535E-09-1.01815230E-13 2 +-9.46834459E+03 1.84373180E+01 5.14987613E+00-1.36709788E-02 4.91800599E-05 3 +-4.84743026E-08 1.66693956E-11-1.02466476E+04-4.64130376E+00 4 +P7B C 1H 4 G 200.000 3500.000 1000.000 1 + 7.48514950E-02 1.33909467E-02-5.73285809E-06 1.22292535E-09-1.01815230E-13 2 +-9.46834459E+03 1.84373180E+01 5.14987613E+00-1.36709788E-02 4.91800599E-05 3 +-4.84743026E-08 1.66693956E-11-1.02466476E+04-4.64130376E+00 4 END REACTIONS @@ -149,4 +163,13 @@ CHEB / -2.3294000e-01 4.0190000e-01 -2.6073000e-02 -5.0486000e-03 / CHEB / -2.9366000e-01 2.8568000e-01 -9.3373000e-03 -4.0102000e-03 / CHEB / -2.2621000e-01 1.6919000e-01 4.8581000e-03 -2.3803000e-03 / CHEB / -1.4322000e-01 7.7111000e-02 1.2708000e-02 -6.4154000e-04 / + +! Bimolecular CHEB with local quantity units +R7+H=P7A+P7B 1.0E0 0.0 0.0 +units / molec / +tcheb / 300.0 2000.0 / pcheb / 0.009869232667160128 98.69232667160128 / +cheb / 4 4 8.2883000e+00 -1.1397000e+00 -1.2059000e-01 1.6034000e-02 / +cheb / 1.9764000e+00 1.0037000e+00 7.2865000e-03 -3.0432000e-02 / +cheb / 3.1770000e-01 2.6889000e-01 9.4806000e-02 -7.6385000e-03 / +cheb / -3.1285000e-02 -3.9412000e-02 4.4375000e-02 1.4458000e-02 / END diff --git a/test/data/pdep-test.xml b/test/data/pdep-test.xml index a113993a2..ffb01960c 100644 --- a/test/data/pdep-test.xml +++ b/test/data/pdep-test.xml @@ -6,10 +6,8 @@ H C - H R1A R1B P1 R2 - P2A P2B R3 P3A P3B - R4 P4 R5 P5A P5B - R6 P6A P6B + H R1A R1B P1 R2 P2A P2B R3 P3A P3B R4 + P4 R5 P5A P5B R6 P6A P6B R7 P7A P7B 300.0 @@ -328,6 +326,57 @@ + + + + C:2 H:7 + + + + 5.149876130E+00, -1.367097880E-02, 4.918005990E-05, -4.847430260E-08, + 1.666939560E-11, -1.024664760E+04, -4.641303760E+00 + + + + 7.485149500E-02, 1.339094670E-02, -5.732858090E-06, 1.222925350E-09, + -1.018152300E-13, -9.468344590E+03, 1.843731800E+01 + + + + + + + C:1 H:4 + + + + 5.149876130E+00, -1.367097880E-02, 4.918005990E-05, -4.847430260E-08, + 1.666939560E-11, -1.024664760E+04, -4.641303760E+00 + + + + 7.485149500E-02, 1.339094670E-02, -5.732858090E-06, 1.222925350E-09, + -1.018152300E-13, -9.468344590E+03, 1.843731800E+01 + + + + + + + C:1 H:4 + + + + 5.149876130E+00, -1.367097880E-02, 4.918005990E-05, -4.847430260E-08, + 1.666939560E-11, -1.024664760E+04, -4.641303760E+00 + + + + 7.485149500E-02, 1.339094670E-02, -5.732858090E-06, 1.222925350E-09, + -1.018152300E-13, -9.468344590E+03, 1.843731800E+01 + + + @@ -518,5 +567,23 @@ R6:1.0 P6A:1.0 P6B:1 + + + + R7 + H [=] P7A + P7B + + 300.0 + 2000.0 + 0.009869232667160128 + 98.69232667160128 + + 2.90681e+01, -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, + -3.12850e-02, -3.94120e-02, 4.43750e-02, 1.44580e-02 + + R7:1.0 H:1 + P7A:1.0 P7B:1 + diff --git a/test/kinetics/pdep.cpp b/test/kinetics/pdep.cpp index c813a1fb9..1f1317a6a 100644 --- a/test/kinetics/pdep.cpp +++ b/test/kinetics/pdep.cpp @@ -57,14 +57,14 @@ Kinetics* PdepTest::kin_ = NULL; TEST_F(PdepTest, reactionCounts) { - EXPECT_EQ((size_t) 6, kin_->nReactions()); + EXPECT_EQ((size_t) 7, kin_->nReactions()); } TEST_F(PdepTest, PlogLowPressure) { // Test that P-log reactions have the right low-pressure limit set_TP(500.0, 1e-7); - vector_fp kf(6); + vector_fp kf(7); kin_->getFwdRateConstants(&kf[0]); // Pre-exponential factor decreases by 10^3 for second-order reaction @@ -84,7 +84,7 @@ TEST_F(PdepTest, PlogHighPressure) { // Test that P-log reactions have the right high-pressure limit set_TP(500.0, 1e10); - vector_fp kf(6); + vector_fp kf(7); kin_->getFwdRateConstants(&kf[0]); // Pre-exponential factor decreases by 10^3 for second-order reaction @@ -100,7 +100,7 @@ TEST_F(PdepTest, PlogDuplicatePressures) { // Test that multiple rate expressions are combined when necessary set_TP(500.0, 1e10); - vector_fp kf(6); + vector_fp kf(7); kin_->getFwdRateConstants(&kf[0]); double kf1 = k(1.3700e+14, -0.79, 17603.0) + k(1.2800e+03, 1.71, 9774.0); @@ -115,7 +115,7 @@ TEST_F(PdepTest, PlogCornerCases) // Test rate evaluation at the corner cases where the pressure // is exactly of the specified interpolation values set_TP(500.0, 101325); - vector_fp kf(6); + vector_fp kf(7); kin_->getFwdRateConstants(&kf[0]); double kf0 = k(4.910800e+28, -4.8507, 24772.8); @@ -130,7 +130,7 @@ TEST_F(PdepTest, PlogCornerCases) TEST_F(PdepTest, PlogIntermediatePressure1) { set_TP(1100.0, 20*101325); - vector_fp ropf(6); + vector_fp ropf(7); kin_->getFwdRatesOfProgress(&ropf[0]); // Expected rates computed using Chemkin @@ -144,7 +144,7 @@ TEST_F(PdepTest, PlogIntermediatePressure1) TEST_F(PdepTest, PlogIntermediatePressure2) { thermo_->setState_TP(1100.0, 0.5*101325); - vector_fp ropf(6); + vector_fp ropf(7); kin_->getFwdRatesOfProgress(&ropf[0]); EXPECT_NEAR(5.244649e+02, ropf[0], 5e-2); @@ -156,7 +156,7 @@ TEST_F(PdepTest, PlogIntermediatePressure2) TEST_F(PdepTest, PlogIntermediatePressure3) { thermo_->setState_TP(800.0, 70*101325); - vector_fp ropf(6); + vector_fp ropf(7); kin_->getFwdRatesOfProgress(&ropf[0]); EXPECT_NEAR(2.274501e+04, ropf[0], 1e+1); @@ -168,31 +168,37 @@ TEST_F(PdepTest, PlogIntermediatePressure3) TEST_F(PdepTest, ChebyshevIntermediate1) { // Test Chebyshev rates in the normal interpolation region - vector_fp kf(6); + vector_fp kf(7); 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); + + // Rate for a reaction specified as "molec" instead of "mol" should + // be higher by a factor of the Avogadro constant (in mol, not kmol). + // Accuracy is limited by the low precision used by ck2cti + EXPECT_NEAR(kf[4], kf[6]/(Avogadro*1e-3), 5e2); } TEST_F(PdepTest, ChebyshevIntermediate2) { // Test Chebyshev rates in the normal interpolation region - vector_fp kf(6); + vector_fp kf(7); 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); + EXPECT_NEAR(kf[4], kf[6]/(Avogadro*1e-3), 1e2); } TEST_F(PdepTest, ChebyshevIntermediateROP) { set_TP(1100.0, 30 * 101325); - vector_fp ropf(6); + vector_fp ropf(7); // Expected rates computed using Chemkin kin_->getFwdRatesOfProgress(&ropf[0]); EXPECT_NEAR(4.552930e+03, ropf[4], 1e-1); @@ -201,7 +207,7 @@ TEST_F(PdepTest, ChebyshevIntermediateROP) TEST_F(PdepTest, ChebyshevEdgeCases) { - vector_fp kf(6); + vector_fp kf(7); // Minimum P set_TP(500.0, 1000.0);