diff --git a/interfaces/python/ck2cti.py b/interfaces/python/ck2cti.py index 073cb5f74..e12d67b31 100755 --- a/interfaces/python/ck2cti.py +++ b/interfaces/python/ck2cti.py @@ -34,6 +34,7 @@ import logging import types import os.path import numpy as np +import re UNIT_OPTIONS = {'CAL/': 'cal/mol', 'CAL/MOL': 'cal/mol', @@ -228,6 +229,7 @@ class Reaction(object): self.reversible = reversible self.duplicate = duplicate self.fwdOrders = fwdOrders if fwdOrders is not None else {} + self.thirdBody = None def _coeff_string(self, coeffs): L = [] @@ -236,8 +238,11 @@ class Reaction(object): L.append('{0} {1}'.format(stoichiometry, species)) else: L.append(str(species)) + expression = ' + '.join(L) + if self.thirdBody: + expression += ' (+ {0})'.format(self.thirdBody) - return ' + '.join(L) + return expression @property def reactantString(self): @@ -483,7 +488,7 @@ class Chebyshev(KineticsModel): return True def to_cti(self, reactantstr, arrow, productstr, indent=0): - rxnstr = reactantstr + ' (+ M)' + arrow + productstr + ' (+ M)' + rxnstr = reactantstr + arrow + productstr prefix = ' '*(indent+19) lines = ['chebyshev_reaction({0!r},'.format(rxnstr), prefix + 'Tmin={0.Tmin}, Tmax={0.Tmax},'.format(self), @@ -588,7 +593,7 @@ class Lindemann(ThirdBody): self.arrheniusLow = arrheniusLow def to_cti(self, reactantstr, arrow, productstr, indent=0): - rxnstr = reactantstr + ' (+ M)' + arrow + productstr + ' (+ M)' + rxnstr = reactantstr + arrow + productstr prefix = ' '*(indent + 17) lines = ['falloff_reaction({0!r},'.format(rxnstr)] lines.append(prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr())) @@ -662,7 +667,7 @@ class Troe(Lindemann): self.T2 = T2 def to_cti(self, reactantstr, arrow, productstr, indent=0): - rxnstr = reactantstr + ' (+ M)' + arrow + productstr + ' (+ M)' + rxnstr = reactantstr + arrow + productstr prefix = ' '*17 lines = ['falloff_reaction({0!r},'.format(rxnstr), prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()), @@ -713,7 +718,7 @@ class Sri(Lindemann): self.E = E def to_cti(self, reactantstr, arrow, productstr, indent=0): - rxnstr = reactantstr + ' (+ M)' + arrow + productstr + ' (+ M)' + rxnstr = reactantstr + arrow + productstr prefix = ' '*17 lines = ['falloff_reaction({0!r},'.format(rxnstr), prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()), @@ -942,16 +947,26 @@ def readKineticsEntry(entry, speciesDict, energyUnits, moleculeUnits): else: raise InputParseError("Failed to find reactant/product delimiter in reaction string.") - reactants = reactants.replace('(+M)','') - reactants = reactants.replace('(+m)','') - products = products.replace('(+M)','') - products = products.replace('(+m)','') - # Create a new Reaction object for this reaction reaction = Reaction(reactants=[], products=[], reversible=reversible) def parseExpression(expression, dest): - thirdBody = False + falloff3b = None + thirdBody = False # simple third body reaction (non-falloff) + + # Look for third-body species for falloff reactions + if re.search(r'\(\+[Mm]\)', expression): + falloff3b = 'M' + expression = re.sub(r'(\(\+[Mm]\))', '', expression) + elif re.search(r'\(\+.*\)', expression): + # See if it matches a known species + for species in speciesDict: + if re.search(r'\(\+%s\)' % re.escape(species), expression): + falloff3b = species + expression = re.sub(r'(\(\+%s\))' % re.escape(species), + '', expression) + break + for term in expression.split('+'): term = term.strip() if not term[0].isalpha(): @@ -970,14 +985,20 @@ def readKineticsEntry(entry, speciesDict, energyUnits, moleculeUnits): if species == 'M' or species == 'm': thirdBody = True elif species not in speciesDict: - raise InputParseError('Unexpected species "{0}" in reaction {1}.'.format(species, reaction)) + raise InputParseError('Unexpected species "{0}" in reaction expression "{1}".'.format(species, expression)) else: dest.append((stoichiometry, speciesDict[species])) - return thirdBody + return falloff3b, thirdBody - thirdBody = parseExpression(reactants, reaction.reactants) - parseExpression(products, reaction.products) + falloff_3b_r, thirdBody = parseExpression(reactants, reaction.reactants) + falloff_3b_p, thirdBody = 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)) + + reaction.thirdBody = falloff_3b_r # Determine the appropriate units for k(T) and k(T,P) based on the number of reactants # This assumes elementary kinetics for all reactions diff --git a/test/data/explicit-third-bodies.inp b/test/data/explicit-third-bodies.inp new file mode 100644 index 000000000..a554bf720 --- /dev/null +++ b/test/data/explicit-third-bodies.inp @@ -0,0 +1,38 @@ +ELEMENTS +H C AR +END + +SPECIES +H +R1A R1B P1 +R2 SP)X +END + +THERMO + 300.000 1000.000 5000.000 +SP)X C 1H 3 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 + +R1A+R1B+m = P1+H+M 3.0E19 -2.0 1900 + DUPLICATE + +R1A+R1B(+ M ) = P1+H(+m) 1.0E18 -2.0 1000 + LOW/4.0E25 -3.0 0/ + R2/0.0/ SP)X/0/ + DUPLICATE + +R1A+R1B(+R2) = P1+H(+ R2) 2.0E18 -3.0 1200 + LOW/3.0E25 -2.0 0/ + DUPLICATE + +! A really pathological species name +R1A+R1B(+SP)X) = P1+H(+SP)X) 3.0E18 -1.0 1400 + LOW/1.0E25 -1.0 0/ + DUPLICATE + +END diff --git a/test/data/explicit-third-bodies.xml b/test/data/explicit-third-bodies.xml new file mode 100644 index 000000000..ed91bf366 --- /dev/null +++ b/test/data/explicit-third-bodies.xml @@ -0,0 +1,205 @@ + + + + + + + H C Ar + + H R1A R1B P1 R2 + SP)X + + + 300.0 + 101325.0 + + + + + + + + + + + + H:1 + + + + 2.500000000E+00, 7.053328190E-13, -1.995919640E-15, 2.300816320E-18, + -9.277323320E-22, 2.547365990E+04, -4.466828530E-01 + + + + 2.500000010E+00, -2.308429730E-11, 1.615619480E-14, -4.735152350E-18, + 4.981973570E-22, 2.547365990E+04, -4.466829140E-01 + + + + + + + H:4 C:1 + + + + 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 + + + + + + + H:4 C:1 + + + + 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 + + + + + + + H:7 C:2 + + + + 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 + + + + + + + H:7 C:2 + + + + 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 + + + + + + + H:3 C:1 + + + + 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 + + + + + + + + + R1A + R1B + M [=] P1 + H + M + + + 3.000000E+13 + -2.0 + 1900.000000 + + + R1B:1 R1A:1.0 + H:1 P1:1.0 + + + + + R1A + R1B (+ M) [=] P1 + H (+ M) + + + 1.000000E+15 + -2.0 + 1000.000000 + + + 4.000000E+19 + -3.0 + 0.000000 + + SP)X:0.0 R2:0.0 + + + R1B:1 R1A:1.0 + H:1 P1:1.0 + + + + + R1A + R1B (+ R2) [=] P1 + H (+ R2) + + + 2.000000E+15 + -3.0 + 1200.000000 + + + 3.000000E+19 + -2.0 + 0.000000 + + R2:1.0 + + + R1B:1 R1A:1.0 + H:1 P1:1.0 + + + + + R1A + R1B (+ SP)X) [=] P1 + H (+ SP)X) + + + 3.000000E+15 + -1.0 + 1400.000000 + + + 1.000000E+19 + -1.0 + 0.000000 + + SP)X:1.0 + + + R1B:1 R1A:1.0 + H:1 P1:1.0 + + + diff --git a/test/python/testConvert.py b/test/python/testConvert.py index d5e5ced97..0cc11f946 100644 --- a/test/python/testConvert.py +++ b/test/python/testConvert.py @@ -124,6 +124,16 @@ class chemkinConverterTest(utilities.CanteraTest): 'sri-falloff.cti') self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) + def test_explicit_third_bodies(self): + if os.path.exists('explicit-third-bodies.cti'): + os.path.remove('explicit-third-bodies.cti') + ck2cti.convertMech('../data/explicit-third-bodies.inp', + thermoFile='../data/dummy-thermo.dat', quiet=True) + + ref, gas = self.checkConversion('explicit-third-bodies.cti', + '../data/explicit-third-bodies.xml') + self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) + def test_explicit_reverse_rate(self): if os.path.exists('explicit-reverse-rate.cti'): os.remove('explicit-reverse-rate.cti')