From 26918dc005ddcece979c21617bab3d2bf34ab0cf Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 14 Mar 2009 03:28:04 +0000 Subject: [PATCH] Fixed a problem that caused the test suite to differ between optimized and debug versions. Added cropping for test suite. --- Cantera/src/equil/ChemEquil.cpp | 16 +-- Cantera/src/equil/MultiPhaseEquil.cpp | 146 ++++++++++++++------------ 2 files changed, 88 insertions(+), 74 deletions(-) diff --git a/Cantera/src/equil/ChemEquil.cpp b/Cantera/src/equil/ChemEquil.cpp index 66dc9b42e..d78714575 100755 --- a/Cantera/src/equil/ChemEquil.cpp +++ b/Cantera/src/equil/ChemEquil.cpp @@ -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 } diff --git a/Cantera/src/equil/MultiPhaseEquil.cpp b/Cantera/src/equil/MultiPhaseEquil.cpp index 4e50dd237..85dd13516 100644 --- a/Cantera/src/equil/MultiPhaseEquil.cpp +++ b/Cantera/src/equil/MultiPhaseEquil.cpp @@ -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