From 754d098884f1afb450048463cb3ce07f4341bf10 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 9 Sep 2003 16:22:26 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/python/Cantera/Flow.py | 726 ------------------------- Cantera/python/Cantera/FlowBoundary.py | 140 ----- Cantera/python/Cantera/FlowPlotter.py | 65 --- Cantera/python/Cantera/OneDim.py | 476 ---------------- Cantera/python/Cantera/Thermo.py | 227 -------- Cantera/python/Cantera/boundaries1D.py | 56 -- Cantera/python/Cantera/ck2ctml.py | 21 - Cantera/python/Cantera/flame.py | 483 ---------------- Cantera/python/Cantera/pip.py | 31 -- Cantera/python/Cantera/refine.py | 295 ---------- 10 files changed, 2520 deletions(-) delete mode 100755 Cantera/python/Cantera/Flow.py delete mode 100755 Cantera/python/Cantera/FlowBoundary.py delete mode 100755 Cantera/python/Cantera/FlowPlotter.py delete mode 100755 Cantera/python/Cantera/OneDim.py delete mode 100755 Cantera/python/Cantera/Thermo.py delete mode 100755 Cantera/python/Cantera/boundaries1D.py delete mode 100644 Cantera/python/Cantera/ck2ctml.py delete mode 100755 Cantera/python/Cantera/flame.py delete mode 100644 Cantera/python/Cantera/pip.py delete mode 100755 Cantera/python/Cantera/refine.py diff --git a/Cantera/python/Cantera/Flow.py b/Cantera/python/Cantera/Flow.py deleted file mode 100755 index fdeac6e38..000000000 --- a/Cantera/python/Cantera/Flow.py +++ /dev/null @@ -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 -""" - diff --git a/Cantera/python/Cantera/FlowBoundary.py b/Cantera/python/Cantera/FlowBoundary.py deleted file mode 100755 index 27199d932..000000000 --- a/Cantera/python/Cantera/FlowBoundary.py +++ /dev/null @@ -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.") - diff --git a/Cantera/python/Cantera/FlowPlotter.py b/Cantera/python/Cantera/FlowPlotter.py deleted file mode 100755 index 673b36a3f..000000000 --- a/Cantera/python/Cantera/FlowPlotter.py +++ /dev/null @@ -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 - diff --git a/Cantera/python/Cantera/OneDim.py b/Cantera/python/Cantera/OneDim.py deleted file mode 100755 index 973d1ce56..000000000 --- a/Cantera/python/Cantera/OneDim.py +++ /dev/null @@ -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() - - diff --git a/Cantera/python/Cantera/Thermo.py b/Cantera/python/Cantera/Thermo.py deleted file mode 100755 index 3859ddd84..000000000 --- a/Cantera/python/Cantera/Thermo.py +++ /dev/null @@ -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 - - - - - - - - - - - diff --git a/Cantera/python/Cantera/boundaries1D.py b/Cantera/python/Cantera/boundaries1D.py deleted file mode 100755 index c75335383..000000000 --- a/Cantera/python/Cantera/boundaries1D.py +++ /dev/null @@ -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) - - diff --git a/Cantera/python/Cantera/ck2ctml.py b/Cantera/python/Cantera/ck2ctml.py deleted file mode 100644 index 17cb88c8d..000000000 --- a/Cantera/python/Cantera/ck2ctml.py +++ /dev/null @@ -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) - - - diff --git a/Cantera/python/Cantera/flame.py b/Cantera/python/Cantera/flame.py deleted file mode 100755 index 2365496c0..000000000 --- a/Cantera/python/Cantera/flame.py +++ /dev/null @@ -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() - - - - - - - - - - diff --git a/Cantera/python/Cantera/pip.py b/Cantera/python/Cantera/pip.py deleted file mode 100644 index 9146a7c56..000000000 --- a/Cantera/python/Cantera/pip.py +++ /dev/null @@ -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) - - - diff --git a/Cantera/python/Cantera/refine.py b/Cantera/python/Cantera/refine.py deleted file mode 100755 index f0aad3b89..000000000 --- a/Cantera/python/Cantera/refine.py +++ /dev/null @@ -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 - - - -