From 5926d2db7c0d4919b75ee50828b0adab4e691a51 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Mon, 5 Oct 2015 14:59:35 -0400
Subject: [PATCH] [ck2cti] Handle reactions containing 'hv' pseudospecies

Convert to an irreversible reaction and ignore the photon.
---
 interfaces/cython/cantera/ck2cti.py           | 22 ++++++--
 .../cython/cantera/test/test_convert.py       | 10 ++++
 test/data/photo-reaction.inp                  |  9 +++
 test/data/photo-reaction.xml                  | 55 +++++++++++++++++++
 4 files changed, 92 insertions(+), 4 deletions(-)
 create mode 100644 test/data/photo-reaction.inp
 create mode 100644 test/data/photo-reaction.xml

diff --git a/interfaces/cython/cantera/ck2cti.py b/interfaces/cython/cantera/ck2cti.py
index 411dbc70b..9d6a5e027 100644
--- a/interfaces/cython/cantera/ck2cti.py
+++ b/interfaces/cython/cantera/ck2cti.py
@@ -1056,7 +1056,8 @@ class Parser(object):
             self.species_tokens.update(k + next_char for k in self.speciesDict)
         self.other_tokens = {'M': 'third-body', 'm': 'third-body',
                              '(+M)': 'falloff3b', '(+m)': 'falloff3b',
-                             '<=>': 'equal', '=>': 'equal', '=': 'equal'}
+                             '<=>': 'equal', '=>': 'equal', '=': 'equal',
+                             'HV': 'photon', 'hv': 'photon'}
         self.other_tokens.update(('(+%s)' % k, 'falloff3b: %s' % k) for k in self.speciesDict)
         self.Slen = max(map(len, self.other_tokens))
 
@@ -1160,6 +1161,7 @@ class Parser(object):
         def parseExpression(expression, dest):
             falloff3b = None
             thirdBody = False  # simple third body reaction (non-falloff)
+            photon = False
             for stoichiometry,species,kind in expression:
                 if kind == 'third-body':
                     thirdBody = True
@@ -1167,18 +1169,30 @@ class Parser(object):
                     falloff3b = 'M'
                 elif kind.startswith('falloff3b:'):
                     falloff3b = kind.split()[1]
+                elif kind == 'photon':
+                    photon = True
                 else:
                     dest.append((stoichiometry, self.speciesDict[species]))
 
-            return falloff3b, thirdBody
+            return falloff3b, thirdBody, photon
 
-        falloff_3b_r, thirdBody = parseExpression(reactants, reaction.reactants)
-        falloff_3b_p, thirdBody = parseExpression(products, reaction.products)
+        falloff_3b_r, thirdBody, photon_r = parseExpression(reactants, reaction.reactants)
+        falloff_3b_p, thirdBody, photon_p = parseExpression(products, reaction.products)
 
         if falloff_3b_r != falloff_3b_p:
             raise InputParseError('Third bodies do not match: "{0}" and "{1}" in'
                 ' reaction entry:\n\n{2}'.format(falloff_3b_r, falloff_3b_p, entry))
 
+        if photon_r:
+            raise InputParseError('Reactant photon not supported. '
+                                  'Found in reaction:\n{0}'.format(entry.strip()))
+        if photon_p and reversible:
+            self.warn('Found reversible reaction containing a product photon:'
+                '\n{0}\nIf the "--permissive" option was specified, this will '
+                'be converted to an irreversible reaction with the photon '
+                'removed.'.format(entry.strip()))
+            reaction.reversible = False
+
         reaction.thirdBody = falloff_3b_r
 
         # Determine the appropriate units for k(T) and k(T,P) based on the number of reactants
diff --git a/interfaces/cython/cantera/test/test_convert.py b/interfaces/cython/cantera/test/test_convert.py
index 55f812490..bb5b7a55e 100644
--- a/interfaces/cython/cantera/test/test_convert.py
+++ b/interfaces/cython/cantera/test/test_convert.py
@@ -274,6 +274,16 @@ class chemkinConverterTest(utilities.CanteraTest):
         self.assertArrayNear(R[:,1], [1, 0, 0, 1])
         self.assertArrayNear(P[:,1], [0, 0.33, 1.67, 0])
 
+    def test_photon(self):
+        convertMech('../data/photo-reaction.inp',
+                    thermoFile='../data/dummy-thermo.dat',
+                    outName='photo-reaction.cti', quiet=True,
+                    permissive=True)
+
+        ref, gas = self.checkConversion('../data/photo-reaction.xml',
+                                        'photo-reaction.cti')
+        self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6])
+
     def test_transport_normal(self):
         convertMech('../../data/inputs/h2o2.inp',
                     transportFile='../../data/transport/gri30_tran.dat',
diff --git a/test/data/photo-reaction.inp b/test/data/photo-reaction.inp
new file mode 100644
index 000000000..4d93f2430
--- /dev/null
+++ b/test/data/photo-reaction.inp
@@ -0,0 +1,9 @@
+ELEM  H  C END
+
+SPECIES
+R1A
+END
+
+REACTIONS
+R1A = R1A + HV          1.0E18   -2.0   1000
+END
diff --git a/test/data/photo-reaction.xml b/test/data/photo-reaction.xml
new file mode 100644
index 000000000..9ba4a6a27
--- /dev/null
+++ b/test/data/photo-reaction.xml
@@ -0,0 +1,55 @@
+<?xml version="1.0"?>
+<ctml>
+  <validate reactions="yes" species="yes"/>
+
+  <!-- phase gas     -->
+  <phase dim="3" id="gas">
+    <elementArray datasrc="elements.xml">H C</elementArray>
+    <speciesArray datasrc="#species_data">R1A</speciesArray>
+    <reactionArray datasrc="#reaction_data"/>
+    <state>
+      <temperature units="K">300.0</temperature>
+      <pressure units="Pa">101325.0</pressure>
+    </state>
+    <thermo model="IdealGas"/>
+    <kinetics model="GasKinetics"/>
+    <transport model="None"/>
+  </phase>
+
+  <!-- species definitions     -->
+  <speciesData id="species_data">
+
+    <!-- u'    species R1A    ' -->
+    <species name="R1A">
+      <atomArray>H:4 C:1 </atomArray>
+      <thermo>
+        <NASA Tmax="1000.0" Tmin="200.0" P0="100000.0">
+           <floatArray name="coeffs" size="7">
+             5.149876130E+00,  -1.367097880E-02,   4.918005990E-05,  -4.847430260E-08, 
+             1.666939560E-11,  -1.024664760E+04,  -4.641303760E+00</floatArray>
+        </NASA>
+        <NASA Tmax="3500.0" Tmin="1000.0" P0="100000.0">
+           <floatArray name="coeffs" size="7">
+             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">
+
+    <!-- reaction 0001    -->
+    <reaction reversible="no" id="0001">
+      <equation>R1A =] R1A</equation>
+      <rateCoeff>
+        <Arrhenius>
+           <A>1.000000E+18</A>
+           <b>-2.0</b>
+           <E units="cal/mol">1000.000000</E>
+        </Arrhenius>
+      </rateCoeff>
+      <reactants>R1A:1.0</reactants>
+      <products>R1A:1.0</products>
+    </reaction>
+  </reactionData>
+</ctml>