*** empty log message ***

This commit is contained in:
Dave Goodwin 2004-08-05 14:26:37 +00:00
parent d7e6c19067
commit d23e5ce699
27 changed files with 4540 additions and 2764 deletions

View file

@ -12,8 +12,22 @@ module cantera_thermo
integer :: nsp
end type phase_t
! these definitions are for use with the equilibrate function.
integer, parameter :: TV = 100
integer, parameter :: HP = 101
integer, parameter :: SP = 102
integer, parameter :: PV = 103
integer, parameter :: TP = 104
integer, parameter :: UV = 105
integer, parameter :: SV = 107
integer, parameter :: VT = -100
integer, parameter :: PH = -101
integer, parameter :: PS = -102
integer, parameter :: VP = -103
integer, parameter :: PT = -104
integer, parameter :: VU = -105
integer, parameter :: VS = -107
contains

View file

@ -0,0 +1,430 @@
module fct
interface
subroutine cantera_error(proc, msg)
character*(*), intent(in) :: proc
character*(*), intent(in) :: msg
end subroutine cantera_error
integer function phase_nelements(n)
integer, intent(in) :: n
end function phase_nelements
integer function phase_nspecies(n)
integer, intent(in) :: n
end function phase_nspecies
double precision function phase_temperature(n)
integer, intent(in) :: n
end function phase_temperature
integer function phase_settemperature(n, t)
integer, intent(in) :: n
double precision, intent(in) :: t
end function phase_settemperature
double precision function phase_density(n)
integer, intent(in) :: n
end function phase_density
integer function phase_setdensity(n, rho)
integer, intent(in) :: n
double precision, intent(in) :: rho
end function phase_setdensity
double precision function phase_molardensity(n)
integer, intent(in) :: n
end function phase_molardensity
double precision function phase_meanmolecularweight(n)
integer, intent(in) :: n
end function phase_meanmolecularweight
integer function phase_elementindex(n, nm)
integer, intent(in) :: n
character*(*), intent(in) :: nm
end function phase_elementindex
integer function phase_speciesindex(n, nm)
integer, intent(in) :: n
character*(*), intent(in) :: nm
end function phase_speciesindex
integer function phase_getmolefractions(n, x)
integer, intent(in) :: n
double precision, intent(out) :: x(*)
end function phase_getmolefractions
double precision function phase_molefraction(n, k)
integer, intent(in) :: n
integer, intent(in) :: k
end function phase_molefraction
integer function phase_getmassfractions(n, y)
integer, intent(in) :: n
double precision, intent(out) :: y(*)
end function phase_getmassfractions
double precision function phase_massfraction(n, k)
integer, intent(in) :: n
integer, intent(in) :: k
end function phase_massfraction
integer function phase_setmolefractions(n, x, norm)
integer, intent(in) :: n
double precision, intent(in) :: x(*)
integer, intent(in) :: norm
end function phase_setmolefractions
integer function phase_setmolefractionsbyname(n, x)
integer, intent(in) :: n
character*(*), intent(in) :: x
end function phase_setmolefractionsbyname
integer function phase_setmassfractions(n, y, norm)
integer, intent(in) :: n
double precision, intent(in) :: y(*)
integer, intent(in) :: norm
end function phase_setmassfractions
integer function phase_setmassfractionsbyname(n, y)
integer, intent(in) :: n
character*(*), intent(in) :: y
end function phase_setmassfractionsbyname
integer function phase_getatomicweights(n, atw)
integer, intent(in) :: n
double precision, intent(out) :: atw(*)
end function phase_getatomicweights
integer function phase_getmolecularweights(n, mw)
integer, intent(in) :: n
double precision, intent(out) :: mw(*)
end function phase_getmolecularweights
integer function phase_getspeciesname(n, k, nm)
integer, intent(in) :: n
integer, intent(in) :: k
character*(*), intent(out) :: nm
end function phase_getspeciesname
integer function phase_getelementname(n, m, nm)
integer, intent(in) :: n
integer, intent(in) :: m
character*(*), intent(out) :: nm
end function phase_getelementname
double precision function phase_natoms(n, k, m)
integer, intent(in) :: n
integer, intent(in) :: k
integer, intent(in) :: m
end function phase_natoms
integer function newthermofromxml(mxml)
integer, intent(in) :: mxml
end function newthermofromxml
integer function th_eostype(n)
integer, intent(in) :: n
end function th_eostype
double precision function th_enthalpy_mole(n)
integer, intent(in) :: n
end function th_enthalpy_mole
double precision function th_intenergy_mole(n)
integer, intent(in) :: n
end function th_intenergy_mole
double precision function th_entropy_mole(n)
integer, intent(in) :: n
end function th_entropy_mole
double precision function th_gibbs_mole(n)
integer, intent(in) :: n
end function th_gibbs_mole
double precision function th_cp_mole(n)
integer, intent(in) :: n
end function th_cp_mole
double precision function th_cv_mole(n)
integer, intent(in) :: n
end function th_cv_mole
double precision function th_pressure(n)
integer, intent(in) :: n
end function th_pressure
double precision function th_enthalpy_mass(n)
integer, intent(in) :: n
end function th_enthalpy_mass
double precision function th_intEnergy_mass(n)
integer, intent(in) :: n
end function th_intEnergy_mass
double precision function th_entropy_mass(n)
integer, intent(in) :: n
end function th_entropy_mass
double precision function th_gibbs_mass(n)
integer, intent(in) :: n
end function th_gibbs_mass
double precision function th_cp_mass(n)
integer, intent(in) :: n
end function th_cp_mass
double precision function th_cv_mass(n)
integer, intent(in) :: n
end function th_cv_mass
integer function th_chempotentials(n, murt)
integer, intent(in) :: n
double precision, intent(out) :: murt(*)
end function th_chempotentials
integer function th_setpressure(n, p)
integer, intent(in) :: n
double precision, intent(in) :: p
end function th_setpressure
integer function th_set_hp(n, v1, v2)
integer, intent(in) :: n
double precision, intent(in) :: v1
double precision, intent(in) :: v2
end function th_set_hp
integer function th_set_uv(n, v1, v2)
integer, intent(in) :: n
double precision, intent(in) :: v1
double precision, intent(in) :: v2
end function th_set_uv
integer function th_set_sv(n, v1, v2)
integer, intent(in) :: n
double precision, intent(in) :: v1
double precision, intent(in) :: v2
end function th_set_sv
integer function th_set_sp(n, v1, v2)
integer, intent(in) :: n
double precision, intent(in) :: v1
double precision, intent(in) :: v2
end function th_set_sp
integer function th_equil(n, XY)
integer, intent(in) :: n
integer, intent(in) :: XY
end function th_equil
double precision function th_refpressure(n)
integer, intent(in) :: n
end function th_refpressure
double precision function th_mintemp(n, k)
integer, intent(in) :: n
integer, intent(in) :: k
end function th_mintemp
double precision function th_maxtemp(n, k)
integer, intent(in) :: n
integer, intent(in) :: k
end function th_maxtemp
integer function th_getenthalpies_rt(n, h_rt)
integer, intent(in) :: n
double precision, intent(out) :: h_rt(*)
end function th_getenthalpies_rt
integer function th_getentropies_r(n, s_r)
integer, intent(in) :: n
double precision, intent(out) :: s_r(*)
end function th_getentropies_r
integer function th_getcp_r(n, lenm, cp_r)
integer, intent(in) :: n
integer, intent(out) :: lenm
double precision, intent(out) :: cp_r(*)
end function th_getcp_r
integer function newkineticsfromxml(mxml, iphase, neighbor1, neighbor2, neighbor3, neighbor4)
integer, intent(in) :: mxml
integer, intent(in) :: iphase
integer, intent(in) :: neighbor1
integer, intent(in) :: neighbor2
integer, intent(in) :: neighbor3
integer, intent(in) :: neighbor4
end function newkineticsfromxml
integer function kin_type(n)
integer, intent(in) :: n
end function kin_type
integer function kin_start(n, p)
integer, intent(in) :: n
integer, intent(in) :: p
end function kin_start
integer function kin_speciesindex(n, nm, ph)
integer, intent(in) :: n
character*(*), intent(in) :: nm
character*(*), intent(in) :: ph
end function kin_speciesindex
integer function kin_ntotalspecies(n)
integer, intent(in) :: n
end function kin_ntotalspecies
integer function kin_nreactions(n)
integer, intent(in) :: n
end function kin_nreactions
double precision function kin_reactantstoichcoeff(n, k, i)
integer, intent(in) :: n
integer, intent(in) :: k
integer, intent(in) :: i
end function kin_reactantstoichcoeff
double precision function kin_productstoichcoeff(n, k, i)
integer, intent(in) :: n
integer, intent(in) :: k
integer, intent(in) :: i
end function kin_productstoichcoeff
integer function kin_reactiontype(n, i)
integer, intent(in) :: n
integer, intent(in) :: i
end function kin_reactiontype
integer function kin_getfwdratesofprogress(n, fwdROP)
integer, intent(in) :: n
double precision, intent(out) :: fwdROP(*)
end function kin_getfwdratesofprogress
integer function kin_getrevratesofprogress(n, revROP)
integer, intent(in) :: n
double precision, intent(out) :: revROP(*)
end function kin_getrevratesofprogress
integer function kin_isreversible(n, i)
integer, intent(in) :: n
integer, intent(in) :: i
end function kin_isreversible
integer function kin_getnetratesofprogress(n, netROP)
integer, intent(in) :: n
double precision, intent(out) :: netROP(*)
end function kin_getnetratesofprogress
integer function kin_getcreationrates(n, cdot)
integer, intent(in) :: n
double precision, intent(out) :: cdot(*)
end function kin_getcreationrates
integer function kin_getdestructionrates(n, ddot)
integer, intent(in) :: n
double precision, intent(out) :: ddot(*)
end function kin_getdestructionrates
integer function kin_getnetproductionrates(n, wdot)
integer, intent(in) :: n
double precision, intent(out) :: wdot(*)
end function kin_getnetproductionrates
double precision function kin_multiplier(n, i)
integer, intent(in) :: n
integer, intent(in) :: i
end function kin_multiplier
integer function kin_getequilibriumconstants(n, kc)
integer, intent(in) :: n
double precision, intent(out) :: kc(*)
end function kin_getequilibriumconstants
integer function kin_getreactionstring(n, i, buf)
integer, intent(in) :: n
integer, intent(in) :: i
character*(*), intent(out) :: buf
end function kin_getreactionstring
integer function kin_setmultiplier(n, i, v)
integer, intent(in) :: n
integer, intent(in) :: i
double precision, intent(out) :: v
end function kin_setmultiplier
integer function kin_advancecoverages(n, tstep)
integer, intent(in) :: n
double precision, intent(in) :: tstep
end function kin_advancecoverages
integer function newtransport(model, ith, loglevel)
character*(*), intent(in) :: model
integer, intent(in) :: ith
integer, intent(in) :: loglevel
end function newtransport
double precision function trans_viscosity(n)
integer, intent(in) :: n
end function trans_viscosity
double precision function trans_thermalConductivity(n)
integer, intent(in) :: n
end function trans_thermalConductivity
integer function trans_getThermalDiffCoeffs(n, dt)
integer, intent(in) :: n
double precision, intent(out) :: dt(*)
end function trans_getThermalDiffCoeffs
integer function trans_getMixDiffCoeffs(n, d)
integer, intent(in) :: n
double precision, intent(out) :: d(*)
end function trans_getMixDiffCoeffs
integer function trans_getBinDiffCoeffs(n, ld, d)
integer, intent(in) :: n
integer, intent(in) :: ld
double precision, intent(out) :: d(*)
end function trans_getBinDiffCoeffs
integer function trans_getMultiDiffCoeffs(n, ld, d)
integer, intent(in) :: n
integer, intent(in) :: ld
double precision, intent(out) :: d(*)
end function trans_getMultiDiffCoeffs
integer function trans_setParameters(n, type, k, d)
integer, intent(in) :: n
integer, intent(in) :: type
integer, intent(in) :: k
double precision, intent(in) :: d(*)
end function trans_setParameters
integer function ctphase_report(nth, buf, show_thermo)
integer, intent(in) :: nth
character*(*), intent(out) :: buf
integer, intent(in) :: show_thermo
end function ctphase_report
integer function ctgetCanteraError(buf)
character*(*), intent(out) :: buf
end function ctgetCanteraError
integer function ctaddCanteraDirectory(buflen, buf)
integer, intent(in) :: buflen
character*(*), intent(in) :: buf
end function ctaddCanteraDirectory
integer function ctbuildSolutionFromXML(src, ixml, id, ith, ikin)
character*(*), intent(in) :: src
integer, intent(in) :: ixml
character*(*), intent(in) :: id
integer, intent(in) :: ith
integer, intent(in) :: ikin
end function ctbuildSolutionFromXML
end interface
end module fct

View file

@ -1,16 +1,16 @@
import sys
bindir = 'c:/cantera/bin'
libdir = 'd:/dgg/dv/sf/cantera/build/lib/i686-pc-win32'
incdir = 'd:/dgg/dv/sf/cantera/build/include'
dflibdir = 'D:\Program Files\Microsoft Visual Studio\DF98\LIB'
bindir = '/Applications/Cantera/bin'
libdir = '/Users/dgg/dv/sf/cantera/build/lib/powerpc-apple-darwin7.4.0'
incdir = '/Users/dgg/dv/sf/cantera/build/include'
dflibdir = ''
bllibstr = "-lctlapack -lctblas"
bllibs = bllibstr.replace('-l',' ')
bllist = bllibs.split()
bldir = "d:/dgg/dv/sf/cantera/build/lib/i686-pc-win32"
bldir = "/Users/dgg/dv/sf/cantera/build/lib/powerpc-apple-darwin7.4.0"
libs = ['clib', 'oneD', 'zeroD', 'transport', 'cantera', 'recipes',
'cvode', 'ctmath', 'tpx']

View file

@ -157,44 +157,10 @@ class XML_Node:
s += line
return s
## def getRef(self):
## if not self["idRef"]: return self
## return find_XML(src = self["src"], root = self.root(),
## id = self["idRef"])
def clear_XML():
_cantera.xml_clear()
## def find_XML(src = "", root = None, id = "", loc = "", name=""):
## doc = None
## r = None
## if src:
## ihash = string.find(src,'#')
## if ihash < 0:
## fname = src
## else:
## fname, idnew = string.split(src,'#')
## if idnew: id = idnew
## if fname:
## doc = XML_Node(name="doc", src=fname)
## root = None
## elif root:
## doc = root
## elif root:
## doc = root
## else:
## raise exceptions.CanteraError("either root or src must be specified.")
## ## try:
## if loc or id or name:
## r = doc.child(loc=loc, id=id, name=name)
## else:
## r = doc
## return r
def getFloatArray(node, convert_units=0):
sz = int(node['size'])
return _cantera.ctml_getFloatArray(node._xml_id, convert_units, sz)

View file

@ -294,9 +294,9 @@ namespace Cantera {
}
}
if (j < m_mm)
return -1;
//throw CanteraError("estimateElementPotentials",
// "too few species (" + int2str(j) + ").");
// return -1;
throw CanteraError("estimateElementPotentials",
"too few species (" + int2str(j) + ").");
for (m = 0; m < m_mm; m++) {
for (n = 0; n < m_mm; n++) {

View file

@ -1,56 +0,0 @@
#ifndef CT_FTNODESYS_H
#define CT_FTNODESYS_H
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
namespace Cantera {
typedef void (*Ftn_RHS_Func) (integer* n, doublereal* t, doublereal* y,
doublereal* ydot);
/**
* A class to integrate a system of ODE's defined by a Fortran function.
* Not currently used.
*/
class FtnODESys : public FuncEval {
public:
FtnODESys(int n, doublereal* y0, Ftn_RHS_Func f) {
m_func = f;
m_n = n;
m_y0.resize(n);
int i;
for (i = 0; i < n; i++) m_y0[i] = y0[i];
}
virtual ~FtnODESys(){}
virtual void eval(double t, double* y, double* ydot) {
m_func(&m_n, &t, y, ydot);
}
virtual void getInitialConditions(double t0, double* y) {
int i;
for (i = 0; i < m_n; i++) y[i] = m_y0[i];
}
virtual int neq() { return m_n; }
protected:
int m_n;
vector_fp m_y0;
Ftn_RHS_Func m_func;
private:
};
}
#endif

View file

@ -1,118 +0,0 @@
/**
* @file ImplicitChem.h
*
* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_IMPCHEM_H
#define CT_IMPCHEM_H
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
#include "FuncEval.h"
#include "CVode.h"
#include "Kinetics.h"
#include "ThermoPhase.h"
namespace Cantera {
/**
* Advances the composition of an associated phase object in time
* by implicitly integrating
* \f[
* \dot Y_k = \frac{\omega_k}{\rho}
* \f]
*/
class ImplicitChem : public FuncEval {
public:
/**
* Constructor.
*/
ImplicitChem(Kinetics& kin, ThermoPhase& therm);
/**
* Destructor. Deletes the integrator.
*/
virtual ~ImplicitChem(){ delete m_integ; }
/**
* Overloads the virtual function
* declared in FuncEval.
*/
virtual void initialize(doublereal t0 = 0.0);
void adiabatic() {
m_energy = true;
}
void isothermal() {
m_energy = false;
}
/**
* Integrate from t0 to t1. The integrator is reinitialized
* first.
*/
void integrate(doublereal t0, doublereal t1) {
m_integ->reinitialize(t0, *this);
m_integ->setMaxStepSize(t1 - t0);
m_rho = m_thermo->density();
m_integ->integrate(t1);
updateState(m_integ->solution());
}
/**
* Integrate from t0 to t1 without reinitializing the
* integrator.
*/
void integrate0(doublereal t0, doublereal t1) {
m_integ->integrate(t1);
updateState(m_integ->solution());
}
// overloaded methods of class FuncEval
virtual int neq() { return m_nsp; }
virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
virtual void getInitialConditions(doublereal t0, size_t leny,
doublereal* y);
protected:
/**
* Set the mixture to a state consistent with solution
* vector y.
*/
void updateState(doublereal* y);
//Kinetics::phase_t* m_mix;
Kinetics* m_kin;
ThermoPhase* m_thermo;
int m_nsp;
Integrator* m_integ; // pointer to integrator
doublereal m_atol, m_rtol; // tolerances
doublereal m_maxstep; // max step size
array_fp m_wt;
doublereal m_rho;
bool m_energy;
doublereal m_h0;
doublereal m_press;
private:
};
}
#endif

View file

@ -1,155 +0,0 @@
deprecated
/**
*
* @file Resid1D.h
*
* >>>>> Under construction! <<<<<
*
* $Author$
* $Date$
* $Revision$
*
* Copyright 2002 California Institute of Technology
*
*/
#ifndef CT_MULTIDOM_H
#define CT_MULTIDOM_H
#include "stringUtils.h"
#include "Resid1D.h"
namespace Cantera {
/**
* Residual function evaluator for a one-dimensional problem.
*/
class MultiDomain {
public:
/**
* Constructor.
* @param nv Number of variables at each grid point.
* @param points Number of grid points.
*/
MultiDomain() : m_bw(0), m_nd(0), m_rdt(0.0), m_jac_ok(false) {}
/// Destructor.
virtual ~MultiDomain(){}
int nDomains() const { return m_nd; }
Resid1D& domain(int i) const { return *m_dom[i]; }
int start(int i) const { return m_start[i]; }
int size() const { return m_size; }
int bandwidth() const { return m_bw; }
virtual void addDomain(Resid1D* d) {
m_dom.push_back(d);
int sz = d->nPoints()
* d->nComponents();
if (m_nd > 0) {
m_start.push_back(m_start.back() + m_states.back());
}
else
m_start.push_back(0);
m_states.push_back(sz);
m_comp.push_back(d->nComponents());
m_points.push_back(d->nPoints());
int bw1, bw2 = 0;
bw1 = 2*d->nComponents() - 1;
if (m_nd > 0) {
bw2 = d->nComponents() + m_dom[m_nd-1]->nComponents() - 1;
}
if (bw1 > m_bw) m_bw = bw1;
if (bw2 > m_bw) m_bw = bw2;
m_nd = m_states.size();
m_size = m_start.back() + m_states.back();
}
void evalDomain(int i, int j, doublereal* x, doublereal* r,
doublereal rdt) {
int jpt;
if (j < 0)
jpt = j;
else
jpt = (j - m_start[i])/m_comp[i];
//cout << "calling Resid1D::eval. jpt, start = " << jpt << " " << m_start[i] << endl;
m_dom[i]->eval(jpt, x + m_start[i], r + m_start[i], rdt);
}
doublereal ssnorm(doublereal* x, doublereal* r) {
//cout << " calling eval to get ss norm " << endl;
eval(-1, x, r, 0.0);
doublereal ss = 0.0;
for (int i = 0; i < m_size; i++)
ss = fmaxx(fabs(r[i]),ss);
return ss;
}
doublereal rdt() const { return m_rdt; }
void initTimeInteg(doublereal dt, doublereal* x) {
doublereal rdt_old = m_rdt;
m_rdt = 1.0/dt;
if (fabs(rdt_old - m_rdt) > Tiny) {
m_jac_ok = false;
}
int i;
for (i = 0; i < m_nd; i++)
m_dom[i]->initTimeInteg(dt, x + m_start[i]);
}
bool transient() const { return (m_rdt != 0.0);}
bool steady() const { return (m_rdt == 0.0); }
void setSteadyMode() {
if (m_rdt > 0)
m_jac_ok = false;
m_rdt = 0.0;
}
void eval(int j, double* x, double* r, doublereal rdt=-1.0) {
int i, jpt;
if (rdt < 0.0) rdt = m_rdt;
if (j < 0) {
for (i = 0; i < m_nd; i++) {
evalDomain(i, j, x, r, rdt);
}
}
else {
for (i = 0; i < m_nd; i++) {
if (j >= m_start[i]) {
//cout << "calling evalDomain with j = " << j << endl;
evalDomain(i, j, x, r, rdt);
jpt = j - m_start[i];
if (jpt == 0 && i > 0)
evalDomain(i-1, j-1, x, r, rdt);
else if (jpt == m_states[i] - 1 && i < m_nd - 1)
evalDomain(i+1, j+1, x, r, rdt);
break;
}
}
}
}
protected:
doublereal m_rdt;
bool m_jac_ok;
int m_nd, m_bw, m_size;
vector_int m_states;
vector_int m_start;
vector_int m_comp, m_points;
vector<Resid1D*> m_dom;
private:
};
}
#endif

View file

@ -1,99 +0,0 @@
/**
* @file PolyThermo.h
*
* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#include "utilities.h"
#ifndef CT_POLYTHERMO_H
#define CT_POLYTHERMO_H
namespace Cantera {
/**
* A polynomial parameterization for one temperature range.
* Seven coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
* \f$ c_p^0(T)\f$, \f$ h^0(T)\f$, and \f$ s^0(T) \f$ as
* polynomials in \f$ T \f$ :
* \f[
* \frac{c_p(T)}{R} = a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4
* \f]
* \f[
* \frac{h^0(T)}{RT} = a_0 + \frac{a_1}{2} T + \frac{a_2}{3} T^2
+ \frac{a_3}{4} T^3 + \frac{a_4}{5} T^4 + a_5.
* \f]
* \f[
* \frac{s^0(T)}{R} = a_0\ln T + a_1 T + \frac{a_2}{2} T^2
+ \frac{a_3}{3} T^3 + \frac{a_4}{4} T^4 + a_6.
* \f]
*
* This class is designed specifically for use by class NasaThermo.
*/
template <int N>
class PolyThermo {
public:
PolyThermo() {}
void setCoefficients(const vector_fp& coeffs) {
m_c.resize(N + 1);
m_h.resize(N + 1);
m_s.resize(N + 1);
//copy(coeffs.begin(), coeffs.begin() + N + 3, m_c.begin());
m_c[0] = coeffs[2];
m_h[0] = m_c[0];
m_s[0] = coeffs[1];
for (int i = 1; i <= N; i++) {
m_c[i] = coeffs[i+2];
m_h[i] = m_c[i]/(i+1);
m_s[i] = m_c[i]/i;
}
m_h0 = coeffs[0];
m_s0 = coeffs[2];
}
virtual ~PolyThermo(){}
// these inline functions are defined only for readability
static doublereal tlog(const doublereal* tt) { return tt[0]; }
static doublereal trecip(const doublereal* tt) { return tt[1]; }
doublereal h0_RT(const doublereal* tt) const {
return m_c[0] * trecip(tt) + m_c[2]; }
doublereal s0_R(const doublereal* tt) const {
return m_c[2] * tlog(tt) + m_c[1]; }
void updateProperties(const doublereal* tt,
doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
doublereal cpt, ht, st;
cpt = m_c[0] + dot4(tt+2,m_c.begin()+1);
ht = m_h0*tt[1] + m_h[0] + dot4(tt+2,m_h.begin()+1);
st = m_s0*tt[0] + m_s[0] + dot4(tt+2,m_s.begin()+1); //s0_R(tt);
*cp_R = cpt;
*h_RT = ht;
*s_R = st;
}
protected:
doublereal m_h0, m_s0;
vector_fp m_c, m_h, m_s;
};
}
#endif

View file

@ -1,96 +0,0 @@
#ifndef CT_POLYTHERMMGR_H
#define CT_POLYTHERMOMGR_H
#include "PolyThermo.h"
#include "ctexceptions.h"
#include "polyfit.h"
#include "SpeciesThermo.h"
namespace Cantera {
/**
* A species thermodynamic property manager for polynomial
* parameterizations.
*/
template<int N>
class PolyThermoMgr : public SpeciesThermo {
public:
PolyThermoMgr() :
m_nsp(0),
m_minTemp(0.0),
m_maxTemp(1.e30),
m_p0(-1.0)
{ m_t.resize(N + 1); }
virtual ~PolyThermoMgr() {}
virtual void install(int index, int type, const vector_fp& coeffs,
doublereal minTemp = 0.0,
doublereal maxTemp = 1.e30,
doublereal refPressure = OneAtm) {
m_thermo.push_back(PolyThermo<N>());
PolyThermo<N>& th = m_thermo.back();
if (coeffs.size() != N + 3)
throw CanteraError("PolyThermoMgr::install",
string("wrong number of coefficients ")
+int2str(coeffs.size()));
th.setCoefficients(coeffs);
if (m_nsp > 0 && refPressure != m_p0) {
throw CanteraError("PolyThermoMgr::install",
"reference pressure mismatch");
}
else if (m_nsp == 0)
m_p0 = refPressure;
m_index.push_back(index);
m_nsp++;
if (minTemp > m_minTemp) m_minTemp = minTemp;
if (maxTemp < m_maxTemp) m_maxTemp = maxTemp;
}
virtual void update(doublereal t, vector_fp& cp_R,
vector_fp& h_RT, vector_fp& s_R) const {
int i;
m_t[0] = log(t);
m_t[1] = 1.0/t;
m_t[2] = t;
// int i;
const int nmax = N + 1;
for (i = 2; i < nmax; i++) m_t[i+1] = m_t[i]*t;
doublereal* bcp = cp_R.begin();
doublereal* bh = h_RT.begin();
doublereal* bs = s_R.begin();
doublereal* tt = m_t.begin();
size_t k;
for (k = 0; k < m_nsp; k++, bcp++, bh++, bs++) {
m_thermo[k].updateProperties(tt, bcp, bh, bs);
}
}
virtual doublereal minTemp(int k=-1) const {return m_minTemp;}
virtual doublereal maxTemp(int k=-1) const {return m_maxTemp;}
virtual doublereal refPressure() const {return m_p0;}
virtual int reportType(int index) const { return POLYNOMIAL_4;}
protected:
size_t m_nsp;
vector<PolyThermo<N> > m_thermo;
vector<int> m_index;
doublereal m_minTemp;
doublereal m_maxTemp;
doublereal m_p0;
mutable vector_fp m_t;
};
}
#endif

View file

@ -1,114 +0,0 @@
/**
* @file FuncEval.h
*
* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_RESIDEVAL_H
#define CT_RESIDEVAL_H
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif
namespace Cantera {
/**
*
* Virtual base class for DAE residual function evaluators.
*
*/
class ResidEval {
public:
/**
* Evaluate the residual function.
* @param t time (input, do not modify)
* @param y solution vector (input, do not modify)
* @param ydot rate of change of solution vector. (input, do
* not modify)
*/
virtual void evalResid(double t, const double deltaT,
const double* y,
const double* ydot,
double* resid)=0;
/** Number of equations. */
virtual int neq()=0;
// virtual bool init(double t0)=0;
virtual double* solution()=0;
virtual double* solution_dot()=0;
/**
* Fill the solution vector with the initial conditions
* at initial time t0.
*/
virtual void getInitialConditions(double t0, size_t leny, double* y)=0;
/**
* Report the analytical result if available and describe
* how close the numerical solution is to the analytical
* solution
*/
virtual void analytical_solution(double, double *)=0;
virtual void analytical_solution_dot(double, double *)=0;
/**
* Write out to a file or to standard output the current solution
* ievent is a description of the event that caused this
* function to be called.
*/
virtual void writeSolution(int ievent, double t,
const double *y, const double *ydot) = 0;
/**
* Apply a filter
*/
virtual void applyFilter(const double time, const double deltaT,
const double *ydot, const double *y,
double *delta_y) = 0;
/**
* Write out to a file or to standard output the current solution
* ievent is a description of the event that caused this
* function to be called.
*/
virtual void writeSolutionFilter(int ievent, double t,
const double * const y,
const double * const ydot,
const double * const delta_y,
const char * const fname) = 0;
/**
* This function returns a value of the delta t constraint
* that may exist in the problem.
*
*/
virtual double delta_t_constraint(const double t, const double *y,
const double *ydot) = 0;
virtual void adjustAtol(double *) = 0;
/**
* This function may be used to create output at various points in the
* execution of an application.
*
*/
virtual void user_out(const int ifunc, const double t, const double *y,
const double *ydot) = 0;
protected:
private:
};
}
#endif

View file

@ -1,116 +0,0 @@
/**
* @file TigerPolynomial.h
*
* $Author$
* $Revision$
* $Date$
*/
WARNING: May be out of date. Test before using!
// Copyright 2001 California Institute of Technology
#ifndef CT_TigerPOLY_H
#define CT_TigerPOLY_H
#include "TempFuncMgr.h"
namespace Cantera {
/**
* The Tiger polynomial parameterization.
*/
class TigerPolynomial {
public:
TigerPolynomial() : m_lowT(0.0), m_highT(0.0), m_Pref(0.0),
m_coeff(vector_fp(11)) {}
/**
* construct a Tiger polynomial parameterization from the
* values in vector coeffs.
* @param coeffs(0) minimum temperature (K)
* @param coeffs(1) maximum temperature (K)
* @param coeffs(2) standard-state pressure (Pa)
* @param coeffs(3)-coeffs(11) polynomial coefficients
*/
TigerPolynomial(int n, const vector_fp& coeffs) :
m_lowT (coeffs[0]),
m_highT (coeffs[1]),
m_Pref (coeffs[2]),
m_coeff (vector_fp(8)) {
m_coeff[0] = coeffs[3];
m_coeff[1] = 1.e-3 * coeffs[4];
m_coeff[2] = 1.e-6 * coeffs[5];
m_coeff[3] = 1.e-9 * coeffs[6];
m_coeff[4] = 1.e3 * coeffs[7];
m_coeff[5] = 1.e6 * coeffs[8];
m_coeff[6] = 1.e9 * coeffs[9];
m_coeff[7] = 1.e9 * coeffs[10] /
(GasConst_cal_mol_K * 298.15);
if (Debug::on) {
Debug::require(coeffs.size(), 11, Debug::NotLessThan);
}
}
doublereal minTemp() const { return m_lowT; }
doublereal maxTemp() const { return m_highT; }
doublereal refPressure() const { return m_Pref; }
void updateProperties(const TempFuncMgr& tt,
vector_fp& cp_R, vector_fp& h_RT, vector_fp& s_R,
vector_fp& g_RT) const {
doublereal ct0 = m_coeff[0];
doublereal ct1 = m_coeff[1]*tt.value(0);
doublereal ct2 = m_coeff[2]*tt.value(T2_INDEX);
doublereal ct3 = m_coeff[3]*tt.value(T3_INDEX);
doublereal ctm1 = m_coeff[4]*tt.value(TRECIP_INDEX);
doublereal ctm2 = m_coeff[5]/tt.value(T2_INDEX);
doublereal ctm3 = m_coeff[6]/tt.value(T3_INDEX);
cp_R[m_index] = ct0 + ct1 + ct2 + ct3
+ ctm1 + ctm2 + ctm3;
h_RT[m_index] = ct0 + 0.5 * ct1 + ct2/3.0
+ 0.25 * ct3
+ m_coeff[4] * tt.value(TLOG_INDEX) / tt.value(0)
- ctm2 - 0.5 * ctm3
+ m_coeff[7]*tt.value(TRECIP_INDEX);
s_R[m_index] = m_coeff[0] * tt.value(TLOG_INDEX) + ct1
+ 0.5 * ct2 + ct3/3.0
- ctm1 - 0.5 * ctm2 - ctm3/3.0
+ m_coeff[8] / GasConst_cal_mol_K ;
g_RT[m_index] = h_RT[m_index] - s_R[m_index];
}
protected:
doublereal m_lowT, m_highT, m_Pref;
vector_fp m_coeff;
int m_index;
private:
};
}
#endif

View file

@ -1,209 +0,0 @@
/**
*
*/
#ifndef CT_FLOW_BOUNDS
#define CT_FLOW_BOUNDS
#include <vector>
#include "ct_defs.h"
//#include "surfKinetics.h"
using namespace Cantera;
namespace FlowBdry {
class Boundary {
public:
Boundary(int nsp = 0)
: rtau(1.e5), m_lr(1), m_nsp(nsp), m_nv(nsp+4),
m_bv(0), m_bp(0), m_temp(0.0), m_mdot(0.0), m_V(0.0) {
m_y.resize(nsp, 0.1);
}
virtual ~Boundary(){}
doublereal rtau;
doublereal T() { return m_temp; }
doublereal mdot() { return m_mdot; }
doublereal V() { return m_V; }
doublereal* Y() { return m_y.begin(); }
int orient() { return m_lr; }
int faceLeft() { m_lr = -1; return m_lr; }
int faceRight() { m_lr = 1; return m_lr; }
void setNSpecies(int n) { m_nsp = n; }
void set_mdot(doublereal mdot) { m_mdot = mdot; }
void set_V(doublereal V) { m_V = V; }
void set_T(doublereal T) { m_temp = T; }
void set_Y(const doublereal* y) {
copy(y, y+m_nsp, m_y.begin());
}
virtual void evalInt(int j, doublereal* x, doublereal* r) {
int i;
for (i = 0; i < m_bv; i++) {
r[i] = x[i];
}
}
virtual void eval(doublereal* x0, doublereal density,
doublereal* data, doublereal* r) {
throw CanteraError("Boundary::eval",
" ERROR!! Base class eval called!");
}
virtual int nIntVars() { return m_bv; }
virtual int nIntPoints() { return m_bp; }
protected:
int offset(int n) { return m_lr*n*m_nv; }
int m_lr, m_nsp, m_nv;
int m_bv, m_bp;
doublereal m_temp, m_mdot, m_V;
vector_fp m_y;
};
/**
* zero axial velocity, zero gradients for everything else.
*/
class SymmPlane : public Boundary {
public:
SymmPlane(int nsp) : Boundary(nsp) {}
virtual ~SymmPlane() {}
virtual void eval(doublereal* x0, doublereal density,
doublereal* data, doublereal* r) {
for (int n = 0; n < m_nv; n++) {
r[n] = x0[n] - x0[offset(1) + n];
}
r[0] = x0[0];
}
protected:
int m_iloc;
};
/**
* Zero gradients for all components, except V, which is zero.
*/
class Outlet : public Boundary {
public:
Outlet(int nsp) : Boundary(nsp) {}
virtual ~Outlet() {}
virtual void eval(doublereal* x0, doublereal density,
doublereal* data, doublereal* r) {
for (int n = 0; n < m_nv; n++) {
r[n] = -rtau*(x0[n] - x0[offset(1) + n]);
}
r[1] = -rtau*x0[1];
}
};
#ifdef INCL_SURF
class Surface : public Boundary {
public:
Surface(int nsp, SurfKinetics* kin=0)
: Boundary(nsp), m_kin(kin), m_dt(1.e3) {
if (m_kin)
m_sdot.resize(m_kin->nTotal());
else
m_sdot.resize(nsp);
}
virtual ~Surface(){}
virtual void eval(doublereal* x0, doublereal density,
doublereal* data, doublereal* r) {
doublereal sum = 0.0, mdot = 0.0;
// updateSurfaceRates(x0[2], density, x0+4);
for (int k = 0; k < m_nsp; k++) {
// r[4+k] = m_sdot[k] - m_lr*data[k];
//mdot += m_sdot[k];
//sum += x0[4+k];
r[4+k] = - m_lr*data[k];
sum += x0[4+k];
}
r[4] = 1.0 - sum;
// match the Stefan velocity at the surface
r[0] = density*x0[0] - m_lr*mdot;
r[1] = x0[1]; // no slip
r[2] = x0[2] - m_temp; // specified T
r[3] = x0[3];
}
/**
* Update the mass species production rates, given
* the species mass fractions at the surface.
* @param y
*/
virtual void updateSurfaceRates(doublereal t,
doublereal rho, doublereal* y) {
int k;
if (m_kin) {
for (k = 0; k < m_nsp; k++) {
m_sdot[k] = fmaxx(Tiny, y[k]);
}
m_kin->bulkPhase(0)->setState_TR(t, rho);
m_kin->bulkPhase(0)->setMassFractions_NoNorm(m_sdot.begin());
//m_kin->integrate(m_dt);
m_kin->getNetProductionRates(m_sdot.begin());
}
else
for (k = 0; k < m_nsp; k++) m_sdot[k] = 0.0;
}
protected:
const doublereal* m_wt;
vector_fp m_sdot;
SurfKinetics* m_kin;
doublereal m_dt;
};
#endif
class Inlet : public Boundary {
public:
Inlet(int nsp, doublereal* wt=0)
: Boundary(nsp) {}
virtual ~Inlet(){}
virtual void eval(doublereal* x0,
doublereal density, doublereal* flux, doublereal* r) {
int k;
for (k = 0; k < m_nsp; k++) {
r[4+k] = rtau*(m_y[k] - x0[k+4] - flux[k]/m_mdot);
}
r[0] = -rtau*(density*x0[0] - m_lr*m_mdot);
r[1] = -rtau*(x0[1] - V());
r[2] = -rtau*(x0[2] - T()); // specified T
r[3] = -rtau*x0[3];
}
};
}
#endif

View file

@ -1,11 +0,0 @@
#ifndef CT_FUNCS_H
#define CT_FUNCS_H
#include "ct_defs.h"
namespace Cantera {
doublereal linearInterp(doublereal x, const vector_fp& xpts,
const vector_fp& fpts);
}
#endif

View file

@ -1,29 +0,0 @@
/**
* @file gases.h
*
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_GASES_H
#define CT_GASES_H
#include "IdealGasMix.h"
#include "TransportFactory.h"
namespace Cantera {
template<class T>
class GasWithTransport : public IdealGasMix, public T {
public:
GasWithTransport(map<string, string>& params,
TransportFactory* f = 0) {
if (f == 0)
f = TransportFactory::factory();
f->initTransport(*self, dbase);
}
};
}
#endif

View file

@ -1,41 +0,0 @@
/**
* @file import.h
*
* functions to import objects from files.
*
* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_IMPORTERS_H
#define CT_IMPORTERS_H
#include "importCK.h"
//#include "importXML.h"
#include "importCTML.h"
namespace Cantera {
/**
* import the specifications for a phase from a file, including
* elements, species, and reactions. The file formats currently
* supported are
*
* - CKML "Chemical Kinetics Markup Language." This XML-based
* markup language has been developed by M. Aivazis and R. Muller
* at Caltech and is used in their Fuego software package.
*
* - CK This is the name we give to the widely-used format
* developed by Kee, Miller, and Rupley for the Chemkin-II
* software package.
*
*/
bool importFromFile(Thermo* th, Kinetics* k, map<string, string>& params);
}
#endif

View file

@ -1,25 +0,0 @@
/**
* @file importXML.h
*
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_IMPORTXML_H
#define CT_IMPORTXML_H
#include <string>
#include "GasKinetics.h"
namespace Cantera {
bool isXMLFile(string infile);
bool importXML(const string& infile, phase_t& gas, Kinetics& kin);
}
#endif

View file

@ -1,95 +0,0 @@
#ifndef CT_MIX_UTILS_H
#define CT_MIX_UTILS_H
#include "ctexceptions.h"
namespace Cantera {
doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y);
// /** Set the temperature (K), pressure (Pa), and mole fractions. */
// template<class S>
// void set_TPX(S& s, doublereal t, doublereal p, doublereal* x) {
// s.setMoleFractions(x); s.setTemperature(t); s.setPressure(p);
// }
// template<class S>
// void set_HP(S& s, doublereal h, doublereal p, doublereal tol) {
// doublereal dt;
// s.setPressure(p);
// for (int n = 0; n < 20; n++) {
// dt = (h - s.enthalpy_mass())/s.cp_mass();
// if (dt > 100.0) dt = 100.0;
// else if (dt < -100.0) dt = -100.0;
// s.setState_TP(s._temp() + dt, p);
// if (fabs(dt) < tol) {
// return;
// }
// }
// throw CanteraError("set_HP","no convergence. dt = " + fp2str(dt));
// }
// template<class S>
// void set_UV(S& s, doublereal u, doublereal v,
// doublereal tol) {
// doublereal dt;
// s.setDensity(1.0/v);
// for (int n = 0; n < 20; n++) {
// dt = (u - s.intEnergy_mass())/s.cv_mass();
// if (dt > 100.0) dt = 100.0;
// else if (dt < -100.0) dt = -100.0;
// s.setTemperature(s._temp() + dt);
// if (fabs(dt) < tol) return;
// }
// throw CanteraError("set_UV","no convergence. dt = " + fp2str(dt));
// }
// template<class S>
// void set_SP(S& s, doublereal entropy, doublereal p,
// doublereal tol) {
// doublereal dt;
// s.setPressure(p);
// for (int n = 0; n < 20; n++) {
// dt = (entropy - s.entropy_mass())*s.temperature()/s.cp_mass();
// if (dt > 100.0) dt = 100.0;
// else if (dt < -100.0) dt = -100.0;
// s.setState_TP(s._temp() + dt, p);
// if (fabs(dt) < tol) return;
// }
// throw CanteraError("set_SP","no convergence. dt = " + fp2str(dt));
// }
// template<class S>
// void set_SV(S& s, doublereal entropy, doublereal v,
// doublereal tol) {
// doublereal dt;
// s.setDensity(1.0/v);
// for (int n = 0; n < 20; n++) {
// dt = (entropy - s.entropy_mass())*s.temperature()/s.cv_mass();
// if (dt > 100.0) dt = 100.0;
// else if (dt < -100.0) dt = -100.0;
// s.setTemperature(s._temp() + dt);
// if (fabs(dt) < tol) return;
// }
// throw CanteraError("set_SV","no convergence. dt = " + fp2str(dt));
// }
template<class S1, class S2>
void mapSpeciesData(const S1& s1, const S2& s2, const doublereal* data1,
doublereal* data2) {
int n1 = s1.nSpecies();
int n2 = s2.nSpecies();
int n, m2;
// zero out the destination array
for (n = 0; n < n2; n++) data2[n] = 0.0;
// copy
for (n = 0; n < n1; n++) {
m2 = s2.speciesIndex(s1.speciesName(n));
if (m2 >= 0) data2[m2] = data1[n];
}
}
}
#endif

View file

@ -83,7 +83,7 @@ namespace Cantera {
doublereal mout = 0.0;
for (int i = 0; i < nout; i++)
mout += m_outlet[i]->massFlowRate();
return volume()/mout;
return mass()/mout;
}
FlowDevice& ReactorBase::inlet(int n) { return *m_inlet[n]; }

View file

@ -104,6 +104,9 @@ tools-install:
hdr-install:
@INSTALL@ -d @ct_incdir@
cp -r -f build/include/cantera @ct_incroot@
ifeq ($(build_f90),1)
cp -f build/include/cantera/*.mod @ct_incroot@
endif
# collect scattered header files and build the include directory
hdr-collect:

5424
config/configure vendored

File diff suppressed because it is too large Load diff

View file

@ -441,7 +441,8 @@ if test "x$BUILD_F90_INTERFACE" = "xy"; then
fi
AC_SUBST(BUILD_F90)
AC_SUBST(F90)
AC_SUBST(F77FLAGS)
AC_SUBST(F90FLAGS)
# AC_SUBST(F77FLAGS)
BUILD_CLIB=1
#
@ -566,7 +567,7 @@ fi
dnl Checks for libraries.
AC_F77_LIBRARY_LDFLAGS()
case $ac_sys_system in
Darwin*) FLIBS='-lg2c -lgcc'; SHARED_CTLIB=0;;
Darwin*) FLIBS='-lstdc++ -lg2c'; SHARED_CTLIB=0;;
esac
@ -685,6 +686,7 @@ AC_OUTPUT(../Cantera/Makefile \
../tools/src/finish_install.py \
../tools/src/package4mac \
../tools/templates/f77/demo.mak \
../tools/templates/f90/demo.mak \
../tools/templates/cxx/demo.mak \
../tools/testtools/Makefile \
../data/inputs/Makefile \

16
configure vendored
View file

@ -116,7 +116,13 @@ BUILD_MATLAB_TOOLBOX=${BUILD_MATLAB_TOOLBOX:="y"}
# A Fortran 90/95 compiler is required.
# UNDER DEVELOPMENT - UNSTABLE!!
BUILD_F90_INTERFACE=${BUILD_F90_INTERFACE:="n"}
BUILD_F90_INTERFACE=${BUILD_F90_INTERFACE:="y"}
# The Fortran 90/95 compiler
F90=${F90:=f95}
# Compiler option flags for the Fortran 90/95 compiler
F90FLAGS=${F90FLAGS:='-O2'}
#----------------------------------------------------------------------
@ -258,9 +264,8 @@ SHARED=${SHARED:="-shared"}
# compiler if you are building the Fortran 90 interface.
F77=${F77:=g77}
F90=${F90:=f95}
# Fortran compiler flags
# Fortran 77 compiler flags. Note that these may be specified separately
# from the options for the Fortran 90/95 compiler (if any).
FFLAGS=${FFLAGS:='-O2'}
# additional flags required to compile fixed-format F77 procedures
@ -324,7 +329,7 @@ export RANLIB
export BLAS_LAPACK_LIBS
export BLAS_LAPACK_DIR
export BUILD_F90
export BUILD_FORTRAN_90_INTERFACE
export BUILD_F90_INTERFACE
export PYTHON_PACKAGE
export BUILD_MATLAB_TOOLBOX
#export MATLAB_CMD
@ -351,6 +356,7 @@ export F90
export F90_EXT
export FFLAGS
export F77FLAGS
export F90FLAGS
export LAPACK_NAMES
export LCXX_FLAGS
export LCXX_END_LIBS

View file

@ -47,6 +47,9 @@ templates-install:
@INSTALL@ -m 644 templates/f77/*.cpp @ct_dir@/templates/f77
@INSTALL@ -m 644 templates/f77/*.f @ct_dir@/templates/f77
@INSTALL@ -m 644 templates/f77/*.mak @ct_dir@/templates/f77
@INSTALL@ -d @ct_dir@/templates/f90
@INSTALL@ -m 644 templates/f90/*.f90 @ct_dir@/templates/f90
@INSTALL@ -m 644 templates/f90/*.mak @ct_dir@/templates/f90
@INSTALL@ -d @ct_dir@/templates/cxx
@INSTALL@ -m 644 templates/cxx/*.cpp @ct_dir@/templates/cxx
@INSTALL@ -m 644 templates/cxx/*.mak @ct_dir@/templates/cxx

View file

@ -62,6 +62,8 @@ f = open(bindir+'/ctnew','w')
f.write("""#!/bin/sh
if test "x$1" = "x-f77"; then
cp """+templdir+"""/f77/*.* .
elif test "x$1" = "x-f90"; then
cp """+templdir+"""/f90/*.* .
else
cp """+templdir+"""/cxx/*.* .
fi

View file

@ -0,0 +1,96 @@
!
! Replace this sample main program with your program
!
! This program uses functions defined in demo_ftnlib.cpp to create
! an ideal gas mixture and print some its properties.
!
! For a C++ version of this program, see ../cxx/demo.cpp.
!
program main
call demo(100, 500)
stop
end
subroutine demo(maxsp, maxrxns)
use cantera
implicit none
integer, intent(in) :: maxsp
integer, intent(in) :: maxrxns
type(phase_t) gas
double precision q(MAXRXNS), qf(MAXRXNS), qr(MAXRXNS)
double precision diff(MAXSP)
character*80 eq
character*20 name
double precision :: t, p, dnu, dlam
integer :: i, irxns, nsp, k
write(*,*)
write(*,*) '******** Fortran 90 Test Program ********'
gas = importPhase('h2o2.cti','ohmech')
t = 1200.0
p = 101325.0
call setState_TPX(gas, t, p, 'H2:1, O2:1, AR:2')
write(*,*) 'Initial state properties:'
write(*,10) temperature(gas), pressure(gas), density(gas), &
enthalpy_mole(gas), entropy_mole(gas), cp_mole(gas)
! compute the equilibrium state holding the specifi!
! enthalpy and pressure constant
call equilibrate(gas, HP)
write(*,*) 'Equilibrium state properties:'
write(*,10) temperature(gas), pressure(gas), density(gas), &
enthalpy_mole(gas), entropy_mole(gas), cp_mole(gas)
10 format(//'Temperature: ',g14.5,' K'/ &
'Pressure: ',g14.5,' Pa'/ &
'Density: ',g14.5,' kg/m3'/ &
'Molar Enthalpy:',g14.5,' J/kmol'/ &
'Molar Entropy: ',g14.5,' J/kmol-K'/ &
'Molar cp: ',g14.5,' J/kmol-K'//)
! Reaction information
irxns = nReactions(gas)
! forward and reverse rates of progress should be equal
! in equilibrium states
call getFwdRatesOfProgress(gas, qf)
call getRevRatesOfProgress(gas, qr)
! net rates of progress should be zero in equilibrium states
call getNetRatesOfProgress(gas, q)
! for each reaction, print the equation and the rates of progress
do i = 1,irxns
call getReactionString(gas, i,eq)
write(*,20) eq,qf(i),qr(i),q(i)
20 format(a27,3e14.5,' kmol/m3/s')
end do
! Transport properties
dnu = viscosity(gas)
dlam = thermalConductivity(gas)
call getMixDiffCoeffs(gas, diff)
write(*,30) dnu, dlam
30 format(//'Viscosity: ',g14.5,' Pa-s'/ &
'Thermal conductivity: ',g14.5,' W/m/K'/)
write(*,*) 'Species Diffusion Coefficient'
nsp = nSpecies(gas)
do k = 1, nsp
call getSpeciesName(gas, k, name)
write(*,40) name, diff(k)
40 format(' ',a20,e14.5,' m2/s')
end do
stop
end subroutine demo

View file

@ -0,0 +1,94 @@
#!/bin/sh
# This Makefile builds a Fortran 90 application that uses Cantera. By
# default, the main program file is 'demo.f90,' which prints out some
# properties of a reacting gas mixture.
# To build program 'demo', simply type 'make', or 'make -f <this
# file>' if this file is named something other than 'Makefile.'
# Once you have verified that the demo runs, edit this file to replace
# object file 'demo.o' with your own object file or files.
#------------------------ edit this block ---------------------------------
# the name of the executable program to be created
PROG_NAME = demo
# the object files to be linked together.
OBJS = demo.o
# additional flags to be passed to the linker. If your program
# requires other external libraries, put them here
LINK_OPTIONS = @LCXX_FLAGS@
#---------------------------------------------------------------------------
# You probably don't need to edit anything below.
# the C++ compiler
CXX = @CXX@
# C++ compile flags
CXX_FLAGS = @CXXFLAGS@
# external libraries
EXT_LIBS = -lfct -lclib @LOCAL_LIBS@ -lctcxx -lstdc++ @FLIBS@
# the Fortran 90/95 compiler
F90 = @F90@
# Fortran compile flags
FORT_FLAGS = @F90FLAGS@
# the directory where the Cantera libraries are located
CANTERA_LIBDIR=@ct_libdir@
# the directory where Cantera include files may be found.
CANTERA_INCDIR=@ct_incroot@
# the directory where Cantera Fortran 90 modules may be found.
CANTERA_MODULE_DIR=@ct_incroot@/cantera
# flags passed to the C++ compiler/linker for the linking step
LCXXFLAGS = -L$(CANTERA_LIBDIR)
# how to compile C++ source files to object files
%.o : %.cpp
$(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS)
# how to compile Fortran 90/95 source files to object files
%.o : %.f90
$(F90) -c $< -I$(CANTERA_MODULE_DIR) $(FORT_FLAGS)
PROGRAM = $(PROG_NAME)$(EXE_EXT)
DEPENDS = $(OBJS:.o=.d)
all: $(PROGRAM)
$(PROGRAM): $(OBJS)
$(F90) -o $(PROGRAM) $(OBJS) $(LCXXFLAGS) $(CANTERA_LIBS) $(LINK_OPTIONS) $(EXT_LIBS) @LIBS@
%.d : %.cpp
g++ -MM -I$(CANTERA_INCDIR) $*.cpp > $*.d
clean:
$(RM) $(OBJS) $(PROGRAM)
depends: $(DEPENDS)
cat *.d > .depends
$(RM) $(DEPENDS)
TAGS:
etags *.h *.cpp
ifeq ($(wildcard .depends), .depends)
include .depends
endif