[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

@@ -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(

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:

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   
+          species="""H    R1A  R1B  P1   R2   P2A  P2B  R3   P3A  P3B  R4
-                     P2A  P2B  R3   P3A  P3B  
+                     P4   R5   P5A  P5B  R6   P6A  P6B  R7   P7A  P7B""",
-                     R4   P4   R5   P5A  P5B  
-                     R6   P6A  P6B""",
           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]])

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

test/data/pdep-test.xml

@@ -6,10 +6,8 @@
   <phase dim="3" id="gas">
     <elementArray datasrc="elements.xml">H C</elementArray>
     <speciesArray datasrc="#species_data">
-      H    R1A  R1B  P1   R2   
+      H    R1A  R1B  P1   R2   P2A  P2B  R3   P3A  P3B  R4
-      P2A  P2B  R3   P3A  P3B  
+      P4   R5   P5A  P5B  R6   P6A  P6B  R7   P7A  P7B</speciesArray>
-      R4   P4   R5   P5A  P5B  
-      R6   P6A  P6B</speciesArray>
     <reactionArray datasrc="#reaction_data"/>
     <state>
       <temperature units="K">300.0</temperature>
@@ -328,6 +326,57 @@
         </NASA>
       </thermo>
     </species>
+
+    <!-- species R7    -->
+    <species name="R7">
+      <atomArray>C:2 H:7 </atomArray>
+      <thermo>
+        <NASA Tmin="200.0" Tmax="1000.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             5.149876130E+00,  -1.367097880E-02,   4.918005990E-05,  -4.847430260E-08,
+             1.666939560E-11,  -1.024664760E+04,  -4.641303760E+00</floatArray>
+        </NASA>
+        <NASA Tmin="1000.0" Tmax="3500.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             7.485149500E-02,   1.339094670E-02,  -5.732858090E-06,   1.222925350E-09,
+             -1.018152300E-13,  -9.468344590E+03,   1.843731800E+01</floatArray>
+        </NASA>
+      </thermo>
+    </species>
+
+    <!-- species P7A    -->
+    <species name="P7A">
+      <atomArray>C:1 H:4 </atomArray>
+      <thermo>
+        <NASA Tmin="200.0" Tmax="1000.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             5.149876130E+00,  -1.367097880E-02,   4.918005990E-05,  -4.847430260E-08,
+             1.666939560E-11,  -1.024664760E+04,  -4.641303760E+00</floatArray>
+        </NASA>
+        <NASA Tmin="1000.0" Tmax="3500.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             7.485149500E-02,   1.339094670E-02,  -5.732858090E-06,   1.222925350E-09,
+             -1.018152300E-13,  -9.468344590E+03,   1.843731800E+01</floatArray>
+        </NASA>
+      </thermo>
+    </species>
+
+    <!-- species P7B    -->
+    <species name="P7B">
+      <atomArray>C:1 H:4 </atomArray>
+      <thermo>
+        <NASA Tmin="200.0" Tmax="1000.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             5.149876130E+00,  -1.367097880E-02,   4.918005990E-05,  -4.847430260E-08,
+             1.666939560E-11,  -1.024664760E+04,  -4.641303760E+00</floatArray>
+        </NASA>
+        <NASA Tmin="1000.0" Tmax="3500.0" P0="100000.0">
+           <floatArray size="7" name="coeffs">
+             7.485149500E-02,   1.339094670E-02,  -5.732858090E-06,   1.222925350E-09,
+             -1.018152300E-13,  -9.468344590E+03,   1.843731800E+01</floatArray>
+        </NASA>
+      </thermo>
+    </species>
   </speciesData>
   <reactionData id="reaction_data">
 
@@ -518,5 +567,23 @@
       <reactants>R6:1.0</reactants>
       <products>P6A:1.0 P6B:1</products>
     </reaction>
+
+    <!-- reaction 0007    -->
+    <reaction id="0007" reversible="yes" type="chebyshev">
+      <equation>R7 + H [=] P7A + P7B</equation>
+      <rateCoeff>
+        <Tmin>300.0</Tmin>
+        <Tmax>2000.0</Tmax>
+        <Pmin units="atm">0.009869232667160128</Pmin>
+        <Pmax units="atm">98.69232667160128</Pmax>
+        <floatArray name="coeffs" degreeT="4" degreeP="4">
+          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</floatArray>
+      </rateCoeff>
+      <reactants>R7:1.0 H:1</reactants>
+      <products>P7A:1.0 P7B:1</products>
+    </reaction>
   </reactionData>
 </ctml>

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);