*** empty log message ***

This commit is contained in:
Dave Goodwin
2003-09-06 14:38:49 +00:00
parent 92e9bee752
commit 3f27dcc438
33 changed files with 494 additions and 219 deletions

View File

@@ -123,10 +123,9 @@ extern "C" {
catch (CanteraError) { return -1; }
}
int DLL_EXPORT domain_setBounds(int i, int nl, double* lower,
int nu, double* upper) {
int DLL_EXPORT domain_setBounds(int i, int n, double lower, double upper) {
try {
_domain(i)->setBounds(nl, lower, nu, upper);
_domain(i)->setBounds(n, lower, upper);
return 0;
}
catch (CanteraError) { return -1; }
@@ -146,10 +145,10 @@ extern "C" {
catch (CanteraError) { return -1.0; }
}
int DLL_EXPORT domain_setTolerances(int i, int nr, double* rtol,
int na, double* atol, int itime) {
int DLL_EXPORT domain_setTolerances(int i, int n, double rtol,
double atol, int itime) {
try {
_domain(i)->setTolerances(nr, rtol, na, atol, itime);
_domain(i)->setTolerances(n, rtol, atol, itime);
return 0;
}
catch (CanteraError) { return -1; }

View File

@@ -13,12 +13,12 @@ extern "C" {
int DLL_IMPORT domain_nPoints(int i);
int DLL_IMPORT domain_componentName(int i, int n, int sz, char* nameout);
int DLL_IMPORT domain_componentIndex(int i, char* name);
int DLL_IMPORT domain_setBounds(int i, int nl, double* lower,
int nu, double* upper);
int DLL_IMPORT domain_setBounds(int i, int n, double lower,
double upper);
double DLL_EXPORT domain_lowerBound(int i, int n);
double DLL_EXPORT domain_upperBound(int i, int n);
int DLL_IMPORT domain_setTolerances(int i, int nr, double* rtol,
int na, double* atol, int itime);
int DLL_IMPORT domain_setTolerances(int i, int n, double rtol,
double atol, int itime);
double DLL_IMPORT domain_rtol(int i, int n);
double DLL_IMPORT domain_atol(int i, int n);
int DLL_IMPORT domain_setupGrid(int i, int npts, double* grid);

View File

@@ -52,9 +52,9 @@ while length(property_argin) >= 2,
if sz == nComponents(a)
setTolerances(a, val(1,:), val(2,:));
elseif length(val) == 2
rt = val(1)*ones(1,nComponents(a));
at = val(2)*ones(1,nComponents(a));
setTolerances(a, rt, at);
%rt = val(1)*ones(1,nComponents(a));
%at = val(2)*ones(1,nComponents(a));
setTolerances(a, 'default', val(1), val(2));
else
error('wrong array size for error tolerances');
end
@@ -63,9 +63,9 @@ while length(property_argin) >= 2,
if sz == nComponents(a)
setTolerances(a, val(1,:), val(2,:));
elseif length(val) == 2
rt = val(1)*ones(1,nComponents(a));
at = val(2)*ones(1,nComponents(a));
setTolerances(a, rt, at, 'ts');
rt = val(1) %*ones(1,nComponents(a));
at = val(2) %*ones(1,nComponents(a));
setTolerances(a, 'default', rt, at, 'ts');
else
error('wrong array size for error tolerances');
end

View File

@@ -1,4 +1,5 @@
function d = setBounds(d, lower, upper)
function d = setBounds(d, component, lower, upper)
% SETBOUNDS -
%
domain_methods(d.dom_id, 51, lower, upper);
n = componentIndex(d,component);
domain_methods(d.dom_id, 51, n, lower, upper);

View File

@@ -1,17 +1,35 @@
function d = setTolerances(d, rtol, atol, typ)
function d = setTolerances(d, component, rtol, atol, typ)
% SETTOLERANCES -
%
ityp = 0;
if nargin == 4
if nargin == 5
switch typ
case 'ts'
itype = -1;
ityp = -1;
case 'time'
itype = -1;
ityp = -1;
case 'ss'
itype = 1;
ityp = 1;
case 'steady'
itype = 1;
ityp = 1;
end
end
domain_methods(d.dom_id, 52, rtol, atol, ityp);
if strcmp(component,'default')
nc = nComponents(d)
for ii = 1:nc
domain_methods(d.dom_id, 52, ii, rtol, atol, ityp);
end
return
end
if iscell(component)
nc = length(component);
for ii = 1:nc
n = componentIndex(d, component{ii});
domain_methods(d.dom_id, 52, n, rtol, atol, ityp);
end
else
n = componentIndex(d, component);
domain_methods(d.dom_id, 52, n, rtol, atol, ityp);
end

View File

@@ -5,5 +5,10 @@ mex private/ctmethods.cpp private/ctfunctions.cpp ...
private/thermomethods.cpp private/kineticsmethods.cpp ...
private/transportmethods.cpp private/reactormethods.cpp ...
private/wallmethods.cpp private/flowdevicemethods.cpp ...
../../../../../lib/cantera14.lib
private/onedimmethods.cpp private/surfmethods.cpp private/write.cpp ...
../../../../../lib/clib.lib ../../../../../lib/oneD.lib ...
../../../../../lib/zeroD.lib ../../../../../lib/transport.lib ...
../../../../../lib/cantera.lib ../../../../../lib/recipes.lib ...
../../../../../lib/cvode.lib ../../../../../lib/ctlapack.lib ...
../../../../../lib/ctmath.lib ../../../../../lib/ctblas.lib
disp('done.');

View File

@@ -52,7 +52,7 @@ teq = temperature(gas);
yeq = massFractions(gas);
rhoeq = density(gas);
z1 = 0.5;
z1 = 0.2;
mdot0 = massFlux(left);
mdot1 = massFlux(right);
t0 = temperature(left);

View File

@@ -83,10 +83,8 @@ fl = flame(gas, burner, f, s);
% solution, uncomment this line and edit the file name and solution id.
%restore(fl,'h2flame2.xml', 'energy')
solve(fl, 1, refine_grid);
%%%%%%%%%%%% enable the energy equation %%%%%%%%%%%%%%%%%%%%%
%
% The energy equation will now be solved to compute the

View File

@@ -208,30 +208,28 @@ void onedimmethods( int nlhs, mxArray *plhs[],
else {
int iok = -1;
double *lower, *upper, *rtol, *atol, *grid, *pos, *values,
double lower, upper, rtol, atol, *grid, *pos, *values,
mdot, t, p, val, *temp, ratio, slope, curve, tstep, *dts,
rdt, prune;
int nlower, nupper, nr, na, npts, np, comp, localPoint, idom,
int nlower, nupper, nr, na, npts, np, nv, comp, localPoint, idom,
loglevel, refine_grid, n, flag, itime, ns, *nsteps, icount,
onoff, ss_age, ts_age;
char *xstr, *fname, *id, *desc, *name;
switch (job) {
case 51:
checkNArgs(5, nrhs);
lower = mxGetPr(prhs[3]);
nlower = mxGetM(prhs[3]) * mxGetN(prhs[3]);
upper = mxGetPr(prhs[4]);
nupper = mxGetM(prhs[4]) * mxGetN(prhs[4]);
iok = domain_setBounds(dom, nlower, lower, nupper, upper);
checkNArgs(6, nrhs);
n = getInt(prhs[3]) - 1;
lower = getDouble(prhs[4]);
upper = getDouble(prhs[5]);
iok = domain_setBounds(dom, n, lower, upper);
break;
case 52:
checkNArgs(6, nrhs);
rtol = mxGetPr(prhs[3]);
nr = mxGetM(prhs[3]) * mxGetN(prhs[3]);
atol = mxGetPr(prhs[4]);
na = mxGetM(prhs[4]) * mxGetN(prhs[4]);
itime = getInt(prhs[5]);
iok = domain_setTolerances(dom, nr, rtol, na, atol, itime);
checkNArgs(7, nrhs);
n = getInt(prhs[3]) - 1;
rtol = getDouble(prhs[4]);
atol = getDouble(prhs[5]);
itime = getInt(prhs[6]);
iok = domain_setTolerances(dom, n, rtol, atol, itime);
break;
case 53:
checkNArgs(4, nrhs);
@@ -269,7 +267,8 @@ void onedimmethods( int nlhs, mxArray *plhs[],
pos = mxGetPr(prhs[3]);
temp = mxGetPr(prhs[4]);
n = mxGetM(prhs[3])*mxGetN(prhs[3]);
iok = stflow_setFixedTempProfile(dom, n, pos, temp);
m = mxGetM(prhs[4])*mxGetN(prhs[4]);
iok = stflow_setFixedTempProfile(dom, n, pos, m, temp);
break;
case 65:
checkNArgs(4, nrhs);
@@ -297,7 +296,8 @@ void onedimmethods( int nlhs, mxArray *plhs[],
pos = mxGetPr(prhs[5]);
values = mxGetPr(prhs[6]);
np = mxGetM(prhs[5])*mxGetN(prhs[5]);
iok = sim1D_setProfile(dom, idom, comp, np, pos, values);
nv = mxGetM(prhs[6])*mxGetN(prhs[6]);
iok = sim1D_setProfile(dom, idom, comp, np, pos, nv, values);
break;
case 102:
checkNArgs(6, nrhs);

View File

@@ -2,5 +2,4 @@
The classes in this package implement one-dimensional reacting flow problems.
"""
from onedim import *
import onedim

View File

@@ -1,14 +1,97 @@
from onedim import *
import Numeric
class BurnerFlame(Stack):
"""A burner-stabilized flat flame."""
def __init__(self, gas = None, type = 'burner'):
self.type = type
self.left = Inlet('burner')
self.right = Outlet('outlet')
self.flow = AxisymmetricFlow('flow',gas = gas)
Stack.__init__(self, [self.left, self.flow, self.right])
self.showSolution()
def __init__(self, gas = None, burner = None, outlet = None, grid = None):
if burner:
self.burner = burner
else:
self.burner = Inlet('burner')
self.gas = gas
self.burner.set(temperature = gas.temperature())
if outlet:
self.outlet = outlet
else:
self.outlet = Outlet('outlet')
self.pressure = gas.pressure()
self.flame = AxisymmetricFlow('flame',gas = gas)
self.flame.setupGrid(grid)
Stack.__init__(self, [self.burner, self.flame, self.outlet])
self.setRefineCriteria()
self._initialized = 0
def init(self):
"""Set the initial guess for the solution."""
gas = self.gas
nsp = gas.nSpecies()
yin = Numeric.zeros(nsp, 'd')
for k in range(nsp):
yin[k] = self.burner.massFraction(k)
gas.setState_TPY(self.burner.temperature(), self.pressure, yin)
u0 = self.burner.mdot()/gas.density()
t0 = self.burner.temperature()
# get adiabatic flame temperature and composition
gas.equilibrate('HP')
teq = gas.temperature()
yeq = gas.massFractions()
u1 = self.burner.mdot()/gas.density()
z1 = 0.2
locs = Numeric.array([0.0, z1, 1.0],'d')
self.setProfile('u', locs, [u0, u1, u1])
self.setProfile('T', locs, [t0, teq, teq])
for n in range(nsp):
self.setProfile(gas.speciesName(n), locs, [yin[n], yeq[n], yeq[n]])
self._initialized = 1
def solve(self, loglevel = 1, refine_grid = 1):
if not self._initialized: self.init()
Stack.solve(self, loglevel = loglevel, refine_grid = refine_grid)
def setRefineCriteria(self, ratio = 10.0, slope = 0.8, curve = 0.8, prune = 0.0):
Stack.setRefineCriteria(self, domain = self.flame, ratio = ratio, slope = slope, curve = curve,
prune = prune)
def setProfile(self, component, locs, vals):
self._initialized = 1
Stack.setProfile(self, self.flame, component, locs, vals)
def set(self, tol = None, energy = '', tol_time = None):
if tol:
self.flame.setTolerances(default = tol)
if tol_time:
self.flame.setTolerances(default = tol_time, time = 1)
if energy:
self.flame.set(energy = energy)
def T(self, point = -1):
return self.solution('T', point)
def u(self, point = -1):
return self.solution('u', point)
def V(self, point = -1):
return self.solution('V', point)
def solution(self, component = '', point = -1):
if point >= 0: return self.value(self.flame, component, point)
else: return self.profile(self.flame, component)
def setGasState(self, j):
nsp = self.gas.nSpecies()
y = Numeric.zeros(nsp, 'd')
for n in range(nsp):
nm = self.gas.speciesName(n)
y[n] = self.solution(nm, j)
self.gas.setState_TPY(self.T(j), self.flame.pressure(), y)

View File

@@ -2,6 +2,8 @@ from Cantera import *
from Cantera import _cantera
import Numeric
_onoff = {'on':1, 'yes':1, 'off':0, 'no':0}
class Domain1D:
"""One-dimensional domains."""
@@ -16,11 +18,11 @@ class Domain1D:
return self._hndl
def type(self):
"""Domain type."""
"""Domain type. Integer."""
return _cantera.domain_type(self._hndl)
def index(self):
"""Index of this domain in a stack."""
"""Index of this domain in a stack. Returns -1 if this domain is not part of a stack."""
return _cantera.domain_index(self._hndl)
def nComponents(self):
@@ -36,18 +38,43 @@ class Domain1D:
return _cantera.domain_componentName(self._hndl, n)
def componentIndex(self, name):
"""Index of the component named 'name'"""
"""Index of the component with name 'name'"""
return _cantera.domain_componentIndex(self._hndl, name)
def setBounds(self, lower, upper):
"""Set the lower and upper bounds on the solution."""
return _cantera.domain_setBounds(self._hndl,
Numeric.asarray(lower),
Numeric.asarray(upper))
def setBounds(self, **bounds):
"""Set the lower and upper bounds on the solution.
The argument list should consist of keyword/value pairs, with
component names as keywords and (lower_bound, upper_bound)
tuples as the values. The keyword 'default' may be used to
specify default bounds for all unspecified components. The
keyword 'Y' can be used to stand for all species mass
fractions in flow domains. """
d = {}
if bounds.has_key('default'):
for n in range(self.nComponents()):
d[self.componentName(n)] = bounds['default']
del bounds['default']
for b in bounds.keys():
if b == 'Y':
if self.type >= 50:
nc = self.nComponents()
for n in range(4, nc):
d[self.componentName(n)] = bounds[b]
else:
raise CanteraError('Y can only be specified in flow domains.')
else:
d[b] = bounds[b]
for b in d.keys():
n = self.componentIndex(b)
_cantera.domain_setBounds(self._hndl, n, d[b][0], d[b][1])
def bounds(self, component):
ic = self.componentIndex(component)
lower = _cantera.domain_lowerBound(self._hndl, ic)
upper = _cantera.domain_upperBound(self._hndl, ic)
upper = _cantera.domain_upperBound(self._hndl, ic)
return (lower, upper)
def tolerances(self, component):
@@ -56,21 +83,48 @@ class Domain1D:
a = _cantera.domain_atol(self._hndl, ic)
return (r, a)
def setTolerances(self, rtol, atol, time=0):
def setTolerances(self, **tol):
"""Set the error tolerances. If 'time' is present and
non-zero, then the values entered will apply to the transient
problem. Otherwise, they will apply to the steady-state
problem. """
return _cantera.domain_setTolerances(self._hndl,
Numeric.asarray(rtol),
Numeric.asarray(atol), itime)
problem.
The argument list should consist of keyword/value pairs, with
component names as keywords and (rtol, atol) tuples as the
values. The keyword 'default' may be used to specify default
bounds for all unspecified components. The keyword 'Y' can be
used to stand for all species mass fractions in flow domains.
"""
d = {}
if tol.has_key('default'):
for n in range(self.nComponents()):
d[self.componentName(n)] = tol['default']
del tol['default']
itime = 0
for b in tol.keys():
if b == 'time': itime = -1
elif b == 'steady': itime = 1
elif b == 'Y':
if self.type >= 50:
nc = self.nComponents()
for n in range(4, nc):
d[self.componentName(n)] = tol[b]
else:
raise CanteraError('Y can only be specified in flow domains.')
else:
d[b] = tol[b]
for b in d.keys():
n = self.componentIndex(b)
_cantera.domain_setTolerances(self._hndl, n, d[b][0], d[b][1], itime)
def setupGrid(self, grid):
"""Specify the grid."""
return _cantera.domain_setupGrid(self._hndl, Numeric.asarray(grid))
def setID(self, id):
print 'id = ',id
return _cantera.domain_setID(self._hndl, id)
def setDesc(self, desc):
@@ -99,26 +153,25 @@ class Domain1D:
else:
raise CanteraError('unknown attribute: '+opt)
def _dict2array(self, d):
a = zeros(self.nComponents(),'d')
def _dict2arrays(self, d = None, array1 = None, array2 = None):
nc = self.nComponents()
if d.has_key('default'):
a += d['default']
for k in d.keys():
a[self.componentIndex(k)] = d[k]
print a
return a
def _dict2arrays(self, d):
a1 = zeros(self.nComponents(),'d')
a2 = zeros(self.nComponents(),'d')
if d.has_key('default'):
a1 += d['default'][0]
a2 += d['default'][1]
a1 = zeros(nc,'d') + d['default'][0]
a2 = zeros(nc,'d') + d['default'][1]
del d['default']
else:
if array1: a1 = array(array1)
else: a1 = zeros(nc,'d')
if array2: a2 = array(array2)
else: a2 = zeros(nc,'d')
for k in d.keys():
a1[self.componentIndex(k)] = d[k][0]
a2[self.componentIndex(k)] = d[k][1]
print a1, a2
c = self.componentIndex(k)
if c >= 0:
a1[self.componentIndex(k)] = d[k][0]
a2[self.componentIndex(k)] = d[k][1]
else:
raise CanteraError('unknown component '+k)
return (a1, a2)
@@ -149,14 +202,16 @@ class Bdry1D(Domain1D):
def set(self, **options):
for opt in options.keys():
v = options[opt]
if opt == 'mdot':
if opt == 'mdot' or opt == 'massflux':
self.setMdot(v)
elif opt == 'temperature':
del options[opt]
elif opt == 'temperature' or opt == 'T':
self.setTemperature(v)
elif opt == 'mole_fractions':
del options[opt]
elif opt == 'mole_fractions' or opt == 'X':
self.setMoleFractions(v)
else:
self._set(options)
del options[opt]
self._set(options)
class Inlet(Bdry1D):
@@ -201,35 +256,83 @@ class Surface(Bdry1D):
class AxisymmetricFlow(Domain1D):
"""An axisymmetric flow"""
def __init__(self, id = 'axisymmetric_flow', gas = 'None'):
"""An axisymmetric flow domain.
In an axisymmetric flow domain, the equations solved are the
similarity equations for the flow in a finite-height gap of
infinite radial extent. The solution variables are
u -- axial velocity
V -- radial velocity divided by radius
T -- temperature
lambda -- (1/r)(dP/dr)
Y_k -- species mass fractions
It may be shown that if the boundary conditions on these variables
are independent of radius, then a similarity solution to the exact
governing equations exists in which these variables are all
independent of radius. This solution holds only in in
low-Mach-number limit, in which case (dP/dz) = 0, and lambda is a
constant. (Lambda is treated as a spatially-varying solution
variable for numerical reasons, but in the final solution it is
always independent of z.) As implemented here, the governing
equations assume an ideal gas mixture. Arbitrary chemistry is
allowed, as well as arbitrary variation of the transport
properties.
"""
def __init__(self, id = 'axisymmetric_flow', gas = None):
Domain1D.__init__(self)
iph = gas.thermo_hndl()
ikin = gas.kinetics_hndl()
itr = gas.transport_hndl()
self._hndl = _cantera.stflow_new(iph, ikin, itr)
if id: self.setID(id)
self._p = -1.0
self.setPressure(gas.pressure())
self.solveEnergyEqn()
def setPressure(self, p):
"""Set the pressure [Pa]. The pressure is a constant, since
the governing equations are those appropriate for the
low-Mach-number limit."""
the governing equations are those for the low-Mach-number limit."""
_cantera.stflow_setPressure(self._hndl, p)
self._p = p
def pressure(self):
return self._p
def setFixedTempProfile(self, temp):
"""Set the fixed temperature profile.
This profile is used whenever the energy equation is disabled.
"""
"""Set the fixed temperature profile. This profile is used
whenever the energy equation is disabled. """
return _cantera.stflow_setFixedTempProfile(self._hndl, temp)
def solveSpeciesEqs(self, flag):
def solveSpeciesEqs(self, flag = 1):
"""Enable or disable solving the species equations. If invoked
with no arguments or with a non-zero argument, the species
equations will be solved. If invoked with a zero argument,
they will not be, and instead the species profiles will be
held at their initial values. Default: species equations
enabled."""
return _cantera.stflow_solveSpeciesEqs(self._hndl, flag)
def solveEnergyEqn(self, flag):
def solveEnergyEqn(self, flag = 1):
"""Enable or disable solving the energy equation. If invoked
with no arguments or with a non-zero argument, the energy
equations will be solved. If invoked with a zero argument,
it will not be, and instead the temperature profiles will be
held to the one specified by the call to setFixedTempProfile.
Default: energy equation enabled."""
return _cantera.stflow_solveEnergyEqn(self._hndl, flag)
def set(self, **opt):
for o in opt.keys():
v = opt[o]
if o == 'P' or o == 'pressure':
self.setPressure(v)
del opt[o]
elif o == 'energy':
self.solveEnergyEqn(flag = _onoff[v])
else:
self._set(opt)
class Stack:
@@ -245,6 +348,7 @@ class Stack:
for n in range(nd):
hndls[n] = domains[n].domain_hndl()
self._hndl = _cantera.sim1D_new(hndls)
self._domains = domains
def __del__(self):
_cantera.sim1D_del(self._hndl)
@@ -278,27 +382,35 @@ class Stack:
def refine(self, loglevel=1):
return _cantera.sim1D_refine(self._hndl, loglevel)
def setRefineCriteria(self, dom, ratio = 10.0, slope = 0.8,
def setRefineCriteria(self, domain = None, ratio = 10.0, slope = 0.8,
curve = 0.8, prune = 0.05):
idom = dom.index()
idom = domain.index()
return _cantera.sim1D_setRefineCriteria(self._hndl,
idom, ratio, slope, curve, prune)
def save(self, fname, id, desc):
return _cantera.sim1D_save(self._hndl, fname, id, desc)
def save(self, file = 'soln.xml', id = 'solution', desc = 'none'):
return _cantera.sim1D_save(self._hndl, file, id, desc)
def restore(self, fname, id):
return _cantera.sim1D_restore(self._hndl, fname, id)
def restore(self, file = 'soln.xml', id = 'solution'):
return _cantera.sim1D_restore(self._hndl, file, id)
def writeStats(self):
def showStats(self):
return _cantera.sim1D_writeStats(self._hndl)
def domainIndex(self, name):
return _cantera.sim1D_domainIndex(self._hndl, name)
def value(self, dom, icomp, localPoint):
idom = dom.index()
def value(self, domain, component, localPoint):
icomp = domain.componentIndex(component)
idom = domain.index()
return _cantera.sim1D_value(self._hndl, idom, icomp, localPoint)
def profile(self, domain, component):
np = domain.nPoints()
x = zeros(np,'d')
for n in range(np):
x[n] = self.value(domain, component, n)
return x
def workValue(self, dom, icomp, localPoint):
idom = dom.index()
return _cantera.sim1D_workValue(self._hndl, idom, icomp, localPoint)

View File

@@ -120,22 +120,13 @@ py_domain_setBounds(PyObject *self, PyObject *args)
{
int _val;
int i;
PyObject* lower;
PyObject* upper;
if (!PyArg_ParseTuple(args, "iOO:domain_setBounds", &i, &lower, &upper))
int n;
double lower;
double upper;
if (!PyArg_ParseTuple(args, "iidd:domain_setBounds", &i, &n, &lower, &upper))
return NULL;
PyArrayObject* lower_array = (PyArrayObject*)lower;
double* lower_data = (double*)lower_array->data;
int lower_len = lower_array->dimensions[0];
PyArrayObject* upper_array = (PyArrayObject*)upper;
double* upper_data = (double*)upper_array->data;
int upper_len = upper_array->dimensions[0];
_val = domain_setBounds(i,lower_len,lower_data,upper_len,upper_data);
_val = domain_setBounds(i,n,lower,upper);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("i",_val);
}
@@ -176,23 +167,13 @@ py_domain_setTolerances(PyObject *self, PyObject *args)
{
int _val;
int i;
PyObject* rtol;
PyObject* atol;
int n;
double rtol, atol;
int itime;
if (!PyArg_ParseTuple(args, "iOOi:domain_setTolerances", &i, &rtol, &atol, &itime))
if (!PyArg_ParseTuple(args, "iiddi:domain_setTolerances", &i, &n, &rtol, &atol, &itime))
return NULL;
PyArrayObject* rtol_array = (PyArrayObject*)rtol;
double* rtol_data = (double*)rtol_array->data;
int rtol_len = rtol_array->dimensions[0];
PyArrayObject* atol_array = (PyArrayObject*)atol;
double* atol_data = (double*)atol_array->data;
int atol_len = atol_array->dimensions[0];
_val = domain_setTolerances(i,rtol_len,rtol_data,atol_len,atol_data,itime);
_val = domain_setTolerances(i,n,rtol,atol,itime);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("i",_val);
}

View File

@@ -17,7 +17,7 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
# Temporarily removed 'filter.o', since it was causing a compile error on Mac OS X.
OBJS = atomicWeightDB.o CKParser.o CKReader.o Reaction.o ckr_utils.o \
thermoFunctions.o writelog.o ck2ct.o ck2ctml.o
thermoFunctions.o writelog.o ck2ct.o
CXX_INCLUDES = -I. -I..
CONV_LIB = @buildlib@/libconverters.a

View File

@@ -12,7 +12,12 @@
#include <iostream>
#include <string>
// on cygwin, #include <sstream> doesn't work
#if defined(__CYGWIN__)
#undef USE_STRINGSTREAM
#else
#define USE_STRINGSTREAM
#endif
#ifdef USE_STRINGSTREAM
#include <sstream>

View File

@@ -44,7 +44,11 @@ namespace ctml {
f << "from Cantera import *\n";
f.close();
int ierr = 0;
string cmd = pypath() + " " + path + " > " + tmpDir() + "/log";
#ifdef WIN32
string cmd = "cmd /C "+pypath() + " " + path + ">> log 2>&1";
#else
string cmd = pypath() + " " + path + " &> " + tmpDir() + "/log";
#endif
try {
ierr = system(cmd.c_str());
if (ierr != 0) {
@@ -85,8 +89,12 @@ namespace ctml {
<< "execfile(file)\n"
<< "write()\n";
f.close();
string cmd = pypath() + " " + path + " > ct2ctml.log";
int ierr;
#ifdef WIN32
string cmd = "cmd /C "+pypath() + " " + path + ">> ct2ctml.log 2>&1";
#else
string cmd = pypath() + " " + path + " &> ct2ctml.log";
#endif
int ierr = 0;
try {
ierr = system(cmd.c_str());
}

View File

@@ -75,8 +75,7 @@ namespace ctml {
#ifdef CTML_VERSION_1_4
XML_Node& f = node.addChild("float",val,fmt);
f.addAttribute("title",title);
#endif
#ifdef CTML_VERSION_1_4_1
#else
XML_Node& f = node.addChild(title,val,fmt);
#endif
if (type != "") f.addAttribute("type",type);

View File

@@ -152,6 +152,11 @@ namespace Cantera {
copy(lower, lower + m_nv, m_min.begin());
}
void setBounds(int n, doublereal lower, doublereal upper) {
m_min[n] = lower;
m_max[n] = upper;
}
void setTolerances(int nr, const doublereal* rtol,
int na, const doublereal* atol, int ts = 0) {
if (nr < m_nv || na < m_nv)
@@ -168,6 +173,17 @@ namespace Cantera {
}
}
void setTolerances(int n, doublereal rtol, doublereal atol, int ts = 0) {
if (ts >= 0) {
m_rtol_ss[n] = rtol;
m_atol_ss[n] = atol;
}
if (ts <= 0) {
m_rtol_ts[n] = rtol;
m_atol_ts[n] = atol;
}
}
/// Relative tolerance of the nth component.
doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); }

View File

@@ -144,9 +144,9 @@ namespace Cantera {
virtual void showSolution(const doublereal* x) {
char buf[80];
sprintf(buf, " Mass Flux: %10.4g kg/m^2/s \n", x[0]);
sprintf(buf, " Mass Flux: %10.4g kg/m^2/s \n", m_mdot);
writelog(buf);
sprintf(buf, " Temperature: %10.4g K \n", x[1]);
sprintf(buf, " Temperature: %10.4g K \n", m_temp);
writelog(buf);
if (m_flow) {
writelog(" Mass Fractions: \n");
@@ -297,7 +297,7 @@ namespace Cantera {
virtual void showSolution(const doublereal* x) {
char buf[80];
sprintf(buf, " Temperature: %10.4g K \n", x[0]);
sprintf(buf, " Temperature: %10.4g K \n", m_temp);
writelog(buf);
writelog("\n");
}
@@ -354,7 +354,7 @@ namespace Cantera {
virtual void showSolution(const doublereal* x) {
char buf[80];
sprintf(buf, " Temperature: %10.4g K \n", x[0]);
sprintf(buf, " Temperature: %10.4g K \n", m_temp);
writelog(buf);
writelog(" Coverages: \n");
for (int k = 0; k < m_nsp; k++) {

View File

@@ -1201,9 +1201,10 @@ namespace Cantera {
writelog(nd["title"]+": "+nd.value()+"\n");
}
map<string, double> params;
getFloats(dom, params);
setPressure(params["pressure"]);
//map<string, double> params;
double pp = -1.0;
pp = getFloat(dom, "pressure", "pressure");
setPressure(pp);
vector<XML_Node*> d;

View File

@@ -23,8 +23,9 @@ namespace Cantera {
m_left_loc(0), m_right_loc(0),
m_left_points(0), m_nv(0),
m_left_nsp(0), m_right_nsp(0),
m_sp_left(0), m_sp_right(0),
m_start_left(0), m_start_right(0),
m_phase_left(0), m_phase_right(0), m_mdot(0.0) {
m_phase_left(0), m_phase_right(0), m_temp(0.0), m_mdot(0.0) {
m_type = cConnectorType;
}
@@ -220,16 +221,16 @@ namespace Cantera {
inlt.addAttribute("type","inlet");
inlt.addAttribute("components",nComponents());
for (int k = 0; k < nComponents(); k++) {
ctml::addFloat(inlt, componentName(k), s[k], "", "",0.0, 1.0);
ctml::addFloat(inlt, componentName(k), s[k], "", "",lowerBound(k), upperBound(k));
}
}
void Inlet1D::
restore(XML_Node& dom, doublereal* soln) {
map<string, double> x;
getFloats(dom, x);
soln[0] = x["mdot"];
soln[1] = x["temperature"];
//map<string, double> x;
//getFloats(dom, x);
soln[0] = getFloat(dom, "mdot", "massflowrate"); // x["mdot"];
soln[1] = getFloat(dom, "temperature", "temperature"); // x["temperature"];
resize(2,1);
}

View File

@@ -43,10 +43,10 @@ namespace Cantera {
for (j = 0; j < np; j++) {
val = x[index(m,j)];
if (loglevel > 0) {
if (val > above + Tiny || val < below - Tiny) {
if (val > above + 1.0e-12 || val < below - 1.0e-12) {
sprintf(buf, "domain %d: %20s(%d) = %10.3e (%10.3e, %10.3e)\n",
r.domainIndex(), r.componentName(m).c_str(), j, val, below, above);
writelog(string("ERROR: solution out of bounds.\n")+buf);
writelog(string("\nERROR: solution out of bounds.\n")+buf);
}
}

View File

@@ -47,6 +47,12 @@ kernel-install:
cd Cantera/clib/src; @MAKE@ install
ranlib @prefix@/lib/cantera/*.a
win-kernel-install:
@INSTALL@ -d @prefix@/lib/cantera
rm -f @prefix@/lib/cantera/*
@INSTALL@ -m 644 @buildlib@/*.lib @prefix@/lib/cantera
cd Cantera/clib/src; @MAKE@ install
data-install:
@INSTALL@ -d @prefix@/cantera/data
@INSTALL@ -m 644 data/inputs/elements.xml @prefix@/cantera/data

View File

@@ -311,15 +311,19 @@ fi
SHARED_CTLIB=0
OS_IS_DARWIN=0
OS_IS_WIN=0
mex_ext=mexglx
dnl Checks for libraries.
AC_F77_LIBRARY_LDFLAGS()
case $ac_sys_system in
Darwin*) FLIBS='-lg2c -lgcc'; SHARED_CTLIB=0; OS_IS_DARWIN=1; mex_ext=mexmac;
Darwin*) FLIBS='-lg2c -lgcc'; SHARED_CTLIB=0; OS_IS_DARWIN=1; mex_ext=mexmac;;
CYGWIN*) OS_IS_WIN=1; mex_ext=dll;;
esac
AC_SUBST(FLIBS)
AC_SUBST(OS_IS_DARWIN)
esac
AC_SUBST(OS_IS_WIN)
AC_SUBST(SHARED_CTLIB)
AC_SUBST(mex_ext)

View File

@@ -14,15 +14,7 @@ namespace Cantera {
{
public:
GRI30() : m_ok(false), m_r(0) {
string path = findInputFile("gri30.xml");
ifstream fin(path.c_str());
if (!fin) {
throw CanteraError("GRI30","could not open "
+path+" for reading.");
}
m_r = new XML_Node("-");
m_r->build(fin);
m_r = get_XML_File("gri30.xml");
m_ok = buildSolutionFromXML(*m_r, "gri30", "phase", this, this);
if (!m_ok) throw CanteraError("GRI30",
"buildSolutionFromXML returned false");

View File

@@ -15,18 +15,11 @@ namespace Cantera {
public:
IdealGasMix(string infile, string id="") : m_ok(false), m_r(0) {
//string path = findInputFile(infile);
// ifstream fin(path.c_str());
// if (!fin) {
// throw CanteraError("IdealGasMix","could not open "
// +path+" for reading.");
// }
m_r = new XML_Node("-");
get_CTML_Tree(m_r, infile); // build(fin);
m_ok = buildSolutionFromXML(*m_r, id, "phase", this, this);
if (!m_ok) throw CanteraError("IdealGasMix",
"buildSolutionFromXML returned false");
m_r = get_XML_File(infile);
if (id == "-") id = "";
m_ok = buildSolutionFromXML(*m_r, id, "phase", this, this);
if (!m_ok) throw CanteraError("IdealGasMix",
"buildSolutionFromXML returned false");
}

View File

@@ -8,6 +8,10 @@ if prefix == '-':
if prefix == '/usr/local':
localinst = 0
icyg = prefix.find('/cygdrive/')
if icyg == 0:
prefix = prefix[10]+':'+prefix[11:]
bindir = prefix+'/bin'
libdir = prefix+'/lib/cantera'
hdrdir = prefix+'/include/cantera'

View File

@@ -12,3 +12,7 @@ make a copy with a different name, open the copy in Visual Studio,
delete file demo.f and replace it with your Fortran program. If you
need to change demo_ftnlib.cpp, then also replace it with your
modified version.
Note that the include and library directories are specified relative
to the location of the project file. If you move it to another
directory, you will need to update these locations under 'Settings'.

View File

@@ -1,7 +1,8 @@
c
c Replace this sample main program with your program
c
c This program uses functions defined in demo_ftnlib.cpp.
c This program uses functions defined in demo_ftnlib.cpp to create
c an ideal gas mixture and print some its properties.
c
program demo
implicit double precision (a-h,o-z)
@@ -10,18 +11,27 @@ c
double precision x(MAXSP), y(MAXSP), wdot(MAXSP)
character*80 eq
c
call newIdealGasMix('h2o2.xml','')
write(*,*) '**** Fortran 77 Test Program ****'
call newIdealGasMix('h2o2.cti','ohmech')
t = 1200.0
p = 101325.0
call setState_TPX_String(t, p,
$ 'H2:2, O2:1, OH:0.01, H:0.01, O:0.01')
call equilibrate('HP')
write(*,*) '**** Fortran 77 Test Program ****'
c
write(*,*) 'Initial state properties:'
write(*,10) temperature(), pressure(), density(),
$ enthalpy_mole(), entropy_mole(), cp_mole()
c compute the equilibrium state holding the specific
c enthalpy and pressure constant
call equilibrate('HP')
write(*,*) 'Equilibrium state properties:'
write(*,10) temperature(), pressure(), density(),
$ enthalpy_mole(), entropy_mole(), cp_mole()
10 format(//'Temperature: ',g14.5,' K'/
$ 'Pressure: ',g14.5,' Pa'/
$ 'Density: ',g14.5,' kg/m3'/
@@ -34,14 +44,22 @@ c
c Reaction information
c
irxns = nReactions()
c forward and reverse rates of progress should be equal
c in equilibrium states
call getFwdRatesOfProgress(qf)
call getRevRatesOfProgress(qr)
c net rates of progress should be zero in equilibrium states
call getNetRatesOfProgress(q)
c for each reaction, print the equation and the rates of progress
do i = 1,irxns
call getReactionEqn(i,eq)
write(*,20) eq,qf(i),qr(i),q(i)
20 format(a20,3g14.5,' kmol/m3/s')
20 format(a27,3e14.5,' kmol/m3/s')
end do
stop
end

View File

@@ -30,6 +30,11 @@ static IdealGasMix* _gas = 0;
// provides access to the pointer for functions in other libraries
IdealGasMix* _gasptr() { return _gas; }
// error handler
void handleError() {
showErrors(cout);
exit(-1);
}
// extern "C" turns off C++ name-mangling, so that the procedure names
// in the object file are exactly as shown here.
@@ -41,19 +46,23 @@ extern "C" {
/**
* Read in a reaction mechanism file and create an IdealGasMix
* object. The file may be in Chemkin-compatible format or in
* CTML. The name of a thermodynamic database may be supplied as a
* second argument. If none is required, enter an empty string as
* the second argument.
* object. The file may be in Cantera input format or in CTML. (If
* you have a file in Chemkin-compatible format, use utility
* program ck2cti first to convert it into Cantera format.)
*/
void newidealgasmix_(char* file, char* thermo,
ftnlen lenfile, ftnlen lenthermo) {
string fin = string(file, lenfile);
string fth = string(thermo, lenthermo);
if (_gas) delete _gas;
_gas = new IdealGasMix(fin, fth);
void newidealgasmix_(char* file, char* id,
ftnlen lenfile, ftnlen lenid) {
try {
string fin = string(file, lenfile);
string fth = string(id, lenid);
if (_gas) delete _gas;
_gas = new IdealGasMix(fin, fth);
}
catch (CanteraError) {
handleError();
}
}
/// integer function nElements()
integer nelements_() { return _gas->nElements(); }
@@ -68,17 +77,26 @@ extern "C" {
// subroutine setState_TPX(T, P, X)
void setstate_tpx_(doublereal* T, doublereal* P, doublereal* X) {
_gas->setState_TPX(*T, *P, X);
try {
_gas->setState_TPX(*T, *P, X);
}
catch (CanteraError) { handleError(); }
}
/// subroutine setState_TPX_String(T, P, X)
void setstate_tpx_string_(doublereal* T, doublereal* P,
char* X, ftnlen lenx) {
_gas->setState_TPX(*T, *P, string(X, lenx));
try {
_gas->setState_TPX(*T, *P, string(X, lenx));
}
catch (CanteraError) { handleError(); }
}
void setstate_try_(doublereal* T, doublereal* rho, doublereal* Y) {
_gas->setState_TRY(*T, *rho, Y);
try {
_gas->setState_TRY(*T, *rho, Y);
}
catch (CanteraError) { handleError(); }
}
//-------------- thermodynamic properties ----------------------
@@ -148,13 +166,24 @@ extern "C" {
return _gas->gibbs_mass();
}
void gotmolefractions_(doublereal* x) {
_gas->getMoleFractions(x);
}
void gotmassfractions_(doublereal* y) {
_gas->getMassFractions(y);
}
void equilibrate_(char* opt, ftnlen lenopt) {
if (lenopt != 2) {
throw CanteraError("equilibrate",
"two-character string required.");
try {
if (lenopt != 2) {
throw CanteraError("equilibrate",
"two-character string required.");
}
string optstr = string(opt, 2);
equilibrate(*_gas, optstr.c_str());
}
string optstr = string(opt, 2);
equilibrate(*_gas, optstr.c_str());
catch (CanteraError) { handleError(); }
}
@@ -205,11 +234,11 @@ int main() {
}
catch (CanteraError) {
showErrors(cerr);
return -1;
exit(-1);
}
catch (...) {
cout << "An exception was trapped. Program terminating." << endl;
return -1;
exit(-1);
}
}

View File

@@ -20,11 +20,10 @@ DEPENDS = $(OBJS:.o=.d)
.cpp.o:
@CXX@ -c $< @DEFS@ -I$(INCDIR) @CXXFLAGS@ $(CXX_FLAGS) -O0
all: csvdiff
cp csvdiff ../../bin
all: $(BINDIR)/csvdiff
csvdiff: mdp_allo.o csvdiff.o tok_input_util.o
@CXX@ -o csvdiff mdp_allo.o csvdiff.o tok_input_util.o \
$(BINDIR)/csvdiff: mdp_allo.o csvdiff.o tok_input_util.o
@CXX@ -o $(BINDIR)/csvdiff mdp_allo.o csvdiff.o tok_input_util.o \
$(LCXX_FLAGS) $(LCXX_END_LIBS)
clean:

View File

@@ -38,9 +38,9 @@
#include <limits.h>
#include <unistd.h>
//#ifndef CYGWIN
//#include <getopt.h>
//#endif
#if defined(__CYGWIN__)
#include <getopt.h>
#endif
#include "mdp_allo.h"
#include "tok_input_util.h"

View File

@@ -333,7 +333,7 @@ static double *smalloc(size_t n)
#endif
double *pntr;
if (n < 0) {
Fprintf(stderr, "smalloc ERROR: Non-positive argument. (%d)\n", n);
Fprintf(stderr, "smalloc ERROR: Non-positive argument. (%d)\n", int(n));
return NULL;
}
else if (n == 0) pntr = NULL;
@@ -344,7 +344,7 @@ static double *smalloc(size_t n)
}
if (pntr == NULL && n != 0) {
Fprintf(stderr, "smalloc : Out of space - number of bytes "
"requested = %d\n", n);
"requested = %d\n", int(n));
}
#ifdef MDP_MEMDEBUG
if (firsttime) {