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

cleanup

Browse source
This commit is contained in:
Dave Goodwin 2005-06-25 05:27:30 +00:00
parent af048f0c7b
commit 0ce92af1c5
16 changed files with 134 additions and 133 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
Cantera/clib/src/ctmultiphase.cpp
View file

@ -74,7 +74,9 @@ namespace Cantera {
extern "C" {
int DLL_EXPORT mix_new() {
mix_t* m = new MultiPhase();
cout << "in mix_new" << endl;
mix_t* m = new MultiPhase;
cout << "did it" << endl;
return Cabinet<mix_t>::cabinet()->add(m);
}
@ -98,7 +100,7 @@ extern "C" {
int DLL_EXPORT mix_nElements(int i) {
return _mix(i)->nElements();
}
}
int DLL_EXPORT mix_elementIndex(int i, char* name) {
return _mix(i)->elementIndex(string(name));
@ -195,6 +197,7 @@ extern "C" {
err, maxsteps, maxiter, loglevel);
}
catch (CanteraError) {
write_logfile("equil_err.html");
return DERR;
}
}

8
Cantera/python/Cantera/OneD/__init__.py
View file

@ -2,9 +2,9 @@
The classes in this package implement one-dimensional reacting flow problems.
"""
from onedim import *
#from BurnerFlame import BurnerFlame
#from BurnerDiffFlame import BurnerDiffFlame
#from CounterFlame import CounterFlame
#from StagnationFlow import StagnationFlow
from BurnerFlame import BurnerFlame
from BurnerDiffFlame import BurnerDiffFlame
from CounterFlame import CounterFlame
from StagnationFlow import StagnationFlow

1
Cantera/python/examples/catcomb.py
View file

@ -14,6 +14,7 @@
from Cantera import *
from Cantera.OneD import *
#from Cantera.OneD.StagnationFlow import StagnationFlow
import math
###############################################################

32
Cantera/python/examples/flames/ohn.cti
View file

@ -1,35 +1,3 @@
#
# see http://reaflow.iwr.uni-heidelberg.de/~Olaf.Deutschmann/ for
# more about this mechanism
#
#---------------------------------------------------------------------!
#***********************************************************************
#**** *
#**** CH4-O2 SURFACE MECHANISM ON PT *
#**** *
#**** Version 1.2 November 1995 *
#**** *
#**** O. Deutschmann, IWR, Heidelberg University, Germany *
#**** *
#**** Kinetic data: *
#**** k = A * T**b * exp (-Ea/RT) A b Ea *
#**** (cm,mol,s) - J/mol *
#**** *
#**** *
#***********************************************************************
#
# Ref:- 1.) Deutschman et al., 26th Symp. (Intl.) on Combustion,1996
# pp. 1747-1754
#----------------------------------------------------------------------
#
# Converted to Cantera format
# by ck2cti on Thu Aug 21 07:58:45 2003
#
#----------------------------------------------------------------------
units(length = "cm", time = "s", quantity = "mol", act_energy = "J/mol")
#
# This definition extracts the O/H/N submechanism from GRI-Mech 3.0

36
Cantera/src/MultiPhase.cpp
View file

@ -7,6 +7,14 @@
namespace Cantera {
MultiPhase::MultiPhase() : m_temp(0.0), m_press(0.0),
m_nel(0), m_nsp(0), m_init(false), m_eloc(-1),
m_equil(0), m_Tmin(1.0), m_Tmax(100000.0) {
}
void MultiPhase::
addPhase(phase_t* p, doublereal moles) {
@ -54,6 +62,15 @@ namespace Cantera {
m_temp = p->temperature();
m_press = p->pressure();
}
//cout << "min, max = " << m_Tmin << " " << m_Tmax << endl;
if (p->nSpecies() > 1) {
double t = p->minTemp();
if (t > m_Tmin) m_Tmin = t;
t = p->maxTemp();
if (t < m_Tmax) m_Tmax = t;
//cout << p->name() << " " << t << " " << m_Tmax << endl;
}
}
@ -184,12 +201,16 @@ namespace Cantera {
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
void MultiPhase::getValidChemPotentials(doublereal not_mu,
doublereal* mu) {
doublereal* mu, bool standard) {
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
if (tempOK(i) || m_phase[i]->nSpecies() > 1)
m_phase[i]->getChemPotentials(mu + loc);
if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
if (!standard)
m_phase[i]->getChemPotentials(mu + loc);
else
m_phase[i]->getStandardChemPotentials(mu + loc);
}
else
fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
loc += m_phase[i]->nSpecies();
@ -359,6 +380,7 @@ namespace Cantera {
if (loglevel > 0) {
addLogEntry("problem type","fixed T,P");
}
// create an equilibrium manager
MultiPhaseEquil e(this);
error = e.equilibrate(XY, err, maxsteps, loglevel-1);
if (loglevel > 0) e.printInfo();
@ -368,18 +390,20 @@ namespace Cantera {
dt = 1.0e2;
h0 = enthalpy();
start = true;
Tlow = 1.0;
Thigh = 1.0e4;
Tlow = m_Tmin; // lower bound on T
Thigh = m_Tmax; // upper bound on T
hlow = 0.0;
hhigh = 0.0;
once = true;
if (loglevel > 0) {
addLogEntry("problem type","fixed H,P");
addLogEntry("H target",fp2str(h0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
}
ferr = 0.1;
for (n = 0; n < maxiter; n++) {
MultiPhaseEquil e(this, start);
MultiPhaseEquil e(this, strt);
start = false;
if (loglevel > 1) {
beginLogGroup("iteration "+int2str(n));

19
Cantera/src/MultiPhase.h
View file

@ -22,21 +22,20 @@ namespace Cantera {
public:
typedef size_t index_t;
typedef ThermoPhase phase_t;
typedef DenseMatrix array_t;
/// Constructor. The constructor takes no arguments, since
/// phases are added using method addPhase.
MultiPhase() : m_temp(0.0), m_press(0.0),
m_nel(0), m_nsp(0), m_init(false), m_eloc(-1),
m_equil(0) {}
MultiPhase();
/// Destructor. Does nothing. Class MultiPhase does not take
/// "ownership" (i.e. responsibility for destroying) the
/// phase objects.
virtual ~MultiPhase() {}
typedef size_t index_t;
typedef ThermoPhase phase_t;
typedef DenseMatrix array_t;
/// Add a phase to the mixture.
/// @param p pointer to the phase object
/// @param moles total number of moles of all species in this phase
@ -88,6 +87,8 @@ namespace Cantera {
return m_spstart[p] + k;
}
doublereal minTemp();
doublereal maxTemp();
doublereal charge();
doublereal phaseCharge(index_t p);
@ -104,7 +105,8 @@ namespace Cantera {
/// chemical potentials of all species with thermo data valid
/// for the current temperature [J/kmol]. For other species,
/// set the chemical potential to the value \c not_mu.
void getValidChemPotentials(doublereal not_mu, doublereal* mu);
void getValidChemPotentials(doublereal not_mu, doublereal* mu,
bool standard = false);
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
@ -195,6 +197,7 @@ namespace Cantera {
int m_eloc;
vector<bool> m_temp_OK;
MultiPhaseEquil* m_equil;
doublereal m_Tmin, m_Tmax;
};
inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) {

117
Cantera/src/MultiPhaseEquil.cpp
View file

@ -105,7 +105,6 @@ namespace Cantera {
"not valid at this temperature, but it has "
"non-zero moles in the initial state.");
}
//cout << "excluding species " << m_mix->speciesName(k) << endl;
}
}
for (k = 0; k < m_nsp_mix; k++) {
@ -118,6 +117,7 @@ namespace Cantera {
// some work arrays for internal use
m_work.resize(m_nsp);
m_work2.resize(m_nsp);
m_work3.resize(m_nsp_mix);
m_mu.resize(m_nsp_mix);
// number of moles of each species
@ -128,9 +128,6 @@ namespace Cantera {
index_t ik;
for (ik = 0; ik < m_nsp; ik++) {
m_moles[ik] = m_mix->speciesMoles(m_species[ik]);
//if (ISNAN(m_moles[ik])) {
// writelog("moles "+int2str(ik)+" initialized to nan \n");
//}
}
// Delta G / RT for each reaction
@ -147,21 +144,12 @@ namespace Cantera {
setInitialMoles();
computeN();
// make sure the components are non-zero
//for (k = 0; k < m_nel; k++) {
// if (m_moles[m_order[k]] <= 0.0) {
// m_moles[m_order[k]] = 1.0e-17;
// }
//}
vector_fp dxi(m_nsp - m_nel, 1.0e-20);
multiply(m_N, dxi.begin(), m_work.begin());
unsort(m_work);
for (k = 0; k < m_nsp; k++) {
m_moles[k] += m_work[k];
//if (ISNAN(m_moles[k])) {
// writelog("moles "+int2str(k)+" is nan 2. \n");
// }
m_lastmoles[k] = m_moles[k];
if (m_mix->solutionSpecies(m_species[k]))
m_dsoln.push_back(1);
@ -170,7 +158,7 @@ namespace Cantera {
}
m_force = false;
setMoles();
}
}
doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
int maxsteps, int loglevel) {
@ -187,17 +175,18 @@ namespace Cantera {
endLogGroup();
}
printInfo();
if (error() == 0.0) {
write_logfile("equil_err.html");
Cantera::error("stopping");
}
//if (error() == 0.0) {
// write_logfile("equil_err.html");
// Cantera::error("stopping");
//}
if (error() < err) break;
}
}
if (i >= maxsteps) {
if (loglevel > 0) {
addLogEntry("Error","no convergence in "+int2str(maxsteps)
+" iterations");
printInfo();
endLogGroup();
}
throw CanteraError("MultiPhaseEquil::equilibrate",
"no convergence in " + int2str(maxsteps) +
@ -209,18 +198,34 @@ namespace Cantera {
addLogEntry("error",fp2str(error()));
endLogGroup();
}
finish();
return error();
}
void MultiPhaseEquil::setMoles() {
vector_fp n(m_nsp_mix, 0.0);
//vector_fp n(m_nsp_mix, 0.0);
fill(m_work3.begin(), m_work3.end(), 0.0);
index_t k;
for (k = 0; k < m_nsp; k++) {
n[m_species[k]] = m_moles[k];
m_work3[m_species[k]] = m_moles[k];
}
m_mix->setMoles(n.begin());
m_mix->setMoles(m_work3.begin());
}
/// Clean up the composition by setting species with negative mole
/// numbers to zero. The solution algorithm can leave some species
/// in stoichiometric condensed phases with very small negative
/// mole numbers. This method simply sets these to zero.
void MultiPhaseEquil::finish() {
fill(m_work3.begin(), m_work3.end(), 0.0);
index_t k;
for (k = 0; k < m_nsp; k++) {
m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0);
}
m_mix->setMoles(m_work3.begin());
}
/**
* Estimate the initial mole fractions. Uses the Simplex method
* to estimate the initial number of moles of each species. The
@ -244,7 +249,9 @@ namespace Cantera {
// get the array of non-dimensional Gibbs functions for the pure
// species
m_mix->getStandardChemPotentials(m_mu.begin());
//m_mix->getStandardChemPotentials(m_mu.begin());
double not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.begin(), true);
int kpp = 0;
index_t k, q;
@ -280,9 +287,6 @@ namespace Cantera {
for (int k = 0; k < int(m_nsp); k++) {
if (ip == ksp) {
m_moles[k] = aa(n+1, 0);
//if (ISNAN(m_moles[k])) {
// writelog("moles "+int2str(k)+" is nan 3. \n");
//}
}
ksp++;
}
@ -486,7 +490,6 @@ namespace Cantera {
for (m = 0; m < m_nel; m++) {
ik = m_order[m];
k = m_species[ik];
addLogEntry("m, k, ik",int2str(m)+int2str(k)+int2str(ik));
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup();
@ -494,7 +497,6 @@ namespace Cantera {
for (m = m_nel; m < m_nsp; m++) {
ik = m_order[m];
k = m_species[ik];
addLogEntry("m, k, ik",int2str(m)+int2str(k)+int2str(ik));
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup();
@ -538,15 +540,6 @@ namespace Cantera {
k = m_order[ik];
m_lastmoles[k] = m_moles[k];
m_moles[k] += omega * deltaN[k];
//if (ISNAN(omega)) {
// writelog("omega is nan\n");
//}
//if (ISNAN(deltaN[k])) {
// writelog("deltaN["+int2str(k)=" is nan\n");
//}
//if (ISNAN(m_moles[k])) {
// writelog("moles "+int2str(k)+" is nan 4. \n");
//}
}
for (ik = m_nel; ik < m_nsp; ik++) {
@ -554,16 +547,9 @@ namespace Cantera {
m_lastmoles[k] = m_moles[k];
if (m_majorsp[k]) {
m_moles[k] += omega * deltaN[k];
//if (ISNAN(m_moles[k])) {
// writelog("moles "+int2str(k)+" is nan 5. \n");
//}
//if (m_moles[k] < 0.0) m_moles[k] = 0.0;
}
else {
m_moles[k] = fabs(m_moles[k])*fminn(10.0, exp(-m_deltaG_RT[ik - m_nel]));
//if (ISNAN(m_moles[k])) {
// writelog("moles "+int2str(k)+" is nan 6. \n");
//}
}
}
setMoles();
@ -580,14 +566,7 @@ namespace Cantera {
index_t ik, j, k = 0;
doublereal grad0 = computeReactionSteps(m_dxi);
//if (grad0 > 0.0) {
//cout << *m_mix << endl;
// cout << "gradient = " << grad0 << endl;
// throw CanteraError("stepComposition", "positive gradient!");
//}
// compute mole the fraction changes.
// compute the mole fraction changes.
//multiply(m_N, dxi.begin(), m_work.begin());
for (ik = 0; ik < m_nsp; ik++) {
m_work[ik] = 0.0;
@ -680,7 +659,7 @@ namespace Cantera {
// current direction. If it is positive, then we have overshot
// the minimum. In this case, interpolate back.
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.begin());
m_mix->getValidChemPotentials(not_mu, m_mu.begin());
doublereal grad1 = 0.0;
for (k = 0; k < m_nsp; k++) {
grad1 += m_work[k] * m_mu[m_species[k]];
@ -780,10 +759,10 @@ namespace Cantera {
fctr = 1.0;
else
fctr = 1.0/(term1 + csum + sum);
if (fctr < -999.0 || fctr > 999.0) {
cout << "fctr, term1, csum, sum = " << fctr << " " << term1 << " " << csum << " " << sum << endl;
cout << reactionString(j) << endl;
}
//if (fctr < -999.0 || fctr > 999.0) {
// cout << "fctr, term1, csum, sum = " << fctr << " " << term1 << " " << csum << " " << sum << endl;
// cout << reactionString(j) << endl;
//}
}
dxi[j] = -fctr*dg_rt;
index_t m;
@ -838,24 +817,30 @@ namespace Cantera {
}
doublereal MultiPhaseEquil::error() {
index_t j, ik, k, maxj;
bool exists = false;
index_t j, ik, k;
doublereal err, maxerr = 0.0;
// examine every reaction
for (j = 0; j < m_nsp - m_nel; j++) {
ik = j + m_nel;
k = m_order[ik];
if (m_dsoln[k] > 0 && fabs(m_moles[k]) <= SmallNumber) err = 0.0;
else if (m_dsoln[k] == 0 && m_moles[k] <= 0.0) {
if (m_deltaG_RT[j] >= 0.0) err = 0.0;
else err = fabs(m_deltaG_RT[j]);//1.0;
}
// don't require formation reactions for solution species
// present in trace amounts to be equilibrated
if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber)
err = 0.0;
// for stoichiometric phase species, no error if not present and
// delta G for the formation reaction is positive
else if (isStoichPhase(ik) && moles(ik) <= 0.0 &&
m_deltaG_RT[j] >= 0.0) err = 0.0;
//else err = fabs(m_deltaG_RT[j]);
//}
else {
exists = true;
err = fabs(m_deltaG_RT[j]);
}
if (err > maxerr) {
maxerr = err;
maxj = j;
}
}

16
Cantera/src/MultiPhaseEquil.h
View file

@ -42,6 +42,7 @@ namespace Cantera {
string reactionString(index_t j);
doublereal error();
void printInfo();
void finish();
protected:
@ -55,6 +56,12 @@ namespace Cantera {
doublereal computeReactionSteps(vector_fp& dxi);
void setMoles();
// moles of the species with sorted index ns
double moles(int ns) const { return m_moles[m_order[ns]]; }
double& moles(int ns) { return m_moles[m_order[ns]]; }
int solutionSpecies(int n) const { return m_dsoln[m_order[n]]; }
bool isStoichPhase(int n) const { return (m_dsoln[m_order[n]] == 0); }
index_t m_nel_mix, m_nsp_mix, m_np;
index_t m_nel, m_nsp;
index_t m_eloc;
@ -63,7 +70,7 @@ namespace Cantera {
doublereal m_press, m_temp;
vector_int m_order;
matrix_t m_N, m_A;
vector_fp m_work, m_work2;
vector_fp m_work, m_work2, m_work3;
vector_fp m_moles, m_lastmoles, m_dxi;
vector_fp m_deltaG_RT, m_mu;
vector<bool> m_majorsp;
@ -71,7 +78,12 @@ namespace Cantera {
vector_int m_lastsort;
vector_int m_dsoln;
vector_int m_incl_element, m_incl_species;
vector_int m_species, m_element;
// Vector of indices for species that are included in the
// calculation. This is used to exclude pure-phase species
// with invalid thermo data
vector_int m_species;
vector_int m_element;
vector<bool> m_solnrxn;
bool m_force;
};

1
Cantera/src/StoichManager.h
View file

@ -304,7 +304,6 @@ namespace Cantera {
m_ic[n] = ic[n];
m_order[n] = order[n];
m_stoich[n] = stoich[n];
cout << "n, stoich[n] = " << n << " " << stoich[n] << endl;
}
}

3
Cantera/src/ThermoPhase.h
View file

@ -579,9 +579,6 @@ namespace Cantera {
* Specific entropy. Units: J/kg/K.
*/
doublereal entropy_mass() const {
//cout << "entropy_mass. " << endl;
//cout << "entropy_mole = " << entropy_mole() << endl;
//cout << "meanMolecularWeight = " << meanMolecularWeight()<< endl;
return entropy_mole()/meanMolecularWeight();
}

6
Cantera/src/oneD/Domain1D.h
View file

@ -87,8 +87,9 @@ namespace Cantera {
}
/**
* Initialize. Base class method does nothing, but may be
* overloaded.
* Initialize. This method is called by OneDim::init() for
* each domain once at the beginning of a simulation. Base
* class method does nothing, but may be overloaded.
*/
virtual void init(){ }
@ -425,6 +426,7 @@ namespace Cantera {
doublereal grid(int point) { return m_z[point]; }
virtual void setupGrid(int n, const doublereal* z) {}
void setGrid(int n, const doublereal* z);
/**
* Writes some or all initial solution values into the global

9
Cantera/src/oneD/OneDim.cpp
View file

@ -447,4 +447,13 @@ namespace Cantera {
s.close();
writelog("Solution saved to file "+fname+" as solution "+id+".\n");
}
void Domain1D::setGrid(int n, const doublereal* z) {
m_z.resize(n);
m_points = n;
int j;
for (j = 0; j < m_points; j++) m_z[j] = z[j];
}
}

4
Cantera/src/oneD/StFlow.cpp
View file

@ -201,7 +201,7 @@ namespace Cantera {
/**
* Change the grid size. Called after grid refinement.
*/
void StFlow::resize(int points) {
void StFlow::resize(int points) {
Domain1D::resize(m_nv, points);
m_rho.resize(m_points, 0.0);
@ -1044,8 +1044,6 @@ namespace Cantera {
writelog("Grid contains "+int2str(np)+
" points.\n");
readgrid = true;
// note that setupGrid also resizes the domain.
setupGrid(np, x.begin());
}
}

4
tools/src/package4mac.in
View file

@ -5,8 +5,8 @@ INSTALL='@INSTALL@'
PYVERSION=2.3
PKGDIR=$HOME/Packages
ctname=Cantera2
mixname=MixMaster2
ctname=Cantera
mixname=MixMaster
CTDIR=$PKGDIR/Cantera/root_dir/Applications/$ctname
PYDIR=$PKGDIR/Cantera/root_dir/Library/Python/$PYVERSION

2
tools/templates/f77/demo.mak.in
View file

@ -25,7 +25,7 @@ OBJS = demo.o demo_ftnlib.o
# additional flags to be passed to the linker. If your program
# requires other external libraries, put them here
LINK_OPTIONS = @LCXX_FLAGS@
LINK_OPTIONS = @LCXX_FLAGS@ @EXTRA_LINK@
#---------------------------------------------------------------------------

2
tools/templates/f90/demo.mak.in
View file

@ -21,7 +21,7 @@ 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@
LINK_OPTIONS = @LCXX_FLAGS@ @EXTRA_LINK@
#---------------------------------------------------------------------------
# You probably don't need to edit anything below.