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

improvements to equilibrate

Browse source
This commit is contained in:
Dave Goodwin 2005-07-14 18:44:55 +00:00
parent 791a00ed1e
commit 673513edae
3 changed files with 360 additions and 129 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
Cantera/src/ChemEquil.h
View file

@ -159,6 +159,7 @@ namespace Cantera {
if (solver > 0) {
MultiPhase mix;
mix.addPhase(&s, 1.0);
mix.init();
mix.setTemperature(s.temperature());
mix.setPressure(s.pressure());
equilibrate(mix, XY, rtol, maxsteps, loglevel);

479
Cantera/src/MultiPhase.cpp
View file

@ -1,11 +1,3 @@
/**
* @file MultiPhase.cpp
*
* $Author$
* $Date$
* $Revision$
*/
#include "MultiPhase.h"
#include "MultiPhaseEquil.h"
@ -15,16 +7,24 @@
namespace Cantera {
/// Constructor.
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) {
m_Tmin(1.0), m_Tmax(100000.0) {
}
void MultiPhase::
addPhase(phase_t* p, doublereal moles) {
addPhases(phase_list& phases, const vector_fp& phaseMoles) {
index_t np = phases.size();
index_t n;
for (n = 0; n < np; n++) {
addPhase(phases[n], phaseMoles[n]);
}
init();
}
void MultiPhase::
addPhase(phase_t* p, doublereal moles) {
if (m_init) {
throw CanteraError("addPhase",
"phases cannot be added after init() has been called.");
@ -56,33 +56,47 @@ namespace Cantera {
// it is a new element. In this case, add the name
// to the list of names, increment the element count,
// and add an entry to the name->(index+1) map.
if (m_enamemap[ename] == 0) {
if (m_enamemap.find(ename) == m_enamemap.end()) {
m_enamemap[ename] = m_nel + 1;
m_enames.push_back(ename);
m_atomicNumber.push_back(p->atomicNumber(m));
// Element 'E' (or 'e') is special. Note its location.
if (ename == "E" || ename == "e") m_eloc = m_nel;
m_nel++;
}
}
// If the mixture temperature hasn't been set, then set the
// temperature and pressure to the values for the phase being
// added.
if (m_temp == 0.0 && p->temperature() > 0.0) {
m_temp = p->temperature();
m_press = p->pressure();
}
//cout << "min, max = " << m_Tmin << " " << m_Tmax << endl;
// If this is a solution phase, update the minimum and maximum
// mixture temperatures. Stoichiometric phases are excluded,
// since a mixture may define multiple stoichiometric phases,
// each of which has thermo data valid only over a limited
// range. For example, a mixture might be defined to contain a
// phase representing water ice and one representing liquid
// water, only one of which should be present if the mixture
// represents an equilibrium state.
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;
}
}
/// Process phases and build atomic composition array. After
/// init() has been called, no more phases may be added.
/// Process phases and build atomic composition array. This method
/// must be called after all phases are added, before doing
/// anything else with the mixture. After init() has been called,
/// no more phases may be added.
void MultiPhase::init() {
if (m_init) return;
index_t ip, kp, k = 0, nsp, m;
@ -149,7 +163,6 @@ namespace Cantera {
/// Moles of species \c k.
doublereal MultiPhase::speciesMoles(index_t k) {
if (!m_init) init();
index_t ip = m_spphase[k];
return m_moles[ip]*m_moleFractions[k];
}
@ -159,7 +172,6 @@ namespace Cantera {
doublereal MultiPhase::elementMoles(index_t m) {
doublereal sum = 0.0, phasesum;
index_t i, k = 0, ik, nsp;
if (!m_init) init();
for (i = 0; i < m_np; i++) {
phasesum = 0.0;
nsp = m_phase[i]->nSpecies();
@ -182,7 +194,10 @@ namespace Cantera {
return sum;
}
/// Charge of one phase
/// Net charge of one phase (Coulombs). The net charge is computed as
/// \f[ Q_p = N_p \sum_k F z_k X_k \f]
/// where the sum runs only over species in phase \a p.
/// @param p index of the phase for which the charge is desired.
doublereal MultiPhase::phaseCharge(index_t p) {
doublereal phasesum = 0.0;
int ik, k, nsp = m_phase[p]->nSpecies();
@ -194,6 +209,7 @@ namespace Cantera {
}
/// Get the chemical potentials of all species in all phases.
void MultiPhase::getChemPotentials(doublereal* mu) {
index_t i, loc = 0;
updatePhases();
@ -203,10 +219,39 @@ namespace Cantera {
}
}
/// Get chemical potentials of species with valid thermo
/// data. This method is designed for use in computing chemical
/// equilibrium by Gibbs minimization. For solution phases (more
/// than one species), this does the same thing as
/// getChemPotentials. But for stoichiometric phases, this writes
/// into array \a mu the user-specified value \a not_mu instead of
/// the chemical potential if the temperature is outside the range
/// for which the thermo data for the one species in the phase are
/// valid. The need for this arises since many condensed phases
/// have thermo data fit only for the temperature range for which
/// they are stable. For example, in the NASA database, the fits
/// for H2O(s) are only done up to 0 C, the fits for H2O(L) are
/// only done from 0 C to 100 C, etc. Using the polynomial fits outside
/// the range for which the fits were done can result in spurious
/// chemical potentials, and can lead to condensed phases
/// appearing when in fact they should be absent.
///
/// By setting \a not_mu to a large positive value, it is possible
/// to force routines which seek to minimize the Gibbs free energy
/// of the mixture to zero out any phases outside the temperature
/// range for which their thermo data are valid.
///
/// If this method is called with \a standard set to true, then
/// the composition-independent standard chemical potentials are
/// returned instead of the composition-dependent chemical
/// potentials.
///
void MultiPhase::getValidChemPotentials(doublereal not_mu,
doublereal* mu, bool standard) {
index_t i, loc = 0;
updatePhases();
updatePhases();
// iterate over the phases
for (i = 0; i < m_np; i++) {
if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
if (!standard)
@ -220,7 +265,7 @@ namespace Cantera {
}
}
/// True if species \a k belongs to a solution phase.
bool MultiPhase::solutionSpecies(index_t k) {
if (m_phase[m_spphase[k]]->nSpecies() > 1)
return true;
@ -228,6 +273,7 @@ namespace Cantera {
return false;
}
/// The Gibbs free energy of the mixture (J).
doublereal MultiPhase::gibbs() const {
index_t i;
doublereal sum = 0.0;
@ -237,6 +283,7 @@ namespace Cantera {
return sum;
}
/// The enthalpy of the mixture (J).
doublereal MultiPhase::enthalpy() const {
index_t i;
doublereal sum = 0.0;
@ -246,6 +293,7 @@ namespace Cantera {
return sum;
}
/// The entropy of the mixture (J/K).
doublereal MultiPhase::entropy() const {
index_t i;
doublereal sum = 0.0;
@ -255,6 +303,9 @@ namespace Cantera {
return sum;
}
/// The specific heat at constant pressure and composition (J/K).
/// Note that this does not account for changes in composition of
/// the mixture with temperature.
doublereal MultiPhase::cp() const {
index_t i;
doublereal sum = 0.0;
@ -264,45 +315,52 @@ namespace Cantera {
return sum;
}
void MultiPhase::updateMoleFractions() {
if (!m_init) init();
// save the current mole fractions for each phase
index_t ip, loc = 0;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
p->getMoleFractions(m_moleFractions.begin() + loc);
loc += p->nSpecies();
}
}
/// Set the mole fractions of phase \a n to the values in
/// array \a x.
void MultiPhase::setPhaseMoleFractions(index_t n, doublereal* x) {
phase_t* p = m_phase[n];
p->setState_TPX(m_temp, m_press, x);
}
/// Set the species moles using a map. The map \a xMap maps
/// species name strings to mole numbers. Mole numbers that are
/// less than or equal to zero will be set to zero.
void MultiPhase::setMolesByName(compositionMap& xMap) {
if (!m_init) init();
int kk = nSpecies();
doublereal x;
vector_fp mf(kk, 0.0);
vector_fp moles(kk, 0.0);
for (int k = 0; k < kk; k++) {
x = xMap[speciesName(k)];
if (x > 0.0) mf[k] = x;
if (x > 0.0) moles[k] = x;
}
setMoles(mf.begin());
setMoles(moles.begin());
}
/// Set the species moles using a string. Unspecified species are
/// set to zero.
void MultiPhase::setMolesByName(const string& x) {
compositionMap xx;
if (!m_init) init();
// add an entry in the map for every species, with value -1.0.
// Function parseCompString (stringUtils.cpp) uses the names
// in the map to specify the allowed species.
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
// build the composition map from the string, and then set the
// moles.
parseCompString(x, xx);
setMolesByName(xx);
}
/// Set the species moles to the values in array \a n. The state
/// of each phase object is also updated to have the specified
/// composition and the mixture temperature and pressure.
void MultiPhase::setMoles(doublereal* n) {
if (!m_init) init();
index_t ip, loc = 0;
@ -328,6 +386,7 @@ namespace Cantera {
}
}
/// The total mixture volume [m^3].
doublereal MultiPhase::volume() {
int i;
doublereal sum = 0;
@ -337,30 +396,6 @@ namespace Cantera {
return sum;
}
/**
* updatePhases() const:
*
* This routine take the state information storred in
* this current object, and sets the state in the underlying
* ThermoPhase objects.
* Currently, this state information is the Temperature
* Pressure, and the mole fraction vector. Temperature and
* pressure are assumed uniform for all phases.
*/
void MultiPhase::updatePhases() const {
//if (!m_init) init();
index_t p, nsp, loc = 0;
for (p = 0; p < m_np; p++) {
nsp = m_phase[p]->nSpecies();
const doublereal* x = m_moleFractions.begin() + loc;
loc += nsp;
m_phase[p]->setState_TPX(m_temp, m_press, x);
m_temp_OK[p] = true;
if (m_temp < m_phase[p]->minTemp()
|| m_temp > m_phase[p]->maxTemp()) m_temp_OK[p] = false;
}
}
doublereal MultiPhase::equilibrate(int XY, doublereal err,
int maxsteps, int maxiter, int loglevel) {
doublereal error;
@ -370,9 +405,11 @@ namespace Cantera {
int n;
bool start, once;
doublereal ferr, hnow, herr = 1.0;
doublereal snow, serr = 1.0, s0;
doublereal Tlow = -1.0, Thigh = -1.0;
doublereal hlow = 0.0, hhigh = 0.0, tnew;
doublereal dta, dtmax;
doublereal Hlow = Undef, Hhigh = Undef, tnew;
doublereal dta, dtmax, cpb;
MultiPhaseEquil* e = 0;
if (!m_init) init();
if (loglevel > 0) {
@ -383,59 +420,177 @@ namespace Cantera {
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();
e = new MultiPhaseEquil(this);
try {
error = e->equilibrate(XY, err, maxsteps, loglevel-1);
}
catch (CanteraError err) {
if (loglevel > 0) {
endLogGroup();
//write_logfile("equil_err.html");
}
delete e;
e = 0;
throw err;
}
goto done;
}
else if (XY == HP) {
h0 = enthalpy();
start = true;
Tlow = m_Tmin; // lower bound on T
Thigh = m_Tmax; // upper bound on T
hlow = 0.0;
hhigh = 0.0;
once = true;
Tlow = 0.5*m_Tmin; // lower bound on T
Thigh = 2.0*m_Tmax; // upper bound on T
if (loglevel > 0) {
addLogEntry("problem type","fixed H,P");
addLogEntry("H target",fp2str(h0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
//addLogEntry("min T",fp2str(Tlow));
//addLogEntry("max T",fp2str(Thigh));
}
for (n = 0; n < maxiter; n++) {
MultiPhaseEquil e(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
start = false;
// if 'strt' is false, the current composition will be used as
// the starting estimate; otherwise it will be estimated
e = new MultiPhaseEquil(this, strt);
// start with a loose error tolerance, but tighten it as we get
// close to the final temperature
if (loglevel > 0) {
beginLogGroup("iteration "+int2str(n));
}
try {
error = e.equilibrate(TP, err, maxsteps, loglevel-1);
error = e->equilibrate(TP, err, maxsteps, loglevel-1);
hnow = enthalpy();
// the equilibrium enthalpy monotonically increases with T;
// if the current value is below the target, the we know the
// current temperature is too low. Set
if (hnow < h0) {
if (m_temp > Tlow) Tlow = m_temp;
if (m_temp > Tlow) {
Tlow = m_temp;
Hlow = hnow;
}
}
// the current enthalpy is greater than the target; therefore the
// current temperature is too high.
else {
if (m_temp < Thigh) {
Thigh = m_temp;
Hhigh = hnow;
}
}
if (Hlow != Undef && Hhigh != Undef) {
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
dt = (h0 - hnow)/cpb;
dta = fabs(dt);
dtmax = 0.5*fabs(Thigh - Tlow);
if (dta > dtmax) dt *= dtmax/dta;
}
else {
if (m_temp < Thigh) Thigh = m_temp;
tnew = sqrt(Tlow*Thigh);
dt = tnew - m_temp;
//cpb = cp();
}
herr = fabs((h0 - hnow)/h0);
if (loglevel > 0) {
addLogEntry("T",fp2str(temperature()));
addLogEntry("H",fp2str(hnow));
addLogEntry("H rel error",fp2str(herr));
addLogEntry("lower T bound",fp2str(Tlow));
addLogEntry("upper T bound",fp2str(Thigh));
endLogGroup();
}
dt = (h0 - hnow)/cp();
dtmax = 0.5*(Thigh - Tlow);
if (herr < err) { // || dta < 1.0e-4) {
if (loglevel > 0) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("H rel error",fp2str(herr));
}
goto done;
}
tnew = m_temp + dt;
if (tnew < 0.0) tnew = 0.5*m_temp;
//dta = fabs(tnew - m_temp);
setTemperature(tnew);
// if the size of Delta T is not too large, use
// the current composition as the starting estimate
if (dta < 100.0) strt = false;
}
catch (CanteraError err) {
if (!strt) {
if (loglevel > 0)
addLogEntry("no convergence",
"try estimating starting composition");
strt = true;
}
else {
tnew = 0.5*(Tlow + Thigh);
setTemperature(tnew);
if (loglevel > 0)
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
endLogGroup();
}
}
delete e;
e = 0;
if (loglevel > 0) {
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
}
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
else if (XY == SP) {
writelog("SP\n");
s0 = entropy();
start = true;
Tlow = 1.0; // m_Tmin; // lower bound on T
Thigh = 1.0e6; // m_Tmax; // upper bound on T
if (loglevel > 0) {
addLogEntry("problem type","fixed S,P");
addLogEntry("S target",fp2str(s0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
}
for (n = 0; n < maxiter; n++) {
e = new MultiPhaseEquil(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
//start = false;
if (loglevel > 0) {
beginLogGroup("iteration "+int2str(n));
}
try {
error = e->equilibrate(TP, err, maxsteps, loglevel-1);
snow = entropy();
if (snow < s0) {
if (m_temp > Tlow) Tlow = m_temp;
}
else {
if (m_temp < Thigh) Thigh = m_temp;
}
serr = fabs((s0 - snow)/s0);
if (loglevel > 0) {
addLogEntry("T",fp2str(temperature()));
addLogEntry("S",fp2str(snow));
addLogEntry("S rel error",fp2str(serr));
endLogGroup();
}
dt = (s0 - snow)*m_temp/cp();
dtmax = 0.5*fabs(Thigh - Tlow);
dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
dta = fabs(dt);
if (dta > dtmax) dt *= dtmax/dta;
if (herr < err || dta < 1.0e-4) {
if (loglevel > 0) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("H rel error",fp2str(herr));
addLogEntry("S rel error",fp2str(serr));
}
goto done;
}
@ -446,7 +601,8 @@ namespace Cantera {
// the current composition as the starting estimate
if (dta < 100.0) strt = false;
}
catch (CanteraError e) {
catch (CanteraError err) {
if (!strt) {
if (loglevel > 0)
addLogEntry("no convergence",
@ -461,58 +617,89 @@ namespace Cantera {
"trying T = "+fp2str(m_temp));
}
endLogGroup();
}
}
delete e;
e = 0;
if (loglevel > 0) {
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
}
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
else if (XY == SP) {
if (loglevel > 0) {
addLogEntry("problem type","fixed S,P");
}
doublereal dt = 1.0e3;
doublereal s0 = entropy();
int n;
bool start = true;
doublereal ferr, snow, serr, tnew;
for (n = 0; n < maxiter; n++) {
MultiPhaseEquil e(this, start);
ferr = 0.1;
start = false;
if (fabs(dt) < 1.0) ferr = err;
if (loglevel > 1) {
beginLogGroup("iteration "+int2str(n));
}
error = e.equilibrate(TP, ferr, maxsteps, loglevel-1);
snow = entropy();
tnew = exp(0.5*(s0 - snow)/cp())*temperature();
serr = fabs((s0 - snow)/s0);
if (loglevel > 1) {
addLogEntry("T",fp2str(temperature()));
addLogEntry("S rel error",fp2str(serr));
endLogGroup();
}
if (serr < err) {
if (loglevel > 0) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("S rel error",fp2str(serr));
}
goto done;
}
setTemperature(tnew);
}
}
// else if (XY == SP) {
// if (loglevel > 0) {
// addLogEntry("problem type","fixed S,P");
// }
// doublereal dt = 1.0e3;
// doublereal s0 = entropy();
// int n;
// bool start = true;
// doublereal ferr, snow, serr, tnew;
// for (n = 0; n < maxiter; n++) {
// e = new MultiPhaseEquil(this, start);
// ferr = 0.1;
// start = false;
// if (fabs(dt) < 1.0) ferr = err;
// if (loglevel > 1) {
// beginLogGroup("iteration "+int2str(n));
// }
// try {
// error = e->equilibrate(TP, ferr, maxsteps, loglevel-1);
// snow = entropy();
// tnew = exp(0.5*(s0 - snow)/cp())*temperature();
// serr = fabs((s0 - snow)/s0);
// if (loglevel > 1) {
// addLogEntry("T",fp2str(temperature()));
// addLogEntry("S rel error",fp2str(serr));
// endLogGroup();
// }
// if (serr < err) {
// if (loglevel > 0) {
// addLogEntry("T iterations",int2str(n));
// addLogEntry("Final T",fp2str(temperature()));
// addLogEntry("S rel error",fp2str(serr));
// }
// goto done;
// }
// setTemperature(tnew);
// }
// catch (CanteraError err) {
// delete e;
// if (!strt) {
// if (loglevel > 0)
// addLogEntry("no convergence",
// "setting strt to True");
// strt = true;
// }
// else {
// tnew = 0.5*(m_temp + Thigh);
// setTemperature(tnew);
// if (loglevel > 0)
// addLogEntry("no convergence",
// "trying T = "+fp2str(m_temp));
// }
// }
// endLogGroup();
// }
// if (loglevel > 0) write_logfile("equil_err.html");
// throw CanteraError("MultiPhase::equilibrate",
// "No convergence for T");
// }
else if (XY == TV) {
if (loglevel > 0) {
addLogEntry("problem type","fixed T, V");
}
//doublereal dt = 1.0e3;
doublereal dt = 1.0e3;
doublereal v0 = volume();
doublereal dVdP;
int n;
bool start = true;
doublereal error, vnow, pnow, verr;
doublereal error, ferr, vnow, pnow, verr, tnew;
for (n = 0; n < maxiter; n++) {
pnow = pressure();
MultiPhaseEquil e(this, start);
@ -549,10 +736,50 @@ namespace Cantera {
}
return -1.0;
done:
delete e;
e = 0;
if (loglevel > 0) {
endLogGroup();
}
return err;
}
//-------------------------------------------------------------
//
// protected methods
//
//-------------------------------------------------------------
/// Update the locally-stored species mole fractions.
void MultiPhase::updateMoleFractions() {
index_t ip, loc = 0;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
p->getMoleFractions(m_moleFractions.begin() + loc);
loc += p->nSpecies();
}
}
/// synchronize the phase objects with the mixture state. This
/// method sets each phase to the mixture temperature and
/// pressure, and sets the phase mole fractions based on the
/// mixture mole numbers.
void MultiPhase::updatePhases() const {
index_t p, nsp, loc = 0;
for (p = 0; p < m_np; p++) {
nsp = m_phase[p]->nSpecies();
const doublereal* x = m_moleFractions.begin() + loc;
loc += nsp;
m_phase[p]->setState_TPX(m_temp, m_press, x);
m_temp_OK[p] = true;
if (m_temp < m_phase[p]->minTemp()
|| m_temp > m_phase[p]->maxTemp()) m_temp_OK[p] = false;
}
}
}

9
Cantera/src/MultiPhase.h
View file

@ -7,8 +7,6 @@
namespace Cantera {
class MultiPhaseEquil;
/// A class for multiphase mixtures. The mixture can contain any
/// number of phases of any type. All phases have the same
/// temperature and pressure, and a specified number of moles.
@ -22,9 +20,11 @@ namespace Cantera {
public:
// some typedefs for convenience
typedef size_t index_t;
typedef ThermoPhase phase_t;
typedef DenseMatrix array_t;
typedef vector<phase_t*> phase_list;
/// Constructor. The constructor takes no arguments, since
/// phases are added using method addPhase.
@ -35,6 +35,8 @@ namespace Cantera {
/// phase objects.
virtual ~MultiPhase() {}
void addPhases(phase_list& phases, const vector_fp& phaseMoles);
/// 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
@ -215,6 +217,8 @@ namespace Cantera {
protected:
// These methods are meant for internal use.
/// update the locally-stored composition to match the current
/// compositions of the phase objects.
void updateMoleFractions();
@ -242,7 +246,6 @@ namespace Cantera {
bool m_init;
int m_eloc;
mutable vector<bool> m_temp_OK;
MultiPhaseEquil* m_equil;
doublereal m_Tmin, m_Tmax;
};