#!/usr/bin/env python # encoding: utf-8 ################################################################################ # # Copyright (c) 2009-2011 by the RMG Team (rmg_dev@mit.edu) # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the 'Software'), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. # ################################################################################ """ This module contains functions for converting Chemkin-format input files to Cantera input files (CTI). """ import logging import types import os.path import numpy as np import re UNIT_OPTIONS = {'CAL/': 'cal/mol', 'CAL/MOL': 'cal/mol', 'CAL/MOLE': 'cal/mol', 'EVOL': 'eV', 'EVOLTS': 'eV', 'JOUL': 'J/mol', 'JOULES/MOL': 'J/mol', 'JOULES/MOLE': 'J/mol', 'KCAL': 'kcal/mol', 'KCAL/MOL': 'kcal/mol', 'KCAL/MOLE': 'kcal/mol', 'KELV': 'K', 'KELVIN': 'K', 'KELVINS': 'K', 'KJOU': 'kJ/mol', 'KJOULES/MOL': 'kJ/mol', 'KJOULES/MOLE': 'kJ/mol', 'MOL': 'mol', 'MOLE': 'mol', 'MOLES': 'mol', 'MOLEC': 'molec', 'MOLECULES': 'molec'} class InputParseError(Exception): """ An exception class for exceptional behavior involving Chemkin-format mechanism files. Pass a string describing the circumstances that caused the exceptional behavior. """ pass class Species(object): def __init__(self, label): self.label = label self.thermo = None self.transport = None self.note = None self.composition = None def __str__(self): return self.label def to_cti(self, indent=0): lines = [] atoms = ' '.join('{0}:{1}'.format(*a) for a in self.composition.iteritems()) prefix = ' '*(indent+8) lines.append('species(name={0!r},'.format(self.label)) lines.append(prefix + 'atoms={0!r},'.format(atoms)) if self.thermo: lines.append(prefix + 'thermo={0},'.format(self.thermo.to_cti(15+indent))) if self.transport: lines.append(prefix + 'transport={0},'.format(self.transport.to_cti(14+indent))) if self.note: lines.append(prefix + 'note={0!r},'.format(self.note)) lines[-1] = lines[-1][:-1] + ')' lines.append('') return '\n'.join(lines) class ThermoModel(object): """ A base class for thermodynamics models, containing several attributes common to all models: =============== =================== ======================================== Attribute Type Description =============== =================== ======================================== `Tmin` ``float`` The minimum temperature at which the model is valid, or ``None`` if unknown or undefined `Tmax` ``float`` The maximum temperature at which the model is valid, or ``None`` if unknown or undefined `comment` ``str`` Information about the model (e.g. its source) =============== =================== ======================================== """ def __init__(self, Tmin=None, Tmax=None, comment=''): if Tmin is not None: self.Tmin = Tmin else: self.Tmin = None if Tmax is not None: self.Tmax = Tmax else: self.Tmax = None self.comment = comment class NASA(ThermoModel): """ A single NASA polynomial for thermodynamic data. The `coeffs` attribute stores the seven or nine polynomial coefficients :math:`\\mathbf{a} = \\left[a_{-2}\\ a_{-1}\\ a_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6 \\right]` from which the relevant thermodynamic parameters are evaluated via the expressions .. math:: \\frac{C_\\mathrm{p}(T)}{R} = a_{-2} T^{-2} + a_{-1} T^{-1} + a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4 .. math:: \\frac{H(T)}{RT} = - a_{-2} T^{-2} + a_{-1} T^{-1} \\ln T + a_0 + \\frac{1}{2} a_1 T + \\frac{1}{3} a_2 T^2 + \\frac{1}{4} a_3 T^3 + \\frac{1}{5} a_4 T^4 + \\frac{a_5}{T} .. math:: \\frac{S(T)}{R} = -\\frac{1}{2} a_{-2} T^{-2} - a_{-1} T^{-1} + a_0 \\ln T + a_1 T + \\frac{1}{2} a_2 T^2 + \\frac{1}{3} a_3 T^3 + \\frac{1}{4} a_4 T^4 + a_6 For the 7 coefficient form, the first two coefficients are taken to be zero. """ def __init__(self, coeffs, **kwargs): ThermoModel.__init__(self, **kwargs) if len(coeffs) not in (7,9): raise InputParseError('Invalid number of NASA polynomial coefficients; ' 'should be 7 or 9.') self.coeffs = coeffs def to_cti(self, indent=0): prefix = ' '*indent vals = ['{0: 15.8E}'.format(i) for i in self.coeffs] if len(self.coeffs) == 7: lines = ['NASA([{0:.2f}, {1:.2f}],'.format(self.Tmin[0], self.Tmax[0]), prefix+' [{0}, {1}, {2},'.format(*vals[0:3]), prefix+' {0}, {1}, {2},'.format(*vals[3:6]), prefix+' {0}]),'.format(vals[6])] else: lines = ['NASA9([{0:.2f}, {1:.2f}],'.format(self.Tmin[0], self.Tmax[0]), prefix+' [{0}, {1}, {2},'.format(*vals[0:3]), prefix+' {0}, {1}, {2},'.format(*vals[3:6]), prefix+' {0}, {1}, {2}]),'.format(*vals[6:9])] return '\n'.join(lines) class MultiNASA(ThermoModel): """ A set of thermodynamic parameters given by NASA polynomials. This class stores a list of :class:`NASA` objects in the `polynomials` attribute. When evaluating a thermodynamic quantity, a polynomial that contains the desired temperature within its valid range will be used. """ def __init__(self, polynomials=None, **kwargs): ThermoModel.__init__(self, **kwargs) self.polynomials = polynomials or [] def to_cti(self, indent=0): prefix = ' '*indent lines = [] for i,p in enumerate(self.polynomials): if i == 0: lines.append('({0}'.format(p.to_cti(indent+1))) elif i != len(self.polynomials)-1: lines.append(prefix + ' {0}'.format(p.to_cti(indent+1))) else: lines.append(prefix + ' {0})'.format(p.to_cti(indent+1)[:-1])) return '\n'.join(lines) class Reaction(object): """ A chemical reaction. The attributes are: =================== =========================== ============================ Attribute Type Description =================== =========================== ============================ `index` :class:`int` A unique nonnegative integer index `reactants` :class:`list` The reactant species (as :class:`Species` objects) `products` :class:`list` The product species (as :class:`Species` objects) `kinetics` :class:`KineticsModel` The kinetics model to use for the reaction `reversible` ``bool`` ``True`` if the reaction is reversible, ``False`` if not `duplicate` ``bool`` ``True`` if the reaction is known to be a duplicate, ``False`` if not `fwdOrders` ``dict`` Reaction order (value) for each specified species (key) =================== =========================== ============================ """ def __init__(self, index=-1, reactants=None, products=None, kinetics=None, reversible=True, duplicate=False, fwdOrders=None): self.index = index self.reactants = reactants # list of (stoichiometry, species) tuples self.products = products # list of (stoichiometry, specis) tuples self.kinetics = kinetics self.reversible = reversible self.duplicate = duplicate self.fwdOrders = fwdOrders if fwdOrders is not None else {} self.thirdBody = None def _coeff_string(self, coeffs): L = [] for stoichiometry,species in coeffs: if stoichiometry != 1: L.append('{0} {1}'.format(stoichiometry, species)) else: L.append(str(species)) expression = ' + '.join(L) if self.thirdBody: expression += ' (+ {0})'.format(self.thirdBody) return expression @property def reactantString(self): return self._coeff_string(self.reactants) @property def productString(self): return self._coeff_string(self.products) def __str__(self): """ Return a string representation of the reaction, in the form 'A + B <=> C + D'. """ arrow = ' <=> ' if self.reversible else ' -> ' return arrow.join([self.reactantString, self.productString]) def to_cti(self, indent=0): arrow = ' <=> ' if self.reversible else ' => ' kinstr = self.kinetics.to_cti(self.reactantString, arrow, self.productString, indent) k_indent = ' ' * (kinstr.find('(') + 1) if self.duplicate: kinstr = kinstr[:-1] + ",\n{0}options='duplicate')".format(k_indent) if self.fwdOrders: order = ' '.join('{0}:{1}'.format(k,v) for (k,v) in self.fwdOrders.items()) kinstr = kinstr[:-1] + ",\n{0}order='{1}')".format(k_indent, order) return kinstr class KineticsModel(object): """ A base class for kinetics models, containing several attributes common to all models: =============== =================== ======================================== Attribute Type Description =============== =================== ======================================== `Tmin` :class:`Quantity` The minimum absolute temperature in K at which the model is valid `Tmax` :class:`Quantity` The maximum absolute temperature in K at which the model is valid `Pmin` :class:`Quantity` The minimum absolute pressure in Pa at which the model is valid `Pmax` :class:`Quantity` The maximum absolute pressure in Pa at which the model is valid `comment` :class:`str` A string containing information about the model (e.g. its source) =============== =================== ======================================== """ def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''): self.Tmin = Tmin self.Tmax = Tmax self.Pmin = Pmin self.Pmax = Pmax self.comment = comment def isPressureDependent(self): """ Return ``True`` if the kinetics are pressure-dependent or ``False`` if they are pressure-independent. This method must be overloaded in the derived class. """ raise InputParseError('Unexpected call to KineticsModel.isPressureDependent();' ' you should be using a class derived from KineticsModel.') def to_cti(self, reactantstr, arrow, productstr): raise InputParseError('to_cti is not implemented for objects of class {0}'.format(self.__class__.__name__)) def efficiencyString(self): if hasattr(self, 'efficiencies'): return ' '.join('{0}:{1}'.format(mol, eff) for mol,eff in self.efficiencies.iteritems()) else: return '' class KineticsData(KineticsModel): """ A kinetics model based around a set of discrete (high-pressure limit) rate coefficients at various temperatures. The attributes are: =========== =================== ============================================ Attribute Type Description =========== =================== ============================================ `Tdata` :class:`Quantity` The temperatures at which the heat capacity data is provided `kdata` :class:`Quantity` The rate coefficients in SI units at each temperature in `Tdata` =========== =================== ============================================ """ def __init__(self, Tdata=None, kdata=None, **kwargs): KineticsModel.__init__(self, **kwargs) self.Tdata = Tdata self.kdata = kdata def isPressureDependent(self): """ Returns ``False`` since KineticsData kinetics are not pressure-dependent. """ return False class Arrhenius(KineticsModel): """ Represent a set of modified Arrhenius kinetics. The kinetic expression has the form .. math:: k(T) = A \\left( \\frac{T}{T_0} \\right)^n \\exp \\left( - \\frac{E_\\mathrm{a}}{RT} \\right) where :math:`A`, :math:`n`, :math:`E_\\mathrm{a}`, and :math:`T_0` are the parameters to be set, :math:`T` is absolute temperature, and :math:`R` is the gas law constant. The attributes are: =============== =================== ======================================== Attribute Type Description =============== =================== ======================================== `A` :class:`Quantity` The preexponential factor in s^-1, m^3/mol*s, etc. `T0` :class:`Quantity` The reference temperature in K `n` :class:`Quantity` The temperature exponent `Ea` :class:`Quantity` The activation energy in J/mol =============== =================== ======================================== """ def __init__(self, A=0.0, n=0.0, Ea=0.0, T0=1.0, **kwargs): KineticsModel.__init__(self, **kwargs) self.A = A self.T0 = T0 self.n = n self.Ea = Ea def isPressureDependent(self): """ Returns ``False`` since Arrhenius kinetics are not pressure-dependent. """ return False def rateStr(self): return '[{0.A[0]:e}, {0.n}, {0.Ea[0]}]'.format(self) def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstring = reactantstr + arrow + productstr return 'reaction({0!r}, {1})'.format(rxnstring, self.rateStr()) class PDepArrhenius(KineticsModel): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the expression .. math:: k(T,P) = A(P) T^{n(P)} \\exp \\left[ \\frac{-E_\\mathrm{a}(P)}{RT} \\right] where the modified Arrhenius parameters are stored at a variety of pressures and interpolated between on a logarithmic scale. The attributes are: =============== ================== ============================================ Attribute Type Description =============== ================== ============================================ `pressures` :class:`list` The list of pressures in Pa `arrhenius` :class:`list` The list of :class:`Arrhenius` objects at each pressure `highPlimit` :class:`Arrhenius` The high (infinite) pressure limiting :class:`Arrhenius` expression =============== ================== ============================================ Note that `highPlimit` is not used in evaluating k(T,P). """ def __init__(self, pressures=None, arrhenius=None, highPlimit=None, **kwargs): KineticsModel.__init__(self, **kwargs) self.pressures = pressures self.arrhenius = arrhenius or [] self.highPlimit = highPlimit or None def isPressureDependent(self): """ Returns ``True`` since PDepArrhenius kinetics are pressure-dependent. """ return True def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstring = reactantstr + arrow + productstr lines = ['pdep_arrhenius({0!r},'.format(rxnstring)] prefix = ' '*(indent+15) template = '[({0}, {1!r}), {2.A[0]:e}, {2.n}, {2.Ea[0]}],' for pressure,arrhenius in zip(self.pressures[0], self.arrhenius): lines.append(prefix + template.format(pressure, self.pressures[1], arrhenius)) lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) class Chebyshev(KineticsModel): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the expression .. math:: \\log k(T,P) = \\sum_{t=1}^{N_T} \\sum_{p=1}^{N_P} \\alpha_{tp} \\phi_t(\\tilde{T}) \\phi_p(\\tilde{P}) where :math:`\\alpha_{tp}` is a constant, :math:`\\phi_n(x)` is the Chebyshev polynomial of degree :math:`n` evaluated at :math:`x`, and .. math:: \\tilde{T} \\equiv \\frac{2T^{-1} - T_\\mathrm{min}^{-1} - T_\\mathrm{max}^{-1}}{T_\\mathrm{max}^{-1} - T_\\mathrm{min}^{-1}} .. math:: \\tilde{P} \\equiv \\frac{2 \\log P - \\log P_\\mathrm{min} - \\log P_\\mathrm{max}}{\\log P_\\mathrm{max} - \\log P_\\mathrm{min}} are reduced temperature and reduced pressures designed to map the ranges :math:`(T_\\mathrm{min}, T_\\mathrm{max})` and :math:`(P_\\mathrm{min}, P_\\mathrm{max})` to :math:`(-1, 1)`. The attributes are: =============== =============== ============================================ Attribute Type Description =============== =============== ============================================ `coeffs` :class:`list` Matrix of Chebyshev coefficients `kunits` ``str`` The units of the generated k(T, P) values `degreeT` :class:`int` The number of terms in the inverse temperature direction `degreeP` :class:`int` The number of terms in the log pressure direction =============== =============== ============================================ """ def __init__(self, coeffs=None, kunits='', **kwargs): KineticsModel.__init__(self, **kwargs) if coeffs is not None: self.coeffs = np.array(coeffs, np.float64) self.degreeT = self.coeffs.shape[0] self.degreeP = self.coeffs.shape[1] else: self.coeffs = None self.degreeT = 0 self.degreeP = 0 self.kunits = kunits def isPressureDependent(self): """ Returns ``True`` since Chebyshev polynomial kinetics are pressure-dependent. """ return True def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstr = reactantstr + arrow + productstr prefix = ' '*(indent+19) lines = ['chebyshev_reaction({0!r},'.format(rxnstr), prefix + 'Tmin={0.Tmin}, Tmax={0.Tmax},'.format(self), prefix + 'Pmin={0.Pmin}, Pmax={0.Pmax},'.format(self)] for i in range(self.degreeT): coeffline = ', '.join('{0: 12.5e}'.format(self.coeffs[i,j]) for j in range(self.degreeP)) if i == 0: lines.append(prefix + 'coeffs=[[{0}],'.format(coeffline)) else: lines.append(prefix + ' [{0}],'.format(coeffline)) lines[-1] = lines[-1][:-1] + '])' return '\n'.join(lines) class ThirdBody(KineticsModel): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the expression .. math:: k(T,P) = k(T) [\\ce{M}] where :math:`k(T)` is an Arrhenius expression and :math:`[\\ce{M}] \\approx P/RT` is the concentration of the third body (i.e. the bath gas). A collision efficiency can be used to further correct the value of :math:`k(T,P)`. The attributes are: =============== ======================= ==================================== Attribute Type Description =============== ======================= ==================================== `arrheniusHigh` :class:`Arrhenius` The Arrhenius kinetics `efficiencies` ``dict`` A mapping of species to collider efficiencies =============== ======================= ==================================== """ def __init__(self, arrheniusHigh=None, efficiencies=None, **kwargs): KineticsModel.__init__(self, **kwargs) self.arrheniusHigh = arrheniusHigh self.efficiencies = {} if efficiencies is not None: for mol, eff in efficiencies.iteritems(): self.efficiencies[mol] = eff def isPressureDependent(self): """ Returns ``True`` since third-body kinetics are pressure-dependent. """ return True def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstr = reactantstr + ' + M' + arrow + productstr + ' + M' prefix = ' '*(indent + 20) lines = ['three_body_reaction({0!r}, {1},'.format(rxnstr, self.arrheniusHigh.rateStr())] if self.efficiencies: lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString())) lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) class Lindemann(ThirdBody): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the expression .. math:: k(T,P) = k_\\infty(T) \\left[ \\frac{P_\\mathrm{r}}{1 + P_\\mathrm{r}} \\right] F where .. math:: P_\\mathrm{r} &= \\frac{k_0(T)}{k_\\infty(T)} [\\ce{M}] k_0(T) &= A_0 T^{n_0} \\exp \\left( - \\frac{E_0}{RT} \\right) k_\\infty(T) &= A_\\infty T^{n_\\infty} \\exp \\left( - \\frac{E_\\infty}{RT} \\right) and :math:`[\\ce{M}] \\approx P/RT` is the concentration of the bath gas. The Arrhenius expressions :math:`k_0(T)` and :math:`k_\\infty(T)` represent the low-pressure and high-pressure limit kinetics, respectively. The former is necessarily one reaction order higher than the latter. For the Lindemann model, :math:`F = 1`. A collision efficiency can be used to further correct the value of :math:`k(T,P)`. The attributes are: =============== ======================= ==================================== Attribute Type Description =============== ======================= ==================================== `arrheniusLow` :class:`Arrhenius` The Arrhenius kinetics at the low-pressure limit `arrheniusHigh` :class:`Arrhenius` The Arrhenius kinetics at the high-pressure limit `efficiencies` ``dict`` A mapping of species to collider efficiencies =============== ======================= ==================================== """ def __init__(self, arrheniusLow=None, **kwargs): ThirdBody.__init__(self, **kwargs) self.arrheniusLow = arrheniusLow def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstr = reactantstr + arrow + productstr prefix = ' '*(indent + 17) lines = ['falloff_reaction({0!r},'.format(rxnstr)] lines.append(prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr())) lines.append(prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr())) if self.efficiencies: lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString())) lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) class Troe(Lindemann): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the expression .. math:: k(T,P) = k_\\infty(T) \\left[ \\frac{P_\\mathrm{r}}{1 + P_\\mathrm{r}} \\right] F where .. math:: P_\\mathrm{r} &= \\frac{k_0(T)}{k_\\infty(T)} [\\ce{M}] k_0(T) &= A_0 T^{n_0} \\exp \\left( - \\frac{E_0}{RT} \\right) k_\\infty(T) &= A_\\infty T^{n_\\infty} \\exp \\left( - \\frac{E_\\infty}{RT} \\right) and :math:`[\\ce{M}] \\approx P/RT` is the concentration of the bath gas. The Arrhenius expressions :math:`k_0(T)` and :math:`k_\\infty(T)` represent the low-pressure and high-pressure limit kinetics, respectively. The former is necessarily one reaction order higher than the latter. A collision efficiency can be used to further correct the value of :math:`k(T,P)`. For the Troe model the parameter :math:`F` is computed via .. math:: \\log F &= \\left\\{1 + \\left[ \\frac{\\log P_\\mathrm{r} + c}{n - d (\\log P_\\mathrm{r} + c)} \\right]^2 \\right\\}^{-1} \\log F_\\mathrm{cent} c &= -0.4 - 0.67 \\log F_\\mathrm{cent} n &= 0.75 - 1.27 \\log F_\\mathrm{cent} d &= 0.14 F_\\mathrm{cent} &= (1 - \\alpha) \\exp \\left( -T/T_3 \\right) + \\alpha \\exp \\left( -T/T_1 \\right) + \\exp \\left( -T_2/T \\right) The attributes are: =============== ======================= ==================================== Attribute Type Description =============== ======================= ==================================== `arrheniusLow` :class:`Arrhenius` The Arrhenius kinetics at the low-pressure limit `arrheniusHigh` :class:`Arrhenius` The Arrhenius kinetics at the high-pressure limit `efficiencies` ``dict`` A mapping of species to collider efficiencies `alpha` :class:`Quantity` The :math:`\\alpha` parameter `T1` :class:`Quantity` The :math:`T_1` parameter `T2` :class:`Quantity` The :math:`T_2` parameter `T3` :class:`Quantity` The :math:`T_3` parameter =============== ======================= ==================================== """ def __init__(self, alpha=0.0, T3=0.0, T1=0.0, T2=None, **kwargs): Lindemann.__init__(self, **kwargs) self.alpha = alpha self.T3 = T3 self.T1 = T1 self.T2 = T2 def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstr = reactantstr + arrow + productstr prefix = ' '*17 lines = ['falloff_reaction({0!r},'.format(rxnstr), prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()), prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr())] if self.T2: troeArgs = 'A={0.alpha[0]}, T3={0.T3[0]}, T1={0.T1[0]}, T2={0.T2[0]}'.format(self) else: troeArgs = 'A={0.alpha[0]}, T3={0.T3[0]}, T1={0.T1[0]}'.format(self) lines.append(prefix + 'falloff=Troe({0}),'.format(troeArgs)) if self.efficiencies: lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString())) # replace trailing comma lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) class Sri(Lindemann): """ A kinetic model of a phenomenological rate coefficient k(T, P) using the "SRI" formulation of the blending function :math:`F` using either 3 or 5 parameters. See :ref:`sec-sri-falloff`. The attributes are: =============== ======================= ==================================== Attribute Type Description =============== ======================= ==================================== `arrheniusLow` :class:`Arrhenius` The Arrhenius kinetics at the low-pressure limit `arrheniusHigh` :class:`Arrhenius` The Arrhenius kinetics at the high-pressure limit `efficiencies` ``dict`` A mapping of species to collider efficiencies `A` ``float`` The :math:`a` parameter `B` ``float`` The :math:`b` parameter `C` ``float`` The :math:`c` parameter `D` ``float`` The :math:`d` parameter `E` ``float`` The :math:`e` parameter =============== ======================= ==================================== """ def __init__(self, A=0.0, B=0.0, C=0.0, D=1.0, E=0.0, **kwargs): Lindemann.__init__(self, **kwargs) self.A = A self.B = B self.C = C self.D = D self.E = E def to_cti(self, reactantstr, arrow, productstr, indent=0): rxnstr = reactantstr + arrow + productstr prefix = ' '*17 lines = ['falloff_reaction({0!r},'.format(rxnstr), prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()), prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr())] if self.D == 1.0 and self.E == 0.0: sriArgs = 'A={0.A}, B={0.B}, C={0.C}'.format(self) else: sriArgs = 'A={0.A}, B={0.B}, C={0.C}, D={0.D}, E={0.E}'.format(self) lines.append(prefix + 'falloff=SRI({0}),'.format(sriArgs)) if self.efficiencies: lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString())) # replace trailing comma lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) class TransportData(object): geometryFlags = ['atom', 'linear', 'nonlinear'] def __init__(self, label, geometry, wellDepth, collisionDiameter, dipoleMoment, polarizability, zRot, comment=None): if not isinstance(label, types.StringTypes): raise InputParseError("Bad label for transport data: " + repr(label)) if int(geometry) not in (0,1,2): raise InputParseError("Bad geometry flag '{0}' for species '{1}'".format(geometry, label)) self.label = label self.geometry = self.geometryFlags[int(geometry)] self.wellDepth = float(wellDepth) self.collisionDiameter = float(collisionDiameter) self.dipoleMoment = float(dipoleMoment) self.polarizability = float(polarizability) self.zRot = float(zRot) self.comment = comment or '' # @todo: include this in the CTI def __repr__(self): return ('TransportData({label!r}, {geometry!r}, {wellDepth!r}, ' '{collisionDiameter!r}, {dipoleMoment!r}, {polarizability!r}, ' '{zRot!r}, {comment!r})').format(**self.__dict__) def to_cti(self, indent=0): prefix = ' '*(indent+18) lines = ['gas_transport(geom={0!r},'.format(self.geometry), prefix+'diam={0},'.format(self.collisionDiameter), prefix+'well_depth={0},'.format(self.wellDepth)] if self.dipoleMoment: lines.append(prefix+'dipole={0},'.format(self.dipoleMoment)) if self.polarizability: lines.append(prefix+'polar={0},'.format(self.polarizability)) if self.zRot: lines.append(prefix+'rot_relax={0},'.format(self.zRot)) lines[-1] = lines[-1][:-1] + ')' return '\n'.join(lines) def fortFloat(s): """ Convert a string representation of a floating point value to a float, allowing for some of the peculiarities of allowable Fortran representations. """ s = s.strip() s = s.replace('D', 'E').replace('d', 'e') s = s.replace('E ', 'E+').replace('e ', 'e+') return float(s) 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 warn(self, message): if self.warning_as_error: raise InputParseError(message) else: logging.warning(message) 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() 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 = 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 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 self.speciesDict: raise InputParseError('Unexpected species "{0}" in reaction expression "{1}".'.format(species, expression)) else: dest.append((stoichiometry, self.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, 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 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]) 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)) 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 #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 = '' line = f.readline() while line != '' and 'END' not in line: lineStartsWithComment = line.startswith('!') line, comment = removeCommentFromLine(line) line = line.strip() comment = comment.strip() 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) 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] == '': kineticsList.pop(0) if len(kineticsList) != len(commentsList): commentsList = ['' for kinetics in kineticsList] 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) 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)) index = 0 for reaction in self.reactions: index += 1 reaction.index = index if transportLines: self.parseTransportData(transportLines) 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. """ 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] 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)) def writeCTI(self, header=None, name='gas', transportModel='Mix', outName='mech.cti'): 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)) 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) speciesNames = '\n'.join(speciesNames).strip() lines = [] if header: lines.extend(header) # 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('') # Write the individual species data lines.append(delimiterLine) lines.append('# Species data') lines.append(delimiterLine) lines.append('') for s in self.speciesList: lines.append(s.to_cti()) # Write the reactions lines.append(delimiterLine) lines.append('# Reaction data') lines.append(delimiterLine) for i,r in enumerate(self.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(self): print """ ck2cti.py: Convert Chemkin-format mechanisms to Cantera input files (.cti) Usage: ck2cti --input= [--thermo=] [--transport=] [--id=] [--output=] [--permissive] [-d | --debug] Example: ck2cti --input=chem.inp --thermo=therm.dat --transport=tran.dat If the output file name is not given, an output file with the same name as the input file, with the extension changed to '.cti'. The '--permissive' option allows certain recoverable parsing errors (e.g. 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) if permissive is not None: self.warning_as_error = not permissive # Read input mechanism files self.loadChemkinFile(inputFile) if thermoFile: self.loadChemkinFile(thermoFile) if transportFile: lines = open(transportFile).readlines() self.parseTransportData(lines) # 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)) if not outName: outName = os.path.splitext(inputFile)[0] + '.cti' # 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__': import getopt import sys longOptions = ['input=', 'thermo=', 'transport=', 'id=', 'output=', 'permissive', 'help', 'debug'] try: optlist, args = getopt.getopt(sys.argv[1:], 'dh', longOptions) options = dict() for o,a in optlist: options[o] = a if args: raise getopt.GetoptError('Unexpected command line option: ' + repr(' '.join(args))) except getopt.GetoptError as e: print 'ck2cti.py: Error parsing arguments:' print e 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: parser.showHelp() sys.exit(0) if '--input' in options: inputFile = options['--input'] else: print 'Error: no mechanism input file specified' sys.exit(1) if '--output' in options: outName = options['--output'] if not outName.endswith('.cti'): outName += '.cti' else: outName = None permissive = '--permissive' in options thermoFile = options.get('--thermo') transportFile = options.get('--transport') phaseName = options.get('--id', 'gas') parser.convertMech(inputFile, thermoFile, transportFile, phaseName, outName, permissive=permissive)