research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
cantera/Cantera/src/equil/MultiPhase.cpp
Ray Speth 4babd5898b Fixing compiler warnings, part 1
2012-01-17 04:10:08 +00:00

963 lines
26 KiB
C++
Raw Blame History

/**
* @file MultiPhase.cpp
* Definitions for the \link Cantera::MultiPhase MultiPhase\endlink
* object that is used to set up multiphase equilibrium problems (see \ref equilfunctions).
*/
#include "MultiPhase.h"
#include "MultiPhaseEquil.h"
#include "ThermoPhase.h"
#include "DenseMatrix.h"
#include "stringUtils.h"
#include "global.h"
using namespace std;
namespace Cantera {
/// Constructor.
MultiPhase::MultiPhase() :
m_np(0),
m_temp(0.0),
m_press(0.0),
m_nel(0),
m_nsp(0),
m_init(false),
m_eloc(-1),
m_Tmin(1.0),
m_Tmax(100000.0)
{
}
void MultiPhase::
addPhases(MultiPhase& mix) {
index_t n;
for (n = 0; n < mix.m_np; n++) {
addPhase(mix.m_phase[n], mix.m_moles[n]);
}
}
void MultiPhase::
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.");
}
// save the pointer to the phase object
m_phase.push_back(p);
// store its number of moles
m_moles.push_back(moles);
m_temp_OK.push_back(true);
// update the number of phases and the total number of
// species
m_np = m_phase.size();
m_nsp += p->nSpecies();
// determine if this phase has new elements
// for each new element, add an entry in the map
// from names to index number + 1:
string ename;
// iterate over the elements in this phase
index_t m, nel = p->nElements();
for (m = 0; m < nel; m++) {
ename = p->elementName(m);
// if no entry is found for this element name, then
// 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.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();
}
// 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;
}
}
// 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;
int mlocal;
string sym;
// allocate space for the atomic composition matrix
m_atoms.resize(m_nel, m_nsp, 0.0);
m_moleFractions.resize(m_nsp, 0.0);
m_elemAbundances.resize(m_nel, 0.0);
// iterate over the elements
// -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
// m_sptart(ip)
for (m = 0; m < m_nel; m++) {
sym = m_enames[m];
k = 0;
// iterate over the phases
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
nsp = p->nSpecies();
mlocal = p->elementIndex(sym);
for (kp = 0; kp < nsp; kp++) {
if (mlocal >= 0) {
m_atoms(m, k) = p->nAtoms(kp, mlocal);
}
if (m == 0) {
m_snames.push_back(p->speciesName(kp));
if (kp == 0) {
m_spstart.push_back(m_spphase.size());
}
m_spphase.push_back(ip);
}
k++;
}
}
}
if (m_eloc >= 0) {
doublereal esum;
for (k = 0; k < m_nsp; k++) {
esum = 0.0;
for (m = 0; m < m_nel; m++) {
if (int(m) != m_eloc)
esum += m_atoms(m,k) * m_atomicNumber[m];
}
//m_atoms(m_eloc, k) += esum;
}
}
/// set the initial composition within each phase to the
/// mole fractions stored in the phase objects
m_init = true;
updateMoleFractions();
}
// Return a reference to phase n. The state of phase n is
// also updated to match the state stored locally in the
// mixture object.
MultiPhase::phase_t& MultiPhase::phase(index_t n) {
if (!m_init) init();
m_phase[n]->setTemperature(m_temp);
m_phase[n]->setMoleFractions_NoNorm(DATA_PTR(m_moleFractions) + m_spstart[n]);
m_phase[n]->setPressure(m_press);
return *m_phase[n];
}
/// Moles of species \c k.
doublereal MultiPhase::speciesMoles(index_t k) const {
index_t ip = m_spphase[k];
return m_moles[ip]*m_moleFractions[k];
}
/// Total moles of element m, summed over all
/// phases
doublereal MultiPhase::elementMoles(index_t m) const {
doublereal sum = 0.0, phasesum;
index_t i, k = 0, ik, nsp;
for (i = 0; i < m_np; i++) {
phasesum = 0.0;
nsp = m_phase[i]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, i);
phasesum += m_atoms(m,k)*m_moleFractions[k];
}
sum += phasesum * m_moles[i];
}
return sum;
}
/// Total charge, summed over all phases
doublereal MultiPhase::charge() const {
doublereal sum = 0.0;
index_t i;
for (i = 0; i < m_np; i++) {
sum += phaseCharge(i);
}
return sum;
}
int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) {
int p = phaseIndex(phaseName);
if (p < 0) {
throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
}
int k = m_phase[p]->speciesIndex(speciesName);
if (k < 0) {
throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
}
return m_spstart[p] + k;
}
/// 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) const {
doublereal phasesum = 0.0;
size_t ik, k, nsp = m_phase[p]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, p);
phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
}
return Faraday*phasesum*m_moles[p];
}
/// Get the chemical potentials of all species in all phases.
void MultiPhase::getChemPotentials(doublereal* mu) const {
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
m_phase[i]->getChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
}
}
// 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) const {
index_t i, loc = 0;
updatePhases();
// iterate over the phases
for (i = 0; i < m_np; i++) {
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();
}
}
/// True if species \a k belongs to a solution phase.
bool MultiPhase::solutionSpecies(index_t k) const {
if (m_phase[m_spphase[k]]->nSpecies() > 1)
return true;
else
return false;
}
/// The Gibbs free energy of the mixture (J).
doublereal MultiPhase::gibbs() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->gibbs_mole() * m_moles[i];
return sum;
}
/// The enthalpy of the mixture (J).
doublereal MultiPhase::enthalpy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->enthalpy_mole() * m_moles[i];
return sum;
}
/// The internal energy of the mixture (J).
doublereal MultiPhase::IntEnergy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->intEnergy_mole() * m_moles[i];
return sum;
}
/// The entropy of the mixture (J/K).
doublereal MultiPhase::entropy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->entropy_mole() * m_moles[i];
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;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->cp_mole() * m_moles[i];
return sum;
}
/// Set the mole fractions of phase \a n to the values in
/// array \a x.
void MultiPhase::setPhaseMoleFractions(const index_t n, const doublereal* const x) {
phase_t* p = m_phase[n];
p->setState_TPX(m_temp, m_press, x);
int nsp = p->nSpecies();
int istart = m_spstart[n];
for (int k = 0; k < nsp; k++) {
m_moleFractions[istart+k] = x[k];
}
}
// 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) {
int kk = nSpecies();
doublereal x;
vector_fp moles(kk, 0.0);
for (int k = 0; k < kk; k++) {
x = xMap[speciesName(k)];
if (x > 0.0) moles[k] = x;
}
setMoles(DATA_PTR(moles));
}
// Set the species moles using a string. Unspecified species are
// set to zero.
void MultiPhase::setMolesByName(const std::string& x) {
compositionMap xx;
// 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);
}
// Get the mole numbers of all species in the multiphase
// object
void MultiPhase::getMoles(doublereal * molNum) const {
/*
* First copy in the mole fractions
*/
copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
index_t ik;
doublereal *dtmp = molNum;
for (index_t ip = 0; ip < m_np; ip++) {
doublereal phasemoles = m_moles[ip];
phase_t* p = m_phase[ip];
index_t nsp = p->nSpecies();
for (ik = 0; ik < nsp; ik++) {
*(dtmp++) *= phasemoles;
}
}
}
/// 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(const doublereal* n) {
if (!m_init) init();
index_t ip, loc = 0;
index_t ik, k = 0, nsp;
doublereal phasemoles;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
nsp = p->nSpecies();
phasemoles = 0.0;
for (ik = 0; ik < nsp; ik++) {
phasemoles += n[k];
k++;
}
m_moles[ip] = phasemoles;
if (nsp > 1) {
if (phasemoles > 0.0) {
p->setState_TPX(m_temp, m_press, n + loc);
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
} else {
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
}
}
else {
m_moleFractions[loc] = 1.0;
}
loc += nsp;
}
}
void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) {
vector_fp tmpMoles(m_nsp, 0.0);
getMoles(DATA_PTR(tmpMoles));
tmpMoles[indexS] += addedMoles;
if (tmpMoles[indexS] < 0.0) {
tmpMoles[indexS] = 0.0;
}
setMoles(DATA_PTR(tmpMoles));
}
void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) {
if (!m_init) init();
m_temp = T;
m_press = Pres;
updatePhases();
}
void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres,
const doublereal *n) {
m_temp = T;
m_press = Pres;
setMoles(n);
}
void MultiPhase::getElemAbundances(doublereal *elemAbundances) const {
index_t eGlobal;
calcElemAbundances();
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
}
}
// Internal routine to calculate the element abundance vector
void MultiPhase::calcElemAbundances() const {
index_t loc = 0;
index_t eGlobal;
index_t ik, kGlobal;
doublereal spMoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] = 0.0;
}
for (index_t ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
int nspPhase = p->nSpecies();
doublereal phasemoles = m_moles[ip];
for (ik = 0; ik < nspPhase; ik++) {
kGlobal = loc + ik;
spMoles = m_moleFractions[kGlobal] * phasemoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
}
}
loc += nspPhase;
}
}
/// The total mixture volume [m^3].
doublereal MultiPhase::volume() const {
int i;
doublereal sum = 0;
for (i = 0; i < int(m_np); i++) {
sum += m_moles[i]/m_phase[i]->molarDensity();
}
return sum;
}
doublereal MultiPhase::equilibrate(int XY, doublereal err,
int maxsteps, int maxiter, int loglevel) {
doublereal error;
bool strt = false;
doublereal dt;
doublereal h0;
int n;
bool start;
doublereal ferr, hnow, herr = 1.0;
doublereal snow, serr = 1.0, s0;
doublereal Tlow = -1.0, Thigh = -1.0;
doublereal Hlow = Undef, Hhigh = Undef, tnew;
doublereal dta=0.0, dtmax, cpb;
MultiPhaseEquil* e = 0;
if (!m_init) init();
if (loglevel > 0)
beginLogGroup("MultiPhase::equilibrate", loglevel);
if (XY == TP) {
if (loglevel > 0) {
addLogEntry("problem type","fixed T,P");
addLogEntry("Temperature",temperature());
addLogEntry("Pressure", pressure());
}
// create an equilibrium manager
e = new MultiPhaseEquil(this);
try {
error = e->equilibrate(XY, err, maxsteps, loglevel);
}
catch (CanteraError &err) {
if (loglevel > 0)
endLogGroup();
delete e;
e = 0;
throw err;
}
goto done;
}
else if (XY == HP) {
h0 = enthalpy();
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));
}
for (n = 0; n < maxiter; n++) {
// if 'strt' is false, the current composition will be used as
// the starting estimate; otherwise it will be estimated
// if (e) {
// cout << "e should be zero, but it is not!" << endl;
// delete e;
// }
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);
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;
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 {
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(); // iteration
}
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*(m_temp + Thigh);
if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0;
setTemperature(tnew);
if (loglevel > 0)
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
if (loglevel > 0)
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) {
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++) {
if (e) delete e;
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);
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("S rel error",fp2str(serr));
}
goto done;
}
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 err) {
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));
}
}
if (loglevel > 0)
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 == TV) {
addLogEntry("problem type","fixed T, V");
// doublereal dt = 1.0e3;
doublereal v0 = volume();
doublereal dVdP;
int n;
bool start = true;
doublereal error, vnow, pnow, verr;
for (n = 0; n < maxiter; n++) {
pnow = pressure();
MultiPhaseEquil e(this, start);
start = false;
beginLogGroup("iteration "+int2str(n));
error = e.equilibrate(TP, err, maxsteps, loglevel);
vnow = volume();
verr = fabs((v0 - vnow)/v0);
addLogEntry("P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
endLogGroup();
if (verr < err) {
addLogEntry("P iterations",int2str(n));
addLogEntry("Final P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
goto done;
}
// find dV/dP
setPressure(pnow*1.01);
dVdP = (volume() - vnow)/(0.01*pnow);
setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
}
}
else {
if (loglevel > 0)
endLogGroup();
throw CanteraError("MultiPhase::equilibrate","unknown option");
}
return -1.0;
done:
delete e;
e = 0;
if (loglevel > 0)
endLogGroup();
return err;
}
#ifdef MULTIPHASE_DEVEL
void importFromXML(string infile, string id) {
XML_Node* root = get_XML_File(infile);
if (id == "-") id = "";
XML_Node* x = get_XML_Node(string("#")+id, root);
if (x.name() != "multiphase")
throw CanteraError("MultiPhase::importFromXML",
"Current XML_Node is not a multiphase element.");
vector<XML_Node*> phases;
x.getChildren("phase",phases);
int np = phases.size();
int n;
ThermoPhase* p;
for (n = 0; n < np; n++) {
XML_Node& ph = *phases[n];
srcfile = infile;
if (ph.hasAttrib("src")) srcfile = ph["src"];
idstr = ph["id"];
p = newPhase(srcfile, idstr);
if (p) {
addPhase(p, ph.value());
}
}
}
#endif
void MultiPhase::setTemperature(const doublereal T) {
if (!m_init) init();
m_temp = T;
updatePhases();
}
// Name of element \a m.
std::string MultiPhase::elementName(size_t m) const {
return m_enames[m];
}
// Index of element with name \a name.
size_t MultiPhase::elementIndex(std::string name) const {
for (size_t e = 0; e < m_nel; e++) {
if (m_enames[e] == name) {
return e;
}
}
return -1;
}
// Name of species with global index \a k.
std::string MultiPhase::speciesName(const size_t k) const {
return m_snames[k];
}
doublereal MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const {
return m_atoms(mGlob, kGlob);
}
void MultiPhase::getMoleFractions(doublereal* const x) const {
std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
}
std::string MultiPhase::phaseName(const index_t iph) const {
const phase_t *tptr = m_phase[iph];
return tptr->id();
}
int MultiPhase::phaseIndex(const std::string &pName) const {
std::string tmp;
for (int iph = 0; iph < (int) m_np; iph++) {
const phase_t *tptr = m_phase[iph];
tmp = tptr->id();
if (tmp == pName) {
return iph;
}
}
return -1;
}
doublereal MultiPhase::phaseMoles(const index_t n) const {
return m_moles[n];
}
void MultiPhase::setPhaseMoles(const index_t n, const doublereal moles) {
m_moles[n] = moles;
}
int MultiPhase::speciesPhaseIndex(const index_t kGlob) const {
return m_spphase[kGlob];
}
doublereal MultiPhase::moleFraction(const index_t kGlob) const{
return m_moleFractions[kGlob];
}
bool MultiPhase::tempOK(const index_t p) const {
return m_temp_OK[p];
}
/// Update the locally-stored species mole fractions.
void MultiPhase::updateMoleFractions() {
uploadMoleFractionsFromPhases();
}
/// Update the locally-stored species mole fractions.
void MultiPhase::uploadMoleFractionsFromPhases() {
index_t ip, loc = 0;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
loc += p->nSpecies();
}
calcElemAbundances();
}
//-------------------------------------------------------------
//
// protected methods
//
//-------------------------------------------------------------
/// 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 = DATA_PTR(m_moleFractions) + 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;
}
}
}
}
Reference in a new issue View git blame Copy permalink