From 183af676d5ba298e61db5fe8539a8780555d1265 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 11 Jan 2013 22:56:25 +0000 Subject: [PATCH] [ck2cti] Added support for the UNITS keyword --- interfaces/python/ck2cti.py | 135 +++++++++++++++++++++++++++--------- test/data/units-custom.inp | 31 +++++++++ test/data/units-default.inp | 29 ++++++++ test/python/testConvert.py | 25 +++++-- 4 files changed, 182 insertions(+), 38 deletions(-) create mode 100644 test/data/units-custom.inp create mode 100644 test/data/units-default.inp diff --git a/interfaces/python/ck2cti.py b/interfaces/python/ck2cti.py index 118ea35d5..df0b4a974 100755 --- a/interfaces/python/ck2cti.py +++ b/interfaces/python/ck2cti.py @@ -36,7 +36,13 @@ import os.path import numpy as np import re -UNIT_OPTIONS = {'CAL/': 'cal/mol', +QUANTITY_UNITS = {'MOL': 'mol', + 'MOLE': 'mol', + 'MOLES': 'mol', + 'MOLEC': 'molec', + 'MOLECULES': 'molec'} + +ENERGY_UNITS = {'CAL/': 'cal/mol', 'CAL/MOL': 'cal/mol', 'CAL/MOLE': 'cal/mol', 'EVOL': 'eV', @@ -52,12 +58,16 @@ UNIT_OPTIONS = {'CAL/': 'cal/mol', 'KELVINS': 'K', 'KJOU': 'kJ/mol', 'KJOULES/MOL': 'kJ/mol', - 'KJOULES/MOLE': 'kJ/mol', - 'MOL': 'mol', - 'MOLE': 'mol', - 'MOLES': 'mol', - 'MOLEC': 'molec', - 'MOLECULES': 'molec'} + 'KJOULES/MOLE': 'kJ/mol'} + + +def compatible_quantities(quantity_basis, units): + if quantity_basis == 'mol': + return 'molec' not in units + elif quantity_basis == 'molec': + return 'molec' in units or 'mol' not in units + else: + raise Exception('Unknown quantity basis: "{0}"'.format(quantity_basis)) class InputParseError(Exception): @@ -291,12 +301,14 @@ class KineticsModel(object): """ - def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''): + def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment='', + parser=None): self.Tmin = Tmin self.Tmax = Tmax self.Pmin = Pmin self.Pmax = Pmax self.comment = comment + self.parser = parser def isPressureDependent(self): """ @@ -381,7 +393,17 @@ class Arrhenius(KineticsModel): return False def rateStr(self): - return '[{0.A[0]:e}, {0.n}, {0.Ea[0]}]'.format(self) + if compatible_quantities(self.parser.quantity_units, self.A[1]): + A = '{0:e}'.format(self.A[0]) + else: + A = "({0:e}, '{1}')".format(*self.A) + + if self.Ea[1] == self.parser.energy_units: + Ea = str(self.Ea[0]) + else: + Ea = "({0}, '{1}')".format(*self.Ea) + + return '[{0}, {1}, {2}]'.format(A, self.n, Ea) def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstring = reactantstr + arrow + productstr @@ -823,6 +845,26 @@ class Parser(object): pass return composition + def getRateConstantUnits(self, length_dims, length_units, quantity_dims, + quantity_units, time_dims=1, time_units='s'): + + units = '' + if length_dims: + units += length_units + if length_dims > 1: + units += str(length_dims) + if quantity_dims: + units += '/' + quantity_units + if quantity_dims > 1: + units += str(quantity_dims) + if time_dims: + units += '/' + time_units + if time_dims > 1: + units += str(time_dims) + if units.startswith('/'): + units = '1' + units + return units + def readThermoEntry(self, entry, TintDefault): """ Read a thermodynamics `entry` for one species in a Chemkin-format file @@ -934,6 +976,28 @@ class Parser(object): reaction and its associated kinetics. """ + # Handle non-default units which apply to this entry + energy_units = self.energy_units + quantity_units = self.quantity_units + if 'units' in entry.lower(): + for units in sorted(QUANTITY_UNITS, key=lambda k: -len(k)): + m = re.search(r'units *\/ *%s *\/' % re.escape(units), + entry, re.IGNORECASE) + if m: + entry = re.sub(r'units *\/ *%s *\/' % re.escape(units), '', + entry, flags=re.IGNORECASE) + quantity_units = QUANTITY_UNITS[units] + break + + for units in sorted(ENERGY_UNITS, key=lambda k: -len(k)): + m = re.search(r'units *\/ *%s *\/' % re.escape(units), + entry, re.IGNORECASE) + if m: + entry = re.sub(r'units *\/ *%s *\/' % re.escape(units), '', + entry, flags=re.IGNORECASE) + energy_units = ENERGY_UNITS[units] + break + lines = entry.strip().splitlines() # The first line contains the reaction equation and a set of modified Arrhenius parameters @@ -1013,26 +1077,25 @@ class Parser(object): # Determine the appropriate units for k(T) and k(T,P) based on the number of reactants # This assumes elementary kinetics for all reactions rStoich = sum(r[0] for r in reaction.reactants) + (1 if thirdBody else 0) - if rStoich == 3: - kunits = "cm^6/(mol^2*s)" - klow_units = "cm^9/(mol^3*s)" - elif rStoich == 2: - kunits = "cm^3/(mol*s)" - klow_units = "cm^6/(mol^2*s)" - elif rStoich == 1: - kunits = "s^-1" - klow_units = "cm^3/(mol*s)" - else: + if rStoich > 3 or rStoich < 1: raise InputParseError('Invalid number of reactant species ({0}) for reaction {1}.'.format(rStoich, reaction)) + length_dim = 3 * (rStoich - 1) + quantity_dim = rStoich - 1 + kunits = self.getRateConstantUnits(length_dim, 'cm', + quantity_dim, quantity_units) + klow_units = self.getRateConstantUnits(length_dim + 3, 'cm', + quantity_dim + 1, quantity_units) + # The rest of the first line contains the high-P limit Arrhenius parameters (if available) #tokens = lines[0][52:].split() tokens = lines[0].split()[1:] arrheniusHigh = Arrhenius( A=(A,kunits), n=n, - Ea=(Ea, self.energy_units), + Ea=(Ea, energy_units), T0=(1,"K"), + parser=self ) if len(lines) == 1: @@ -1061,8 +1124,9 @@ class Parser(object): arrheniusLow = Arrhenius( A=(float(tokens[0].strip()),klow_units), n=float(tokens[1].strip()), - Ea=(float(tokens[2].strip()),"kcal/mol"), + Ea=(float(tokens[2].strip()),energy_units), T0=(1,"K"), + parser=self ) elif 'rev' in line.lower(): @@ -1076,8 +1140,9 @@ class Parser(object): revReaction.kinetics = Arrhenius( A=(float(tokens[0].strip()),klow_units), n=float(tokens[1].strip()), - Ea=(float(tokens[2].strip()),"kcal/mol"), + Ea=(float(tokens[2].strip()),energy_units), T0=(1,"K"), + parser=self ) elif 'ford' in line.lower(): @@ -1100,6 +1165,7 @@ class Parser(object): T3=(T3,"K"), T1=(T1,"K"), T2=(T2,"K") if T2 is not None else None, + parser=self ) elif 'sri' in line.lower(): # SRI falloff parameters @@ -1115,9 +1181,9 @@ class Parser(object): E = None if D is None or E is None: - sri = Sri(A=A, B=B, C=C) + sri = Sri(A=A, B=B, C=C, parser=self) else: - sri = Sri(A=A, B=B, C=C, D=D, E=E) + sri = Sri(A=A, B=B, C=C, D=D, E=E, parser=self) elif 'cheb' in line.lower(): # Chebyshev parameters @@ -1153,10 +1219,10 @@ class Parser(object): pdepArrhenius.append([float(tokens[0].strip()), Arrhenius( A=(float(tokens[1].strip()),kunits), n=float(tokens[2].strip()), - Ea=(float(tokens[3].strip()),"kcal/mol"), + Ea=(float(tokens[3].strip()),energy_units), T0=(1,"K"), + parser=self )]) - else: # Assume a list of collider efficiencies for collider, efficiency in zip(tokens[0::2], tokens[1::2]): @@ -1179,6 +1245,7 @@ class Parser(object): reaction.kinetics = PDepArrhenius( pressures=([P for P, arrh in pdepArrhenius],"atm"), arrhenius=[arrh for P, arrh in pdepArrhenius], + parser=self ) elif troe is not None: troe.arrheniusHigh = arrheniusHigh @@ -1191,10 +1258,13 @@ class Parser(object): sri.efficiencies = efficiencies reaction.kinetics = sri elif arrheniusLow is not None: - reaction.kinetics = Lindemann(arrheniusHigh=arrheniusHigh, arrheniusLow=arrheniusLow) + reaction.kinetics = Lindemann(arrheniusHigh=arrheniusHigh, + arrheniusLow=arrheniusLow, + parser=self) reaction.kinetics.efficiencies = efficiencies elif thirdBody: - reaction.kinetics = ThirdBody(arrheniusHigh=arrheniusHigh) + reaction.kinetics = ThirdBody(arrheniusHigh=arrheniusHigh, + parser=self) reaction.kinetics.efficiencies = efficiencies else: reaction.kinetics = arrheniusHigh @@ -1332,14 +1402,13 @@ class Parser(object): except IndexError: pass - #global PROCESSED_UNITS, ENERGY_UNITS, QUANTITY_UNITS if not self.processed_units: self.processed_units = True - self.energy_units = UNIT_OPTIONS[energyUnits] - self.quantity_units = UNIT_OPTIONS[moleculeUnits] + self.energy_units = ENERGY_UNITS[energyUnits] + self.quantity_units = QUANTITY_UNITS[moleculeUnits] else: - if (self.energy_units != UNIT_OPTIONS[energyUnits] or - self.quantity_units != UNIT_OPTIONS[moleculeUnits]): + if (self.energy_units != ENERGY_UNITS[energyUnits] or + self.quantity_units != QUANTITY_UNITS[moleculeUnits]): raise InputParseError("Multiple REACTIONS sections with " "different units are not supported.") diff --git a/test/data/units-custom.inp b/test/data/units-custom.inp new file mode 100644 index 000000000..fddce3f1a --- /dev/null +++ b/test/data/units-custom.inp @@ -0,0 +1,31 @@ +ELEMENTS +H C AR +END + +SPECIES +H +R1A R1B P1 +R2 P2A P2B +R3 P3A P3B +R4 P4 +R6 P6A P6B +END + +REACTIONS KELVIN MOLECULES +R1A+R1B = P1+H 1.660538921E-6 -2.0 503.21956 +R2+H = P2A+P2B 3.8192395183E-9 -1.0 500 + UNITS / CAL/MOL / +R3+H+M = P3A+P3B+M 4.0E21 0.0 1207.726956 + H/2/ + UNITS / MOL / +R4 = P4 1.0E12 2.0 7531.2 + UNITS / JOULES/MOL / +R6(+M) = P6A + P6B(+M) 2.0E3 0.0 754.829347 + LOW / 4.98161676E0 -0.400 452.897608/ + TROE / 0.650 7050. 123.0/ + DUPLICATE +P6A + P6B(+M) = R6(+M) 6.64215568E-12 2.0 754.829347 + LOW / 8.2721685E-12 -1.400 452.897608/ + TROE / 0.650 7050. 123.0/ + DUPLICATE +END diff --git a/test/data/units-default.inp b/test/data/units-default.inp new file mode 100644 index 000000000..5c54b5b95 --- /dev/null +++ b/test/data/units-default.inp @@ -0,0 +1,29 @@ +ELEMENTS +H C AR +END + +SPECIES +H +R1A R1B P1 +R2 P2A P2B +R3 P3A P3B +R4 P4 +R6 P6A P6B +END + +REACTIONS + +R1A+R1B = P1+H 1.0E18 -2.0 1000 +R2+H = P2A+P2B 2.3E15 -1.0 500 +R3+H+M = P3A+P3B+M 4.0E21 0.0 2400 + H/2/ +R4 = P4 1.0E12 2.0 1800 +R6(+M) = P6A + P6B(+M) 2.0E3 0.0 1500 + LOW / 3.0E24 -0.400 900.00/ + TROE / 0.650 7050. 123.0/ + DUPLICATE +P6A + P6B(+M) = R6(+M) 4.0E12 2.0 1500 + LOW / 3.0E36 -1.400 900.00/ + TROE / 0.650 7050. 123.0/ + DUPLICATE +END diff --git a/test/python/testConvert.py b/test/python/testConvert.py index 9b568760a..059cb2af0 100644 --- a/test/python/testConvert.py +++ b/test/python/testConvert.py @@ -22,7 +22,10 @@ class chemkinConverterTest(utilities.CanteraTest): self.assertEqual(ref.elementNames(), gas.elementNames()) self.assertEqual(ref.speciesNames(), gas.speciesNames()) - self.assertTrue((ref.reactantStoichCoeffs() == gas.reactantStoichCoeffs()).all()) + coeffs_ref = ref.reactantStoichCoeffs() + coeffs_gas = gas.reactantStoichCoeffs() + self.assertEqual(coeffs_gas.shape, coeffs_ref.shape) + self.assertTrue((coeffs_gas == coeffs_ref).all()) compositionA = [[ref.nAtoms(i,j) for j in range(ref.nElements())] for i in range(ref.nSpecies())] @@ -47,7 +50,7 @@ class chemkinConverterTest(utilities.CanteraTest): self.assertNear(ref_h[i], gas_h[i], 1e-7) self.assertNear(ref_s[i], gas_s[i], 1e-7) - def checkKinetics(self, ref, gas, temperatures, pressures): + def checkKinetics(self, ref, gas, temperatures, pressures, tol=1e-8): for T,P in itertools.product(temperatures, pressures): ref.set(T=T, P=P) gas.set(T=T, P=P) @@ -56,8 +59,8 @@ class chemkinConverterTest(utilities.CanteraTest): gas_kf = gas.fwdRateConstants() gas_kr = gas.revRateConstants() for i in range(gas.nReactions()): - self.assertNear(ref_kf[i], gas_kf[i]) - self.assertNear(ref_kr[i], gas_kr[i]) + self.assertNear(ref_kf[i], gas_kf[i], rtol=tol) + self.assertNear(ref_kr[i], gas_kr[i], rtol=tol) def test_gri30(self): convertMech('../../data/inputs/gri30.inp', @@ -148,8 +151,20 @@ class chemkinConverterTest(utilities.CanteraTest): 'explicit-forward-order.cti') self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) - def test_transport_normal(self): + def test_reaction_units(self): + convertMech('../data/units-default.inp', + thermoFile='../data/dummy-thermo.dat', + outName='units-default.cti', quiet=True) + convertMech('../data/units-custom.inp', + thermoFile='../data/dummy-thermo.dat', + outName='units-custom.cti', quiet=True) + default, custom = self.checkConversion('units-default.cti', + 'units-custom.cti') + self.checkKinetics(default, custom, + [300, 800, 1450, 2800], [5e3, 1e5, 2e6], 1e-7) + + def test_transport_normal(self): convertMech('../../data/inputs/h2o2.inp', transportFile='../../data/transport/gri30_tran.dat', outName='h2o2_transport_normal.cti', quiet=True)