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-27 14:51:05 +00:00
parent dfa29b0f70
commit 88d45066c7
4 changed files with 174 additions and 212 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

50
Cantera/src/MultiPhase.cpp
View file

@ -8,13 +8,12 @@
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) {
}
void MultiPhase::
addPhase(phase_t* p, doublereal moles) {
@ -187,8 +186,6 @@ namespace Cantera {
}
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
void MultiPhase::getChemPotentials(doublereal* mu) {
index_t i, loc = 0;
updatePhases();
@ -198,8 +195,6 @@ 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, bool standard) {
index_t i, loc = 0;
@ -218,17 +213,6 @@ namespace Cantera {
}
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
void MultiPhase::getStandardChemPotentials(doublereal* mu) {
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
m_phase[i]->getStandardChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
}
}
bool MultiPhase::solutionSpecies(index_t k) {
if (m_phase[m_spphase[k]]->nSpecies() > 1)
return true;
@ -236,7 +220,7 @@ namespace Cantera {
return false;
}
doublereal MultiPhase::gibbs() {
doublereal MultiPhase::gibbs() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
@ -245,7 +229,7 @@ namespace Cantera {
return sum;
}
doublereal MultiPhase::enthalpy() {
doublereal MultiPhase::enthalpy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
@ -254,7 +238,7 @@ namespace Cantera {
return sum;
}
doublereal MultiPhase::entropy() {
doublereal MultiPhase::entropy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
@ -263,7 +247,7 @@ namespace Cantera {
return sum;
}
doublereal MultiPhase::cp() {
doublereal MultiPhase::cp() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
@ -345,12 +329,12 @@ namespace Cantera {
return sum;
}
void MultiPhase::updatePhases() {
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();
doublereal* x = m_moleFractions.begin() + loc;
const doublereal* x = m_moleFractions.begin() + loc;
loc += nsp;
m_phase[p]->setState_TPX(m_temp, m_press, x);
m_temp_OK[p] = true;
@ -383,11 +367,10 @@ namespace Cantera {
// create an equilibrium manager
MultiPhaseEquil e(this);
error = e.equilibrate(XY, err, maxsteps, loglevel-1);
if (loglevel > 0) e.printInfo();
// if (loglevel > 0) e.printInfo();
goto done;
}
else if (XY == HP) {
dt = 1.0e2;
h0 = enthalpy();
start = true;
Tlow = m_Tmin; // lower bound on T
@ -401,9 +384,10 @@ namespace Cantera {
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
}
ferr = 0.1;
for (n = 0; n < maxiter; n++) {
MultiPhaseEquil e(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
start = false;
if (loglevel > 0) {
beginLogGroup("iteration "+int2str(n));
@ -413,14 +397,10 @@ namespace Cantera {
hnow = enthalpy();
if (hnow < h0) {
if (m_temp > Tlow) {
Tlow = m_temp;
}
if (m_temp > Tlow) Tlow = m_temp;
}
else {
if (m_temp < Thigh) {
Thigh = m_temp;
}
if (m_temp < Thigh) Thigh = m_temp;
}
herr = fabs((h0 - hnow)/h0);
if (loglevel > 0) {
@ -443,12 +423,16 @@ namespace Cantera {
}
tnew = m_temp + dt;
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 e) {
if (!strt) {
if (loglevel > 0)
addLogEntry("no convergence","setting strt to True");
addLogEntry("no convergence",
"setting strt to True");
strt = true;
}
else {

98
Cantera/src/MultiPhase.h
View file

@ -22,6 +22,10 @@ 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();
@ -31,22 +35,27 @@ namespace Cantera {
/// 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
void addPhase(phase_t* p, doublereal moles);
/// Number of elements.
int nElements() { return int(m_nel); }
/// Name of element \a m.
string elementName(int m) { return m_enames[m]; }
/// Index of element with name \a name.
int elementIndex(string name) { return m_enamemap[name] - 1;}
/// Number of species, summed over all phases.
int nSpecies() { return int(m_nsp); }
/// Name of species with index \a k.
string speciesName(int k) { return m_snames[k]; }
/// Number of atoms of element \a m in species \a k.
doublereal nAtoms(int k, int m) {
if (!m_init) init();
return m_atoms(m,k);
@ -87,36 +96,56 @@ namespace Cantera {
return m_spstart[p] + k;
}
/// Minimum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal minTemp();
/// Maximum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal maxTemp();
/// Total charge (Coulombs).
doublereal charge();
/// Charge (Coulombs) of phase with index \a p.
doublereal phaseCharge(index_t p);
/// Total moles of element m, summed over all
/// phases
/// Total moles of element \a m, summed over all phases.
doublereal elementMoles(index_t m);
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
/// Chemical potentials. Write into array \a mu the chemical
/// potentials of all species [J/kmol]. The chemical
/// potentials are related to the activities by
/// \f[ \mu_k = \mu_k^0(T, P) + RT \ln a_k. \f].
void getChemPotentials(doublereal* mu);
/// Valid chemical potentials. Write into array \c mu the
/// Valid chemical potentials. Write into array \a mu the
/// 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.
/// set the chemical potential to the value \a not_mu. If \a
/// standard is set to true, then the values returned are
/// standard chemical potentials.
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].
void getStandardChemPotentials(doublereal* mu);
/// Temperature [K].
doublereal temperature() {
return m_temp;
}
doublereal temperature() { return m_temp; }
/// Set the mixture to a state of chemical equilibrium.
/// @param XY Integer flag specifying properties to hold fixed.
/// @param err Error tolerance for \f$\Delta \mu/RT \f$ for
/// all reactions. Also used as the relative error tolerance
/// for the outer loop.
/// @param maxsteps Maximum number of steps to take in solving
/// the fixed TP problem.
/// @param maxiter Maximum number of "outer" iterations for
/// problems holding fixed something other than (T,P).
/// @param loglevel Level of diagnostic output, written to a
/// file in HTML format.
doublereal equilibrate(int XY, doublereal err = 1.0e-9,
int maxsteps = 1000, int maxiter = 200, int loglevel = 0);
@ -127,26 +156,39 @@ namespace Cantera {
updatePhases();
}
/// Pressure [Pa].
doublereal pressure() {
return m_press;
}
/// Volume [m^3].
doublereal volume();
/// Set the pressure [Pa].
void setPressure(doublereal P) {
m_press = P;
updatePhases();
}
doublereal enthalpy();
doublereal entropy();
doublereal gibbs();
doublereal cp();
/// Enthalpy [J].
doublereal enthalpy() const;
/// Entropy [J/K].
doublereal entropy() const;
/// Gibbs function [J].
doublereal gibbs() const;
/// Heat capacity at constant pressure [J/K].
doublereal cp() const;
/// Number of phases.
index_t nPhases() {
return m_np;
}
/// Return true is species \a k is a species in a
/// multicomponent solution phase.
bool solutionSpecies(index_t k);
index_t speciesPhaseIndex(index_t k) {
@ -157,8 +199,6 @@ namespace Cantera {
return m_moleFractions[k];
}
void updateMoleFractions();
void setPhaseMoleFractions(index_t n, doublereal* x);
void setMolesByName(compositionMap& xMap);
@ -167,16 +207,22 @@ namespace Cantera {
void setMoles(doublereal* n);
/// Return true if the phase \a p has valid thermo data for
/// the current temperature.
bool tempOK(index_t p) {
return m_temp_OK[p];
}
protected:
/// update the locally-stored composition to match the current
/// compositions of the phase objects.
void updateMoleFractions();
/// Set the states of the phase objects to the locally-stored
/// state. Note that if individual phases have T and P different
/// than that stored locally, the phase T and P will be modified.
void updatePhases();
void updatePhases() const;
vector_fp m_moles;
vector<phase_t*> m_phase;
@ -195,7 +241,7 @@ namespace Cantera {
index_t m_nsp;
bool m_init;
int m_eloc;
vector<bool> m_temp_OK;
mutable vector<bool> m_temp_OK;
MultiPhaseEquil* m_equil;
doublereal m_Tmin, m_Tmax;
};

228
Cantera/src/MultiPhaseEquil.cpp
View file

@ -1,7 +1,6 @@
#include "MultiPhaseEquil.h"
#include "MultiPhase.h"
#include "sort.h"
#include "recipes.h"
#include "global.h"
#include <math.h>
@ -19,8 +18,6 @@ using namespace std;
#endif
#endif
#undef DEBUG_MULTIPHASE_EQUIL
namespace Cantera {
const doublereal TINY = 1.0e-20;
@ -140,8 +137,9 @@ namespace Cantera {
m_N.resize(m_nsp, m_nsp - m_nel);
m_order.resize(m_nsp, 0);
if (start)
if (start) {
setInitialMoles();
}
computeN();
vector_fp dxi(m_nsp - m_nel, 1.0e-20);
@ -157,9 +155,10 @@ namespace Cantera {
m_dsoln.push_back(0);
}
m_force = false;
setMoles();
updateMixMoles();
}
doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
int maxsteps, int loglevel) {
int i;
@ -175,10 +174,6 @@ namespace Cantera {
endLogGroup();
}
if (loglevel > 2) printInfo();
//if (error() == 0.0) {
// write_logfile("equil_err.html");
// Cantera::error("stopping");
//}
if (error() < err) break;
}
if (i >= maxsteps) {
@ -202,8 +197,7 @@ namespace Cantera {
return error();
}
void MultiPhaseEquil::setMoles() {
//vector_fp n(m_nsp_mix, 0.0);
void MultiPhaseEquil::updateMixMoles() {
fill(m_work3.begin(), m_work3.end(), 0.0);
index_t k;
for (k = 0; k < m_nsp; k++) {
@ -212,10 +206,10 @@ namespace Cantera {
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.
/// Clean up the composition. 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;
@ -226,73 +220,69 @@ namespace Cantera {
}
/**
* Estimate the initial mole fractions. Uses the Simplex method
* to estimate the initial number of moles of each species. The
* linear Gibbs minimization problem is solved, neglecting the
* free energy of mixing terms. This procedure produces a good
* estimate of the low-temperature equilibrium composition.
*
* @param s phase object
* @param elementMoles vector of elemental moles
*/
/// Extimate the initial mole numbers. This is done by running
/// each reaction as far forward or backward as possible, subject
/// to the constraint that all mole numbers remain
/// non-negative. Reactions for which \f$ \Delta \mu^0 \f$ are
/// positive are run in reverse, and ones for which it is negative
/// are run in the forward direction. The end result is equivalent
/// to solving the linear programming problem of minimizing the
/// linear Gibbs function subject to the element and
/// non-negativity constraints.
int MultiPhaseEquil::setInitialMoles() {
index_t m, n;
doublereal lp = log(m_press/OneAtm);
DenseMatrix aa(m_nel+2, m_nsp+1, 0.0);
// first column contains fixed element moles
for (m = 0; m < m_nel; m++) {
aa(m+1,0) = m_mix->elementMoles(m_element[m]);
}
index_t m, n, ik, j;
// get the array of non-dimensional Gibbs functions for the pure
// species
//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;
doublereal rt = GasConstant * m_temp;
for (k = 0; k < m_nsp; k++) {
kpp++;
aa(0, kpp) = -m_mu[m_species[k]]/rt;
aa(0, kpp) -= m_dsoln[k]*lp; // ideal gas
for (q = 0; q < m_nel; q++)
aa(q+1, kpp) = -m_mix->nAtoms(m_species[k], m_element[q]);
}
doublereal dg_rt;
integer mp = m_nel+2; // parameters for SIMPLX
integer np = m_nsp+1;
integer m1 = 0;
integer m2 = 0;
integer m3 = m_nel;
integer icase=0;
integer nel = m_nel;
integer nsp = m_nsp;
vector_int iposv(m_nel);
vector_int izrov(m_nsp);
// solve the linear programming problem
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
bool redo = true;
int iter = 0;
while (redo) {
simplx_(&aa(0,0), &nel, &nsp, &mp, &np, &m1, &m2, &m3,
&icase, izrov.begin(), iposv.begin());
fill(m_moles.begin(), m_moles.end(), 0.0);
for (n = 0; n < m_nel; n++) {
int ksp = 0;
int ip = iposv[n] - 1;
for (int k = 0; k < int(m_nsp); k++) {
if (ip == ksp) {
m_moles[k] = aa(n+1, 0);
// choose a set of components based on the current
// composition
computeN();
redo = false;
iter++;
if (iter > 4) break;
// loop over all reactions
for (j = 0; j < m_nsp - m_nel; j++) {
dg_rt = 0.0;
dxi_min = 1.0e10;
for (ik = 0; ik < m_nsp; ik++) {
dg_rt += mu(ik) * m_N(ik,j);
}
// fwd or rev direction
idir = (dg_rt < 0.0 ? 1 : -1);
for (ik = 0; ik < m_nsp; ik++) {
nu = m_N(ik, j);
// set max change in progress variable by
// non-negativity requirement
if (nu*idir < 0) {
delta_xi = fabs(moles(ik)/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (delta_xi < SmallNumber && ik < m_nel) redo = true;
if (delta_xi < dxi_min) dxi_min = delta_xi;
}
}
// step the composition by dxi_min
for (ik = 0; ik < m_nsp; ik++) {
moles(ik) += m_N(ik, j) * idir*dxi_min;
}
ksp++;
}
// set the moles of the phase objects to match
updateMixMoles();
}
setMoles();
return icase;
return 0;
}
@ -307,30 +297,20 @@ namespace Cantera {
/// The constituent species are taken to be the first M species
/// in array 'species' that have linearly-independent compositions.
///
/// Arguments:
/// @param order On entry, vector \a order should contain species
/// index numbers in the order of decreasing desirability as a
/// constituent. For example, if it is desired to choose the
/// constituents from among the major species, this array might
/// list species index numbers in decreasing order of mole
/// fraction. If array 'species' does not have length =
/// nSpecies(), then the species will be considered as candidates
/// to be constituents in declaration order, beginning with the
/// first phase added.
///
/// On entry, vector species shold contain species index numbers
/// in the order of decreasing desirability as a constituent. For
/// example, if it is desired to choose the constituents from
/// among the major species, this array might list species index
/// numbers in decreasing order of mole fraction. If array
/// 'species' does not have length = nSpecies(), then the species
/// will be considered as candidates to be constituents in
/// declaration order, beginning with the first phase added.
///
/// On return, the first M entries of array 'species' contain the index
/// numbers of the constituent species.
///
/// Matrix nu is an output array that contains the stoichiometric
/// coefficents for a set of K - M formation reactions for the
/// non-constituent species, such that nu(k,i) is the net
/// stoichiometric coefficent of species k in reaction i. Matrix
/// nu will be resized to (K, K-M) and its initial values, if
/// any, will be erased.
void MultiPhaseEquil::getComponents(const vector_int& order) {
index_t m, k, j;
int n;
// if the input species array has the wrong size, ignore it
// and consider the species for constituents in declarationi order.
if (order.size() != m_nsp) {
@ -346,13 +326,6 @@ namespace Cantera {
index_t nColumns = m_nsp;
doublereal fctr;
#ifdef DEBUG_MULTIPHASE_EQUIL
cout << "most abundant:" << endl;
for (m = 0; m < nRows; m++) {
cout << m_mix->speciesName(m_species[m_order[m]]) << " " << m_moles[m_order[m]] << endl;
}
#endif
// set up the atomic composition matrix
for (m = 0; m < nRows; m++) {
for (k = 0; k < nColumns; k++) {
@ -380,17 +353,10 @@ namespace Cantera {
m_A(n, kmax) = tmp;
}
// exchange the species labels on the columns
#ifdef DEBUG_MULTIPHASE_EQUIL
cout << "in row " << m << ", pivot is zero" << endl;
cout << "exchanging " << m_mix->speciesName(m_species[m_order[m]]) << " for " << m_mix->speciesName(m_species[m_order[kmax]]) << endl;
#endif
itmp = m_order[m];
m_order[m] = m_order[kmax];
m_order[kmax] = itmp;
// throw an exception if the entire row is zero
// if (k >= m_nsp)
// throw CanteraError("getComponents","all zeros!");
}
// scale row m so that the diagonal element is unity
@ -423,23 +389,6 @@ namespace Cantera {
}
}
#ifdef DEBUG_MULTIPHASE_EQUIL
// check
bool ok = true;
for (m = 0; m < nRows; m++) {
cout << m_mix->speciesName(m_species[m_order[m]]) << " " << m_moles[m_order[m]] << endl;
if (m_A(m,m) != 1.0) ok = false;
for (n = 0; n < nRows; n++) {
if (n != m && fabs(m_A(m,n)) > TINY)
ok = false;
}
}
if (!ok) {
cout << m_A << endl;
throw CanteraError("getComponents","error in A matrix");
}
#endif
// create stoichometric coefficient matrix.
for (n = 0; n < int(m_nsp); n++) {
if (n < int(m_nel))
@ -462,17 +411,6 @@ namespace Cantera {
}
}
/// Re-arrange a vector of species properties in sequential form
/// into sorted (components first) form.
void MultiPhaseEquil::sort(vector_fp& x) {
copy(x.begin(), x.end(), m_work2.begin());
index_t k;
for (k = 0; k < m_nsp; k++) {
x[k] = m_work2[m_order[k]];
}
}
/// Re-arrange a vector of species properties in sorted form
/// (components first) into unsorted, sequential form.
void MultiPhaseEquil::unsort(vector_fp& x) {
@ -549,12 +487,14 @@ namespace Cantera {
m_moles[k] += omega * deltaN[k];
}
else {
m_moles[k] = fabs(m_moles[k])*fminn(10.0, exp(-m_deltaG_RT[ik - m_nel]));
m_moles[k] = fabs(m_moles[k])*fminn(10.0,
exp(-m_deltaG_RT[ik - m_nel]));
}
}
setMoles();
updateMixMoles();
}
/// Take one step in composition, given the gradient of G at the
/// starting point, and a vector of reaction steps dxi.
doublereal MultiPhaseEquil::
@ -567,14 +507,7 @@ namespace Cantera {
doublereal grad0 = computeReactionSteps(m_dxi);
// 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;
k = m_order[ik];
for (j = 0; j < m_nsp - m_nel; j++) {
m_work[ik] += m_N(ik, j) * m_dxi[j];
}
}
multiply(m_N, m_dxi.begin(), m_work.begin());
// change to sequential form
unsort(m_work);
@ -643,7 +576,6 @@ namespace Cantera {
for (k = 0; k < m_nsp; k++) {
grad1 += m_work[k] * m_mu[m_species[k]];
}
// doublereal grad1 = dot(m_work.begin(), m_work.end(), m_work2.begin());
omega = omegamax;
if (grad1 > 0.0) {
@ -739,10 +671,6 @@ 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;
//}
}
dxi[j] = -fctr*dg_rt;
index_t m;

10
Cantera/src/MultiPhaseEquil.h
View file

@ -47,13 +47,14 @@ namespace Cantera {
void getComponents(const vector_int& order);
int setInitialMoles();
int setInitialMoles2();
void computeN();
doublereal stepComposition(int loglevel);
void sort(vector_fp& x);
//void sort(vector_fp& x);
void unsort(vector_fp& x);
void step(doublereal omega, vector_fp& deltaN);
doublereal computeReactionSteps(vector_fp& dxi);
void setMoles();
void updateMixMoles();
void finish();
// moles of the species with sorted index ns
@ -61,7 +62,10 @@ namespace Cantera {
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); }
doublereal mu(int n) const { return m_mu[m_species[m_order[n]]]; }
string speciesName(int n) const { return
m_mix->speciesName(m_species[m_order[n]]); }
index_t m_nel_mix, m_nsp_mix, m_np;
index_t m_nel, m_nsp;
index_t m_eloc;