From 8ae3ca5a17951de267c2e0022232f3e524b13bd5 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 11 Jan 2013 22:56:16 +0000 Subject: [PATCH] [ck2cti] Refactored to add class Parser and remove global variables --- interfaces/python/ck2cti.py | 1460 +++++++++++++++++------------------ test/python/testConvert.py | 98 +-- 2 files changed, 772 insertions(+), 786 deletions(-) diff --git a/interfaces/python/ck2cti.py b/interfaces/python/ck2cti.py index dd1311334..118ea35d5 100755 --- a/interfaces/python/ck2cti.py +++ b/interfaces/python/ck2cti.py @@ -59,11 +59,6 @@ UNIT_OPTIONS = {'CAL/': 'cal/mol', 'MOLEC': 'molec', 'MOLECULES': 'molec'} -PROCESSED_UNITS = False -ENERGY_UNITS = 'cal/mol' -QUANTITY_UNITS = 'mol' -WARNING_AS_ERROR = True - class InputParseError(Exception): """ @@ -74,13 +69,6 @@ class InputParseError(Exception): pass -def warn(message): - if WARNING_AS_ERROR: - raise InputParseError(message) - else: - logging.warning(message) - - class Species(object): def __init__(self, label): self.label = label @@ -798,756 +786,750 @@ def fortFloat(s): return float(s) -def parseComposition(elements, nElements, width): - """ - Parse the elemental composition from a 7 or 9 coefficient NASA polynomial - entry. - """ - composition = {} - for i in range(nElements): - symbol = elements[width*i:width*i+2].strip() - count = elements[width*i+2:width*i+width].strip() - if not symbol: - continue - try: - count = int(float(count)) - if count: - composition[symbol.capitalize()] = count - except ValueError: - pass - return composition +class Parser(object): + def __init__(self): + self.processed_units = False + self.energy_units = 'cal/mol' + self.quantity_units = 'mol' + self.warning_as_error = True + self.elements = [] + self.speciesList = [] + self.speciesDict = {} + self.reactions = [] -def readThermoEntry(entry, TintDefault): - """ - Read a thermodynamics `entry` for one species in a Chemkin-format file - (consisting of two 7-coefficient NASA polynomials). Returns the label of - the species, the thermodynamics model as a :class:`MultiNASA` object, the - elemental composition of the species, and the comment/note associated with - the thermo entry. - """ - lines = entry.splitlines() - identifier = lines[0][0:24].split() - species = identifier[0].strip() - - if len(identifier) > 1: - note = ''.join(identifier[1:]).strip() - else: - note = '' - - # Extract the NASA polynomial coefficients - # Remember that the high-T polynomial comes first! - try: - Tmin = fortFloat(lines[0][45:55]) - Tmax = fortFloat(lines[0][55:65]) - try: - Tint = fortFloat(lines[0][65:75]) - except ValueError: - Tint = TintDefault - - coeffs_high = [fortFloat(lines[i][j:k]) - for i,j,k in [(1,0,15), (1,15,30), (1,30,45), (1,45,60), - (1,60,75), (2,0,15), (2,15,30)]] - coeffs_low = [fortFloat(lines[i][j:k]) - for i,j,k in [(2,30,45), (2,45,60), (2,60,75), (3,0,15), - (3,15,30), (3,30,45), (3,45,60)]] - - except (IndexError, ValueError): - raise InputParseError('Error while reading thermo entry for species {0}'.format(species)) - - composition = parseComposition(lines[0][24:44], 4, 5) - - # Non-standard extended elemental composition data may be located beyond - # column 80 on the first line of the thermo entry - if len(lines[0]) > 80: - elements = lines[0][80:] - composition2 = parseComposition(elements, len(elements)/10, 10) - composition.update(composition2) - - # Construct and return the thermodynamics model - thermo = MultiNASA( - polynomials=[ - NASA(Tmin=(Tmin,"K"), Tmax=(Tint,"K"), coeffs=coeffs_low), - NASA(Tmin=(Tint,"K"), Tmax=(Tmax,"K"), coeffs=coeffs_high) - ], - Tmin=(Tmin,"K"), - Tmax=(Tmax,"K"), - ) - - return species, thermo, composition, note - - -def readNasa9Entry(entry): - """ - Read a thermodynamics `entry` for one species given as one or more - 9-coefficient NASA polynomials, written in the format described in - Appendix A of NASA Reference Publication 1311 (McBride and Gordon, 1996). - Returns the label of the species, the thermodynamics model as a - :class:`MultiNASA` object, the elemental composition of the species, and - the comment/note associated with the thermo entry. - """ - tokens = entry[0].split() - species = tokens[0] - note = ' '.join(tokens[1:]) - N = int(entry[1][:2]) - note2 = entry[1][3:9].strip() - if note and note2: - note = '{0} [{1}]'.format(note, note2) - elif note2: - note = note2 - - composition = parseComposition(entry[1][10:50], 5, 8) - - polys = [] - totalTmin = 1e100 - totalTmax = -1e100 - try: - for i in range(N): - A,B,C = entry[2+3*i:2+3*(i+1)] - Tmin = fortFloat(A[1:11]) - Tmax = fortFloat(A[11:21]) - coeffs = [fortFloat(B[0:16]), fortFloat(B[16:32]), - fortFloat(B[32:48]), fortFloat(B[48:64]), - fortFloat(B[64:80]), fortFloat(C[0:16]), - fortFloat(C[16:32]), fortFloat(C[48:64]), - fortFloat(C[64:80])] - polys.append(NASA(Tmin=(Tmin,"K"), Tmax=(Tmax,"K"), coeffs=coeffs)) - totalTmin = min(Tmin, totalTmin) - totalTmax = max(Tmax, totalTmax) - except (IndexError, ValueError): - raise InputParseError('Error while reading thermo entry for species {0}'.format(species)) - - thermo = MultiNASA(polynomials=polys, - Tmin=(totalTmin,"K"), - Tmax=(totalTmax,"K")) - - return species, thermo, composition, note - - -def readKineticsEntry(entry, speciesDict, energyUnits, moleculeUnits): - """ - Read a kinetics `entry` for a single reaction as loaded from a - Chemkin-format file. The associated mapping of labels to species - `speciesDict` should also be provided. Returns a :class:`Reaction` object - with the reaction and its associated kinetics. - """ - - lines = entry.strip().splitlines() - - # The first line contains the reaction equation and a set of modified Arrhenius parameters - tokens = lines[0].split() - A = float(tokens[-3]) - n = float(tokens[-2]) - Ea = float(tokens[-1]) - reaction = ''.join(tokens[:-3]) - revReaction = None - - # Split the reaction equation into reactants and products - if '<=>' in reaction: - reversible = True - reactants, products = reaction.split('<=>') - elif '=>' in reaction: - reversible = False - reactants, products = reaction.split('=>') - elif '=' in reaction: - reversible = True - reactants, products = reaction.split('=') - else: - raise InputParseError("Failed to find reactant/product delimiter in reaction string.") - - # Create a new Reaction object for this reaction - reaction = Reaction(reactants=[], products=[], reversible=reversible) - - def parseExpression(expression, dest): - 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(): - # This allows for for non-unity stoichiometric coefficients, e.g. - # 2A=B+C or .85A+.15B=>C - j = [i for i,c in enumerate(term) if c.isalpha()][0] - if term[:j].isdigit(): - stoichiometry = int(term[:j]) - else: - stoichiometry = float(term[:j]) - species = term[j:] - else: - species = term - stoichiometry = 1 - - if species == 'M' or species == 'm': - thirdBody = True - elif species not in speciesDict: - raise InputParseError('Unexpected species "{0}" in reaction expression "{1}".'.format(species, expression)) - else: - dest.append((stoichiometry, speciesDict[species])) - - return falloff3b, thirdBody - - 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 - 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: - raise InputParseError('Invalid number of reactant species ({0}) for reaction {1}.'.format(rStoich, reaction)) - - # 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, energyUnits), - T0=(1,"K"), - ) - - if len(lines) == 1: - # If there's only one line then we know to use the high-P limit kinetics as-is - reaction.kinetics = arrheniusHigh - else: - # There's more kinetics information to be read - arrheniusLow = None - troe = None - sri = None - chebyshev = None - pdepArrhenius = None - efficiencies = {} - chebyshevCoeffs = [] - - # Note that the subsequent lines could be in any order - for line in lines[1:]: - tokens = line.split('/') - if 'dup' in line.lower(): - # Duplicate reaction - reaction.duplicate = True - - elif 'low' in line.lower(): - # Low-pressure-limit Arrhenius parameters - tokens = tokens[1].split() - arrheniusLow = Arrhenius( - A=(float(tokens[0].strip()),klow_units), - n=float(tokens[1].strip()), - Ea=(float(tokens[2].strip()),"kcal/mol"), - T0=(1,"K"), - ) - - elif 'rev' in line.lower(): - reaction.reversible = False - - # Create a reaction proceeding in the opposite direction - revReaction = Reaction(reactants=reaction.products, - products=reaction.reactants, - reversible=False) - tokens = tokens[1].split() - revReaction.kinetics = Arrhenius( - A=(float(tokens[0].strip()),klow_units), - n=float(tokens[1].strip()), - Ea=(float(tokens[2].strip()),"kcal/mol"), - T0=(1,"K"), - ) - - elif 'ford' in line.lower(): - tokens = tokens[1].split() - reaction.fwdOrders[tokens[0].strip()] = tokens[1].strip() - - elif 'troe' in line.lower(): - # Troe falloff parameters - tokens = tokens[1].split() - alpha = float(tokens[0].strip()) - T3 = float(tokens[1].strip()) - T1 = float(tokens[2].strip()) - try: - T2 = float(tokens[3].strip()) - except (IndexError, ValueError): - T2 = None - - troe = Troe( - alpha=(alpha,''), - T3=(T3,"K"), - T1=(T1,"K"), - T2=(T2,"K") if T2 is not None else None, - ) - elif 'sri' in line.lower(): - # SRI falloff parameters - tokens = tokens[1].split() - A = float(tokens[0].strip()) - B = float(tokens[1].strip()) - C = float(tokens[2].strip()) - try: - D = float(tokens[3].strip()) - E = float(tokens[4].strip()) - except (IndexError, ValueError): - D = None - E = None - - if D is None or E is None: - sri = Sri(A=A, B=B, C=C) - else: - sri = Sri(A=A, B=B, C=C, D=D, E=E) - - elif 'cheb' in line.lower(): - # Chebyshev parameters - if chebyshev is None: - chebyshev = Chebyshev() - tokens = [t.strip() for t in tokens] - if 'TCHEB' in line: - index = tokens.index('TCHEB') - tokens2 = tokens[index+1].split() - chebyshev.Tmin = float(tokens2[0].strip()) - chebyshev.Tmax = float(tokens2[1].strip()) - if 'PCHEB' in line: - index = tokens.index('PCHEB') - tokens2 = tokens[index+1].split() - chebyshev.Pmin = (float(tokens2[0].strip()), 'atm') - chebyshev.Pmax = (float(tokens2[1].strip()), 'atm') - if 'TCHEB' in line or 'PCHEB' in line: - pass - elif chebyshev.degreeT == 0 or chebyshev.degreeP == 0: - tokens2 = tokens[1].split() - chebyshev.degreeT = int(float(tokens2[0].strip())) - chebyshev.degreeP = int(float(tokens2[1].strip())) - chebyshev.coeffs = np.zeros((chebyshev.degreeT,chebyshev.degreeP), np.float64) - else: - tokens2 = tokens[1].split() - chebyshevCoeffs.extend([float(t.strip()) for t in tokens2]) - - elif 'plog' in line.lower(): - # Pressure-dependent Arrhenius parameters - if pdepArrhenius is None: - pdepArrhenius = [] - tokens = tokens[1].split() - 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"), - T0=(1,"K"), - )]) - - else: - # Assume a list of collider efficiencies - for collider, efficiency in zip(tokens[0::2], tokens[1::2]): - efficiencies[collider.strip()] = float(efficiency.strip()) - - # Decide which kinetics to keep and store them on the reaction object - # Only one of these should be true at a time! - if chebyshev is not None: - if chebyshev.Tmin is None or chebyshev.Tmax is None: - raise InputParseError('Missing TCHEB line for reaction {0}'.format(reaction)) - if chebyshev.Pmin is None or chebyshev.Pmax is None: - raise InputParseError('Missing PCHEB line for reaction {0}'.format(reaction)) - index = 0 - for t in range(chebyshev.degreeT): - for p in range(chebyshev.degreeP): - chebyshev.coeffs[t,p] = chebyshevCoeffs[index] - index += 1 - reaction.kinetics = chebyshev - elif pdepArrhenius is not None: - reaction.kinetics = PDepArrhenius( - pressures=([P for P, arrh in pdepArrhenius],"atm"), - arrhenius=[arrh for P, arrh in pdepArrhenius], - ) - elif troe is not None: - troe.arrheniusHigh = arrheniusHigh - troe.arrheniusLow = arrheniusLow - troe.efficiencies = efficiencies - reaction.kinetics = troe - elif sri is not None: - sri.arrheniusHigh = arrheniusHigh - sri.arrheniusLow = arrheniusLow - sri.efficiencies = efficiencies - reaction.kinetics = sri - elif arrheniusLow is not None: - reaction.kinetics = Lindemann(arrheniusHigh=arrheniusHigh, arrheniusLow=arrheniusLow) - reaction.kinetics.efficiencies = efficiencies - elif thirdBody: - reaction.kinetics = ThirdBody(arrheniusHigh=arrheniusHigh) - reaction.kinetics.efficiencies = efficiencies + def warn(self, message): + if self.warning_as_error: + raise InputParseError(message) else: - reaction.kinetics = arrheniusHigh + logging.warning(message) - return reaction, revReaction + def parseComposition(self, elements, nElements, width): + """ + Parse the elemental composition from a 7 or 9 coefficient NASA polynomial + entry. + """ + composition = {} + for i in range(nElements): + symbol = elements[width*i:width*i+2].strip() + count = elements[width*i+2:width*i+width].strip() + if not symbol: + continue + try: + count = int(float(count)) + if count: + composition[symbol.capitalize()] = count + except ValueError: + pass + return composition + def readThermoEntry(self, entry, TintDefault): + """ + Read a thermodynamics `entry` for one species in a Chemkin-format file + (consisting of two 7-coefficient NASA polynomials). Returns the label of + the species, the thermodynamics model as a :class:`MultiNASA` object, the + elemental composition of the species, and the comment/note associated with + the thermo entry. + """ + lines = entry.splitlines() + identifier = lines[0][0:24].split() + species = identifier[0].strip() -def loadChemkinFile(path, speciesList=None): - """ - Load a Chemkin-format input file to `path` on disk, returning lists of - the species and reactions in the Chemkin file. - """ - elementList = [] - speciesDict = {} - if speciesList is None: - speciesList = [] - else: - for species in speciesList: - speciesDict[species.label] = species - - reactionList = [] - transportLines = [] - - def removeCommentFromLine(line): - if '!' in line: - index = line.index('!') - comment = line[index+1:-1] - line = line[0:index] + '\n' - return line, comment + if len(identifier) > 1: + note = ''.join(identifier[1:]).strip() else: - comment = '' - return line, comment + note = '' - with open(path, 'r') as f: - line = f.readline() - while line != '': - line = removeCommentFromLine(line)[0] - line = line.strip() - tokens = line.split() + # Extract the NASA polynomial coefficients + # Remember that the high-T polynomial comes first! + try: + Tmin = fortFloat(lines[0][45:55]) + Tmax = fortFloat(lines[0][55:65]) + try: + Tint = fortFloat(lines[0][65:75]) + except ValueError: + Tint = TintDefault - if 'ELEMENTS' in line: - index = tokens.index('ELEMENTS') - tokens = tokens[index+1:] - while 'END' not in tokens: - line = f.readline() - line = removeCommentFromLine(line)[0] - line = line.strip() - tokens.extend(line.split()) + coeffs_high = [fortFloat(lines[i][j:k]) + for i,j,k in [(1,0,15), (1,15,30), (1,30,45), (1,45,60), + (1,60,75), (2,0,15), (2,15,30)]] + coeffs_low = [fortFloat(lines[i][j:k]) + for i,j,k in [(2,30,45), (2,45,60), (2,60,75), (3,0,15), + (3,15,30), (3,30,45), (3,45,60)]] - for token in tokens: - if token == 'END': + except (IndexError, ValueError): + raise InputParseError('Error while reading thermo entry for species {0}'.format(species)) + + composition = self.parseComposition(lines[0][24:44], 4, 5) + + # Non-standard extended elemental composition data may be located beyond + # column 80 on the first line of the thermo entry + if len(lines[0]) > 80: + elements = lines[0][80:] + composition2 = self.parseComposition(elements, len(elements)/10, 10) + composition.update(composition2) + + # Construct and return the thermodynamics model + thermo = MultiNASA( + polynomials=[ + NASA(Tmin=(Tmin,"K"), Tmax=(Tint,"K"), coeffs=coeffs_low), + NASA(Tmin=(Tint,"K"), Tmax=(Tmax,"K"), coeffs=coeffs_high) + ], + Tmin=(Tmin,"K"), + Tmax=(Tmax,"K"), + ) + + return species, thermo, composition, note + + def readNasa9Entry(self, entry): + """ + Read a thermodynamics `entry` for one species given as one or more + 9-coefficient NASA polynomials, written in the format described in + Appendix A of NASA Reference Publication 1311 (McBride and Gordon, 1996). + Returns the label of the species, the thermodynamics model as a + :class:`MultiNASA` object, the elemental composition of the species, and + the comment/note associated with the thermo entry. + """ + tokens = entry[0].split() + species = tokens[0] + note = ' '.join(tokens[1:]) + N = int(entry[1][:2]) + note2 = entry[1][3:9].strip() + if note and note2: + note = '{0} [{1}]'.format(note, note2) + elif note2: + note = note2 + + composition = self.parseComposition(entry[1][10:50], 5, 8) + + polys = [] + totalTmin = 1e100 + totalTmax = -1e100 + try: + for i in range(N): + A,B,C = entry[2+3*i:2+3*(i+1)] + Tmin = fortFloat(A[1:11]) + Tmax = fortFloat(A[11:21]) + coeffs = [fortFloat(B[0:16]), fortFloat(B[16:32]), + fortFloat(B[32:48]), fortFloat(B[48:64]), + fortFloat(B[64:80]), fortFloat(C[0:16]), + fortFloat(C[16:32]), fortFloat(C[48:64]), + fortFloat(C[64:80])] + polys.append(NASA(Tmin=(Tmin,"K"), Tmax=(Tmax,"K"), coeffs=coeffs)) + totalTmin = min(Tmin, totalTmin) + totalTmax = max(Tmax, totalTmax) + except (IndexError, ValueError): + raise InputParseError('Error while reading thermo entry for species {0}'.format(species)) + + thermo = MultiNASA(polynomials=polys, + Tmin=(totalTmin,"K"), + Tmax=(totalTmax,"K")) + + return species, thermo, composition, note + + def readKineticsEntry(self, entry): + """ + Read a kinetics `entry` for a single reaction as loaded from a + Chemkin-format file. Returns a :class:`Reaction` object with the + reaction and its associated kinetics. + """ + + lines = entry.strip().splitlines() + + # The first line contains the reaction equation and a set of modified Arrhenius parameters + tokens = lines[0].split() + A = float(tokens[-3]) + n = float(tokens[-2]) + Ea = float(tokens[-1]) + reaction = ''.join(tokens[:-3]) + revReaction = None + + # Split the reaction equation into reactants and products + if '<=>' in reaction: + reversible = True + reactants, products = reaction.split('<=>') + elif '=>' in reaction: + reversible = False + reactants, products = reaction.split('=>') + elif '=' in reaction: + reversible = True + reactants, products = reaction.split('=') + else: + raise InputParseError("Failed to find reactant/product delimiter in reaction string.") + + # Create a new Reaction object for this reaction + reaction = Reaction(reactants=[], products=[], reversible=reversible) + + def parseExpression(expression, dest): + 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 self.speciesDict: + if re.search(r'\(\+%s\)' % re.escape(species), expression): + falloff3b = species + expression = re.sub(r'(\(\+%s\))' % re.escape(species), + '', expression) break - elementList.append(token.capitalize()) - elif 'SPECIES' in line: - # List of species identifiers - index = tokens.index('SPECIES') - tokens = tokens[index+1:] - while 'END' not in tokens: - line = f.readline() - line = removeCommentFromLine(line)[0] - line = line.strip() - tokens.extend(line.split()) - - for token in tokens: - if token == 'END': - break - if token in speciesDict: - species = speciesDict[token] + for term in expression.split('+'): + term = term.strip() + if not term[0].isalpha(): + # This allows for for non-unity stoichiometric coefficients, e.g. + # 2A=B+C or .85A+.15B=>C + j = [i for i,c in enumerate(term) if c.isalpha()][0] + if term[:j].isdigit(): + stoichiometry = int(term[:j]) else: - species = Species(label=token) - speciesDict[token] = species - speciesList.append(species) + stoichiometry = float(term[:j]) + species = term[j:] + else: + species = term + stoichiometry = 1 - elif 'THERM' in line.upper() and 'NASA9' in line: - entryPosition = 0 - entryLength = None - entry = [] - while not line.startswith('END'): - line = f.readline() - line = removeCommentFromLine(line)[0] - if not line: - continue + if species == 'M' or species == 'm': + thirdBody = True + elif species not in self.speciesDict: + raise InputParseError('Unexpected species "{0}" in reaction expression "{1}".'.format(species, expression)) + else: + dest.append((stoichiometry, self.speciesDict[species])) - if entryLength is None: - entryLength = 0 - # special case if (redundant) temperature ranges are - # given as the first line - try: - s = line.split() - float(s[0]), float(s[1]), float(s[2]) + return falloff3b, thirdBody + + 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 + 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: + raise InputParseError('Invalid number of reactant species ({0}) for reaction {1}.'.format(rStoich, reaction)) + + # 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), + T0=(1,"K"), + ) + + if len(lines) == 1: + # If there's only one line then we know to use the high-P limit kinetics as-is + reaction.kinetics = arrheniusHigh + else: + # There's more kinetics information to be read + arrheniusLow = None + troe = None + sri = None + chebyshev = None + pdepArrhenius = None + efficiencies = {} + chebyshevCoeffs = [] + + # Note that the subsequent lines could be in any order + for line in lines[1:]: + tokens = line.split('/') + if 'dup' in line.lower(): + # Duplicate reaction + reaction.duplicate = True + + elif 'low' in line.lower(): + # Low-pressure-limit Arrhenius parameters + tokens = tokens[1].split() + arrheniusLow = Arrhenius( + A=(float(tokens[0].strip()),klow_units), + n=float(tokens[1].strip()), + Ea=(float(tokens[2].strip()),"kcal/mol"), + T0=(1,"K"), + ) + + elif 'rev' in line.lower(): + reaction.reversible = False + + # Create a reaction proceeding in the opposite direction + revReaction = Reaction(reactants=reaction.products, + products=reaction.reactants, + reversible=False) + tokens = tokens[1].split() + revReaction.kinetics = Arrhenius( + A=(float(tokens[0].strip()),klow_units), + n=float(tokens[1].strip()), + Ea=(float(tokens[2].strip()),"kcal/mol"), + T0=(1,"K"), + ) + + elif 'ford' in line.lower(): + tokens = tokens[1].split() + reaction.fwdOrders[tokens[0].strip()] = tokens[1].strip() + + elif 'troe' in line.lower(): + # Troe falloff parameters + tokens = tokens[1].split() + alpha = float(tokens[0].strip()) + T3 = float(tokens[1].strip()) + T1 = float(tokens[2].strip()) + try: + T2 = float(tokens[3].strip()) + except (IndexError, ValueError): + T2 = None + + troe = Troe( + alpha=(alpha,''), + T3=(T3,"K"), + T1=(T1,"K"), + T2=(T2,"K") if T2 is not None else None, + ) + elif 'sri' in line.lower(): + # SRI falloff parameters + tokens = tokens[1].split() + A = float(tokens[0].strip()) + B = float(tokens[1].strip()) + C = float(tokens[2].strip()) + try: + D = float(tokens[3].strip()) + E = float(tokens[4].strip()) + except (IndexError, ValueError): + D = None + E = None + + if D is None or E is None: + sri = Sri(A=A, B=B, C=C) + else: + sri = Sri(A=A, B=B, C=C, D=D, E=E) + + elif 'cheb' in line.lower(): + # Chebyshev parameters + if chebyshev is None: + chebyshev = Chebyshev() + tokens = [t.strip() for t in tokens] + if 'TCHEB' in line: + index = tokens.index('TCHEB') + tokens2 = tokens[index+1].split() + chebyshev.Tmin = float(tokens2[0].strip()) + chebyshev.Tmax = float(tokens2[1].strip()) + if 'PCHEB' in line: + index = tokens.index('PCHEB') + tokens2 = tokens[index+1].split() + chebyshev.Pmin = (float(tokens2[0].strip()), 'atm') + chebyshev.Pmax = (float(tokens2[1].strip()), 'atm') + if 'TCHEB' in line or 'PCHEB' in line: + pass + elif chebyshev.degreeT == 0 or chebyshev.degreeP == 0: + tokens2 = tokens[1].split() + chebyshev.degreeT = int(float(tokens2[0].strip())) + chebyshev.degreeP = int(float(tokens2[1].strip())) + chebyshev.coeffs = np.zeros((chebyshev.degreeT,chebyshev.degreeP), np.float64) + else: + tokens2 = tokens[1].split() + chebyshevCoeffs.extend([float(t.strip()) for t in tokens2]) + + elif 'plog' in line.lower(): + # Pressure-dependent Arrhenius parameters + if pdepArrhenius is None: + pdepArrhenius = [] + tokens = tokens[1].split() + 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"), + T0=(1,"K"), + )]) + + else: + # Assume a list of collider efficiencies + for collider, efficiency in zip(tokens[0::2], tokens[1::2]): + efficiencies[collider.strip()] = float(efficiency.strip()) + + # Decide which kinetics to keep and store them on the reaction object + # Only one of these should be true at a time! + if chebyshev is not None: + if chebyshev.Tmin is None or chebyshev.Tmax is None: + raise InputParseError('Missing TCHEB line for reaction {0}'.format(reaction)) + if chebyshev.Pmin is None or chebyshev.Pmax is None: + raise InputParseError('Missing PCHEB line for reaction {0}'.format(reaction)) + index = 0 + for t in range(chebyshev.degreeT): + for p in range(chebyshev.degreeP): + chebyshev.coeffs[t,p] = chebyshevCoeffs[index] + index += 1 + reaction.kinetics = chebyshev + elif pdepArrhenius is not None: + reaction.kinetics = PDepArrhenius( + pressures=([P for P, arrh in pdepArrhenius],"atm"), + arrhenius=[arrh for P, arrh in pdepArrhenius], + ) + elif troe is not None: + troe.arrheniusHigh = arrheniusHigh + troe.arrheniusLow = arrheniusLow + troe.efficiencies = efficiencies + reaction.kinetics = troe + elif sri is not None: + sri.arrheniusHigh = arrheniusHigh + sri.arrheniusLow = arrheniusLow + sri.efficiencies = efficiencies + reaction.kinetics = sri + elif arrheniusLow is not None: + reaction.kinetics = Lindemann(arrheniusHigh=arrheniusHigh, arrheniusLow=arrheniusLow) + reaction.kinetics.efficiencies = efficiencies + elif thirdBody: + reaction.kinetics = ThirdBody(arrheniusHigh=arrheniusHigh) + reaction.kinetics.efficiencies = efficiencies + else: + reaction.kinetics = arrheniusHigh + + return reaction, revReaction + + def loadChemkinFile(self, path): + """ + Load a Chemkin-format input file to `path` on disk. + """ + + transportLines = [] + + def removeCommentFromLine(line): + if '!' in line: + index = line.index('!') + comment = line[index+1:-1] + line = line[0:index] + '\n' + return line, comment + else: + comment = '' + return line, comment + + with open(path, 'r') as f: + line = f.readline() + while line != '': + line = removeCommentFromLine(line)[0] + line = line.strip() + tokens = line.split() + + if 'ELEMENTS' in line: + index = tokens.index('ELEMENTS') + tokens = tokens[index+1:] + while 'END' not in tokens: + line = f.readline() + line = removeCommentFromLine(line)[0] + line = line.strip() + tokens.extend(line.split()) + + for token in tokens: + if token == 'END': + break + self.elements.append(token.capitalize()) + + elif 'SPECIES' in line: + # List of species identifiers + index = tokens.index('SPECIES') + tokens = tokens[index+1:] + while 'END' not in tokens: + line = f.readline() + line = removeCommentFromLine(line)[0] + line = line.strip() + tokens.extend(line.split()) + + for token in tokens: + if token == 'END': + break + if token in self.speciesDict: + species = self.speciesDict[token] + else: + species = Species(label=token) + self.speciesDict[token] = species + self.speciesList.append(species) + + elif 'THERM' in line.upper() and 'NASA9' in line: + entryPosition = 0 + entryLength = None + entry = [] + while not line.startswith('END'): + line = f.readline() + line = removeCommentFromLine(line)[0] + if not line: continue - except (IndexError, ValueError): - pass - if entryPosition == 0: - entry.append(line) - elif entryPosition == 1: - entryLength = 2 + 3 * int(line.split()[0]) - entry.append(line) - elif entryPosition < entryLength: - entry.append(line) - - if entryPosition == entryLength-1: - label, thermo, comp, note = readNasa9Entry(entry) - try: - speciesDict[label].thermo = thermo - speciesDict[label].composition = comp - speciesDict[label].note = note - except KeyError: - logging.info('Skipping unexpected species "{0}" while reading thermodynamics entry.'.format(label)) - - entryPosition = -1 - entry = [] - - entryPosition += 1 - - elif 'THERM' in line: - # List of thermodynamics (hopefully one per species!) - line = f.readline() - TintDefault = float(line.split()[1]) - thermo = '' - while line != '' and 'END' not in line: - line = removeCommentFromLine(line)[0] - if len(line) >= 80 and line[79] in ['1', '2', '3', '4']: - thermo += line - if line[79] == '4': - label, thermo, comp, note = readThermoEntry(thermo, TintDefault) + if entryLength is None: + entryLength = 0 + # special case if (redundant) temperature ranges are + # given as the first line try: - speciesDict[label].thermo = thermo - speciesDict[label].composition = comp - speciesDict[label].note = note + s = line.split() + float(s[0]), float(s[1]), float(s[2]) + continue + except (IndexError, ValueError): + pass + + if entryPosition == 0: + entry.append(line) + elif entryPosition == 1: + entryLength = 2 + 3 * int(line.split()[0]) + entry.append(line) + elif entryPosition < entryLength: + entry.append(line) + + if entryPosition == entryLength-1: + label, thermo, comp, note = self.readNasa9Entry(entry) + try: + self.speciesDict[label].thermo = thermo + self.speciesDict[label].composition = comp + self.speciesDict[label].note = note except KeyError: logging.info('Skipping unexpected species "{0}" while reading thermodynamics entry.'.format(label)) - thermo = '' + + entryPosition = -1 + entry = [] + + entryPosition += 1 + + elif 'THERM' in line: + # List of thermodynamics (hopefully one per species!) line = f.readline() + TintDefault = float(line.split()[1]) + thermo = '' + while line != '' and 'END' not in line: + line = removeCommentFromLine(line)[0] + if len(line) >= 80 and line[79] in ['1', '2', '3', '4']: + thermo += line + if line[79] == '4': + label, thermo, comp, note = self.readThermoEntry(thermo, TintDefault) + try: + self.speciesDict[label].thermo = thermo + self.speciesDict[label].composition = comp + self.speciesDict[label].note = note + except KeyError: + logging.info('Skipping unexpected species "{0}" while reading thermodynamics entry.'.format(label)) + thermo = '' + line = f.readline() - elif 'REACTIONS' in line: - # Reactions section - energyUnits = 'CAL/MOL' - moleculeUnits = 'MOLES' - try: - energyUnits = tokens[1] - moleculeUnits = tokens[2] - except IndexError: - pass + elif 'REACTIONS' in line: + # Reactions section + energyUnits = 'CAL/MOL' + moleculeUnits = 'MOLES' + try: + energyUnits = tokens[1] + moleculeUnits = tokens[2] + except IndexError: + pass - global PROCESSED_UNITS, ENERGY_UNITS, QUANTITY_UNITS - if not PROCESSED_UNITS: - PROCESSED_UNITS = True - ENERGY_UNITS = UNIT_OPTIONS[energyUnits] - QUANTITY_UNITS = UNIT_OPTIONS[moleculeUnits] - else: - if (ENERGY_UNITS != UNIT_OPTIONS[energyUnits] or - QUANTITY_UNITS != UNIT_OPTIONS[moleculeUnits]): - raise InputParseError("Multiple REACTIONS sections with " - "different units are not supported.") + #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] + else: + if (self.energy_units != UNIT_OPTIONS[energyUnits] or + self.quantity_units != UNIT_OPTIONS[moleculeUnits]): + raise InputParseError("Multiple REACTIONS sections with " + "different units are not supported.") - kineticsList = [] - commentsList = [] - kinetics = '' - comments = '' + kineticsList = [] + commentsList = [] + kinetics = '' + comments = '' - line = f.readline() - while line != '' and 'END' not in line: + line = f.readline() + while line != '' and 'END' not in line: - lineStartsWithComment = line.startswith('!') - line, comment = removeCommentFromLine(line) - line = line.strip() - comment = comment.strip() + lineStartsWithComment = line.startswith('!') + line, comment = removeCommentFromLine(line) + line = line.strip() + comment = comment.strip() - if '=' in line and not lineStartsWithComment: - # Finish previous record + if '=' in line and not lineStartsWithComment: + # Finish previous record + kineticsList.append(kinetics) + commentsList.append(comments) + kinetics = '' + comments = '' + + if line: + kinetics += line + '\n' + if comment: + comments += comment + '\n' + + line = f.readline() + + # Don't forget the last reaction! + if kinetics.strip() != '': kineticsList.append(kinetics) commentsList.append(comments) - kinetics = '' - comments = '' - if line: - kinetics += line + '\n' - if comment: - comments += comment + '\n' - - line = f.readline() - - # Don't forget the last reaction! - if kinetics.strip() != '': - kineticsList.append(kinetics) - commentsList.append(comments) - - if kineticsList[0] == '' and commentsList[-1] == '': - # True for mechanism files generated from RMG-Py - kineticsList.pop(0) - commentsList.pop(-1) - elif kineticsList[0] == '' and commentsList[0] == '': - # True for mechanism files generated from RMG-Java - kineticsList.pop(0) - commentsList.pop(0) - else: - # In reality, comments can occur anywhere in the mechanism - # file (e.g. either or both of before and after the - # reaction equation) - # If we can't tell what semantics we are using, then just - # throw the comments away - # (This is better than failing to load the mechanism file at - # all, which would likely occur otherwise) - if kineticsList[0] == '': + if kineticsList[0] == '' and commentsList[-1] == '': + # True for mechanism files generated from RMG-Py kineticsList.pop(0) - if len(kineticsList) != len(commentsList): - commentsList = ['' for kinetics in kineticsList] + commentsList.pop(-1) + elif kineticsList[0] == '' and commentsList[0] == '': + # True for mechanism files generated from RMG-Java + kineticsList.pop(0) + commentsList.pop(0) + else: + # In reality, comments can occur anywhere in the mechanism + # file (e.g. either or both of before and after the + # reaction equation) + # If we can't tell what semantics we are using, then just + # throw the comments away + # (This is better than failing to load the mechanism file at + # all, which would likely occur otherwise) + if kineticsList[0] == '': + kineticsList.pop(0) + if len(kineticsList) != len(commentsList): + commentsList = ['' for kinetics in kineticsList] - for kinetics, comments in zip(kineticsList, commentsList): - reaction,revReaction = readKineticsEntry(kinetics, speciesDict, energyUnits, moleculeUnits) - reactionList.append(reaction) - if revReaction is not None: - reactionList.append(revReaction) + for kinetics, comments in zip(kineticsList, commentsList): + reaction,revReaction = self.readKineticsEntry(kinetics) + self.reactions.append(reaction) + if revReaction is not None: + self.reactions.append(revReaction) + + elif 'TRAN' in line: + line = f.readline() + while 'END' not in line: + transportLines.append(line) - elif 'TRAN' in line: line = f.readline() - while 'END' not in line: - transportLines.append(line) - line = f.readline() + # Check for marked (and unmarked!) duplicate reactions + # Raise exception for unmarked duplicate reactions + for index1 in range(len(self.reactions)): + reaction1 = self.reactions[index1] + for index2 in range(index1+1, len(self.reactions)): + reaction2 = self.reactions[index2] + if reaction1.reactants == reaction2.reactants and reaction1.products == reaction2.products: + if reaction1.duplicate and reaction2.duplicate: + pass + elif reaction1.kinetics.isPressureDependent() == reaction2.kinetics.isPressureDependent(): + # If both reactions are pressure-independent or both are pressure-dependent, then they need duplicate tags + # pdep and non-pdep reactions are treated as different, so those are okay + raise InputParseError('Encountered unmarked duplicate reaction {0}.'.format(reaction1)) - # Check for marked (and unmarked!) duplicate reactions - # Raise exception for unmarked duplicate reactions - for index1 in range(len(reactionList)): - reaction1 = reactionList[index1] - for index2 in range(index1+1, len(reactionList)): - reaction2 = reactionList[index2] - if reaction1.reactants == reaction2.reactants and reaction1.products == reaction2.products: - if reaction1.duplicate and reaction2.duplicate: - pass - elif reaction1.kinetics.isPressureDependent() == reaction2.kinetics.isPressureDependent(): - # If both reactions are pressure-independent or both are pressure-dependent, then they need duplicate tags - # pdep and non-pdep reactions are treated as different, so those are okay - raise InputParseError('Encountered unmarked duplicate reaction {0}.'.format(reaction1)) + index = 0 + for reaction in self.reactions: + index += 1 + reaction.index = index - index = 0 - for reaction in reactionList: - index += 1 - reaction.index = index + if transportLines: + self.parseTransportData(transportLines) - if transportLines: - parseTransportData(transportLines, speciesList) + def parseTransportData(self, lines): + """ + Parse the Chemkin-format transport data in ``lines`` (a list of strings) + and add that transport data to the previously-loaded species. + """ - return elementList, speciesList, reactionList + for line in lines: + line = line.strip() + if not line or line.startswith('!'): + continue + if line.startswith('END'): + break + data = line.split() + if len(data) < 7: + raise InputParseError('Unable to parse transport data: not enough parameters') + if len(data) >= 8: + # comment may contain spaces. Rejoin into a single field. + comment = ''.join(data[7:]).lstrip('!') + data = data[:7] + [comment] -def parseTransportData(lines, speciesList): - """ - Parse the Chemkin-format transport data in ``lines`` (a list of strings) - and add that transport data to the species in ``speciesList``. - """ - speciesDict = dict((species.label, species) for species in speciesList) - for line in lines: - line = line.strip() - if not line or line.startswith('!'): - continue - if line.startswith('END'): - break + speciesName = data[0] + if speciesName in self.speciesDict: + if self.speciesDict[speciesName].transport is None: + self.speciesDict[speciesName].transport = TransportData(*data) + else: + self.warn('Ignoring duplicate transport data' + ' for species "{0}".'.format(speciesName)) - data = line.split() - if len(data) < 7: - raise InputParseError('Unable to parse transport data: not enough parameters') - if len(data) >= 8: - # comment may contain spaces. Rejoin into a single field. - comment = ''.join(data[7:]).lstrip('!') - data = data[:7] + [comment] + def writeCTI(self, header=None, name='gas', transportModel='Mix', + outName='mech.cti'): - speciesName = data[0] - if speciesName in speciesDict: - if speciesDict[speciesName].transport is None: - speciesDict[speciesName].transport = TransportData(*data) - else: - warn('Ignoring duplicate transport data' - ' for species "{0}".'.format(speciesName)) + delimiterLine = '#' + '-'*79 + haveTransport = True + speciesNameLength = 1 + elementsFromSpecies = set() + for s in self.speciesList: + if not s.transport: + haveTransport = False + if s.composition is None: + raise InputParseError('No thermo data found for species: {0!r}'.format(s.label)) + elementsFromSpecies.update(s.composition) + speciesNameLength = max(speciesNameLength, len(s.label)) + # validate list of elements + missingElements = elementsFromSpecies - set(self.elements) + if missingElements: + raise InputParseError('Undefined elements: ' + str(missingElements)) -def writeCTI(elements, - species, - reactions=None, - header=None, - name='gas', - transportModel='Mix', - outName='mech.cti'): + speciesNames = [''] + for i,s in enumerate(self.speciesList): + if i and not i % 5: + speciesNames.append(' '*21) + speciesNames[-1] += '{0:{1}s}'.format(s.label, speciesNameLength+2) - delimiterLine = '#' + '-'*79 - haveTransport = True - speciesNameLength = 1 - elementsFromSpecies = set() - for s in species: - if not s.transport: - haveTransport = False - if s.composition is None: - raise InputParseError('No thermo data found for species: {0!r}'.format(s.label)) - elementsFromSpecies.update(s.composition) - speciesNameLength = max(speciesNameLength, len(s.label)) + speciesNames = '\n'.join(speciesNames).strip() - # validate list of elements - missingElements = elementsFromSpecies - set(elements) - if missingElements: - raise InputParseError('Undefined elements: ' + str(missingElements)) + lines = [] + if header: + lines.extend(header) - speciesNames = [''] - for i,s in enumerate(species): - if i and not i % 5: - speciesNames.append(' '*21) - speciesNames[-1] += '{0:{1}s}'.format(s.label, speciesNameLength+2) + # Write the gas definition + lines.append("units(length='cm', time='s', quantity={0!r}, act_energy={1!r})".format(self.quantity_units, self.energy_units)) + lines.append('') + lines.append('ideal_gas(name={0!r},'.format(name)) + lines.append(' elements="{0}",'.format(' '.join(self.elements))) + lines.append(' species="""{0}""",'.format(speciesNames)) + if self.reactions: + lines.append(" reactions='all',") + if haveTransport: + lines.append(" transport={0!r},".format(transportModel)) + lines.append(' initial_state=state(temperature=300.0, pressure=OneAtm))') + lines.append('') - speciesNames = '\n'.join(speciesNames).strip() + # Write the individual species data + lines.append(delimiterLine) + lines.append('# Species data') + lines.append(delimiterLine) + lines.append('') - lines = [] - if header: - lines.extend(header) + for s in self.speciesList: + lines.append(s.to_cti()) - # Write the gas definition - lines.append("units(length='cm', time='s', quantity={0!r}, act_energy={1!r})".format(QUANTITY_UNITS, ENERGY_UNITS)) - lines.append('') - lines.append('ideal_gas(name={0!r},'.format(name)) - lines.append(' elements="{0}",'.format(' '.join(elements))) - lines.append(' species="""{0}""",'.format(speciesNames)) - if reactions: - lines.append(" reactions='all',") - if haveTransport: - lines.append(" transport={0!r},".format(transportModel)) - lines.append(' initial_state=state(temperature=300.0, pressure=OneAtm))') - lines.append('') + # Write the reactions + lines.append(delimiterLine) + lines.append('# Reaction data') + lines.append(delimiterLine) - # Write the individual species data - lines.append(delimiterLine) - lines.append('# Species data') - lines.append(delimiterLine) - lines.append('') + for i,r in enumerate(self.reactions): + lines.append('\n# Reaction {0}'.format(i+1)) + lines.append(r.to_cti()) - for s in species: - lines.append(s.to_cti()) + lines.append('') - # Write the reactions - lines.append(delimiterLine) - lines.append('# Reaction data') - lines.append(delimiterLine) + f = open(outName, 'w') + f.write('\n'.join(lines)) - for i,r in enumerate(reactions): - lines.append('\n# Reaction {0}'.format(i+1)) - lines.append(r.to_cti()) - - lines.append('') - - f = open(outName, 'w') - f.write('\n'.join(lines)) - - -def showHelp(): - print """ + def showHelp(self): + print """ ck2cti.py: Convert Chemkin-format mechanisms to Cantera input files (.cti) Usage: @@ -1570,40 +1552,38 @@ duplicate transport data) to be ignored. """ + def convertMech(self, inputFile, thermoFile=None, + transportFile=None, phaseName='gas', + outName=None, quiet=False, permissive=None): + if quiet: + logging.basicConfig(level=logging.ERROR) -def convertMech(inputFile, thermoFile=None, - transportFile=None, phaseName='gas', - outName=None, quiet=False, permissive=None): - if quiet: - logging.basicConfig(level=logging.ERROR) + if permissive is not None: + self.warning_as_error = not permissive - if permissive is not None: - global WARNING_AS_ERROR - WARNING_AS_ERROR = not permissive + # Read input mechanism files + self.loadChemkinFile(inputFile) - # Read input mechanism files - elements, species, reactions = loadChemkinFile(inputFile) + if thermoFile: + self.loadChemkinFile(thermoFile) - if thermoFile: - _, species, _ = loadChemkinFile(thermoFile, species) + if transportFile: + lines = open(transportFile).readlines() + self.parseTransportData(lines) - if transportFile: - lines = open(transportFile).readlines() - parseTransportData(lines, species) + # Transport validation: make sure all species have transport data + for s in self.speciesList: + if s.transport is None: + raise InputParseError("No transport data for species '{0}'.".format(s)) - # Transport validation: make sure all species have transport data - for s in species: - if s.transport is None: - raise InputParseError("No transport data for species '{0}'.".format(s)) + if not outName: + outName = os.path.splitext(inputFile)[0] + '.cti' - if not outName: - outName = os.path.splitext(inputFile)[0] + '.cti' - - # Write output file - writeCTI(elements, species, reactions, name=phaseName, outName=outName) - if not quiet: - print 'Wrote CTI mechanism file to {0!r}.'.format(outName) - print 'Mechanism contains {0} species and {1} reactions.'.format(len(species), len(reactions)) + # Write output file + self.writeCTI(name=phaseName, outName=outName) + if not quiet: + print 'Wrote CTI mechanism file to {0!r}.'.format(outName) + print 'Mechanism contains {0} species and {1} reactions.'.format(len(self.speciesList), len(self.reactions)) if __name__ == '__main__': @@ -1629,8 +1609,10 @@ if __name__ == '__main__': print 'Run "ck2cti.py --help" to see usage help.' sys.exit(1) + parser = Parser() + if not options or '-h' in options or '--help' in options: - showHelp() + parser.showHelp() sys.exit(0) if '--input' in options: @@ -1646,12 +1628,10 @@ if __name__ == '__main__': else: outName = None - if '--permissive' in options: - WARNING_AS_ERROR = False - + permissive = '--permissive' in options thermoFile = options.get('--thermo') transportFile = options.get('--transport') - phaseName = options.get('--id', 'gas') - convertMech(inputFile, thermoFile, transportFile, phaseName, outName) + parser.convertMech(inputFile, thermoFile, transportFile, phaseName, + outName, permissive=permissive) diff --git a/test/python/testConvert.py b/test/python/testConvert.py index 1914bd517..8125180cb 100644 --- a/test/python/testConvert.py +++ b/test/python/testConvert.py @@ -7,6 +7,12 @@ import ck2cti import utilities import Cantera as ct + +def convertMech(*args, **kwargs): + parser = ck2cti.Parser() + parser.convertMech(*args, **kwargs) + + class chemkinConverterTest(utilities.CanteraTest): def checkConversion(self, refFile, testFile): ref = ct.IdealGasMix(refFile) @@ -56,9 +62,9 @@ class chemkinConverterTest(utilities.CanteraTest): if os.path.exists('gri30_test.cti'): os.remove('gri30_test.cti') - ck2cti.convertMech('../../data/inputs/gri30.inp', - transportFile='../../data/transport/gri30_tran.dat', - outName='gri30_test.cti', quiet=True) + convertMech('../../data/inputs/gri30.inp', + transportFile='../../data/transport/gri30_tran.dat', + outName='gri30_test.cti', quiet=True) ref, gas = self.checkConversion('gri30.xml', 'gri30_test.cti') self.checkKinetics(ref, gas, [300, 1500], [5e3, 1e5, 2e6]) @@ -67,9 +73,9 @@ class chemkinConverterTest(utilities.CanteraTest): if os.path.exists('soot_test.cti'): os.remove('soot_test.cti') - ck2cti.convertMech('../data/soot.inp', - thermoFile='../data/soot-therm.dat', - outName='soot_test.cti', quiet=True) + convertMech('../data/soot.inp', + thermoFile='../data/soot-therm.dat', + outName='soot_test.cti', quiet=True) ref, gas = self.checkConversion('../data/soot.xml', 'soot_test.cti') self.checkThermo(ref, gas, [300, 1100]) @@ -79,8 +85,8 @@ class chemkinConverterTest(utilities.CanteraTest): if os.path.exists('pdep_test.cti'): os.remove('pdep_test.cti') - ck2cti.convertMech('../data/pdep-test.inp', - outName='pdep_test.cti', quiet=True) + convertMech('../data/pdep-test.inp', + outName='pdep_test.cti', quiet=True) ref, gas = self.checkConversion('../data/pdep-test.xml', 'pdep_test.cti') self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) @@ -91,23 +97,23 @@ class chemkinConverterTest(utilities.CanteraTest): os.remove('h2o2_missingElement.cti') self.assertRaises(ck2cti.InputParseError, - lambda: ck2cti.convertMech('../data/h2o2_missingElement.inp', - quiet=True)) + lambda: convertMech('../data/h2o2_missingElement.inp', + quiet=True)) def test_missingThermo(self): if os.path.exists('h2o2_missingThermo.cti'): os.remove('h2o2_missingThermo.cti') self.assertRaises(ck2cti.InputParseError, - lambda: ck2cti.convertMech('../data/h2o2_missingThermo.inp', - quiet=True)) + lambda: convertMech('../data/h2o2_missingThermo.inp', + quiet=True)) def test_nasa9(self): if os.path.exists('nasa9_test.cti'): os.remove('nasa9_test.cti') - ck2cti.convertMech('../data/nasa9-test.inp', - thermoFile='../data/nasa9-test-therm.dat', - outName='nasa9_test.cti', quiet=True) + convertMech('../data/nasa9-test.inp', + thermoFile='../data/nasa9-test-therm.dat', + outName='nasa9_test.cti', quiet=True) ref, gas = self.checkConversion('../data/nasa9-test.xml', 'nasa9_test.cti') @@ -116,9 +122,9 @@ class chemkinConverterTest(utilities.CanteraTest): def test_sri_falloff(self): if os.path.exists('sri-falloff.cti'): os.remove('sri-falloff.cti') - ck2cti.convertMech('../data/sri-falloff.inp', - thermoFile='../data/dummy-thermo.dat', - outName='sri-falloff.cti', quiet=True) + convertMech('../data/sri-falloff.inp', + thermoFile='../data/dummy-thermo.dat', + outName='sri-falloff.cti', quiet=True) ref, gas = self.checkConversion('../data/sri-falloff.xml', 'sri-falloff.cti') @@ -127,8 +133,8 @@ class chemkinConverterTest(utilities.CanteraTest): 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) + 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') @@ -137,9 +143,9 @@ class chemkinConverterTest(utilities.CanteraTest): def test_explicit_reverse_rate(self): if os.path.exists('explicit-reverse-rate.cti'): os.remove('explicit-reverse-rate.cti') - ck2cti.convertMech('../data/explicit-reverse-rate.inp', - thermoFile='../data/dummy-thermo.dat', - outName='explicit-reverse-rate.cti', quiet=True) + convertMech('../data/explicit-reverse-rate.inp', + thermoFile='../data/dummy-thermo.dat', + outName='explicit-reverse-rate.cti', quiet=True) ref, gas = self.checkConversion('../data/explicit-reverse-rate.xml', 'explicit-reverse-rate.cti') self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) @@ -157,9 +163,9 @@ class chemkinConverterTest(utilities.CanteraTest): def test_explicit_forward_order(self): if os.path.exists('explicit-forward-order.cti'): os.remove('explicit-forward-order.cti') - ck2cti.convertMech('../data/explicit-forward-order.inp', - thermoFile='../data/dummy-thermo.dat', - outName='explicit-forward-order.cti', quiet=True) + convertMech('../data/explicit-forward-order.inp', + thermoFile='../data/dummy-thermo.dat', + outName='explicit-forward-order.cti', quiet=True) ref, gas = self.checkConversion('../data/explicit-forward-order.xml', 'explicit-forward-order.cti') self.checkKinetics(ref, gas, [300, 800, 1450, 2800], [5e3, 1e5, 2e6]) @@ -168,9 +174,9 @@ class chemkinConverterTest(utilities.CanteraTest): if os.path.exists('h2o2_transport_normal.cti'): os.remove('h2o2_transport_normal.cti') - ck2cti.convertMech('../../data/inputs/h2o2.inp', - transportFile='../../data/transport/gri30_tran.dat', - outName='h2o2_transport_normal.cti', quiet=True) + convertMech('../../data/inputs/h2o2.inp', + transportFile='../../data/transport/gri30_tran.dat', + outName='h2o2_transport_normal.cti', quiet=True) gas = ct.IdealGasMix('h2o2_transport_normal.cti') gas.set(X='H2:1.0, O2:1.0', T=300, P=101325) @@ -181,10 +187,10 @@ class chemkinConverterTest(utilities.CanteraTest): os.remove('h2o2_transport_missing_species.cti') def convert(): - ck2cti.convertMech('../../data/inputs/h2o2.inp', - transportFile='../data/h2o2-missing-species-tran.dat', - outName='h2o2_transport_missing_species.cti', - quiet=True) + convertMech('../../data/inputs/h2o2.inp', + transportFile='../data/h2o2-missing-species-tran.dat', + outName='h2o2_transport_missing_species.cti', + quiet=True) self.assertRaises(ck2cti.InputParseError, convert) @@ -193,20 +199,20 @@ class chemkinConverterTest(utilities.CanteraTest): os.remove('h2o2_transport_duplicate_species.cti') def convert(): - ck2cti.convertMech('../../data/inputs/h2o2.inp', - transportFile='../data/h2o2-duplicate-species-tran.dat', - outName='h2o2_transport_duplicate_species.cti', - quiet=True) + convertMech('../../data/inputs/h2o2.inp', + transportFile='../data/h2o2-duplicate-species-tran.dat', + outName='h2o2_transport_duplicate_species.cti', + quiet=True) # This should fail self.assertRaises(ck2cti.InputParseError, convert) # This should succeed - ck2cti.convertMech('../../data/inputs/h2o2.inp', - transportFile='../data/h2o2-duplicate-species-tran.dat', - outName='h2o2_transport_duplicate_species.cti', - quiet=True, - permissive=True) + convertMech('../../data/inputs/h2o2.inp', + transportFile='../data/h2o2-duplicate-species-tran.dat', + outName='h2o2_transport_duplicate_species.cti', + quiet=True, + permissive=True) def test_transport_bad_geometry(self): @@ -214,9 +220,9 @@ class chemkinConverterTest(utilities.CanteraTest): os.remove('h2o2_transport_bad_geometry.cti') def convert(): - ck2cti.convertMech('../../data/inputs/h2o2.inp', - transportFile='../data/h2o2-bad-geometry-tran.dat', - outName='h2o2_transport_bad_geometry.cti', - quiet=True) + convertMech('../../data/inputs/h2o2.inp', + transportFile='../data/h2o2-bad-geometry-tran.dat', + outName='h2o2_transport_bad_geometry.cti', + quiet=True) self.assertRaises(ck2cti.InputParseError, convert)