research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[ck2cti] Added support for the UNITS keyword

Browse source
This commit is contained in:
Ray Speth 2013-01-11 22:56:25 +00:00
parent 7ca8d1830e
commit 183af676d5
4 changed files with 182 additions and 38 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

135
interfaces/python/ck2cti.py
View file

@ -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.")

31
test/data/units-custom.inp Normal file
View file

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

29
test/data/units-default.inp Normal file
View file

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

25
test/python/testConvert.py
View file

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