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

Fixed a problem that caused the test suite to differ between optimized

Browse source 
and debug versions. Added cropping for test suite.
This commit is contained in:
Harry Moffat 2009-03-14 03:28:04 +00:00
parent 8055ac0e3a
commit 26918dc005
2 changed files with 88 additions and 74 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
Cantera/src/equil/ChemEquil.cpp
View file

@ -287,14 +287,14 @@ namespace Cantera {
#ifdef DEBUG_MODE
if (ChemEquil_print_lvl > 0) {
PrintCtrl pc(std::cout, -10, PrintCtrl::CT_OFF_GLOBALOBEY);
PrintCtrl pc(std::cout, -28, PrintCtrl::CT_OFF_GLOBALOBEY);
writelog("setInitialMoles: Estimated Mole Fractions\n");
writelogf(" Temperature = %g\n", s.temperature());
writelogf(" Pressure = %g\n", s.pressure());
for (int k = 0; k < m_kk; k++) {
string nnn = s.speciesName(k);
double mf = s.moleFraction(k);
mf = pc.cropAbs10(mf, -18);
mf = pc.cropAbs10(mf, -28);
writelogf(" %-12s % -10.5g\n", nnn.c_str(), mf);
}
writelog(" Element_Name ElementGoal ElementMF\n");
@ -386,7 +386,7 @@ namespace Cantera {
#ifdef DEBUG_MODE
if (ChemEquil_print_lvl > 0) {
PrintCtrl pc(std::cout, -18, PrintCtrl::CT_OFF_GLOBALOBEY);
PrintCtrl pc(std::cout, -28, PrintCtrl::CT_OFF_GLOBALOBEY);
for (m = 0; m < m_nComponents; m++) {
int isp = m_component[m];
string nnn = s.speciesName(isp);
@ -399,7 +399,7 @@ namespace Cantera {
writelog(" id Name MF mu/RT \n");
for (n = 0; n < s.nSpecies(); n++) {
string nnn = s.speciesName(n);
double mf = pc.cropAbs10(xMF_est[n], -18);
double mf = pc.cropAbs10(xMF_est[n], -28);
writelogf("%10d %15s %10.5g %10.5g\n",
n, nnn.c_str(), mf, mu_RT[n]);
}
@ -1161,10 +1161,12 @@ namespace Cantera {
#ifdef DEBUG_MODE
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
PrintCtrl pc(std::cout, -14, PrintCtrl::CT_OFF_GLOBALOBEY);
writelog("Residual: ElFracGoal ElFracCurrent Resid\n");
for (n = 0; n < m_mm; n++) {
double rrr = pc.cropAbs10(resid[n], -14);
writelogf(" % -14.7E % -14.7E % -10.5E\n",
elmFracGoal[n], elmFrac[n], resid[n]);
elmFracGoal[n], elmFrac[n], rrr);
}
}
#endif
@ -1183,9 +1185,11 @@ namespace Cantera {
#ifdef DEBUG_MODE
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
PrintCtrl pc(std::cout, -14, PrintCtrl::CT_OFF_GLOBALOBEY);
writelog(" Goal Xvalue Resid\n");
writelogf(" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]);
writelogf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
double rrr = pc.cropAbs10(resid[m_skip], -14);
writelogf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, rrr);
}
#endif
}

146
Cantera/src/equil/MultiPhaseEquil.cpp
View file

@ -169,6 +169,9 @@ namespace Cantera {
m_A.resize(m_nel, m_nsp, 0.0);
m_N.resize(m_nsp, m_nsp - m_nel);
m_order.resize(m_nsp, 0);
for (k = 0; k < m_nsp; k++) {
m_order[k] = k;
}
// if the 'start' flag is set, estimate the initial mole
// numbers by doing a linear Gibbs minimization. In this case,
@ -266,82 +269,89 @@ namespace Cantera {
}
/// 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(int loglevel) {
index_t ik, j;
/// 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(int loglevel) {
index_t ik, j;
double not_mu = 1.0e12;
if (loglevel > 0)
beginLogGroup("MultiPhaseEquil::setInitialMoles");
double not_mu = 1.0e12;
if (loglevel > 0)
beginLogGroup("MultiPhaseEquil::setInitialMoles");
m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu), true);
doublereal dg_rt;
m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu), true);
doublereal dg_rt;
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
bool redo = true;
int iter = 0;
while (redo) {
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
bool redo = true;
int iter = 0;
// choose a set of components based on the current
// composition
computeN();
if (loglevel > 0)
addLogEntry("iteration",iter);
redo = false;
iter++;
if (iter > 4) break;
while (redo) {
// 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);
// choose a set of components based on the current
// composition
computeN();
if (loglevel > 0)
addLogEntry("iteration",iter);
redo = false;
iter++;
if (iter > 4) break;
for (ik = 0; ik < m_nsp; ik++) {
nu = m_N(ik, j);
// 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 (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
if (loglevel > 0)
addLogEntry("component too small",speciesName(ik));
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;
}
}
// set the moles of the phase objects to match
updateMixMoles();
}
for (ik = 0; ik < m_nsp; ik++)
if (moles(ik) != 0.0) addLogEntry(speciesName(ik), moles(ik));
if (loglevel > 0)
endLogGroup("MultiPhaseEquil::setInitialMoles");
return 0;
// set max change in progress variable by
// non-negativity requirement
// -> Note, 0.99 factor is so that difference of 2 numbers
// isn't zero. This causes differences between
// optimized and debug versions of the code
if (nu*idir < 0) {
delta_xi = fabs(0.99*moles(ik)/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
if (loglevel > 0)
addLogEntry("component too small",speciesName(ik));
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;
}
}
// set the moles of the phase objects to match
updateMixMoles();
}
for (ik = 0; ik < m_nsp; ik++)
if (moles(ik) != 0.0) addLogEntry(speciesName(ik), moles(ik));
if (loglevel > 0)
endLogGroup("MultiPhaseEquil::setInitialMoles");
return 0;
}
/// This method finds a set of component species and a complete