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

fixed convergence problem in MultiPhaseEquil

Browse source
This commit is contained in:
Dave Goodwin 2004-12-09 22:10:16 +00:00
parent 80dfdaa8fe
commit d5d13be7d4
3 changed files with 128 additions and 47 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

78
Cantera/src/MultiPhase.h
View file

@ -3,6 +3,7 @@
#include "ThermoPhase.h"
#include "DenseMatrix.h"
#include "stringUtils.h"
#include <iostream>
@ -37,15 +38,13 @@ namespace Cantera {
/// 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) {
if (m_init) {
throw CanteraError("addPhase","phases cannot be added after init() has been called.");
throw CanteraError("addPhase",
"phases cannot be added after init() has been called.");
}
// set this false so that init() will be called to
// recompute the atomic composition array
m_init = false;
// save the pointer to the phase object
m_phase.push_back(p);
@ -104,7 +103,9 @@ namespace Cantera {
copy(m_moleFractions.begin(), m_moleFractions.end(), x);
}
// process phases and build atomic composition array
/// Process phases and build atomic composition array. After
/// init() has been called, no more phases may be added.
void init() {
if (m_init) return;
index_t ip, kp, k = 0, nsp, m;
@ -130,8 +131,9 @@ namespace Cantera {
}
if (m == 0) {
m_snames.push_back(p->speciesName(kp));
if (kp == 0)
if (kp == 0) {
m_spstart.push_back(m_spphase.size());
}
m_spphase.push_back(ip);
}
k++;
@ -150,23 +152,32 @@ namespace Cantera {
return m_moles[n];
}
/// Set the number of moles of phase with index p.
/// Set the number of moles of phase with index n.
void setPhaseMoles(index_t n, doublereal moles) {
m_moles[n] = moles;
}
/// Return a reference to phase n.
/// Return a reference to phase n. The state of phase n is
/// also updated to match the state stored locally in the
/// mixture object.
phase_t& phase(index_t n) {
if (!m_init) init();
m_phase[n]->setState_TPX(m_temp, m_press,
m_moleFractions.begin() + m_spstart[n]);
return *m_phase[n];
}
/// Return a const reference to phase n.
const phase_t& phase(index_t n) const {
if (!m_init) init();
m_phase[n]->setState_TPX(m_temp,
m_press, m_moleFractions.begin() + m_spstart[n]);
return *m_phase[n];
}
/// Moles of species \c k.
doublereal speciesMoles(index_t k) {
if (!m_init) init();
index_t ip = m_spphase[k];
return m_moles[ip]*m_moleFractions[k];
}
@ -197,7 +208,8 @@ namespace Cantera {
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
void getChemPotentials(doublereal* mu) {
index_t i, k = 0, loc = 0;
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
m_phase[i]->getChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
@ -207,7 +219,8 @@ namespace Cantera {
/// Chemical potentials. Write into array \c mu the chemical
/// potentials of all species [J/kmol].
void getStandardChemPotentials(doublereal* mu) {
index_t i, k = 0, loc = 0;
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
m_phase[i]->getStandardChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
@ -237,6 +250,7 @@ namespace Cantera {
doublereal gibbs() {
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;
@ -272,6 +286,32 @@ namespace Cantera {
}
}
void setPhaseMoleFractions(index_t n, doublereal* x) {
phase_t* p = m_phase[n];
p->setState_TPX(m_temp, m_press, x);
}
void setMolesByName(compositionMap& xMap) {
int kk = nSpecies();
doublereal x;
vector_fp mf(kk, 0.0);
for (int k = 0; k < kk; k++) {
x = xMap[speciesName(k)];
if (x > 0.0) mf[k] = x;
}
setMoles(mf.begin());
}
void setMolesByName(const string& x) {
compositionMap xx;
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
parseCompString(x, xx);
setMolesByName(xx);
}
void setMoles(doublereal* n) {
if (!m_init) init();
index_t ip, loc = 0;
@ -330,16 +370,16 @@ namespace Cantera {
bool m_init;
};
inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) {
int ip;
for (ip = 0; ip < x.nPhases(); ip++) {
s << "*************** Phase " << ip << " *****************" << endl;
s << "Moles: " << x.phaseMoles(ip) << endl;
inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) {
size_t ip;
for (ip = 0; ip < x.nPhases(); ip++) {
s << "*************** Phase " << ip << " *****************" << endl;
s << "Moles: " << x.phaseMoles(ip) << endl;
s << report(x.phase(ip)) << endl;
}
return s;
s << report(x.phase(ip)) << endl;
}
return s;
}
}
#endif

91
Cantera/src/MultiPhaseEquil.cpp
View file

@ -121,6 +121,7 @@ namespace Cantera {
else
m_dsoln.push_back(0);
}
m_force = false;
setMoles();
}
@ -144,8 +145,8 @@ namespace Cantera {
* @param elementMoles vector of elemental moles
*/
int MultiPhaseEquil::setInitialMoles() {
int m, n;
double lp = log(m_press/OneAtm);
index_t m, n;
doublereal lp = log(m_press/OneAtm);
DenseMatrix aa(m_nel+2, m_nsp+1, 0.0);
@ -159,12 +160,13 @@ namespace Cantera {
m_mix->getStandardChemPotentials(m_mu.begin());
int kpp = 0;
index_t k, q;
doublereal rt = GasConstant * m_temp;
for (int k = 0; k < m_nsp; k++) {
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 (int q = 0; q < m_nel; q++)
for (q = 0; q < m_nel; q++)
aa(q+1, kpp) = -m_mix->nAtoms(m_species[k], m_element[q]);
}
@ -188,7 +190,7 @@ namespace Cantera {
for (n = 0; n < m_nel; n++) {
int ksp = 0;
int ip = iposv[n] - 1;
for (int k = 0; k < m_nsp; k++) {
for (int k = 0; k < int(m_nsp); k++) {
if (ip == ksp) {
m_moles[k] = aa(n+1, 0);
}
@ -233,8 +235,8 @@ namespace Cantera {
/// any, will be erased.
void MultiPhaseEquil::getComponents(const vector_int& order) {
int m, n, k, j;
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) {
@ -314,6 +316,7 @@ namespace Cantera {
// check
bool ok = true;
for (m = 0; m < nRows; m++) {
cout << m_mix->speciesName(m_species[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)
@ -369,6 +372,26 @@ namespace Cantera {
}
}
void MultiPhaseEquil::printInfo() {
index_t m, ik, k;
cout << "components: " << endl;
for (m = 0; m < m_nel; m++) {
ik = m_order[m];
k = m_species[ik];
cout << m_mix->speciesName(k) << " " << m_moles[ik] << endl;
}
cout << "non-components: " << endl;
for (m = m_nel; m < m_nsp; m++) {
ik = m_order[m];
k = m_species[ik];
cout << m_mix->speciesName(k) << " " << m_moles[ik] << endl;
}
cout << "Error = " << error() << endl;
for (k = 0; k < m_nsp - m_nel; k++) {
cout << reactionString(k) << " " << m_deltaG_RT[k] << endl;
}
}
/// Return a string specifying the jth reaction.
string MultiPhaseEquil::reactionString(index_t j) {
string sr = "", sp = "";
@ -391,8 +414,16 @@ namespace Cantera {
return sr + " <=> " + sp;
}
doublereal MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) {
void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) {
index_t k, ik;
//if (m_iter > 500) {
// for (ik = 0; ik < m_nsp; ik++) {
//k = m_order[ik];
//if (ik < m_nel) cout << "*";
//cout << m_mix->speciesName(m_species[k]) <<
// ": " << m_moles[k] << " += " << omega << " * " << deltaN[k] << endl;
// }
//}
if (omega < 0.0)
throw CanteraError("step","negative omega");
@ -421,11 +452,14 @@ namespace Cantera {
stepComposition() {
m_iter++;
index_t m, ip, ik, nsp, j, k = 0;
index_t ik, j, k = 0;
doublereal grad0 = computeReactionSteps(m_dxi);
if (grad0 > 0.0)
throw CanteraError("stepComposition", "positive gradient!");
//if (grad0 > 0.0) {
//cout << *m_mix << endl;
// cout << "gradient = " << grad0 << endl;
// throw CanteraError("stepComposition", "positive gradient!");
//}
// compute mole the fraction changes.
@ -442,13 +476,17 @@ namespace Cantera {
unsort(m_work);
// scale omega to keep the major species non-negative
const doublereal FCTR = 0.99;
doublereal FCTR = 0.99;
const doublereal MAJOR_THRESHOLD = 1.0e-12;
doublereal omega = 1.0, omax, omegamax = 1.0;
for (ik = 0; ik < m_nsp; ik++) {
k = m_order[ik];
if (ik < m_nel) {
FCTR = 0.99;
if (m_moles[k] < MAJOR_THRESHOLD) m_force = true;
}
else FCTR = 0.9;
// if species k is in a multi-species solution phase, then its
// mole number must remain positive, unless the entire phase
// goes away. First we'll determine an upper bound on omega,
@ -456,11 +494,13 @@ namespace Cantera {
if (m_dsoln[k] == 1) {
if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) {
if (m_moles[k] < MAJOR_THRESHOLD) m_force = true;
omax = m_moles[k]*FCTR/(fabs(m_work[k]) + TINY);
if (m_work[k] < 0.0 && omax < omegamax) {
omegamax = omax;
#ifdef DEBUG_MULTIPHASE_EQUIL
if (omegamax < 1.0e-5) {
m_force = true;
#ifdef DEBUG_MULTIPHASE_EQUIL
cout << m_mix->speciesName(m_species[k]) << " results in "
<< " omega = " << omegamax << endl;
//cout << m_moles[k] << " " << m_work[k] << endl;
@ -469,8 +509,8 @@ namespace Cantera {
for (nk = 0; nk < m_nel; nk++) {
cout << "component " << m_mix->speciesName(m_species[m_order[nk]]) << " " << m_moles[m_order[nk]] << endl;
}
}
#endif
}
}
m_majorsp[k] = true;
}
@ -483,14 +523,15 @@ namespace Cantera {
omax = -m_moles[k]/m_work[k];
if (omax < omegamax) {
omegamax = omax*1.000001;
#ifdef DEBUG_MULTIPHASE_EQUIL
if (omegamax < 1.0e-5) {
m_force = true;
#ifdef DEBUG_MULTIPHASE_EQUIL
cout << m_mix->speciesName(m_species[k]) << " results in "
<< " omega = " << omegamax << endl;
//cout << m_moles[k] << " " << m_work[k] << endl;
if (ik < m_nel) cout << "component" << endl;
}
#endif
}
}
}
m_majorsp[k] = true;
@ -512,7 +553,7 @@ namespace Cantera {
omega = omegamax;
if (grad1 > 0.0) {
omega *= -grad0 / (grad1 - grad0);
omega *= fabs(grad0) / (grad1 + fabs(grad0));
for (k = 0; k < m_nsp; k++) m_moles[k] = m_lastmoles[k];
step(omega, m_work);
}
@ -524,9 +565,8 @@ namespace Cantera {
doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) {
index_t i, j, k, ik, kc, ip;
int inu;
doublereal stoich, nmoles, csum, term1, fctr, dg_rt;
index_t j, k, ik, kc, ip;
doublereal stoich, nmoles, csum, term1, fctr;
vector_fp nu;
const doublereal TINY = 1.0e-20;
doublereal grad = 0.0;
@ -612,9 +652,7 @@ namespace Cantera {
}
void MultiPhaseEquil::computeN() {
index_t m, k, isp;
const doublereal THRESHOLD = 0.01;
index_t m, k;
// get the species moles
@ -643,10 +681,11 @@ namespace Cantera {
}
ok = false;
for (ij = 0; ij < m_nel; ij++) {
if (k == m_order[ij]) ok = true;
if (int(k) == m_order[ij]) ok = true;
}
if (!ok) {
if (!ok || m_force) {
getComponents(m_sortindex);
m_force = true;
break;
}
}

6
Cantera/src/MultiPhaseEquil.h
View file

@ -45,6 +45,7 @@ namespace Cantera {
if (error() < err) break;
}
if (i >= maxsteps) {
printInfo();
throw CanteraError("MultiPhaseEquil::equilibrate",
"no convergence in " + int2str(maxsteps) +
" iterations. Error = " + fp2str(error()));
@ -54,7 +55,7 @@ namespace Cantera {
string reactionString(index_t j);
doublereal error();
void printInfo();
protected:
void getComponents(const vector_int& order);
@ -63,7 +64,7 @@ namespace Cantera {
doublereal stepComposition();
void sort(vector_fp& x);
void unsort(vector_fp& x);
doublereal step(doublereal omega, vector_fp& deltaN);
void step(doublereal omega, vector_fp& deltaN);
doublereal computeReactionSteps(vector_fp& dxi);
void setMoles();
@ -84,6 +85,7 @@ namespace Cantera {
vector_int m_incl_element, m_incl_species;
vector_int m_species, m_element;
vector<bool> m_solnrxn;
bool m_force;
};
//-----------------------------------------------------------