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

*** empty log message ***

Browse source
This commit is contained in:
Dave Goodwin 2003-09-09 16:22:26 +00:00
parent 3132e39e1f
commit 754d098884
10 changed files with 0 additions and 2520 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

726
Cantera/python/Cantera/Flow.py
View file

@ -1,726 +0,0 @@
"""
One-dimensional reacting flows.
"""
from Cantera import getCanteraError
from FlowPlotter import FlowPlotter
from exceptions import *
import refine
import sys, types, copy, tempfile
import _cantera
import interp
import math
from Numeric import array, zeros, ones, transpose, size, sort, shape, asarray
_flows = {'Stagnation':0, 'Stag':0,
'OneDimensional':1, '1D':1, 'OneD':1, 'OneDim':1,
'Free':2}
_geom = {'Axisymmetric':0, 'Axi':0, 'Planar':1}
class Flow1D:
""" One-dimensional reacting flows.
Class Flow1D simulates a one-dimensional flow domain. To use
Flow1D objects, they must be installed in a container, which is an
object of class OneDim. Each Flow1D domain must be terminated by
boundary domains.
Class Flow1D can model several types of steady 'one dimensional'
reacting flows. The flows are one-dimensional in the sense that
the governing equations for the steady-state solution can be cast
in the form of a set of ordinary differential equations in one
axial coordinate (z). For the case of stagnation flows, this
results from a similarity transformation that reduces the
physically two-dimensional problem to one that is mathematically
one-dimensional.
The types of flows that may be simulated are:
- One-dimensional reacting flows, such as burner-stabilized
premixed flames;
- Axisymmetric stagnation-point flows
- Planar stagnation-point flows
"""
# Allowed option keywords are defined here so that only these
# keywords may be added to the _opt dictionary.
_timeint_options = ['ftime', 'min_timestep', 'max_timestep',
'nsteps', 'timestep']
_newton_options = ['max_jac_age', 'rtol', 'atol']
_output_options = ['loglevel', 'plotfile']
_options = _newton_options + _timeint_options + _output_options
def __init__(self,
flow_type = 'Stagnation',
flow_geom = 'Axisymmetric',
gas = None,
grid = None,
pressure = 1.01325e5):
"""Flow1D Constructor.
"""
self.__flow_id = -1
self.domainType = 0
self.loglevel = 1
self.initial = {}
if not grid:
raise CanteraError('Grid not specified!')
if not gas:
raise CanteraError('Gas mixture object not specified!')
self.gas = gas
self.nsp = self.gas.nSpecies()
self.nv = self.nsp + 4
if _flows.has_key(flow_type):
self.type = _flows[flow_type]
else:
raise CanteraError('unsupported flow type: '+flow_type)
if self.type == 0:
if flow_geom == 'Planar': self.type = 3
# Create the kernel object. This is an instance of a subclass of
# C++ Cantera class 'StFlow'.
self.__flow_id = _cantera.Flow(self.type, self.gas.phase_id(), len(grid))
# Set the grid. Must be done _after_ creating the kernel
# object, since it sets the grid there too.
self.setGrid(asarray(grid))
# Set the pressure. Note that the pressure is constant
# throughout the flowfield, due to the assumption of low Mach
# number.
self.setPressure(pressure)
# set the thermo, kinetics, and transport managers to those of
# object self.gas.
self.setThermo(self.gas)
self.setKinetics(self.gas)
self.setTransport(self.gas)
# create a NumPy array to hold the solution. Since NumPy
# arrays are stored by row, while Cantera expects arrays
# stored by column, this array is defined with the grid point
# number as the first index. In this way, all solution
# variables for one grid point will be stored in contiguous
# locations.
self.x = zeros((self.npts, self.nsp + 4), 'd')
self.xnew = zeros((self.npts, self.nsp + 4), 'd')
# Default error tolerances
self.setTolerances(
u = (1.e-8, 1.e-15),
V = (1.e-8, 1.e-15),
T = (1.e-8, 1.e-15),
L = (1.e-8, 1.e-15),
Y = (1.e-8, 1.e-15))
self.energy = 0
self.names = ['u','V','T','L']+ list(self.gas.speciesNames())
# finish setting default parameter values
self.restoreDefaults()
self.time = 0.0
# For refiner, only refine the grid based on species mass
# fraction and velocity profiles until the energy equation is
# enabled.
self.refine_components = range(4,4+self.nsp)
self.refiner = refine.Refiner(components = self.refine_components,
names = self.names)
# Create an object to handle plotting results.
self.plotter = FlowPlotter(self)
#================> end of method '__init__' <===================
def __del__(self):
"""Delete the Flow1D instance."""
if self.__flow_id >= 0:
_cantera.flow_delete(self.__flow_id)
def shape(self):
"""Return (rows, columns) of solution matrix."""
return (len(self.z), self.gas.nSpecies() + 4)
def nPoints(self):
"""Number of grid points."""
return self.npts
def flow_id(self):
"""ID used to accesss the kernel object."""
return self.__flow_id
def option(self, key):
"""Return the value of an option."""
return self._opt[key]
def setGrid(self, z):
"""Set the grid to the values in sequence z.
The values will be sorted, and so the input sequence
does not need to be monotonic.
"""
self.z = array(z)
sort(self.z)
self.npts = len(z)
return _cantera.flow_setupgrid(self.flow_id(), z)
def setThermo(self, th):
"""Set the thermodynamic property manager."""
id = _cantera.flow_setthermo(self.flow_id(), th.phase_id())
return id
def setKinetics(self, kin):
"""Set the kinetics manager."""
id = _cantera.flow_setkinetics(self.__flow_id, kin.ckin)
return id
def setTransport(self, tr, soret=0):
"""Set the transport manager."""
id = _cantera.flow_settransport(self.__flow_id,
tr.transport_id(), soret)
return id
def setPressure(self, p):
"""Set the pressure [Pa].
Since the flow has very nearly the same pressure everywhere,
this pressure value is used in all computations involving the
equation of state.
"""
self.p = p
_cantera.flow_setpressure(self.__flow_id, p)
def holdTemperature(self, points, t0):
"""Hold the temperature at grid points 'points' to t0,
and disable the energy equation.
"""
self.x[points,2] = t0
_cantera.flow_settemperature(self.__flow_id, points, t0)
self.setEnergyEqn('off')
def holdMassFraction(self, j, k, y0):
self.x[j,4+k] = y0
_cantera.flow_setmassfraction(self.__flow_id, j, k, y0)
def setInitialProfiles(self, datatable=None):
"""Set initial velocity, temperature, and/or species profiles.
datatable -- Dictionary mapping variable names to sequences of
(position, value) pairs. The position is specified in
relative terms, as a number in the range [0,1], where the
value zero corresponds to the smallest grid value, and the
value one to the largest.
The keys of datatable must be 'u', 'V', 'T', or a species
name. Velocity and temperature values are entered in SI units,
and species values are entered in arbitrary molar units, and
will be normalized to produce mole fractions.
The profile will be linearly interpolated onto the grid from
the data provided. Each variable may be specified at different
locations.
Example:
data = {}
data['T'] = [(0, 500), (0.3, 2000), (0.8, 2200), (1, 1500)]
data['u'] = [(0, 0.0), (1, 2)]
data['H2'] = [(0, 0.2), (1, 0.3)]
data['O2'] = [(0, 0.8), (1, 0.7)]
flow.setInitialProfiles(data)
"""
if datatable:
self.datatable = datatable
else:
datatable = self.datatable
vars = datatable.keys()
x = zeros((self.npts, self.nsp),'d')
equil = 0
for var in vars:
data = datatable[var]
if not var == 'equil':
data.sort()
zz = []
v = []
for item in data:
zz.append(self.z[0] + item[0]*(self.z[-1] - self.z[0]))
v.append(item[1])
self.initial[var] = (zz, v)
if (var == 'u'):
for j in range(self.npts):
self.x[j,0] = interp.interp(self.z[j],zz,v)
elif (var == 'V'):
for j in range(self.npts):
self.x[j,1] = interp.interp(self.z[j],zz,v)
elif (var == 'T'):
for j in range(self.npts):
self.holdTemperature(j,interp.interp(self.z[j],zz,v))
else:
k = self.gas.speciesIndex(var)
if k < 0:
raise CanteraError('Unknown species name: '+var)
for j in range(self.npts):
x[j,k] = interp.interp(self.z[j],zz,v)
else:
equil = 1
if equil == 1:
xin = self.left.X
x[0,:] = xin
for j in range(1,self.npts):
self.gas.setState_TPX(self.x[j,2],self.p,xin)
try:
self.gas.equilibrate('TP')
except:
pass
x[j,:] = self.gas.moleFractions()
# convert input mole fractions to mass fractions,
# and set the mass fraction profiles
for j in range(self.npts):
self.gas.setMoleFractions(x[j,:])
y = self.gas.massFractions()
for k in range(self.nsp):
self.x[j, k+4] = y[k]
self.holdMassFraction(j,k,y[k])
self.enableSpecies()
def regrid(self, grid):
oldx = self.x
oldgrid = self.z
np, nv = shape(oldx)
v = range(nv)
self.setGrid(grid)
for j in range(self.npts):
for n in v:
self.x[j,n] = interp.interp(self.z[j], oldgrid, oldx)
def __repr__(self):
return self.show()
def show(self, x = None):
fname = tempfile.mktemp('.dat')
x = self.x
_cantera.flow_showsolution(self.__flow_id, fname, x)
f = open(fname,'r')
y = f.readlines()
print
for line in y:
print line,
print
f.close()
def showResid(self):
"""Print the current residual values"""
print transpose(self.resid())
def T(self,j):
"""Temperature at grid point j [K]."""
return self.x[j,2]
def u(self,j):
"""Axial velocity at grid point j [m/s]."""
return self.x[j, 0]
def V(self,j):
"""Radial velocity divided by radius at grid point j [1/s]."""
return self.x[j, 1]
def lamb(self,j):
"""(1/r)(dP/dr) at grid point j.
If the solution has converged, this will be the same at all
grid points.
"""
return self.x[j, 3]
def massFraction(self,sp,j):
"""Mass fraction of species 'sp', which may be referenced by
name or by index number."""
k = self.gas.speciesIndex(sp)
return self.x[j, k + 4]
def density(self, j):
"""Density [kg/m^3]."""
self.setGas(j)
return self.gas.density()
def molWt(self, j):
"""Mean molecular weight [kg/kmol]."""
self.setGas(j)
return self.gas.meanMolecularWeight()
def setGas(self, j):
"""Set the state of the internal gas mixture object to be
consistent with the solution at grid point j."""
self.gas.setTemperature(self.T(j))
y = self.x[j,4:]
self.gas.setMassFractions(y)
self.gas.setPressure(self.p)
def setTolerances(self, u = None, V = None, T = None,
L = None, Y = None):
"""Set error tolerances.
The inputs are tuples of (relative, absolute) tolerances for
each of u, V, T, L, and Y.
"""
default = (1.e-7, 1.e-15)
if u == None: u = default
if V == None: V = default
if T == None: T = default
if L == None: L = default
if Y == None: Y = default
rtol = zeros(self.nsp+4,'d')
atol = zeros(self.nsp+4,'d')
rtol[0] = u[0]
rtol[1] = V[0]
rtol[2] = T[0]
rtol[3] = L[0]
for k in range(self.nsp):
rtol[4+k] = Y[0]
atol[0] = u[1]
atol[1] = V[1]
atol[2] = T[1]
atol[3] = L[1]
for k in range(self.nsp):
atol[4+k] = Y[1]
_cantera.flow_settolerances(self.__flow_id, len(rtol), rtol,
len(atol), atol)
def disableSpecies(self):
off = zeros(self.nsp,'d')
_cantera.flow_solvespecies(self.__flow_id, self.nsp, off)
def enableSpecies(self):
o = ones(self.nsp,'d')
_cantera.flow_solvespecies(self.__flow_id, self.nsp, o)
def setEnergyEqn(self, o, loglevel = 0, pt = -1):
"""Enable or disable the energy equation."""
if o == 'on':
self.energy = 1
_cantera.flow_energy(self.__flow_id, pt, 1)
if loglevel > 0: print '\n%%%%%%%%%%%% Enabling energy equation %%%%%%%%%%%%\n'
elif o == 'off':
self.energy = 0
_cantera.flow_energy(self.__flow_id, pt, 0)
if loglevel > 0: print '\n%%%%%%%%%%%% Disabling energy equation %%%%%%%%%%%%\n'
def setEnergyFactor(self, e):
_cantera.flow_setenergyfactor(self.__flow_id, e)
## def resid(self, point = -1):
## """Return the residual vector.
## If 'point' is specified, the residual is only evaluated at
## this point and those adjacent to it. Otherwise, it is
## evaluated at all grid points.
## """
## r = zeros(shape(self.x),'d')
## if point >= 0 and point < self.npts:
## j = point
## else:
## j = -1
## _cantera.flow_eval(self.__flow_id, j, self.x, r)
## return r
def ssnorm(self):
"""Absolute maximum value of the steady-state residual for any
component at any point."""
a = self.container.ssnorm(self.x,self.xnew)
return a
def integrateChem(self, dt, loglevel=1):
"""Condition the species profiles by integrating the
constant-pressure kinetics rate equations
\[
\dot Y_k = \dot\omega_k M_k / \rho.
\]
at each grid point for time 'dt'.
"""
if loglevel > 0:
print '\nIntegrating the chemical ource terms for %10.4g s...' % dt,
_cantera.flow_integratechem(self.__flow_id, self.x, dt)
print _cantera.readlog()
def restoreDefaults(self):
"""Restore default options."""
self._opt = {}
self.setOptions(
max_jac_age = 5,
timestep = 1.e-6,
min_timestep = 1.e-12,
max_timestep = 0.1,
nsteps = 20,
ftime = 3.0,
plotfile = ""
)
def setOptions(self, **options):
"""Set options."""
for kw in options.keys():
if kw in Flow1D._options:
self._opt[kw] = options[kw]
else:
raise OptionError(kw)
def refine(self, loglevel = 2):
"""Refine the grid.
"""
r = self.refiner
r.components = range(4,self.nsp+4)
if self.energy:
r.components.append(2)
if self.type == 0:
r.components += [0,1]
#dsave = r.delta
#while 1 > 0:
# try:
znew, zadded, xn, ok = r.refine(grid = self.z, solution = self.x)
nin = len(znew) - len(self.z)
# break
# except:
# r.delta = (0.9*r.delta[0], 0.9*r.delta[1])
#r.delta = dsave
if not ok:
self.setGrid(znew)
self.x = array(xn,'d')
self.xnew = zeros(shape(xn),'d')
# update the fixed temperature values if the energy
# equation is not being solved
if self.energy == 0:
for j in range(self.npts):
zz, tt = self.initial['T']
t = interp.interp(self.z[j],zz,tt)
self.holdTemperature(j,t)
else:
for j in range(self.npts):
self.holdTemperature(j,self.x[j,2])
self.setEnergyEqn('on')
if loglevel > 0:
print 'Refine: ',
print 'added',nin,'points.'
print 'Grid size = ',len(self.z)
return nin
def prune(self, loglevel = 2):
"""Prune the grid.
"""
r = self.refiner
r.components = range(4,self.nsp+4)
if self.energy:
r.components.append(2)
znew, xn = r.prune(grid = self.z, solution = self.x)
nout = len(self.z) - len(znew)
if nout > 0:
self.setGrid(znew)
self.x = array(xn,'d')
self.xnew = zeros(shape(xn),'d')
# update the fixed temperature values if the energy
# equation is not being solved
if self.energy == 0:
for j in range(self.npts):
zz, tt = self.initial['T']
t = interp.interp(self.z[j],zz,tt)
self.holdTemperature(j,t)
else:
for j in range(self.npts):
self.holdTemperature(j,self.x[j,2])
self.setEnergyEqn('on')
if loglevel > 0:
print 'Prune: ',
print 'removed',nout,'points.'
print 'Grid size = ',len(self.z)
return nout
def save(self, filename, id, desc="", append=0):
"""Save a solution to a file.
filename -- file name. If the file, does not exist, it will
be created. The save files are xml files, and the
filename should have the extension '.xml'. If it
does not, this extension will be appended to the
name.
id -- the ID tag of the solution. Multiple solutions may
be saved to the same file. Specifying a unique ID
tag allows this solution to selected later by
method 'restore'.
append -- If append > 0, the solution will be appended to the
file. Otherwise, the file will be overwritten if it
exists.
"""
appnd = append
fn = filename
extn = filename[-4:]
if extn <> '.xml' and extn <> '.XML':
fn = filename + '.xml'
#if self.loglevel > 0:
#print 'Solution saved to file',filename,'as solution',`id`
#_cantera.flow_save(self.__flow_id, fn, id, appnd, self.x)
_cantera.flow_save(self.__flow_id, fn, id, desc, self.x)
print _cantera.readlog()
np, nv = shape(self.x)
f = open('ctsoln.dat','w')
for j in range(np):
f.write('%14.6e %14.6e ' % (100.0*self.z[j], self.x[j,2]))
for k in range(4,nv):
f.write('%14.6e ' % (self.x[j,k],))
f.write('\n')
f.close()
def restore(self, filename, id):
"""Restore a previously-saved solution.
filename -- name of a file containing a solution previously
saved by a call to 'save'
id -- the ID tag of the solution
"""
#try:
(z, s) = _cantera.flow_restore(self.__flow_id, 0, filename, id)
#except:
# raise CanteraError()
self.setGrid(z)
self.x = array(s,'d')
self.xnew = array(self.x,'d')
for j in range(self.npts):
self.holdTemperature(j,self.x[j,2])
self.initial['T'] = (self.z, self.x[:,2])
self.setEnergyEqn('off')
if self.loglevel > 0:
print 'Solution ',`id`,'read from file',filename
print _cantera.readlog()
def outputTEC(self, plotfile="", title="", zone="c0", append=0):
"""Write the current solution to a file in TECPLOT format.
plotfile -- file name (required)
title -- plot title
zone -- zone name
append -- if append > 0, the output is appended to the file
"""
self.plotter.plot(fname = plotfile, title = title,
zone = zone, append=append)
def outputCSV(self, plotfile="", append=0):
"""Write the current solution to a file in CSV format.
plotfile -- file name (required)
append -- if append > 0, the output is appended to the file
"""
self.plotter.plot(fname = plotfile, fmt = 'EXCEL',
append=append)
def plot(self, i):
"""Plot solution component i. Requires the scipy package."""
from scipy import gplt
return gplt.plot(self.z, self.x[:,i])
def setBoundaries(self, left = None, right = None):
"""Install the boundary objects.
The type of boundary object determines the boundary conditions.
"""
nleft = 0
nright = 0
if left:
self.left = left
nleft = left.bdry_id()
if right:
self.right = right
nright = right.bdry_id()
_cantera.flow_setboundaries(self.__flow_id, nleft, nright)
"""
$Author$
$Revision$
$Date$
Copyright 2001 California Institute of Technology
"""

140
Cantera/python/Cantera/FlowBoundary.py
View file

@ -1,140 +0,0 @@
from Cantera import CanteraError
import _cantera
from Numeric import *
class FlowBoundary:
def __init__(self, type, phase, kin=None):
self.phase = phase
self.kin = kin
if kin:
self.__bdry_id = _cantera.bdry_new(type, phase.phase_id(), kin.kin_id())
else:
self.__bdry_id = _cantera.bdry_new(type, phase.phase_id(), 0)
FlowBoundary.set(self, mdot = 0.0, V = 0.0, T = phase.temperature(),
X = phase.moleFractions())
def __del__(self):
_cantera.bdry_del(self.__bdry_id)
def bdry_id(self):
return self.__bdry_id
def set(self, mdot = -999.0, V = -999.0, T = -999.0, X = None, Y = None):
if mdot > 0.0:
self.mdot = mdot
_cantera.bdry_set(self.__bdry_id, 1, mdot, None)
if V > 0.0:
self.V = V
_cantera.bdry_set(self.__bdry_id, 2, V, None)
if T > 0.0:
self.T = T
_cantera.bdry_set(self.__bdry_id, 3, T, None)
if X:
self.phase.setMoleFractions(X)
yy = self.phase.massFractions()
self.X = self.phase.moleFractions()
_cantera.bdry_set(self.__bdry_id, 4, 0.0, yy)
if Y:
self.phase.setMassFractions(Y)
yy = self.phase.massFractions()
self.X = self.phase.moleFractions()
_cantera.bdry_set(self.__bdry_id, 4, 0.0, yy)
class Inlet(FlowBoundary):
def __init__(self, phase):
FlowBoundary.__init__(self, 0, phase)
class Outlet(FlowBoundary):
def __init__(self, phase):
FlowBoundary.__init__(self, 1, phase)
def set(self, mdot = -999.0, V = -999.0, T = -999.0, X = None, Y = None):
raise CanteraError("outlet properties cannot be set.")
class Surface(FlowBoundary):
def __init__(self, phase, surf):
self._surf = surf
self.domainType = 1
FlowBoundary.__init__(self, 2, phase, surf)
self.__surf1d_id = _cantera.surf1d_new(self.kin.kin_id())
self.x = zeros((1,self.kin.nSpecies()),'d')
self.x[0,:] = self.kin.coverages()
self.species_on = []
self.species_off = []
def surf_id(self):
return self.__surf1d_id
def setCoverages(self, theta):
self.kin.setCoverages(theta)
self.x[0,:] = self.kin.coverages()
def coverages(self):
return self.x[0] # kin.coverages()
def show(self):
self.kin.setCoverages(self.x)
self.kin.show()
def shape(self):
return (1, self.kin.nSpecies())
def integrate(self, dt):
self.kin.setCoverages(self.x)
self.kin.integrate(dt)
c = self.kin.coverages()
self.x[0,:] = c
for k in range(len(c)):
#print k, c[k]
self.fixSpecies(k,c[k])
self.setSpeciesEqn(on = self.species_on, off = self.species_off)
def fixSpecies(self, k, c):
_cantera.surf1d_fixspecies(self.__surf1d_id, k, c)
def setTemperature(self, t):
_cantera.surf1d_settemperature(self.__surf1d_id, t)
def temperature(self):
return _cantera.surf1d_temperature(self.__surf1d_id)
def setMultiplier(self,f):
for k in range(self.kin.nSpecies()):
_cantera.surf1d_setmultiplier(self.__surf1d_id,k,f)
def setSpeciesEqn(self, on=None, off=None, loglevel=0):
"""Enable or disable surface species equations."""
self.species_on = on
self.species_off = off
on_msg = `on`
off_msg = `off`
if on == 'all':
on = self.kin.speciesNames()
on_msg = 'all species'
if off == 'all':
off = self.kin.speciesNames()
off_msg = 'all species'
a = zeros(self.kin.nSpecies(),'d')
if loglevel > 0:
if on:
print '\n\n--------------------------------------------------\n'
print 'Enabling species equation for',on_msg
print '\n--------------------------------------------------\n\n'
if off:
print '\n\n--------------------------------------------------\n'
print 'Disabling species equation for',off_msg
print '\n--------------------------------------------------\n\n'
if on:
for sp in on:
k = self.kin.speciesIndex(sp)
a[k] = 1.0
if off:
for sp in off:
k = self.kin.speciesIndex(sp)
a[k] = 0.0
_cantera.surf1d_solvespecies(self.__surf1d_id, len(a), a)
class SymmPlane(FlowBoundary):
def __init__(self, phase):
FlowBoundary.__init__(self, 3, phase)
def set(self, mdot = -999.0, V = -999.0, T = -999.0, X = None, Y = None):
raise CanteraError("outlet properties cannot be set.")

65
Cantera/python/Cantera/FlowPlotter.py
View file

@ -1,65 +0,0 @@
"""
Plotting of Flow1D solutions.
"""
from Numeric import shape, zeros, ones, array
class FlowPlotter:
"""FlowPlotter objects handle generating plots of Flow1D solutions.
This class is primarily designed for internal use by class Flow1D.
"""
def __init__(self, flow):
self.flow = flow
self.np = 0
self.nv = 0
self.names = ['z [m]', 'u [m/s]', 'V [1/s]',
'T [K]', 'lambda', 'P [Pa]'] + list(flow.gas.speciesNames())
def value(self, j, n):
if n == 0:
return self.flow.z[j]
elif n < 5:
return self.flow.x[j,n-1]
elif n == 5:
return self.flow.gas.pressure()
else:
return self.mf[j,n-6]
def plot(self, fmt = 'TECPLOT',
fname = 'plot.dat', title = 'plot',
zone = 'zone1',
moles = 1,
append = 0):
self.np, self.nv = shape(self.flow.x)
self.mf = zeros((self.np, self.nv - 4), 'd')
if moles:
for j in range(self.np):
self.flow.setGas(j)
self.mf[j,:] = self.flow.gas.moleFractions()
else:
for j in range(self.np):
self.flow.setGas(j)
self.mf[j,:] = self.flow.gas.massFractions()
if fmt == 'TECPLOT':
from tecplot import write_TECPLOT_zone
data = zeros((self.np,self.flow.gas.nSpecies()+6),'d')
data[:,0] = self.flow.z
data[:,1:5] = self.flow.x[:,0:4]
data[:,5] = ones(self.np)*self.flow.gas.pressure()
data[:,6:] = self.mf
write_TECPLOT_zone(fname, title, zone, self.names,
self.np, self.nv+2, append, data)
elif fmt == 'EXCEL':
from excel import write_CSV_data
write_CSV_data(fname, self.names,
self.np, self.nv+2, append, self)
else:
raise 'unknown format'+fmt

476
Cantera/python/Cantera/OneDim.py
View file

@ -1,476 +0,0 @@
from exceptions import *
import refine
import sys, types, copy, tempfile
import interp
import math
import _cantera
from Numeric import *
def print_heading(msg):
print '\n\n\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n'
print msg
print '\n'
class OneDim:
"""
One-dimensional, multi-domain problems.
Class OneDim allows solving multi-domain, one-dimensional,
steady-state problems implicitly. Each domain has a set of one or
more grid points, and each may have a different number of solution
components.
At ... N(I) algebraic residual
equations are defined for the N(I) solution components. Each
domain may have a different number of components.
The domains are linked in a linear chain, and are of two
types. Standard domains know nothing of their neighbors, and
evaluate their residual functions using only information in their
own domain. 'Connector' domains can also modify the residual
equations of their immediate neighbors, but only at the nearest
grid point. Every standard domain must be attached to a connector
at both ends. Connectors also serve to terminate the ends of the
chain.
"""
_timeint_options = ['ftime', 'min_timestep', 'max_timestep',
'nsteps', 'timestep', 'ts_jac_age']
_newton_options = ['max_jac_age', 'rtol', 'atol']
_output_options = ['loglevel', 'plotfile']
_options = _newton_options + _timeint_options + _output_options
def __init__(self, domains):
"""Create a new one-didmensional model from a list of domains. """
# instance variables
self._size = []
self._start = []
self._end = []
self._domain = [] # all domains
self._flow = [] # extended domains
self._shape = []
self._loc = 0
self._opt = {}
self.time = 0.0
self.x = array([0.0,],'d')
self._surf = []
self.npts = []
# local variables
dtype = [] # list of integer domain types
dlist = [] # list of integer domain ids
# add each domain
for d in domains:
if d.domainType == 0:
self.addFlow(d)
dlist.append(d.flow_id())
self.npts.append(d.nPoints())
elif d.domainType == 1:
self.addSurface(d)
dlist.append(d.surf_id())
self.npts.append(1)
elif d.domainType == 2:
self.addBoundary(d)
dlist.append(d.bndry_id())
self.npts.append(1)
else:
raise 'unknown domain type'
dtype.append(d.domainType)
self.__onedim_id = _cantera.onedim_new(len(dlist),
array(dlist,'i'),
array(dtype,'i'))
self.collect()
self.restoreDefaults();
self.ienergy = 0
self.ts_jac_age = 50
def __del__(self):
"""Delete the kernel object.
This does not delete the individual domains."""
_cantera.onedim_del(self.__onedim_id)
def addFlow(self, flow):
# add the domain to the list of all domains and to the list of
# extended domains
self._domain.append(flow)
self._flow.append(flow)
# set the index of this domain
flow.index = len(self._domain) - 1
np, nv = flow.shape()
self._shape.append((np,nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
def index(self, n, j, i):
np, nv = self._shape[i]
return self._start[n] + nv*j + i
## def resetEnergy(self):
## i = 0
## ilast = self.ienergy
## for f in self._flow:
## i += f.resetEnergy()
## self.ienergy = 0
## return 0# ilast
def finish(self):
"""
Update the solution in each domain based on the global solution.
This method is called by function 'solve' when a converged
solution has been found, just prior to grid refinement.
"""
for i in range(len(self._domain)):
self._domain[i].x = self.solution(i)
def solution(self, i):
""" Return the solution array for domain i.
The returned array has the shape (points,
components) appropriate for domain i. """
x = self.x[self._start[i]:self._end[i]]
dx = reshape(x, self._domain[i].shape())
return dx
def resid(self, i):
"""
The residual matrix for domain i.
The returned array has the shape
(points, components) appropriate for domain i.
"""
self.ssnorm()
x = self.xnew[self._start[i]:self._end[i]]
dx = reshape(x, self._domain[i].shape())
return dx
def addSurface(self, surf):
"""Add a surface domain."""
self._surf.append(surf)
self._domain.append(surf)
surf.index = len(self._domain) - 1
nv = surf.kin.nSpecies()
np = 1
self._shape.append((np,nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
def addBoundary(self, b):
"""Add a boundary domain."""
#self._surf.append(surf)
self._domain.append(b)
b.index = len(self._domain) - 1
nv = 2 # surf.kin.nSpecies()
np = 1
self._shape.append((np,nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
def setNewtonOptions(self, max_jac_age = 5):
_cantera.onedim_setnewtonoptions(self.__onedim_id, max_jac_age)
def newton_solve(self, loglevel = 0):
"""Damped Newton iteration.
This method invokes C++ method 'solve' of kernel class
'OneDim' on the current solution. The solution is only
modified if the damped Newton process leads to a fully
converged solution. Otherwise, an exception is raised.
"""
iok = _cantera.onedim_solve(self.__onedim_id, self.x,
self.xnew, loglevel)
#if loglevel > 0: print _cantera.readlog()
if iok >= 0:
_cantera.copy(size(self.x),self.xnew,self.x)
elif iok > -10:
raise CanteraError()
else:
raise 'iok = '+`iok`
return iok
def collect(self):
"""Collect the state information from each domain to
construct the global solution vector."""
n = 0
strt = [] # list of start locations for each domain
self.npts = []
nd = len(self._domain)
for d in self._domain:
strt.append(n)
n += size(d.x)
self.npts.append(d.shape()[0])
strt.append(n)
self.x = zeros(n,'d')
self.xnew = zeros(n,'d')
# set the portion of the global solution vector corresponding
# to each domain to the flattened solution matrix for that
# domain
for i in range(nd):
self.x[strt[i]:strt[i+1]] = reshape(self._domain[i].x,(-1,))
def ssnorm(self):
"""Max norm of the steady-state residual."""
n = _cantera.onedim_ssnorm(self.__onedim_id, self.x, self.xnew)
return n
def setSteadyMode(self):
"""Prepare to solve the steady-state problem."""
return _cantera.onedim_setsteadymode(self.__onedim_id)
def setTransientMode(self, dt):
"""Prepare for time-stepping with timestep dt.
Must be called before each step."""
return _cantera.onedim_settransientmode(self.__onedim_id, dt, self.x)
def option(self, key):
"""Return the value of an option."""
return self._opt[key]
def restoreDefaults(self):
"""Restore default options."""
self._opt = {}
self.setOptions(
max_jac_age = 20,
ts_jac_age = 30,
timestep = 1.e-6,
min_timestep = 1.e-12,
max_timestep = 0.1,
nsteps = [1,2,4,8,20],
ftime = 3.0,
plotfile = ""
)
def setOptions(self, **options):
"""
Set options.
Time stepping:
nsteps -- number of steps.
min_timestep -- minimum timestep
max_timestep -- maximum timestep
ftime -- factor by which to increase the timestep for next
set of 'nsteps' timesteps
Newton solver:
max_jac_age -- maximum number of times Jacobian will be used
before re-evaluating
rtol -- relative error tolerance
atol -- absolute error tolerance
Output:
loglevel -- controls amount of diagnostic output
plotfile -- file to write plot data for intermediate solutions
"""
for kw in options.keys():
if kw in OneDim._options:
self._opt[kw] = options[kw]
else:
raise OptionError(kw)
if kw in OneDim._newton_options:
self.setNewtonOptions(max_jac_age =
self._opt['max_jac_age'])
def refine(self, loglevel = 2):
"""Refine the grid of every flow domain."""
new_points = 0
for f in self._flow:
new_points += f.refine(loglevel)
if new_points > 0:
self.collect()
_cantera.onedim_resize(self.__onedim_id)
self._shape = []
self._size = []
self._start = []
self._end = []
self._loc = 0
for d in self._domain:
np, nv = d.shape()
self._shape.append((np, nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
return new_points
def prune(self, loglevel = 2):
"""Prune the grid of every flow domain."""
rem_points = 0
for f in self._flow:
rem_points += f.prune(loglevel)
if rem_points > 0:
self.collect()
_cantera.onedim_resize(self.__onedim_id)
self._shape = []
self._size = []
self._start = []
self._end = []
self._loc = 0
for d in self._domain:
np, nv = d.shape()
self._shape.append((np, nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
return rem_points
def setEnergyFactor(self, e):
for f in self._flow:
f.setEnergyFactor(e)
def restore(self, n, file, soln, loglevel = 2):
"""Read the solution for domain n from a file."""
self._domain[n].restore(file, soln)
self.collect()
_cantera.onedim_resize(self.__onedim_id)
self._shape = []
self._size = []
self._start = []
self._end = []
self._loc = 0
for d in self._domain:
np, nv = d.shape()
self._shape.append((np, nv))
self._size.append(np*nv)
self._start.append(self._loc)
self._loc += np*nv
self._end.append(self._loc)
def c_timeStep(self, nsteps, dt, loglevel = 0):
dtnew = _cantera.onedim_timestep(self.__onedim_id, nsteps, dt,
self.x, self.xnew, loglevel)
#print _cantera.readlog()
return dtnew
def py_timeStep(self, nsteps, dt, loglevel = 0):
"""Take time steps using Backward Euler.
nsteps -- number of steps
dt -- initial step size
loglevel -- controls amount of printed diagnostics
"""
self.setNewtonOptions(max_jac_age = self._opt['ts_jac_age'])
if loglevel > 0:
print_heading('Begin time integration.\n\n')
print(' step size (s) log10(ss) ')
print('===============================')
n = 0
maxdt = self._opt['max_timestep']
while n < nsteps:
if loglevel > 0:
ss = self.ssnorm()
str = ' %4d %10.4g %10.4g' % (n,dt,math.log10(ss))
print str,
try:
self.setTransientMode(dt)
m = self.newton_solve(loglevel-1)
self.time += dt
n += 1
if m == 100:
dt *= 1.5
if dt > maxdt: dt = maxdt
if loglevel > 0: print
except CanteraError:
#print self.resid(1)[:,0]
if loglevel > 0: print '...failure.'
dt *= 0.5
if dt < 1.e-16:
self._domain[0].show()
raise CanteraError('Time integration failed.')
self.setSteadyMode()
self.setNewtonOptions(max_jac_age =
self._opt['max_jac_age'])
return dt
def show(self):
for d in self._domain:
d.show()
def showStatistics(self):
_cantera.onedim_writestats(self.__onedim_id)
#print _cantera.readlog()
def save(self, filename, id, desc=""):
"""Save a solution to a file.
filename -- file name. If the file, does not exist, it will
be created. The save files are xml files, and the
filename should have the extension '.xml'. If it
does not, this extension will be appended to the
name.
id -- the ID tag of the solution. Multiple solutions may
be saved to the same file. Specifying a unique ID
tag allows this solution to selected later by
method 'restore'.
"""
fn = filename
extn = filename[-4:]
if extn <> '.xml' and extn <> '.XML':
fn = filename + '.xml'
_cantera.onedim_save(self.__onedim_id, fn, id, desc, self.x)
#print _cantera.readlog()

227
Cantera/python/Cantera/Thermo.py
View file

@ -1,227 +0,0 @@
"""
Module Thermo
"""
DEPRECATED
from Cantera.Phase import Phase
import ctthermo
import ctphase
import types
def thermoIndex(id):
return ctthermo.thermoIndex(id)
class Thermo:
_equilmap = {'TP':104,'TV':100,'HP':101,'SP':102,'SV':107,'UV':105,
'PT':104,'VT':100,'PH':101,'PS':102,'VS':107,'VU':105}
def __init__(self, eostype=1, phase=None, sptherm=0,
root=None, id=None, index=-1):
self.__phase = None
self.cthermo = None
self._owner = 1
self.idtag = ""
#if thermoIndex(id) > 0:
# index = thermoIndex(id)
if index >= 0:
# create a Python wrapper for an existing kernel
# Thermo instance
self.cthermo = index
self.__phase = Phase(index = ctthermo.phase(index))
self._owner = 0
elif root:
# create a new kernel instance from an XML specification
self.cthermo, ph = ctthermo.ThermoFromXML(root._xml_id, id)
self.__phase = Phase(ph)
self.idtag = id
else:
# create a new kernel instance with specified parameters
self.cthermo = ctthermo.Thermo(eostype, phase.phase_id(), sptherm)
self.__phase = phase
def __del__(self):
if self._owner:
ctthermo.delete(self.cthermo)
def importFromXML(self, xml_root, id):
ctthermo.import_xml(self.cthermo, xml_root._xml_id, id)
def thermophase(self):
return self.__phase
def refPressure(self):
"""Reference pressure [Pa].
All standard-state thermodynamic properties are for this pressure.
"""
return ctthermo.refpressure(self.cthermo)
def minTemp(self, sp=-1):
""" Minimum temperature for which the parameterization of
standard-state thermodynamic properties vs. T for species 'sp'
is valid. If no species is specified, the value returned is
the maximum value of minTemp for any one species, and
therefore is the minimum temperature at which mixture
thermodynamic properties are valid."""
return ctthermo.mintemp(self.cthermo, self.speciesIndex(sp))
def maxTemp(self, sp=-1):
""" Maximum temperature for which the parameterization of
standard-state thermodynamic properties vs. T for species 'sp'
is valid. If no species is specified, the value returned is
the minimum value of maxTemp for any one species, and
therefore is the maximum temperature at which mixture
thermodynamic properties are valid."""
return ctthermo.maxtemp(self.cthermo, self.speciesIndex(sp))
def enthalpy_mole(self):
""" The molar enthalpy [J/kmol]."""
return ctthermo.getfp(self.cthermo,1)
def intEnergy_mole(self):
""" The molar internal energy [J/kmol]."""
return ctthermo.getfp(self.cthermo,2)
def entropy_mole(self):
""" The molar entropy [J/kmol/K]."""
return ctthermo.getfp(self.cthermo,3)
def gibbs_mole(self):
""" The molar Gibbs function [J/kmol]."""
return ctthermo.getfp(self.cthermo,4)
def cp_mole(self):
""" The molar heat capacity at constant pressure [J/kmol/K]."""
return ctthermo.getfp(self.cthermo,5)
def cv_mole(self):
""" The molar heat capacity at constant volume [J/kmol/K]."""
return ctthermo.getfp(self.cthermo,6)
def pressure(self):
""" The pressure [Pa]."""
return ctthermo.getfp(self.cthermo,7)
def chemPotentials(self):
"""Species chemical potentials.
This method returns an array containing the species
chemical potentials [J/kmol]. The expressions used to
compute these depend on the model implemented by the
underlying kernel thermo manager."""
return ctthermo.getarray(self.cthermo,20)
def enthalpies_RT(self):
"""Pure species non-dimensional enthalpies.
This method returns an array containing the pure-species
standard-state enthalpies divided by RT. For gaseous species,
these values are ideal gas enthalpies."""
return ctthermo.getarray(self.cthermo,23)
def entropies_R(self):
"""Pure species non-dimensional entropies.
This method returns an array containing the pure-species
standard-state entropies divided by R. For gaseous species,
these values are ideal gas entropies."""
return ctthermo.getarray(self.cthermo,24)
def gibbs_RT(self):
"""Pure species non-dimensional Gibbs free energies.
This method returns an array containing the pure-species
standard-state Gibbs free energies divided by R.
For gaseous species, these are ideal gas values."""
return (ctthermo.getarray(self.cthermo,23)
- ctthermo.getarray(self.cthermo,24))
def cp_R(self):
"""Pure species non-dimensional heat capacities
at constant pressure.
This method returns an array containing the pure-species
standard-state heat capacities divided by R. For gaseous
species, these values are ideal gas heat capacities."""
return ctthermo.getarray(self.cthermo,25)
def setPressure(self, p):
"""Set the pressure [Pa]."""
ctthermo.setfp(self.cthermo,1,p,0.0)
def enthalpy_mass(self):
"""Specific enthalpy [J/kg]."""
return ctthermo.getfp(self.cthermo,8)
def intEnergy_mass(self):
"""Specific internal energy [J/kg]."""
return ctthermo.getfp(self.cthermo,9)
def entropy_mass(self):
"""Specific entropy [J/kg/K]."""
return ctthermo.getfp(self.cthermo,10)
def gibbs_mass(self):
"""Specific Gibbs free energy [J/kg]."""
return ctthermo.getfp(self.cthermo,11)
def cp_mass(self):
"""Specific heat at constant pressure [J/kg/K]."""
return ctthermo.getfp(self.cthermo,12)
def cv_mass(self):
"""Specific heat at constant volume [J/kg/K]."""
return ctthermo.getfp(self.cthermo,13)
def setState_HP(self, h, p):
"""Set the state by specifying the specific enthalpy and
the pressure."""
ctthermo.setfp(self.cthermo, 2, h, p)
def setState_UV(self, u, v):
"""Set the state by specifying the specific internal
energy and the specific volume."""
ctthermo.setfp(self.cthermo, 3, u, v)
def setState_SV(self, s, v):
"""Set the state by specifying the specific entropy
and the specific volume."""
ctthermo.setfp(self.cthermo, 4, s, v)
def setState_SP(self, s, p):
"""Set the state by specifying the specific entropy
energy and the pressure."""
ctthermo.setfp(self.cthermo, 5, s, p)
def equilibrate(self, XY):
"""Set to a state of chemical equilibrium holding property pair
'XY' constant. The pair is specified by a two-letter string,
which must be one of the set
['TP','TV','HP','SP','SV','UV','PT','VT','PH','PS','VS','VU'].
If H, U, S, or V is specified, the value must be the specific
value (per unit mass).
"""
ixy = Thermo._equilmap[XY]
if ixy > 0:
ctthermo.equil(self.cthermo, ixy)
else:
raise 'invalid equilibrium option: '+XY

56
Cantera/python/Cantera/boundaries1D.py
View file

@ -1,56 +0,0 @@
from Cantera import CanteraError
import _cantera
from Numeric import *
import types
class Boundary1D:
def __init__(self, type):
self.domainType = 2
self.__bndry_id = _cantera.bndry_new(type)
self.x = zeros((1,2),'d')
self.x[0,1] = 300.0
def __del__(self):
_cantera.bndry_del(self.__bndry_id)
def shape(self):
return (1,2)
def show(self):
pass
def restore(self, file='', solution=''):
pass
def bndry_id(self):
return self.__bndry_id
def set(self, mdot = -999.0, V = -999.0, T = -999.0, X = None, Y = None):
if mdot > 0.0:
self.mdot = mdot
self.x[0,0] = mdot
_cantera.bndry_setmdot(self.__bndry_id, mdot)
if T > 0.0:
self.T = T
self.x[0,1] = T
_cantera.bndry_settemperature(self.__bndry_id, T)
if V <> -999.0:
self.V = V
_cantera.bndry_setspreadrate(self.__bndry_id, V)
if X:
self.X = X
if type(X) == types.StringType:
_cantera.bndry_setxinbyname(self.__bndry_id, X)
else:
_cantera.bndry_setxin(self.__bndry_id, X)
def Inlet1D():
return Boundary1D(1)
def Symm1D():
return Boundary1D(2)
def Surf1D():
return Boundary1D(3)

21
Cantera/python/Cantera/ck2ctml.py
View file

@ -1,21 +0,0 @@
import _cantera
from exceptions import *
import os
def ck2ctml(infile = '', thermo = '-', transport = '-', outfile = '', id = ''):
if not infile:
raise CanteraError('No input file specified')
fname = os.path.basename(infile)
ff = os.path.splitext(fname)
if len(ff) == 2:
mechname = ff[0]
else:
mechname = ff
if not outfile:
outfile = mechname + '.xml'
if not id:
id = mechname
_cantera.ck2ctml(infile, thermo, transport, outfile, id)

483
Cantera/python/Cantera/flame.py
View file

@ -1,483 +0,0 @@
from Cantera import OneAtm
from Cantera.exceptions import CanteraError
from Cantera.Flow import Flow1D
from Cantera.boundaries1D import Inlet1D, Surf1D, Symm1D
from Numeric import array, zeros, arrayrange
from Cantera.gases import IdealGasMix, GRI30
from Cantera.solve import solve
#from Cantera.esolve import esolve
from Cantera.OneDim import OneDim
from Cantera.FlowBoundary import Inlet, Outlet, SymmPlane
from Cantera import stoich
import math
class BurnerFlame:
"""One-dimensional flat, premixed flames.
flame = BurnerFlame(gas, domain, fuel, oxidizer, inert, grid, pressure)
arguments:
gas --- an object representing the gas mixture
domain --- [zmin, zmax]
fuel --- a string specifying the fuel stream composition
oxidizer --- a string specifying the oxidizer stream composition
inert --- a string specifying the composition of an inert
stream (optional)
grid --- a sequence defining the initial grid. The first point
should be zmin, and the last one zmax. If omitted,
a default grid will be used.
pressure --- the pressure, which is treated as constant.
example:
flame = BurnerFlame(gas = GRI30(),
domain = [0.0, 10.0*units.cm],
fuel = 'CH4:1',
oxidizer = 'O2:1,N2:3.76',
grid = [0.0, 0.01, 0.03, 0.06, 0.1],
pressure = OneAtm)
"""
def __init__(self, gas = None, domain = None,
fuel = '', oxidizer = '', inert = '',
grid = None, pressure = -1.0):
# check that all required inputs have been specified
if not gas or not domain or not fuel or not oxidizer or not pressure:
raise self.__doc__
self.gas = gas
self.p = pressure
dx = (domain[1] - domain[0])
# if no grid specified, use this one that concentrates points
# near the burner
if grid == None:
grid = dx * array([0.0, 0.01, 0.03, 0.1, 0.3, 0.6, 1.0])
self.__flow = Flow1D(flow_type = 'OneDim', gas = gas,
grid = grid, pressure = self.p)
#------ these methods are deprecated, but still needed for now.
self.inlet = Inlet(gas)
self.outlet = Outlet(gas)
self.__flow.setBoundaries(left = self.inlet, right = self.outlet)
#----------------------------------------------------------------
# The container contains only the Flow1D object. Should be
# modified at some point to contain an Inlet1D and an Outlet1D
# object.
self.__container = OneDim([self.__flow])
self.start = 0
# get the compositions of the fuel and oxidizer streams, and
# calculate the fuel/oxidizer ratio for stoichiometric
# combustion
gas.setMoleFractions(fuel)
self._xfuel = gas.moleFractions()
gas.setMoleFractions(oxidizer)
self._xox = gas.moleFractions()
if inert:
gas.setMoleFractions(inert)
self._xinert = gas.moleFractions()
else:
self._xinert = zeros(gas.nSpecies(),'d')
self._stoich_FO = stoich.stoich_fuel_to_oxidizer(gas, fuel, oxidizer)
# TODO: accout for inert stream
def setEquivRatio(self, phi):
"""Set the equivalence ratio."""
f_flow = self._stoich_FO * phi
comp = f_flow * self._xfuel + self._xox
self.gas.setState_PX(self.p, comp)
self.inlet.set(X = self.gas.moleFractions())
def setEquilProducts(self):
"""Generate a starting estimate for the flame state.
The following procedure is used:
1) At the burner, the composition is the specified inlet
composition;
2) The last 80% of the domain has constant composition and
temperature corresponding to the adiabatic equilibrium
solution;
3) In the initial 20%, the composition and temperature vary linearly
from the inlet values to the equilibrium values.
"""
x0 = self.inlet.X
self.gas.setState_TPX(self.inlet.T, self.p, x0)
rho0 = self.gas.density()
mdot = self.inlet.mdot
self.gas.equilibrate('HP')
xp = self.gas.moleFractions()
xinit = {}
z0 = 0.2
teq = self.gas.temperature()
rhoeq = self.gas.density()
xinit['T'] = [(0.0, self.inlet.T), (z0, teq), (1.0, teq)]
xinit['u'] = [(0.0, mdot/rho0), (z0, mdot/rhoeq), (1.0, mdot/rhoeq)]
for k in range(self.gas.nSpecies()):
nm = self.gas.speciesName(k)
x = [(0.0, x0[k]), (z0, xp[k]), (1.0, xp[k])]
xinit[nm] = x
self.__flow.setInitialProfiles(xinit)
def plot(self, plotfile = '', title = '', fmt = 'TECPLOT',
zone = 'c0', append = 0):
"""Plot the current solution."""
self.__flow.plotter.plot(fname = plotfile, title = title,
fmt = fmt, zone = zone, append=append)
def setInitialProfiles(self, **init):
"""Specify estimates for the initial profiles.
"""
self.__flow.setInitialProfiles(init)
self.start = 1
def restore(self, src = '', solution = ''):
"""Start from a previously-saved solution."""
self.__container.restore(0, src, solution)
self.start = 1
def setTolerances(self, V = None, T = None, Y = None):
"""Set tolerances for convergence for velocity, temperature,
and mass fractions."""
self.__flow.setTolerances( V, V, T, Y)
def prune(self, loglevel = 2):
"""Remove unneeded grid points.
This method attempts to remove each grid point one by one, and
calls 'refine' each time to see whether it puts it back. If it does,
the point is not removed, otherwise it is.
"""
self.__container.prune(loglevel)
def refine(self, loglevel = 2):
"""Refine the grid using the current grid refinement parameters."""
self.__container.refine(loglevel)
def show(self):
"""Print a summary of the current solution to the screen."""
self.__flow.show()
def stretch(self, factor):
"""Stretch the grid by 'factor'"""
self.__flow.setGrid(factor*self.__flow.z)
def set(self, **opt):
"""Set options.
The options that may be set are:
energy --- 'on' or 'off'. If 'on', the energy equation is
solved; otherwise, the temperature is held to the specified
profile.
pressure --- the pressure in Pa.
mdot --- the inlet mass flow rate per unit area.
equiv_ratio --- the equivalence ratio
T_burner --- the burner surface temperature [K].
refine --- a triplet specifying the refinement criteria.
See refine.py for more information.
tol --- error tolerances for u, V, T, and Y.
max_jac_age --- the maximum number of times to use a Jacobian
before recomputing it.
timesteps --- number and duration of time steps to take
when Newton iteration fails. The format is
( number_sequence, initial_stepsize )
These parameters can be changed as the solution proceeds."""
# TODO: is this necessary?
if self.__container == None:
self.__container = OneDim([self.__flow,])
for o in opt.keys():
v = opt[o]
if o == 'energy':
self.__flow.setEnergyEqn(v,loglevel=1)
elif o == 'pressure':
self.p = v
self.__flow.setPressure(v)
elif o == 'mdot':
self.inlet.set(mdot = v)
elif o == 'equiv_ratio':
self.setEquivRatio(v)
elif o == 'T_burner':
self.inlet.set(T = v)
elif o == 'refine':
self.__flow.refiner.delta = v
elif o == 'tol':
self.__flow.setTolerances(u = v, V = v, T = v, Y = v)
elif o == 'max_jac_age':
self.__container.setOptions(max_jac_age = v)
elif o == 'timesteps':
self.__container.setOptions(nsteps = v[0], timestep = v[1])
def solve(self, loglevel = 0):
""" Solve the flame equations.
If no starting estimate has been given, setEquilProducts()
is called to generate one.
"""
if not self.start:
self.setEquilProducts()
self.start = 1
solve(self.__container, loglevel = loglevel, refine_grid = 1)
## def esolve(self, loglevel = 0, efactor = 1.0e4):
## if not self.start:
## self.setEquilProducts()
## self.start = 1
## esolve(self.__container, efactor = efactor, loglevel = loglevel, refine_grid = 1)
def save(self, soln, desc, file = 'flame.xml'):
"""Save the current solution.
soln --- string to identify this solution in the file.
desc --- descriptive text string.
file --- file name.
"""
self.__container.save(file, soln, desc)
def showStatistics(self):
"""Show numerical statistics."""
self.__container.showStatistics()
class StagnationFlame:
"""Axisymmetric premixed stagnation-point flames.
flame = StagnationFlame(gas, domain, fuel, oxidizer, inert, grid, pressure)
example:
flame = BurnerFlame(gas = GRI30(),
domain = [0.0, 10.0*units.cm],
fuel = 'CH4:1',
oxidizer = 'O2:1,N2:3.76',
grid = [0.0, 0.01, 0.03, 0.06, 0.1],
pressure = OneAtm)
"""
def __init__(self, gas = None, domain = None,
fuel = '', oxidizer = '', inert = '',
grid = None, pressure = -1.0):
if not gas or not domain or not fuel or not oxidizer or not pressure:
raise self.__doc__
self.gas = gas
self.p = pressure
dx = (domain[1] - domain[0])
self.dx = dx
if grid == None:
grid = dx * array([0.0, 0.01, 0.03, 0.1, 0.3, 0.6, 1.0])
self.__flow = Flow1D(flow_type = 'Stag', gas = gas,
grid = grid, pressure = self.p)
self.__left = Inlet1D()
self.__right = Surf1D()
self.__container = OneDim([self.__left, self.__flow, self.__right])
self.start = 0
# get the compositions of the fuel and oxidizer streams, and
# calculate the fuel/oxidizer ratio for stoichiometric
# combustion
gas.setMoleFractions(fuel)
self._xfuel = gas.moleFractions()
gas.setMoleFractions(oxidizer)
self._xox = gas.moleFractions()
if inert:
gas.setMoleFractions(inert)
self._xinert = gas.moleFractions()
else:
self._xinert = zeros(gas.nSpecies(),'d')
self._stoich_FO = stoich.stoich_fuel_to_oxidizer(gas, fuel, oxidizer)
def nPoints(self):
return len(self.__flow.z)
def setEquivRatio(self, phi):
"""Set the equivalence ratio."""
f_flow = self._stoich_FO * phi
comp = f_flow * self._xfuel + self._xox
self.gas.setState_PX(self.p, comp)
self.__left.set(X = self.gas.moleFractions())
def setEquilProducts(self):
"""Set the flame state to chemical equilibrium.
This is useful to generate a starting estimate.
"""
x0 = self.__left.X
self.gas.setState_TPX(self.__left.T, self.p, x0)
rho0 = self.gas.density()
mdot = self.__left.mdot
self.gas.equilibrate('HP')
xp = self.gas.moleFractions()
xinit = {}
z0 = 0.2
teq = self.gas.temperature()
rhoeq = self.gas.density()
re = self.dx * mdot / self.gas.viscosity()
z1 = 1.0 - 1.0/math.sqrt(re)
tw = self.__right.T
self.gas.setState_TPX(tw, self.p, x0)
self.gas.equilibrate('TP')
x1 = self.gas.moleFractions()
rho1 = self.gas.density()
xinit['T'] = [(0.0, self.__left.T), (z0, teq), (z1, teq),
(1.0, tw)]
xinit['u'] = [(0.0, mdot/rho0), (1.0, 0.0)]
xinit['V'] = [(0.0, 0.0), (z1, mdot/(rhoeq*z1*self.dx)), (1.0, 0.0)]
for k in range(self.gas.nSpecies()):
nm = self.gas.speciesName(k)
x = [(0.0, x0[k]), (z0, xp[k]), (z1, xp[k]), (1.0, x1[k])]
xinit[nm] = x
self.__flow.setInitialProfiles(xinit)
def plot(self, plotfile = '', title = '', fmt = 'TECPLOT',
zone = 'c0', append = 0):
self.__flow.plotter.plot(fname = plotfile, title = title,
fmt = fmt, zone = zone, append=append)
def setInitialProfiles(self, **init):
self.__flow.setInitialProfiles(init)
self.start = 1
def resid(self):
return self.__container.resid(1)
def restore(self, src = '', solution = ''):
self.__container.restore(1,src, solution)
self.start = 1
def setTolerances(self, V = None, T = None, Y = None):
self.__flow.setTolerances( V, V, T, Y)
def show(self):
self.__flow.show()
def stretch(self, factor):
self.__flow.setGrid(factor*self.__flow.z)
def enableEnergy(self, pt):
self.__flow.setEnergyEqn('on',loglevel=1,pt=pt)
def prune(self, loglevel = 2):
self.__container.prune(loglevel)
def refine(self, loglevel = 2):
self.__container.refine(loglevel)
def set(self, **opt):
if self.__container == None:
self.__container = OneDim([self.__flow,])
for o in opt.keys():
v = opt[o]
if o == 'energy':
self.__flow.setEnergyEqn(v,loglevel=1)
elif o == 'pressure':
self.p = v
self.__flow.setPressure(v)
elif o == 'mdot':
self.__left.set(mdot = v)
elif o == 'equiv_ratio':
self.setEquivRatio(v)
elif o == 'T_burner':
self.__left.set(T = v)
elif o == 'spreadingRate':
self.__left.set(V = v)
elif o == 'T_surface':
self.__right.set(T = v)
elif o == 'refine':
self.__flow.refiner.delta = v
elif o == 'efactor':
self.__flow.setEnergyFactor(v)
elif o == 'tol':
self.__flow.setTolerances(u = v, V = v, T = v, Y = v)
elif o == 'max_jac_age':
self.__container.setOptions(max_jac_age = v)
elif o == 'jac_age':
self.__container.setOptions(max_jac_age = v[0])
self.__container.setOptions(ts_jac_age = v[1])
elif o == 'timesteps':
self.__container.setOptions(nsteps = v[0], timestep = v[1])
else:
raise CanteraError("unknown option: "+o)
def solve(self, loglevel = 0):
if not self.start:
self.setEquilProducts()
self.start = 1
solve(self.__container, loglevel = loglevel, refine_grid = 1)
def save(self, soln, desc, file = 'flame.xml'):
self.__container.save(file, soln, desc)
def showStatistics(self):
self.__container.showStatistics()

31
Cantera/python/Cantera/pip.py
View file

@ -1,31 +0,0 @@
import sys, os
from tempfile import mktemp
def process(name):
parts = name.split('.')
base = parts[0]
if len(parts) == 2:
ext = parts[1]
fname = mktemp('.py')
fo = open(fname,'w')
txt = """from Cantera.ctml_writer import *
import sys, os
f = sys.argv[1]
b = sys.argv[2]
try:
os.remove(b+'.xml')
except:
pass
execfile(f)
write()
"""
fo.write(txt)
fo.close()
cmd = sys.executable+' '+fname+' '+name+' '+base
err = os.system(cmd)
os.remove(fname)
if err:
sys.exit(-1)

295
Cantera/python/Cantera/refine.py
View file

@ -1,295 +0,0 @@
"""Grid refinement.
Suppose you have a monotonic NumPy array of grid points 'z', and a
solution array soln[j,n] that contains 3 three solution components
denoted 'a', 'b', and 'c', evaluated at the grid points. To refine the
grid based on components 'a' and 'b' but not 'c', do the following.
>>> from refine import Refiner
>>> r = Refiner([(0, 'a'), (1, 'b')])
>>> new_grid, new_soln = r.refine(grid, soln)
"""
import Numeric
import math
from Cantera import CanteraError
from Cantera import interp
def eps():
"""Return the square root of machine precision."""
e = 1.0
while 1.0 + e <> 1.0: e = 0.5*e
return math.sqrt(e)
def delta(f):
"""Given an array f, return an array of the difference in
adjacent values."""
n = len(f)
d = Numeric.zeros(n-1,'d')
for j in range(n-1):
d[j] = f[j+1] - f[j]
return d
def slope(z, f):
"""Given arrays z and f, return an array of the slopes df/dz in
each interval."""
n = len(z)
s = Numeric.zeros(n-1,'d')
for j in range(n-1):
s[j] = ((f[j+1] - f[j])/(z[j+1] - z[j]))
return Numeric.array(s,'d')
class RefineError(CanteraError):
def __init__(self, msg):
self.msg = 'Grid refinement error!\n'+msg
class Refiner:
"""Grid refiner.
Attributes:
components -- sequence of (number, name) pairs specifying the
components of the solution to use for grid refinement. The number
is used to access the component in the solution array, and the
name is used only for diagnostic messages.
max_delta -- Maximum tolerated difference in solution values
between neighboring grid points, expressed as a fraction between 0
and 1 of the total range of the component over all grid
points. Default: 0.8 (minimal refinement).
max_delta_slope -- Maximum tolerated difference in solution slopes
between neighboring grid intervals, expressed as a fraction between 0
and 1 of the total range of the component over all grid
points. Default: 0.8 (minimal refinement).
"""
def __init__(self, components = [], delta = (2.0, 0.1, 0.2), names = []):
self.components = components
self.delta = delta
self.names = names
self.loglevel = 2
self.eps = eps()
self.min_range = 0.01
self.direction = 1
self.fctr = 1.0
self.ok = 0
def prune(self, grid = None, solution = None, threshold = None):
n0 = len(grid)
g = grid
sol = solution
self.fctr = 1.0
savedir = self.direction
ll = self.loglevel
#self.loglevel = 0
j = 1
while j < len(g)-1:
g0 = g
s0 = sol
pt = g[j]
nn = len(g)
# remove point j
g = Numeric.take(g,range(0,j)+range(j+1,nn))
# remove row j
sol = Numeric.take(sol, range(0,j)+range(j+1,nn))
np = len(g)
self.direction = 1
gnew, gn, snew, ok = self.refine(g, sol, threshold)
if (len(gnew) > np):
g = g0
sol = s0
j += 1
if ll > 0:
print 'cannot remove point at ',pt
else:
if ll > 0:
print 'removed point at ',pt
self.loglevel = ll
self.fctr = 0.2
self.direction = savedir
return (g, sol)
def refine(self, grid = None, solution = None, threshold = None, prune = 1):
self.ok = 0
# grid parameters
n0 = len(grid)
dz0 = grid[-1] - grid[0]
maxpts = self.fctr*n0 + 1
ncomp = Numeric.shape(solution)[1]
if threshold:
self.threshold = threshold
else:
self.threshold = self.eps * Numeric.ones(ncomp, 'd')
if Numeric.shape(solution)[0] <> n0:
raise RefineError('Number of solution points differs from '+
'number of grid points.')
# if the solution components to examine for refinement have
# not been specified, use all components.
nc = Numeric.shape(solution)[1]
if not self.components: self.components = range(nc)
c = {}
p = {}
dz = delta(grid)
for j in range(1,n0-1):
if dz[j] > self.delta[0]*dz[j-1]:
p[j] = 1
c['point '+`j`] = 1
if dz[j] < dz[j-1]/self.delta[0]:
p[j-1] = 1
c['point '+`j-1`] = 1
for i in self.components:
try:
name = self.names[i]
except:
name = 'component '+`i`
# get component i at all points, and compute its slope
v = solution[:,i]
s = slope(grid, v)
# compute the change in value and slope
dv = delta(v)
ds = delta(s)
# find the range of values and slopes
vmin = min(v)
vmax = max(v)
smin = min(s)
smax = max(s)
# max absolute values of v and s
aa = max((abs(vmax), abs(vmin)))
ss = max((abs(smax), abs(smin)))
# refine based on component i only if the range of v is
# greater than a fraction 'min_range' of max |v|. This
# eliminates components that consist of small fluctuations
# on a constant background.
if (vmax - vmin) > self.min_range*aa:
# maximum allowable difference in value between
# adjacent points.
dmax = self.delta[1]*(vmax - vmin) + self.threshold[i]
for j in range(len(dv)):
r = abs(dv[j])/dmax
if r > 1.0:
p[j] = 1
c[name] = 1
# refine based on the slope of component i only if the
# range of s is greater than a fraction 'min_range' of max
# |s|. This eliminates components that consist of small
# fluctuations on a constant slope background.
if (smax - smin) > self.min_range*ss:
# maximum allowable difference in slope between
# adjacent points.
dmax = self.delta[2]*(smax - smin)
for j in range(len(ds)):
r = abs(ds[j]) / (dmax + self.threshold[i]/dz[j])
if r > 1:
c[name] = 1
p[j] = 1
p[j+1] = 1
if len(p) == 0: self.ok = 1
znew = []
nnew = len(p)
nadded = nnew
if self.loglevel > 0:
if nnew > 0:
print '\nRefining grid.'
print 'New points inserted after grid points',
for j in range(n0 - 1):
znew.append(grid[j])
if p.has_key(j):
if self.loglevel > 0: print j,
znew.append(0.5*(grid[j] + grid[j+1]))
if self.loglevel > 0: print
znew.append(grid[-1])
if self.loglevel > 0 and nnew > 0:
print 'to resolve ',
ck = c.keys()
for s in ck:
if s <> ck[-1]:
print s+',',
else:
print s,
print
npts = len(znew)
newsoln = Numeric.zeros((npts, ncomp),'d')
for i in range(ncomp):
for j in range(npts):
newsoln[j,i] = interp.interp(znew[j],grid,solution[:,i])
return (Numeric.array(znew), Numeric.array(znew), newsoln, self.ok)
def refine(grid = None, solution = None, components = [], delta = (0.8, 1.0), threshold = None):
"""Refine a grid and interpolate the solution onto the new grid."""
r = Refiner(components = components, delta = delta)
return r.refine(grid, solution, threshold)
def prune(grid = None, solution = None, components = [], delta = (0.8, 1.0), threshold = None):
"""Remove unneeded points from a grid and solution array."""
r = Refiner(components = components, delta = delta)
return r.prune(grid, solution, threshold)
# test it
if __name__ == '__main__':
grid = Numeric.array([0.0, 0.2, 0.3, 1.0, 4.0])
soln = Numeric.array([[100.0, 0.4, -9.0],
[500.0, 0.0, -89.0],
[700.0, 0.9, 99.0],
[-99.0, 8.0, 77.0],
[567.0, 8.0, 0.0]])
grid_new, soln_new = refine(grid, soln,
components = [0,2],
delta = (0.5, 0.8))
print 'new grid = ',grid_new
print 'new solution = ',soln_new