From e58bd0985910b6b4c6e6d72cd7366ae4c864a53f Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 18 Nov 2015 11:02:30 -0500 Subject: [PATCH] Eliminate debug_verbose option and DEBUG_MODE macro Vebose printing can always be enabled at runtime by setting flags on various classes. --- SConstruct | 7 - doc/sphinx/scons-options.txt | 5 - include/cantera/base/config.h.in | 13 - include/cantera/equil/vcs_solve.h | 8 +- include/cantera/thermo/HMWSoln.h | 8 +- include/cantera/thermo/PDSS_HKFT.h | 4 - include/cantera/thermo/RedlichKisterVPSSTP.h | 2 - src/base/ctml.cpp | 10 - src/equil/BasisOptimize.cpp | 18 +- src/equil/ChemEquil.cpp | 54 ++-- src/equil/vcs_MultiPhaseEquil.cpp | 4 - src/equil/vcs_VolPhase.cpp | 19 +- src/equil/vcs_elem.cpp | 8 +- src/equil/vcs_elem_rearrange.cpp | 4 +- src/equil/vcs_inest.cpp | 22 +- src/equil/vcs_nondim.cpp | 4 +- src/equil/vcs_phaseStability.cpp | 40 +-- src/equil/vcs_prob.cpp | 10 - src/equil/vcs_rxnadj.cpp | 225 +++++--------- src/equil/vcs_setMolesLinProg.cpp | 14 +- src/equil/vcs_solve.cpp | 4 - src/equil/vcs_solve_TP.cpp | 301 ++++++++----------- src/numerics/RootFind.cpp | 50 ++- src/thermo/HMWSoln.cpp | 249 +++++++-------- src/thermo/IonsFromNeutralVPSSTP.cpp | 52 ++-- src/thermo/LatticeSolidPhase.cpp | 12 +- src/thermo/PDSS_HKFT.cpp | 4 - src/thermo/RedlichKisterVPSSTP.cpp | 3 +- src/thermo/RedlichKwongMFTP.cpp | 3 +- src/transport/GasTransport.cpp | 30 +- src/transport/MMCollisionInt.cpp | 22 +- 31 files changed, 483 insertions(+), 726 deletions(-) diff --git a/SConstruct b/SConstruct index b8ec37ddc..40300ddba 100644 --- a/SConstruct +++ b/SConstruct @@ -390,11 +390,6 @@ config_options = [ ('FORTRANFLAGS', 'Compilation options for the Fortran (90) compiler.', '-O3'), - BoolVariable( - 'debug_verbose', - """Enable extra printing to aid in debugging. This code is marked - by the preprocessor macros DEBUG_MODE and DEBUG_MODE_ENABLED.""", - False), BoolVariable( 'coverage', """Enable collection of code coverage information with gcov. @@ -1201,8 +1196,6 @@ def cdefine(definevar, configvar, comp=True, value=1): else: configh[definevar] = None -cdefine('DEBUG_MODE', 'debug_verbose') - # Need to test all of these to see what platform.system() returns configh['SOLARIS'] = 1 if env['OS'] == 'Solaris' else None configh['DARWIN'] = 1 if env['OS'] == 'Darwin' else None diff --git a/doc/sphinx/scons-options.txt b/doc/sphinx/scons-options.txt index 06c24bce7..c20759cad 100644 --- a/doc/sphinx/scons-options.txt +++ b/doc/sphinx/scons-options.txt @@ -132,11 +132,6 @@ running 'scons build'. The format of this file is: Compilation options for the Fortran (90) compiler. - default: '-O3' -* debug_verbose: [ yes | no ] - Enable extra printing to aid in debugging. This code is marked by - the preprocessor macros DEBUG_MODE and DEBUG_MODE_ENABLED. - - default: 'no' - * coverage: [ yes | no ] Enable collection of code coverage information with gcov. Available only when compiling with gcc. diff --git a/include/cantera/base/config.h.in b/include/cantera/base/config.h.in index 2b8beb2f4..2549c17ef 100644 --- a/include/cantera/base/config.h.in +++ b/include/cantera/base/config.h.in @@ -8,19 +8,6 @@ // Just the major + minor version (i.e. 2.2 instead of 2.2.0) %(CANTERA_SHORT_VERSION)s -//------------------------ Development flags ------------------// -// -// Compile in additional debug printing where available. -// Note, the printing may need to be turned on via a switch. -// This just compiles in the code. -%(DEBUG_MODE)s - -#ifdef DEBUG_MODE -#define DEBUG_MODE_ENABLED 1 -#else -#define DEBUG_MODE_ENABLED 0 -#endif - //------------------------ Fortran settings -------------------// // define types doublereal, integer, and ftnlen to match the diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 76e1d250e..7ec9d8c24 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -509,9 +509,8 @@ public: * flags with total moles in each phase. */ double vcs_tmoles(); -#ifdef DEBUG_MODE + void check_tmoles() const; -#endif //! This subroutine calculates reaction free energy changes for //! all noncomponent formation reactions. @@ -1291,11 +1290,8 @@ private: */ double l2normdg(double dg[]) const; -#ifdef DEBUG_MODE - //! Print out and check the elemental abundance vector void prneav() const; -#endif void checkDelta1(double* const ds, double* const delTPhMoles, size_t kspec); @@ -1895,8 +1891,6 @@ public: * * 5 each iteration in solve_TP gets a report with one line per species * * 6 Each decision in solve_TP gets a line per species in addition to 4 * * 10 Additionally Hessian matrix is printed out - * - * Levels of printing above 4 are only accessible when DEBUG_MODE is turned on */ int m_debug_print_lvl; diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h index 419d98b10..5cfe42f1a 100644 --- a/include/cantera/thermo/HMWSoln.h +++ b/include/cantera/thermo/HMWSoln.h @@ -2789,16 +2789,10 @@ private: static int interp_est(const std::string& estString); public: - /*! - * Turn on copious debug printing when this - * is true and DEBUG_MODE is turned on. - */ + //! Turn on copious debug printing when this is true mutable int m_debugCalc; //! Return int specifying the amount of debug printing - /*! - * This will return 0 if DEBUG_MODE is not turned on - */ int debugPrinting(); }; diff --git a/include/cantera/thermo/PDSS_HKFT.h b/include/cantera/thermo/PDSS_HKFT.h index e80339379..95281c69a 100644 --- a/include/cantera/thermo/PDSS_HKFT.h +++ b/include/cantera/thermo/PDSS_HKFT.h @@ -86,7 +86,6 @@ public: virtual doublereal enthalpy_mole() const; -#ifdef DEBUG_MODE //! Return the molar enthalpy in units of J kmol-1 /*! * Returns the species standard state enthalpy in J kmol-1 at the @@ -97,7 +96,6 @@ public: * @return returns the species standard state enthalpy in J kmol-1 */ doublereal enthalpy_mole2() const; -#endif virtual doublereal intEnergy_mole() const; virtual doublereal entropy_mole() const; @@ -221,14 +219,12 @@ private: */ doublereal deltaS() const; -#ifdef DEBUG_MODE //! Routine that actually calculates the enthalpy difference //! between the reference state at Tr, Pr and T,P /*! * This is an extra routine that was added to check the arithmetic */ doublereal deltaH() const; -#endif //! Internal formula for the calculation of a_g() /*! diff --git a/include/cantera/thermo/RedlichKisterVPSSTP.h b/include/cantera/thermo/RedlichKisterVPSSTP.h index 6af271c07..1388c20f4 100644 --- a/include/cantera/thermo/RedlichKisterVPSSTP.h +++ b/include/cantera/thermo/RedlichKisterVPSSTP.h @@ -432,7 +432,6 @@ private: */ void s_update_dlnActCoeff_dX_() const; -#ifdef DEBUG_MODE public: //! Utility routine that calculates a literature expression /*! @@ -442,7 +441,6 @@ public: * nonideal term and mf term */ void Vint(double& VintOut, double& voltsOut); -#endif protected: //! number of binary interaction expressions diff --git a/src/base/ctml.cpp b/src/base/ctml.cpp index 1536b8ba1..37dfa7258 100644 --- a/src/base/ctml.cpp +++ b/src/base/ctml.cpp @@ -214,18 +214,8 @@ doublereal getFloatCurrent(const XML_Node& node, const std::string& type) fctr = toSI(units); } else if (type != "" && units != "") { fctr = toSI(units); -#ifdef DEBUG_MODE writelog("\nWarning: conversion toSI() was done on node value " + node.name() + "but wasn't explicitly requested. Type was \"" + type + "\"\n"); -#endif -#ifdef DEBUG_MODE_MORE - } else if (type == "" && units != "") { - // Note, below currently produces a lot of output due to transport blocks. - // This needs to be addressed. - writelog("\nWarning: XML node " + node.name() + - "has a units attribute, \"" + units + "\"," - "but no conversion was done because the getFloat() command didn't have a type\n"); -#endif } return fctr*x; } diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index dc9785255..a83da399a 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -56,7 +56,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, } } - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) { + if (BasisOptimize_print_lvl >= 1) { writelog(" "); for (i=0; i<77; i++) { writelog("-"); @@ -118,9 +118,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // For debugging purposes keep an unmodified copy of the array. vector_fp molNumBase; - if (DEBUG_MODE_ENABLED) { - molNumBase = molNum; - } + molNumBase = molNum; double molSave = 0.0; size_t jr = 0; @@ -157,9 +155,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Assign a small negative number to the component that we have // just found, in order to take it out of further consideration. -#ifdef DEBUG_MODE molSave = molNum[kk]; -#endif molNum[kk] = USEDBEFORE; // CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES @@ -208,7 +204,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // REARRANGE THE DATA if (jr != k) { - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) { + if (BasisOptimize_print_lvl >= 1) { kk = orderVectorSpecies[k]; writelogf(" --- %-12.12s", mphase->speciesName(kk)); jj = orderVectorSpecies[jr]; @@ -283,7 +279,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne, &ipiv[0], &formRxnMatrix[0], ne, info); - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) { + if (BasisOptimize_print_lvl >= 1) { writelog(" ---\n"); writelogf(" --- Number of Components = %d\n", nComponents); writelog(" --- Formula Matrix:\n"); @@ -379,7 +375,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, size_t nspecies = mphase->nSpecies(); double test = -1.0E10; - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { + if (BasisOptimize_print_lvl > 0) { writelog(" "); for (i=0; i<77; i++) { writelog("-"); @@ -456,7 +452,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, if (k == nelements) { // When we are here, there is an error usually. // We haven't found the number of elements necessary. - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { + if (BasisOptimize_print_lvl > 0) { writelogf("Error exit: returning with nComponents = %d\n", jr); } throw CanteraError("ElemRearrange", "Required number of elements not found."); @@ -517,7 +513,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, } // REARRANGE THE DATA if (jr != k) { - if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { + if (BasisOptimize_print_lvl > 0) { kk = orderVectorElements[k]; writelog(" --- "); writelogf("%-2.2s", mphase->elementName(kk)); diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index b5dacbdcc..bd5110b4b 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -204,7 +204,7 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, // and element abundance vectors kept within the ChemEquil object. update(s); - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog("setInitialMoles: Estimated Mole Fractions\n"); writelogf(" Temperature = %g\n", s.temperature()); writelogf(" Pressure = %g\n", s.pressure()); @@ -259,7 +259,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, doublereal rrt = 1.0/(GasConstant* s.temperature()); scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt); - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { for (size_t m = 0; m < m_nComponents; m++) { size_t isp = m_component[m]; writelogf("isp = %d, %s\n", isp, s.speciesName(isp)); @@ -293,7 +293,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, lambda_RT[m_orderVectorElements[m]] = 0.0; } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog(" id CompSpecies ChemPot EstChemPot Diff\n"); for (size_t m = 0; m < m_nComponents; m++) { size_t isp = m_component[m]; @@ -639,7 +639,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // Compute the Jacobian matrix equilJacobian(s, x, elMolesGoal, jac, xval, yval); - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf("Jacobian matrix %d:\n", iter); for (m = 0; m <= m_mm; m++) { writelog(" [ "); @@ -698,7 +698,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, fctr = std::min(fctr, 0.2/fabs(res_trial[mm])); } } - if (fctr != 1.0 && DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (fctr != 1.0 && ChemEquil_print_lvl > 0) { writelogf("WARNING Soln Damping because of bounds: %g\n", fctr); } @@ -756,7 +756,7 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, for (size_t m = 0; m < x.size(); m++) { x[m] = oldx[m] + damp * step[m]; } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf("Solution Unknowns: damp = %g\n", damp); writelog(" X_new X_old Step\n"); for (size_t m = 0; m < m_mm; m++) { @@ -794,7 +794,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !m_doResPerturb) { + if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { writelog("Residual: ElFracGoal ElFracCurrent Resid\n"); for (size_t n = 0; n < m_mm; n++) { writelogf(" % -14.7E % -14.7E % -10.5E\n", @@ -807,7 +807,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, resid[m_mm] = xx/xval - 1.0; resid[m_skip] = yy/yval - 1.0; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !m_doResPerturb) { + if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { 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]); @@ -957,7 +957,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog("estimateEP_Brinkley::\n\n"); double temp = s.temperature(); double pres = s.pressure(); @@ -991,7 +991,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } x_old[m_mm] = n_t; // Calculate the mole numbers of species - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf("START ITERATION %d:\n", iter); } // Calculate the mole numbers of species and elements. @@ -1002,7 +1002,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); for (k = 0; k < m_kk; k++) { writelogf("%15s: %10.5g %10.5g\n", @@ -1032,7 +1032,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !normalStep) { + if (ChemEquil_print_lvl > 0 && !normalStep) { writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); } if (!normalStep) { @@ -1087,7 +1087,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } nCutoff = 1.0E-9 * n_t_calc; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog(" Lump Sum Elements Calculation: \n"); } for (m = 0; m < m_mm; m++) { @@ -1111,7 +1111,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf(" %5s %3d : %5d %5d\n", s.elementName(m), lumpSum[m], kMSp, kMSp2); } @@ -1164,7 +1164,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, for (m = 0; m < m_mm; m++) { if (a1(m,m) < 1.0E-50) { - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m); } for (n = 0; n < m_mm; n++) { @@ -1183,7 +1183,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, resid[m_mm] = n_t - n_t_calc; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog("Matrix:\n"); for (m = 0; m <= m_mm; m++) { writelog(" ["); @@ -1196,7 +1196,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, tmp = resid[m_mm] /(n_t + 1.0E-15); sum += tmp * tmp; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf("(it %d) Convergence = %g\n", iter, sum); } @@ -1216,7 +1216,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, tmp += fabs(a1(m,n)); } if (m < m_mm && tmp < 1.0E-30) { - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing row %d\n", m); } for (n = 0; n <= m_mm; n++) { @@ -1233,7 +1233,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, resid[m] *= tmp; } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelog("Row Summed Matrix:\n"); for (m = 0; m <= m_mm; m++) { writelog(" ["); @@ -1279,7 +1279,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } if (sameAsRow != npos || lumpSum[m]) { - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { if (lumpSum[m]) { writelogf("Lump summing row %d, due to rank deficiency analysis\n", m); } else if (sameAsRow != npos) { @@ -1296,7 +1296,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && modifiedMatrix) { + if (ChemEquil_print_lvl > 0 && modifiedMatrix) { writelog("Row Summed, MODIFIED Matrix:\n"); for (m = 0; m <= m_mm; m++) { writelog(" ["); @@ -1328,7 +1328,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, beta = std::min(beta, -1.0 / resid[m]); } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && beta != 1.0) { + if (ChemEquil_print_lvl > 0 && beta != 1.0) { writelogf("(it %d) Beta = %g\n", iter, beta); } } @@ -1338,7 +1338,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } n_t *= exp(beta * resid[m_mm]); - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); for (m = 0; m < m_mm; m++) { writelogf(" %5s %10.5g %10.5g %10.5g\n", @@ -1347,7 +1347,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, writelogf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { double temp = s.temperature(); double pres = s.pressure(); @@ -1377,7 +1377,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) size_t maxNegEloc = npos; double maxPosVal = -1.0; double maxNegVal = -1.0; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + if (ChemEquil_print_lvl > 0) { for (k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0 && m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { maxPosVal = m_molefractions[k]; @@ -1407,7 +1407,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) return; } double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) { + if (ChemEquil_print_lvl > 0 && factor < 0.9999999999) { writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", s.speciesName(maxPosEloc), m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); @@ -1419,7 +1419,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) } } else { double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg; - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) { + if (ChemEquil_print_lvl > 0 && factor < 0.9999999999) { writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", s.speciesName(maxNegEloc), m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index fc7679e63..b6b3eddaa 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -579,9 +579,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) double Temp = m_mix->temperature(); double pres = m_mix->pressure(); vector_fp& mf = m_vprob.mf; -#ifdef DEBUG_MODE double* fe = &m_vprob.m_gibbsSpecies[0]; -#endif vector_fp VolPM; vector_fp activity; vector_fp ac; @@ -692,7 +690,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) } } -#ifdef DEBUG_MODE // Check consistency: These should be equal tref.getChemPotentials(fe+istart); for (size_t k = 0; k < nSpecies; k++) { @@ -702,7 +699,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "incompatibility!"); } } -#endif } fclose(FP); diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 684fd5c30..0006d9f54 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -430,9 +430,10 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc, molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_old[0]; } else if (stateCalc == VCS_STATECALC_NEW) { molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_new[0]; - } else if (DEBUG_MODE_ENABLED) { - throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here"); } - } else if (DEBUG_MODE_ENABLED && m_owningSolverObject) { + } else { + throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here"); + } + } else if (m_owningSolverObject) { if (stateCalc == VCS_STATECALC_OLD) { if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_old[0]) { throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here"); @@ -890,20 +891,16 @@ void vcs_VolPhase::setExistence(const int existence) { if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) { if (v_totalMoles != 0.0) { - if (DEBUG_MODE_ENABLED) { - throw CanteraError("vcs_VolPhase::setExistence", - "setting false existence for phase with moles"); - } else { - v_totalMoles = 0.0; - } + throw CanteraError("vcs_VolPhase::setExistence", + "setting false existence for phase with moles"); } - } else if (DEBUG_MODE_ENABLED && m_totalMolesInert == 0.0) { + } else if (m_totalMolesInert == 0.0) { if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) { throw CanteraError("vcs_VolPhase::setExistence", "setting true existence for phase with no moles"); } } - if (DEBUG_MODE_ENABLED && m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) { + if (m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) { throw CanteraError("vcs_VolPhase::setExistence", "Trying to set existence of an electron phase to false"); } diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 518b0bfab..9a3bfaa63 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -94,7 +94,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) { int retn = 0; -#ifdef DEBUG_MODE vector_fp ga_save(m_elemAbundances); if (m_debug_print_lvl >= 2) { plogf(" --- vcsc_elcorr: Element abundances correction routine"); @@ -112,7 +111,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) l2before += x[i] * x[i]; } l2before = sqrt(l2before/m_numElemConstraints); -#endif // Special section to take out single species, single component, // moles. These are species which have non-zero entries in the @@ -190,7 +188,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[i] / atomComp; if (m_molNumSpecies_old[kspec] > maxPermissible) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) { + if (m_debug_print_lvl >= 3) { plogf(" --- vcs_elcorr: Reduced species %s from %g to %g " "due to %s max bounds constraint\n", m_speciesName[kspec], m_molNumSpecies_old[kspec], @@ -205,7 +203,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } else { m_speciesStatus[kspec] = VCS_SPECIES_ACTIVEBUTZERO; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_elcorr: Zeroed species %s and changed " "status to %d due to max bounds constraint\n", m_speciesName[kspec], m_speciesStatus[kspec]); @@ -422,7 +420,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) L_CLEANUP: ; vcs_tmoles(); -#ifdef DEBUG_MODE double l2after = 0.0; for (size_t i = 0; i < m_numElemConstraints; ++i) { l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2); @@ -439,7 +436,6 @@ L_CLEANUP: plogf(" --- Diff_Norm: %20.12E %20.12E\n", l2before, l2after); } -#endif return retn; } diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index ed3226aca..9edcab763 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -21,7 +21,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, double* const sm, double* const ss) { size_t ncomponents = m_numComponents; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" "); for (size_t i=0; i<77; i++) { plogf("-"); @@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, } // REARRANGE THE DATA if (jr != k) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-2.2s", m_elementName[k]); plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]); diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 191983194..819b8078e 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -28,7 +28,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, // molNum(I) .GE. 0.0. Note, both of these programs do this. vcs_setMolesLinProg(); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n", pprefix); plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix); @@ -116,7 +116,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } vcs_deltag(0, true, VCS_STATECALC_NEW); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { for (size_t kspec = 0; kspec < nspecies; ++kspec) { plogf("%s", pprefix); plogf("%-12.12s", m_speciesName[kspec]); @@ -171,7 +171,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { for (size_t kspec = 0; kspec < nspecies; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("%sdirection (", pprefix); @@ -269,7 +269,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, lt = 1; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%s Final Mole Numbers produced by inest:\n", pprefix); plogf("%s SPECIES MOLE_NUMBER\n", pprefix); @@ -290,14 +290,14 @@ int VCS_SOLVE::vcs_inest_TP() // Calculate the elemental abundances vcs_elab(); if (vcs_elabcheck(0)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%s Initial guess passed element abundances on input\n", pprefix); plogf("%s m_doEstimateEquil = 1 so will use the input mole " "numbers as estimates", pprefix); plogendl(); } return retn; - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf("%s Initial guess failed element abundances on input\n", pprefix); plogf("%s m_doEstimateEquil = 1 so will discard input " "mole numbers and find our own estimate", pprefix); @@ -312,7 +312,7 @@ int VCS_SOLVE::vcs_inest_TP() vector_fp aw(m_numSpeciesTot+ m_numElemConstraints, 0.0); // Go get the estimate of the solution - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%sGo find an initial estimate for the equilibrium problem", pprefix); plogendl(); @@ -332,7 +332,7 @@ int VCS_SOLVE::vcs_inest_TP() // matrix. bool rangeCheck = vcs_elabcheck(1); if (!vcs_elabcheck(0)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%sInitial guess failed element abundances\n", pprefix); plogf("%sCall vcs_elcorr to attempt fix", pprefix); plogendl(); @@ -347,7 +347,7 @@ int VCS_SOLVE::vcs_inest_TP() plogendl(); retn = -1; } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (rangeCheck) { plogf("%sInitial guess now satisfies element abundances", pprefix); plogendl(); @@ -361,7 +361,7 @@ int VCS_SOLVE::vcs_inest_TP() } } } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (rangeCheck) { plogf("%sInitial guess satisfies element abundances", pprefix); plogendl(); @@ -375,7 +375,7 @@ int VCS_SOLVE::vcs_inest_TP() } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix, vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_new[0], &m_tPhaseMoles_old[0])); diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index d98e42f57..b1d386ab1 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -105,7 +105,7 @@ void VCS_SOLVE::vcs_nondim_TP() if (m_totalMoleScale != 1.0) { if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale); plogendl(); } @@ -149,7 +149,7 @@ void VCS_SOLVE::vcs_redim_TP() } if (m_totalMoleScale != 1.0) { if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale); plogendl(); } diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index ff6d7a02c..aa5c33b5f 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -205,24 +205,18 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) doublereal FephaseMax = -1.0E30; doublereal Fephase = -1.0E30; -#ifdef DEBUG_MODE char anote[128]; if (m_debug_print_lvl >= 2) { plogf(" --- vcs_popPhaseID() called\n"); plogf(" --- Phase Status F_e MoleNum\n"); plogf(" --------------------------------------------------------------------------\n"); } -#else - char* anote; -#endif for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* Vphase = m_VolPhaseList[iph]; int existence = Vphase->exists(); - if (DEBUG_MODE_ENABLED) { - strcpy(anote, ""); - } + strcpy(anote, ""); if (existence > 0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %18s %5d NA %11.3e\n", Vphase->PhaseName, existence, m_tPhaseMoles_old[iph]); } @@ -234,24 +228,20 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) doublereal deltaGRxn = m_deltaGRxn_old[irxn]; Fephase = exp(-deltaGRxn) - 1.0; if (Fephase > 0.0) { - if (DEBUG_MODE_ENABLED) { - strcpy(anote," (ready to be birthed)"); - } + strcpy(anote," (ready to be birthed)"); if (Fephase > FephaseMax) { iphasePop = iph; FephaseMax = Fephase; - if (DEBUG_MODE_ENABLED) { - strcpy(anote," (chosen to be birthed)"); - } + strcpy(anote," (chosen to be birthed)"); } } - if (DEBUG_MODE_ENABLED && Fephase < 0.0) { + if (Fephase < 0.0) { strcpy(anote," (not stable)"); AssertThrowMsg(m_tPhaseMoles_old[iph] <= 0.0, "VCS_SOLVE::vcs_popPhaseID", "shouldn't be here"); } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %18s %5d %10.3g %10.3g %s\n", Vphase->PhaseName, existence, Fephase, m_tPhaseMoles_old[iph], anote); @@ -268,13 +258,13 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) } else { FephaseMax = std::max(FephaseMax, Fephase); } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %18s %5d %11.3g %11.3g\n", Vphase->PhaseName, existence, Fephase, m_tPhaseMoles_old[iph]); } } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %18s %5d blocked %11.3g\n", Vphase->PhaseName, existence, m_tPhaseMoles_old[iph]); @@ -290,7 +280,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) // Insert logic here to figure out if phase pops are linked together. Only // do one linked pop at a time. - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" ---------------------------------------------------------------------\n"); } return iphasePop; @@ -312,7 +302,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes", "called for a phase that exists!"); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n", Vphase->PhaseName, iphasePop); } @@ -494,13 +484,13 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) if (doSuccessiveSubstitution) { int KP = 0; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_phaseStabilityTest() called\n"); plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]" " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP); plogf(" --------------------------------------------------------------" "--------------------------------------------------------\n"); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { + } else if (m_debug_print_lvl == 1) { plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph); } @@ -679,7 +669,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) fracDelta_new[k] = fracDelta_old[k] + damp * delFrac[k]; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP], delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability); } @@ -703,10 +693,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } else { throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet"); } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" ------------------------------------------------------------" "-------------------------------------------------------------\n"); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { + } else if (m_debug_print_lvl == 1) { if (funcPhaseStability > 0.0) { plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability); } else { diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 06c07e9eb..20ca98750 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -460,16 +460,6 @@ void VCS_PROB::reportCSV(const std::string& reportFile) } } - if (DEBUG_MODE_ENABLED) { - // Check consistency: These should be equal - tp->getChemPotentials(&m_gibbsSpecies[0]+istart); - for (size_t k = 0; k < nSpeciesPhase; k++) { - if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) { - fclose(FP); - throw CanteraError("VCS_PROB::reportCSV", "incompatibility"); - } - } - } iK += nSpeciesPhase; } fclose(FP); diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index f5f1039cd..5a987aec6 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -21,7 +21,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) { size_t iphDel = npos; size_t k = 0; -#ifdef DEBUG_MODE char ANOTE[128]; if (m_debug_print_lvl >= 2) { plogf(" "); @@ -38,9 +37,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) plogf(" --- Species KMoles Rxn_Adjustment DeltaG" " | Comment\n"); } -#else - char* ANOTE = 0; -#endif // We update the matrix dlnActCoeffdmolNumber[][] at the top of the loop, // when necessary @@ -50,16 +46,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // LOOP OVER THE FORMATION REACTIONS for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Normal Calc"); - } + sprintf(ANOTE, "Normal Calc"); size_t kspec = m_indexRxnToSpecies[irxn]; if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) { m_deltaMolNumSpecies[kspec] = 0.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); - } + sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); } else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { // MULTISPECIES PHASE WITH total moles equal to zero @@ -78,37 +70,27 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES; if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { m_deltaMolNumSpecies[kspec] = 0.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); } else { m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); - } - } - } else { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E", + sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); } + } else { + sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); m_deltaMolNumSpecies[kspec] = 0.0; if (Vphase->exists() > 0 && trphmoles > 0.0) { m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, - "MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, + "MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E", + vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); } } } else { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15), - m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15), + m_deltaGRxn_new[irxn]); m_deltaMolNumSpecies[kspec] = 0.0; } } else { @@ -118,14 +100,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // bother if superconvergence has already been achieved in this // mode. if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]); - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); - } + sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]); + if (m_debug_print_lvl >= 2) { + plogf(" --- %-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E %12.4E | %s\n", + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } continue; } @@ -133,14 +113,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // Don't calculate for minor or nonexistent species if their // values are to be decreasing anyway. if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); - } + sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); + if (m_debug_print_lvl >= 2) { + plogf(" --- %-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E %12.4E | %s\n", + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } continue; } @@ -170,10 +148,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) if (m_useActCoeffJac) { double s_old = s; s = vcs_Hessian_diag_adj(irxn, s_old); - if (DEBUG_MODE_ENABLED && s_old != s) { - sprintf(ANOTE, "Normal calc: diag adjusted from %g " - "to %g due to act coeff", s_old, s); - } + sprintf(ANOTE, "Normal calc: diag adjusted from %g " + "to %g due to act coeff", s_old, s); } m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s; @@ -184,18 +160,14 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) double negChangeComp = -stoicC * m_deltaMolNumSpecies[kspec]; if (negChangeComp > m_molNumSpecies_old[j]) { if (m_molNumSpecies_old[j] > 0.0) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Delta damped from %g " - "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); - } + sprintf(ANOTE, "Delta damped from %g " + "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[j] / stoicC; } else { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Delta damped from %g " - "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); - } + sprintf(ANOTE, "Delta damped from %g " + "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str()); m_deltaMolNumSpecies[kspec] = 0.0; } } @@ -204,11 +176,9 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // Implement a damping term that limits m_deltaMolNumSpecies // to the size of the mole number if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Delta damped from %g " - "to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec], - m_speciesName[kspec].c_str()); - } + sprintf(ANOTE, "Delta damped from %g " + "to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec], + m_speciesName[kspec].c_str()); m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; } } else { @@ -261,12 +231,10 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // having all of the mole numbers of that phases. it // seems that we can suggest a zero of the species // and the code will recover. - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); - } + sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec]); plogf(" %12.4E %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], @@ -287,31 +255,29 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) iphDel = m_phaseID[k]; kSpecial = k; - if (DEBUG_MODE_ENABLED) { - if (k != kspec) { - sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel, - m_speciesName[k].c_str()); - } else { - sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel); - } - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); - plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s", - m_speciesName[k]); - plogendl(); - } + if (k != kspec) { + sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel, + m_speciesName[k].c_str()); + } else { + sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel); + } + if (m_debug_print_lvl >= 2) { + plogf(" --- %-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E %12.4E | %s\n", + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); + plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s", + m_speciesName[k]); + plogendl(); } if (k != kspec) { forceComponentCalc = 1; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Force a component recalculation \n"); plogendl(); } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" "); writeline('-', 82); } @@ -319,7 +285,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } } } // End of regular processing - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec]); plogf(" %12.4E %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], @@ -327,7 +293,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) } } // End of loop over m_speciesUnknownType } // End of loop over non-component stoichiometric formation reactions - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" "); writeline('-', 82); } @@ -337,7 +303,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) int VCS_SOLVE::vcs_rxn_adj_cg() { int soldel = 0; -#ifdef DEBUG_MODE char ANOTE[128]; plogf(" "); for (size_t j = 0; j < 77; j++) { @@ -345,17 +310,12 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } plogf("\n --- Subroutine rxn_adj_cg() called\n"); plogf(" --- Species Moles Rxn_Adjustment | Comment\n"); -#else - char* ANOTE = 0; -#endif // Precalculation loop -> we calculate quantities based on loops over the // number of species. We also evaluate whether the matrix is appropriate for // this algorithm. If not, we bail out. for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Normal Calc"); - } + sprintf(ANOTE, "Normal Calc"); size_t kspec = m_indexRxnToSpecies[irxn]; if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) { @@ -364,16 +324,12 @@ int VCS_SOLVE::vcs_rxn_adj_cg() // HKM -> the statment below presupposes units in m_deltaGRxn_new[]. // It probably should be replaced with something more relativistic if (m_deltaGRxn_new[irxn] < -1.0e-4) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]); m_deltaMolNumSpecies[kspec] = 1.0e-10; m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; --m_numRxnMinorZeroed; } else { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]); m_deltaMolNumSpecies[kspec] = 0.0; } } else { @@ -382,26 +338,22 @@ int VCS_SOLVE::vcs_rxn_adj_cg() // First take care of cases where we want to bail out. Don't bother // if superconvergence has already been achieved in this mode. if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); - plogf(" --- "); - plogf("%-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - } + sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); + plogf(" --- "); + plogf("%-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); continue; } // Don't calculate for minor or nonexistent species if their values // are to be decreasing anyway. if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); - plogf(" --- "); - plogf("%-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - } + sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); + plogf(" --- "); + plogf("%-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); continue; } @@ -474,11 +426,9 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } m_molNumSpecies_old[k] = 0.0; m_tPhaseMoles_old[m_phaseID[k]] = 0.0; - if (DEBUG_MODE_ENABLED) { - plogf(" --- vcs_st2 Special section to delete "); - plogf("%-12.12s", m_speciesName[k]); - plogf("\n --- Immediate return - Restart iteration\n"); - } + plogf(" --- vcs_st2 Special section to delete "); + plogf("%-12.12s", m_speciesName[k]); + plogf("\n --- Immediate return - Restart iteration\n"); // We need to immediately recompute the component basis, // because we just zeroed it out. @@ -491,25 +441,21 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } } } // End of regular processing - if (DEBUG_MODE_ENABLED) { - plogf(" --- "); - plogf("%-12.12s", m_speciesName[kspec]); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - } + plogf(" --- "); + plogf("%-12.12s", m_speciesName[kspec]); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); } // End of loop over non-component stoichiometric formation reactions // When we form the Hessian we must be careful to ensure that it is a // symmetric positive definite matrix, still. This means zeroing out columns // when we zero out rows as well. I suggest writing a small program to make // sure of this property. - if (DEBUG_MODE_ENABLED) { - plogf(" "); - for (size_t j = 0; j < 77; j++) { - plogf("-"); - } - plogf("\n"); + plogf(" "); + for (size_t j = 0; j < 77; j++) { + plogf("-"); } + plogf("\n"); return soldel; } @@ -613,7 +559,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* if (deltaGOrig > 0.0) { if (dx_orig > 0.0) { dx = 0.0; - if (DEBUG_MODE_ENABLED && ANOTE) { + if (ANOTE) { sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx>0 dg > 0"); } return dx; @@ -621,7 +567,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* } else if (deltaGOrig < 0.0) { if (dx_orig < 0.0) { dx = 0.0; - if (DEBUG_MODE_ENABLED && ANOTE) { + if (ANOTE) { sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx<0 dg < 0"); } return dx; @@ -699,17 +645,12 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* finalize: vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); if (its >= MAXITS) { -#ifdef DEBUG_MODE sprintf(ANOTE, "Rxn reduced to zero step size from %g to %g (MAXITS)", dx_orig, dx); return dx; -#endif } -#ifdef DEBUG_MODE if (dx != dx_orig) { sprintf(ANOTE, "Line Search reduced step size from %g to %g", dx_orig, dx); } -#endif - return dx; } diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 76aec21e9..721faefc0 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -34,7 +34,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() size_t ik, irxn; double test = -1.0E-10; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- call setInitialMoles\n"); } @@ -59,13 +59,13 @@ int VCS_SOLVE::vcs_setMolesLinProg() } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); } while (redo) { if (!vcs_elabcheck(0)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- seMolesLinProg Mole numbers failing element abundances\n"); plogf(" --- seMolesLinProg Call vcs_elcorr to attempt fix\n"); } @@ -88,7 +88,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() return retn; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("iteration %d\n", iter); } redo = false; @@ -124,7 +124,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() // if a component has nearly zero moles, redo // with a new set of components if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Component too small: %s\n", m_speciesName[jcomp]); } redo = true; @@ -152,12 +152,12 @@ int VCS_SOLVE::vcs_setMolesLinProg() } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { + if (m_debug_print_lvl == 1) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); plogf(" --- setInitialMoles end\n"); } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index cf9ae3d10..d2e0d4027 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -362,11 +362,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) // NPhase = number of phases m_numPhases = nph; -#ifdef DEBUG_MODE m_debug_print_lvl = pub->vcs_debug_print_lvl; -#else - m_debug_print_lvl = std::min(2, pub->vcs_debug_print_lvl); -#endif // FormulaMatrix[] -> Copy the formula matrix over for (size_t i = 0; i < nspecies; i++) { diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 2b9f87c87..6b7e3ca8f 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -58,13 +58,9 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) bool uptodate_minors = true; int forceComponentCalc = 1; -#ifdef DEBUG_MODE char ANOTE[128]; // Set the debug print lvl to the same as the print lvl. m_debug_print_lvl = printDetails; -#else - char* ANOTE = 0; -#endif if (printDetails > 0 && print_lvl == 0) { print_lvl = 1; } @@ -256,7 +252,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) npb = vcs_add_all_deleted(); if (npb > 0) { iti = 0; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) { + if (m_debug_print_lvl >= 1) { plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!"); plogendl(); } @@ -328,7 +324,7 @@ int VCS_SOLVE::solve_tp_component_calc(bool& allMinorZeroedSpecies) // EVALUATE THE ELELEMT ABUNDANCE CHECK if (! vcs_elabcheck(0)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Element Abundance check failed"); plogendl(); } @@ -338,7 +334,7 @@ int VCS_SOLVE::solve_tp_component_calc(bool& allMinorZeroedSpecies) // Update the phase objects with the contents of the soln vector vcs_updateVP(VCS_STATECALC_OLD); vcs_deltag(0, false, VCS_STATECALC_OLD); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf(" --- Element Abundance check passed"); plogendl(); } @@ -384,9 +380,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // Calculate the total moles in each phase -> old solution // -> Needed for numerical stability when phases disappear. // -> the phase moles tend to drift off without this step. -#ifdef DEBUG_MODE check_tmoles(); -#endif vcs_tmoles(); // COPY OLD into NEW and ZERO VECTORS // Copy the old solution into the new solution as an initial guess @@ -411,7 +405,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, int soldel = vcs_popPhaseRxnStepSizes(iphasePop); if (soldel == 3) { iphasePop = npos; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich " "prevented phase %d popping\n"); } @@ -428,7 +422,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // minor species in the future too) kspec = npos; iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf(" --- vcs_RxnStepSizes not called because alternative" "phase creation delta was used instead\n"); } @@ -453,7 +447,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // m_molNumSpecies_old[kspec] m_molNumSpecies_new[kspec] Species Mole Numbers // m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers if (iphaseDelete != npos) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Main Loop Treatment -> Circumvented due to Phase Deletion "); plogendl(); } @@ -490,7 +484,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // only if ITI NE 0 vcs_deltag(0, false, VCS_STATECALC_NEW); } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Main Loop Treatment of each non-component species "); if (iti == 0) { plogf("- Full Calculation:\n"); @@ -505,18 +499,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); size_t iph = m_phaseID[kspec]; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; - if (DEBUG_MODE_ENABLED) { - ANOTE[0] = '\0'; - } + ANOTE[0] = '\0'; double dx; if (iphasePop != npos) { if (iph == iphasePop) { dx = m_deltaMolNumSpecies[kspec]; m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Phase pop"); - } + sprintf(ANOTE, "Phase pop"); } else { dx = 0.0; m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; @@ -529,7 +519,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { // ZEROED OUT SPECIES bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) { + if (m_debug_print_lvl >= 3) { plogf(" --- %s currently zeroed (SpStatus=%-2d):", m_speciesName[kspec], m_speciesStatus[kspec]); plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n", @@ -540,17 +530,15 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; resurrect = false; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]); - if (m_deltaGRxn_new[irxn] < 0.0) { - if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { - sprintf(ANOTE, "Species stays zeroed even though dg neg due to " - "STOICH/PHASEPOP constraint: DG = %11.4E", - m_deltaGRxn_new[irxn]); - } else { - sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed", - m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]); + if (m_deltaGRxn_new[irxn] < 0.0) { + if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { + sprintf(ANOTE, "Species stays zeroed even though dg neg due to " + "STOICH/PHASEPOP constraint: DG = %11.4E", + m_deltaGRxn_new[irxn]); + } else { + sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed", + m_deltaGRxn_new[irxn]); } } } else { @@ -561,11 +549,9 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Species stays zeroed even though dG " - "neg, because of %s elemAbund", - m_elementName[j].c_str()); - } + sprintf(ANOTE, "Species stays zeroed even though dG " + "neg, because of %s elemAbund", + m_elementName[j].c_str()); resurrect = false; break; } @@ -577,14 +563,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // Resurrect the species if (resurrect) { if (Vphase->exists() == VCS_PHASE_EXIST_NO) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Zeroed species changed to major: "); plogf("%-12s\n", m_speciesName[kspec]); } m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; allMinorZeroedSpecies = false; } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Zeroed species changed to minor: "); plogf("%-12s\n", m_speciesName[kspec]); } @@ -598,9 +584,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec]; } m_deltaMolNumSpecies[kspec] = dx; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]); - } + sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]); } else { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; @@ -615,16 +599,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE,"minor species not considered"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); - plogf("%-12s", m_speciesName[kspec]); - plogf("%3d%11.4E%11.4E%11.4E | %s", - m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - plogendl(); - } + sprintf(ANOTE,"minor species not considered"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); + plogf("%-12s", m_speciesName[kspec]); + plogf("%3d%11.4E%11.4E%11.4E | %s", + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + plogendl(); } continue; } @@ -648,7 +630,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, if (soldel_ret) { // DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF // MOLE NUMBER - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Delete minor species in multispec phase: %-12s", m_speciesName[kspec]); plogendl(); @@ -675,9 +657,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } else { // MAJOR SPECIES - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "Normal Major Calc"); - } + sprintf(ANOTE, "Normal Major Calc"); // Check for superconvergence of the formation reaction. Do // nothing if it is superconverged. Skip to the end of the irxn @@ -686,16 +666,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "major species is converged"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); - plogf("%-12s", m_speciesName[kspec]); - plogf("%3d%11.4E%11.4E%11.4E | %s", - m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - plogendl(); - } + sprintf(ANOTE, "major species is converged"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); + plogf("%-12s", m_speciesName[kspec]); + plogf("%3d%11.4E%11.4E%11.4E | %s", + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + plogendl(); } continue; } @@ -712,10 +690,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } else { dx = 0.0; m_deltaMolNumSpecies[kspec] = 0.0; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "dx set to 0, DG flipped sign due to " - "changed initial point"); - } + sprintf(ANOTE, "dx set to 0, DG flipped sign due to " + "changed initial point"); } //Form a tentative value of the new species moles @@ -726,10 +702,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // the outcome, we branch to sections below, or we restart the // entire iteration. if (m_molNumSpecies_new[kspec] <= 0.0) { - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "initial nonpos kmoles= %11.3E", - m_molNumSpecies_new[kspec]); - } + sprintf(ANOTE, "initial nonpos kmoles= %11.3E", + m_molNumSpecies_new[kspec]); // NON-POSITIVE MOLES OF MAJOR SPECIES // // We are here when a tentative value of a mole fraction @@ -768,19 +742,15 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; if (m_molNumSpecies_new[kspec] > 0.0) { m_deltaMolNumSpecies[kspec] = dx; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, - "zeroing SS phase created a neg component species " - "-> reducing step size instead"); - } + sprintf(ANOTE, + "zeroing SS phase created a neg component species " + "-> reducing step size instead"); } else { // We are going to zero the single species phase. // Set the existence flag iph = m_phaseID[kspec]; Vphase = m_VolPhaseList[iph]; - if (DEBUG_MODE_ENABLED) { - sprintf(ANOTE, "zeroing out SS phase: "); - } + sprintf(ANOTE, "zeroing out SS phase: "); // Change the base mole numbers for the iteration. // We need to do this here, because we have decided @@ -789,13 +759,11 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, m_molNumSpecies_new[kspec] = 0.0; doPhaseDeleteIph = iph; - if (DEBUG_MODE_ENABLED) { - doPhaseDeleteKspec = kspec; - if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) { - plogf(" --- SS species changed to zeroedss: "); - plogf("%-12s", m_speciesName[kspec]); - plogendl(); - } + doPhaseDeleteKspec = kspec; + if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) { + plogf(" --- SS species changed to zeroedss: "); + plogf("%-12s", m_speciesName[kspec]); + plogendl(); } m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; ++m_numRxnMinorZeroed; @@ -842,10 +810,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } - if (DEBUG_MODE_ENABLED) { - checkDelta1(&m_deltaMolNumSpecies[0], - &m_deltaPhaseMoles[0], kspec+1); - } + checkDelta1(&m_deltaMolNumSpecies[0], + &m_deltaPhaseMoles[0], kspec+1); // Branch point for returning if (m_debug_print_lvl >= 2) { @@ -860,7 +826,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } if (doPhaseDeleteIph != npos) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-12.12s Main Loop Special Case deleting phase with species: ", m_speciesName[doPhaseDeleteKspec]); @@ -869,7 +835,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, break; } } // END OF MAIN LOOP OVER FORMATION REACTIONS - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { for (size_t k = 0; k < m_numComponents; k++) { plogf(" --- "); plogf("%-12.12s", m_speciesName[k]); @@ -888,7 +854,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria // to limit its magnitude. double par = 0.5; - size_t ll; // only used in DEBUG_MODE + size_t ll; for (size_t k = 0; k < m_numComponents; ++k) { if (m_molNumSpecies_old[k] > 0.0) { double xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k]; @@ -908,7 +874,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, if (par <= 1.01 && par > 0.0) { // Reduce the size of the step by the multiplicative factor, par par *= 0.99; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Reduction in step size due to component "); plogf("%s", m_speciesName[ll]); plogf(" going negative = %11.3E", par); @@ -923,10 +889,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } else { par = 1.0; } - if (DEBUG_MODE_ENABLED) { - checkDelta1(&m_deltaMolNumSpecies[0], - &m_deltaPhaseMoles[0], m_numSpeciesTot); - } + checkDelta1(&m_deltaMolNumSpecies[0], + &m_deltaPhaseMoles[0], m_numSpeciesTot); // Now adjust the wt[kspec]'s so that the reflect the decrease in the // overall length of m_deltaMolNumSpecies[kspec] just calculated. At the @@ -1074,7 +1038,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0], &m_tPhaseMoles_new[0])); plogendl(); - if (DEBUG_MODE_ENABLED && m_VCount->Its > 550) { + if (m_VCount->Its > 550) { plogf(" --- Troublesome solve"); plogendl(); } @@ -1101,7 +1065,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // Increment the iteration counters ++m_VCount->Its; ++it1; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Increment counter increased, step is accepted: %4d", m_VCount->Its); plogendl(); @@ -1117,7 +1081,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, for (size_t iph = 0; iph < m_numPhases; iph++) { if (!m_VolPhaseList[iph]->m_singleSpecies && m_tPhaseMoles_old[iph] != 0.0 && m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) { + if (m_debug_print_lvl >= 1) { plogf(" --- Setting microscopic phase %d to zero", iph); plogendl(); } @@ -1146,12 +1110,12 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } // CHECK FOR ELEMENT ABUNDANCE - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Normal element abundance check"); } vcs_elab(); if (! vcs_elabcheck(0)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" - failed -> redoing element abundances."); plogendl(); } @@ -1160,7 +1124,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); vcs_deltag(0, true, VCS_STATECALC_OLD); uptodate_minors = true; - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf(" - passed"); plogendl(); } @@ -1192,7 +1156,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because "); plogf("%s", m_speciesName[k]); plogf(" is better than comp "); @@ -1206,7 +1170,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Check for an optimum basis passed"); plogendl(); } @@ -1219,7 +1183,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, if (iti != 0) { return; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Reevaluate major-minor status of noncomponents:\n"); } m_numRxnMinorZeroed = 0; @@ -1227,13 +1191,13 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, size_t kspec = m_indexRxnToSpecies[irxn]; int speciesType = vcs_species_type(kspec); if (speciesType < VCS_SPECIES_MINOR) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { + if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { plogf(" --- major/minor species is now zeroed out: %s\n", m_speciesName[kspec]); } ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MINOR) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { + if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) { plogf(" --- Noncomponent turned from major to minor: "); } else if (kspec < m_numComponents) { @@ -1247,7 +1211,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MAJOR) { if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { plogf(" --- Noncomponent turned from minor to major: "); } else if (kspec < m_numComponents) { @@ -1275,7 +1239,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, int& stage, bool& lec) { if (! allMinorZeroedSpecies) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Equilibrium check for major species: "); } for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { @@ -1288,7 +1252,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, // -> we have run out of iterations! stage = RETURN_A; } else { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]); } // Convergence amongst major species has not been achieved. @@ -1299,11 +1263,11 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, return; } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" MAJOR SPECIES CONVERGENCE achieved"); plogendl(); } - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf(" MAJOR SPECIES CONVERGENCE achieved " "(because there are no major species)"); plogendl(); @@ -1320,7 +1284,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, vcs_deltag(1, false, VCS_STATECALC_OLD); uptodate_minors = true; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Equilibrium check for minor species: "); } for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { @@ -1333,7 +1297,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, stage = RETURN_A; return; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]); } @@ -1344,7 +1308,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, return; } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" CONVERGENCE achieved\n"); } } @@ -1357,14 +1321,14 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, // LEC is only true when we are near the end game if (lec) { if (!giveUpOnElemAbund) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Check the Full Element Abundances: "); } // Final element abundance check: If we fail then we need to go back // and correct the element abundances, and then go do a major step if (! vcs_elabcheck(1)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (! vcs_elabcheck(0)) { plogf(" failed\n"); } else { @@ -1375,7 +1339,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies, stage = ELEM_ABUND_CHECK; return; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" passed\n"); } } @@ -1437,7 +1401,7 @@ void VCS_SOLVE::solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec, } else if (ncAfter) { if (!neAfter) { // Probably an unrecoverable range error - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n"); plogf(" --- vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n"); plogf(" --- vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n"); @@ -1477,7 +1441,7 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF; } dg_irxn = std::max(dg_irxn, -200.0); - if (DEBUG_MODE_ENABLED && ANOTE) { + if (ANOTE) { sprintf(ANOTE,"minor species alternative calc"); } if (dg_irxn >= 23.0) { @@ -1542,7 +1506,7 @@ L_ZERO_SPECIES: // Voltage calculation // Need to check the sign -> This is good for electrons dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim; - if (DEBUG_MODE_ENABLED && ANOTE) { + if (ANOTE) { sprintf(ANOTE,"voltage species alternative calc"); } } @@ -1608,7 +1572,7 @@ int VCS_SOLVE::vcs_zero_species(const size_t kspec) double dx = -m_molNumSpecies_old[kspec]; if (dx != 0.0) { retn = delta_species(kspec, &dx); - if (DEBUG_MODE_ENABLED && !retn && m_debug_print_lvl >= 1) { + if (!retn && m_debug_print_lvl >= 1) { plogf("vcs_zero_species: Couldn't zero the species %d, " "did delta of %g. orig conc of %g", kspec, dx, m_molNumSpecies_old[kspec] + dx); @@ -1683,7 +1647,7 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec) void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) { size_t iph = m_phaseID[kspec]; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Add back a deleted species: %-12s\n", m_speciesName[kspec]); } @@ -1738,7 +1702,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) // set the phase existence flag to dead Vphase->setTotalMoles(0.0); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName); } @@ -1753,7 +1717,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) int retn = delta_species(kspec, &dxTent); if (retn != 1) { successful = false; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n", iph, Vphase->PhaseName, m_speciesName[kspec]); plogf(" --- delta attempted: %g achieved: %g " @@ -1781,7 +1745,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0; for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { if (m_phaseID[kcomp] == iph) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n", kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]); } @@ -1818,7 +1782,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) } if (m_molNumSpecies_old[kcomp] != 0.0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n", kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]); plogf(" --- zeroing it \n"); @@ -1843,7 +1807,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) ++m_numRxnRdc; ++m_numSpeciesRdc; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Make %s", m_speciesName[kspec]); plogf(" an active but zeroed species because its phase " "was zeroed\n"); @@ -1868,7 +1832,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) int VCS_SOLVE::vcs_recheck_deleted() { vector_fp& xtcutoff = m_TmpPhase; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Start rechecking deleted species in multispec phases\n"); } if (m_numSpeciesRdc == m_numSpeciesTot) { @@ -2019,20 +1983,20 @@ size_t VCS_SOLVE::vcs_add_all_deleted() double dx = m_molNumSpecies_new[kspec]; retn = delta_species(kspec, &dx); if (retn == 0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl) { + if (m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", m_speciesName[kspec], kspec, dx); } if (dx > 1.0E-50) { dx = 1.0E-50; retn = delta_species(kspec, &dx); - if (DEBUG_MODE_ENABLED && retn == 0 && m_debug_print_lvl) { + if (retn == 0 && m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", m_speciesName[kspec], kspec, dx); } } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (retn != 0) { plogf(" --- add_deleted(): species %s added back in with mol number %g", m_speciesName[kspec], dx); @@ -2057,7 +2021,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() VCS_DELETE_MINORSPECIES_CUTOFF) || (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { retn++; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", m_speciesName[kspec], m_molNumSpecies_old[kspec], m_deltaGRxn_old[irxn]); @@ -2092,13 +2056,13 @@ bool VCS_SOLVE::vcs_globStepDamp() } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1); plogf(" --- subroutine FORCE: End Slope = %g\n", s2); } if (s1 > 0.0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments,"); if (s1 < 1.0E-40) { plogf(" s1 positive but really small"); @@ -2111,7 +2075,7 @@ bool VCS_SOLVE::vcs_globStepDamp() } if (s2 <= 0.0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments, s2 < 0"); plogendl(); } @@ -2124,17 +2088,17 @@ bool VCS_SOLVE::vcs_globStepDamp() al = s1 / (s1 - s2); } if (al >= 0.95 || al < 0.0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments (al = %g)\n", al); } return false; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced a damping factor = %g\n", al); } // ADJUST MOLE NUMBERS, CHEM. POT - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { m_deltaGRxn_tmp = m_deltaGRxn_new; } @@ -2148,7 +2112,7 @@ bool VCS_SOLVE::vcs_globStepDamp() } vcs_updateVP(VCS_STATECALC_NEW); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE adjusted the mole " "numbers, AL = %10.3f\n", al); } @@ -2172,7 +2136,7 @@ bool VCS_SOLVE::vcs_globStepDamp() } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE: Adj End Slope = %g", s2); plogendl(); } @@ -2187,7 +2151,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], size_t jlose = npos; double* scrxn_ptr; clockWC tickTock; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" "); for (size_t i=0; i<77; i++) { plogf("-"); @@ -2363,7 +2327,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], for (size_t i = 0; i < m_numSpeciesTot; ++i) { m_indexRxnToSpecies[i] = ncTrial + i; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Total number of components found = %3d (ne = %d)\n ", ncTrial, m_numElemConstraints); } @@ -2417,7 +2381,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], // REARRANGE THE DATA if (jr != k) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[k]); if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("(Volts = %9.2g)", m_molNumSpecies_old[k]); @@ -2434,7 +2398,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], } vcs_switch_pos(false, jr, k); std::swap(aw[jr], aw[k]); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[k]); if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("(Volts = %9.2g) remains ", m_molNumSpecies_old[k]); @@ -2560,7 +2524,7 @@ L_END_LOOP: m_scSize[i] = szTmp; } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Components:"); for (size_t j = 0; j < ncTrial; j++) { plogf(" %3d", j); @@ -2756,7 +2720,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %s can not be nonzero because" " needed element %s is zero\n", m_speciesName[kspec], m_elementName[j]); @@ -2782,7 +2746,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const double negChangeComp = - stoicC; if (negChangeComp > 0.0) { if (m_molNumSpecies_old[j] < 1.0E-60) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %s is prevented from popping into existence because" " a needed component to be consumed, %s, has a zero mole number\n", m_speciesName[kspec], m_speciesName[j]); @@ -2797,7 +2761,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const size_t jph = m_phaseID[j]; vcs_VolPhase* jVPhase = m_VolPhaseList[jph]; if (jVPhase->exists() <= 0) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- %s is prevented from popping into existence because" " a needed component %s is in a zeroed-phase that would be " "popped into existence at the same time\n", @@ -2954,7 +2918,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, tPhMoles_ptr = &m_tPhaseMoles_new[0]; actCoeff_ptr = &m_actCoeffSpecies_new[0]; molNum = &m_molNumSpecies_new[0]; - } else if (DEBUG_MODE_ENABLED) { + } else { throw CanteraError("VCS_SOLVE::vcs_dfe", "Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc); } @@ -2962,7 +2926,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, AssertThrowMsg(m_unitsState != VCS_DIMENSIONAL_G, "VCS_SOLVE::vcs_dfe", "called with wrong units state"); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { if (ll == 0) { if (lbot != 0) { plogf(" --- Subroutine vcs_dfe called for one species: "); @@ -3239,7 +3203,6 @@ void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const writeline('-', 132); } -#ifdef DEBUG_MODE void VCS_SOLVE::prneav() const { vector_fp eav(m_numElemConstraints, 0.0); @@ -3276,7 +3239,6 @@ void VCS_SOLVE::prneav() const plogf("--------------------------------------------------"); plogendl(); } -#endif double VCS_SOLVE::l2normdg(double dgLocal[]) const { @@ -3320,7 +3282,6 @@ double VCS_SOLVE::vcs_tmoles() return m_totalMolNum; } -#ifdef DEBUG_MODE void VCS_SOLVE::check_tmoles() const { double sum = 0.0; @@ -3341,7 +3302,6 @@ void VCS_SOLVE::check_tmoles() const } } } -#endif void VCS_SOLVE::vcs_updateVP(const int vcsState) { @@ -3355,7 +3315,7 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState) Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, &m_molNumSpecies_new[0], &m_tPhaseMoles_new[0]); - } else if (DEBUG_MODE_ENABLED) { + } else { throw CanteraError("VCS_SOLVE::vcs_updateVP", "wrong stateCalc value: {}", vcsState); } @@ -3365,15 +3325,15 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState) bool VCS_SOLVE::vcs_evaluate_speciesType() { m_numRxnMinorZeroed = 0; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Species Status decision is reevaluated: All species are minor except for:\n"); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) { + } else if (m_debug_print_lvl >= 5) { plogf(" --- Species Status decision is reevaluated"); plogendl(); } for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { m_speciesStatus[kspec] = vcs_species_type(kspec); - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) { + if (m_debug_print_lvl >= 5) { plogf(" --- %-16s: ", m_speciesName[kspec]); if (kspec < m_numComponents) { plogf("(COMP) "); @@ -3383,7 +3343,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() plogf(" %10.3g ", m_molNumSpecies_old[kspec]); const char* sString = vcs_speciesType_string(m_speciesStatus[kspec], 100); plogf("%s\n", sString); - } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + } else if (m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { switch (m_speciesStatus[kspec]) { case VCS_SPECIES_COMPONENT: @@ -3425,7 +3385,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() ++m_numRxnMinorZeroed; } } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" ---"); plogendl(); } @@ -3459,7 +3419,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, throw CanteraError("VCS_SOLVE::vcs_deltag", "bad vcsState"); } - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Subroutine vcs_deltag called for "); if (L < 0) { plogf("major noncomponents\n"); @@ -3755,7 +3715,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, feSpecies = &m_feSpecies_old[0]; deltaGRxn = &m_deltaGRxn_old[0]; actCoeffSpecies = &m_actCoeffSpecies_old[0]; - } else if (DEBUG_MODE_ENABLED) { + } else { throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc"); } @@ -3765,7 +3725,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, } vcs_VolPhase* vPhase = m_VolPhaseList[iphase]; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + if (m_debug_print_lvl >= 2) { plogf(" --- Subroutine vcs_deltag_Phase called for phase %d\n", iphase); } @@ -3858,8 +3818,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k if (k1 == k2) { return; } - if (DEBUG_MODE_ENABLED && (k1 >= m_numSpeciesTot || - k2 >= m_numSpeciesTot)) { + if (k1 >= m_numSpeciesTot || k2 >= m_numSpeciesTot) { plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n", k1, k2); } @@ -3915,7 +3874,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k // Find the Rxn indices corresponding to the two species size_t i1 = k1 - m_numComponents; size_t i2 = k2 - m_numComponents; - if (DEBUG_MODE_ENABLED && (i1 > m_numRxnTot || i2 >= m_numRxnTot)) { + if (i1 > m_numRxnTot || i2 >= m_numRxnTot) { plogf("switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n", i1 , i2); } diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 1e5fe2fcf..ea7cb99a8 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -8,14 +8,6 @@ */ #include "cantera/numerics/RootFind.h" - -// turn on debugging for now -#ifndef DEBUG_MODE -#define DEBUG_MODE -#undef DEBUG_MODE_ENABLED -#define DEBUG_MODE_ENABLED 1 -#endif - #include "cantera/base/utilities.h" #include "cantera/base/stringUtils.h" @@ -201,7 +193,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfT.reasoning = "First Point: "; callNum++; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fp = fopen(fmt::format("RootFind_%d.log", callNum).c_str(), "w"); fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n"); fprintf(fp, "-----------------------------------------------------" @@ -256,7 +248,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun int its = 1; f1 = func(x1); - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { print_funcEval(fp, x1, f1, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1); } @@ -305,7 +297,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun deltaX2 = x2 - x1; its++; f2 = func(x2); - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { print_funcEval(fp, x2, f2, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2); } @@ -359,7 +351,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun while (!converged && its < itmax) { // Find an estimate of the next point, xnew, to try based on a linear // approximation from the last two points. - if (DEBUG_MODE_ENABLED && fabs(x2 - x1) < 1.0E-14) { + if (fabs(x2 - x1) < 1.0E-14) { writelogf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1); } doublereal delXtmp = deltaXControlled(x2, x1); @@ -387,7 +379,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } rfT.reasoning += "Slope is good. "; } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlin = %-11.5E", xnew); } deltaXnew = xnew - x2; @@ -478,13 +470,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xDelMin = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < xDelMin) { xnew = x1 + sign(xnew-x1) * xDelMin; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } if (fabs(xnew - x2) < 0.1 * xDelMin) { xnew = x2 + sign(xnew-x2) * 0.1 * xDelMin; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } @@ -498,7 +490,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } if (fabs(xDelMax) < fabs(xnew - x2)) { xnew = x2 + sign(xnew-x2) * xDelMax; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitsize = %-11.5E", xnew); } } @@ -510,13 +502,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xDelMin = 0.1 * fabs(x2 - x1); if (fabs(xnew - x2) < xDelMin) { xnew = x2 + sign(xnew - x2) * xDelMin; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } if (fabs(xnew - x1) < xDelMin) { xnew = x1 + sign(xnew - x1) * xDelMin; - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } @@ -558,7 +550,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_ && xorig != xnew) { + if (printLvl >= 3 && writeLogAllowed_ && xorig != xnew) { fprintf(fp, " | xstraddle = %-11.5E", xnew); } } @@ -596,7 +588,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = xmax; } } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitmax = %-11.5E", xnew); } } @@ -620,7 +612,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = xmin; } } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitmin = %-11.5E", xnew); } } @@ -628,7 +620,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun its++; fnew = func(xnew); - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp,"\n"); print_funcEval(fp, xnew, fnew, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew); @@ -947,7 +939,7 @@ done: if (printLvl >= 1) { writelogf("RootFind success: convergence achieved\n"); } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | RootFind success in %d its, fnorm = %g\n", its, fnorm); } rfHistory_.push_back(rfT); @@ -972,7 +964,7 @@ done: } rfT.reasoning += "Maximum iterations exceeded without convergence, cause unknown"; } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fprintf(fp, "\nRootFind failure in %d its\n", its); } @@ -982,7 +974,7 @@ done: rfT.fval = f2; rfHistory_.push_back(rfT); } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + if (printLvl >= 3 && writeLogAllowed_) { fclose(fp); } @@ -996,13 +988,9 @@ done: doublereal RootFind::func(doublereal x) { doublereal r; - if (DEBUG_MODE_ENABLED) { - checkFinite(x); - } + checkFinite(x); m_residFunc->evalSS(0.0, &x, &r); - if (DEBUG_MODE_ENABLED) { - checkFinite(r); - } + checkFinite(r); doublereal ff = r - m_funcTargetValue; if (x >= x_maxTried_) { x_maxTried_ = x; diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index cfc8c0f49..9e39f9279 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -1492,7 +1492,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // molalitysum is the sum of the molalities over all solutes, even those // with zero charge. double molalitysumUncropped = 0.0; - debuglog("\n Debugging information from hmw_act \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog("\n Debugging information from hmw_act \n", m_debugCalc); // Make sure the counter variables are setup counterIJ_setup(); @@ -1510,7 +1510,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // Store the ionic molality in the object for reference. m_IionicMolality = Is; sqrtIs = sqrt(Is); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelog(" Step 1: \n"); writelogf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); @@ -1522,11 +1522,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // Step 2: Find the coefficients E-theta and E-thetaprime for all // combinations of positive unlike charges up to 4 - debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 2: \n", m_debugCalc); for (int z1 = 1; z1 <=4; z1++) { for (int z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } @@ -1535,7 +1535,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const debuglog(" Step 3: \n" " Species Species g(x) hfunc(x)\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // calculate g(x) and hfunc(x) for each cation-anion pair MX. In the // original literature, hfunc, was called gprime. However, it's not the @@ -1574,7 +1574,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const m_gfunc_IJ[counterIJ] = 0.0; m_hfunc_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i), speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]); } @@ -1585,7 +1585,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // Agrees with Pitzer, Eq. (49), (51), (55) debuglog(" Step 4: \n" " Species Species BMX BprimeMX BphiMX\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk - 1; i++) { for (size_t j = i+1; j < m_kk; j++) { @@ -1600,7 +1600,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const + m_Beta1MX_ij[counterIJ] * m_gfunc_IJ[counterIJ] + m_Beta2MX_ij[counterIJ] * m_g2func_IJ[counterIJ]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf("%d %g: %g %g %g %g\n", counterIJ, m_BMX_IJ[counterIJ], m_Beta0MX_ij[counterIJ], m_Beta1MX_ij[counterIJ], m_Beta2MX_ij[counterIJ], m_gfunc_IJ[counterIJ]); @@ -1617,7 +1617,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const m_BprimeMX_IJ[counterIJ] = 0.0; m_BphiMX_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n", speciesName(i), speciesName(j), m_BMX_IJ[counterIJ], m_BprimeMX_IJ[counterIJ], m_BphiMX_IJ[counterIJ]); @@ -1628,8 +1628,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // SUBSECTION TO CALCULATE CMX // Agrees with Pitzer, Eq. (53). debuglog(" Step 5: \n" - " Species Species CMX\n", - DEBUG_MODE_ENABLED && m_debugCalc); + " Species Species CMX\n", m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -1644,7 +1643,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } else { m_CMX_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f \n", speciesName(i), speciesName(j), m_CMX_IJ[counterIJ]); } @@ -1655,7 +1654,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // Agrees with Pitzer, Eq. 72, 73, 74 debuglog(" Step 6: \n" " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -1675,7 +1674,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const m_Phiprime_IJ[counterIJ] = 0.0; m_PhiPhi_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n", speciesName(i), speciesName(j), m_Phi_IJ[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ[counterIJ]); @@ -1685,11 +1684,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // SUBSECTION FOR CALCULATION OF F // Agrees with Pitzer Eqn. (65) - debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 7: \n", m_debugCalc); double Aphi = A_Debye_TP() / 3.0; double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" initial value of F = %10.6f \n", F); } for (size_t i = 1; i < m_kk-1; i++) { @@ -1709,13 +1708,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(i)*charge(j) > 0) { F += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ]; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" F = %10.6f \n", F); } } } debuglog(" Step 8: Summing in All Contributions to Activity Coefficients \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk; i++) { @@ -1723,12 +1722,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // equations agree with my notes, Eqn. (118). // Equations agree with Pitzer, eqn.(63) if (charge(i) > 0.0) { - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i)); } // species i is the cation (positive) to calc the actcoeff double zsqF = charge(i)*charge(i)*F; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Unary term: z*z*F = %10.5f\n", zsqF); } double sum1 = 0.0; @@ -1745,7 +1744,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all anions sum1 += molality[j] * (2.0*m_BMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { std::string snj = speciesName(j) + ":"; writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj, molality[j]*2.0*m_BMX_IJ[counterIJ]); @@ -1761,7 +1760,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*m_Psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) { + if (m_debugCalc && m_Psi_ijk[n] != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1775,7 +1774,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all cations if (j != i) { sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) { + if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) { std::string snj = speciesName(j) + ":"; writelogf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj, molality[j]*(2.0*m_Phi_IJ[counterIJ])); @@ -1786,7 +1785,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // two inner sums over anions n = k + j * m_kk + i * m_kk * m_kk; sum2 += molality[j]*molality[k]*m_Psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) { + if (m_debugCalc && m_Psi_ijk[n] != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1797,7 +1796,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const size_t counterIJ2 = m_CounterIJ[n]; sum4 += (fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]); - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) { + if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj, fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]); @@ -1809,7 +1808,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // Handle neutral j species if (charge(j) == 0) { sum5 += molality[j]*2.0*m_Lambda_nj(j,i); - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { + if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { std::string snj = speciesName(j) + ":"; writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj, molality[j]*2.0*m_Lambda_nj(j,i)); @@ -1824,7 +1823,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double zeta = m_Psi_ijk[n]; if (zeta != 0.0) { sum5 += molality[j]*molality[k]*zeta; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1839,7 +1838,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // activity coefficients (molality scale) m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n", speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } @@ -1849,12 +1848,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // equations agree with my notes, Eqn. (119). // Equations agree with Pitzer, eqn.(64) if (charge(i) < 0) { - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i)); } // species i is an anion (negative) double zsqF = charge(i)*charge(i)*F; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Unary term: z*z*F = %10.5f\n", zsqF); } double sum1 = 0.0; @@ -1871,7 +1870,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(j) > 0) { sum1 += molality[j]* (2.0*m_BMX_IJ[counterIJ]+molarcharge*m_CMX_IJ[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { std::string snj = speciesName(j) + ":"; writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj, molality[j]*2.0*m_BMX_IJ[counterIJ]); @@ -1884,7 +1883,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) > 0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*m_Psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) { + if (m_debugCalc && m_Psi_ijk[n] != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1899,7 +1898,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all anions if (j != i) { sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) { + if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) { std::string snj = speciesName(j) + ":"; writelogf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj, molality[j]*(2.0*m_Phi_IJ[counterIJ])); @@ -1910,7 +1909,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // two inner sums over cations n = k + j * m_kk + i * m_kk * m_kk; sum2 += molality[j]*molality[k]*m_Psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) { + if (m_debugCalc && m_Psi_ijk[n] != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1920,7 +1919,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const size_t counterIJ2 = m_CounterIJ[n]; sum4 += fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]; - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) { + if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj, fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]); @@ -1932,7 +1931,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // for Anions, do the neutral species interaction if (charge(j) == 0.0) { sum5 += molality[j]*2.0*m_Lambda_nj(j,i); - if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { + if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { std::string snj = speciesName(j) + ":"; writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj, molality[j]*2.0*m_Lambda_nj(j,i)); @@ -1947,7 +1946,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double zeta = m_Psi_ijk[n]; if (zeta != 0.0) { sum5 += molality[j]*molality[k]*zeta; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1959,7 +1958,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n", speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } @@ -1969,14 +1968,14 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // equations agree with my notes, // Equations agree with Pitzer, if (charge(i) == 0.0) { - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i)); } double sum1 = 0.0; double sum3 = 0.0; for (size_t j = 1; j < m_kk; j++) { sum1 += molality[j]*2.0*m_Lambda_nj(i,j); - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Lambda_nj(i,j) != 0.0) { + if (m_debugCalc && m_Lambda_nj(i,j) != 0.0) { std::string snj = speciesName(j) + ":"; writelogf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj, molality[j]*2.0*m_Lambda_nj(i,j)); @@ -1987,7 +1986,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { size_t n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*m_Psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) { + if (m_debugCalc && m_Psi_ijk[n] != 0.0) { std::string snj = speciesName(j) + "," + speciesName(k) + ":"; writelogf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj, molality[j]*molality[k]*m_Psi_ijk[n]); @@ -1997,19 +1996,19 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } } double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]; - if (DEBUG_MODE_ENABLED && m_debugCalc && m_Mu_nnn[i] != 0.0) { + if (m_debugCalc && m_Mu_nnn[i] != 0.0) { writelogf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n", 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]); } m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n", speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } } } - debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 9: \n", m_debugCalc); // SUBSECTION FOR CALCULATING THE OSMOTIC COEFF // equations agree with my notes, Eqn. (117). @@ -2133,7 +2132,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } else { osmotic_coef = 1.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", term1, sum1, sum2, sum3, sum4, sum5); @@ -2152,7 +2151,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double xmolSolvent = moleFraction(m_indexSolvent); double xx = std::max(m_xmolSolventMIN, xmolSolvent); m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { double wateract = exp(lnwateract); writelogf(" Weight of Solvent = %16.7g\n", m_weightSolvent); writelogf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped); @@ -2195,9 +2194,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // immediately preceding the calling of this routine. Therefore, some // quantities do not need to be recalculated in this routine. -#ifdef DEBUG_MODE - m_debugCalc = 0; -#endif // HKM -> Assumption is made that the solvent is species 0. if (m_indexSolvent != 0) { throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT", @@ -2222,7 +2218,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const double molalitysum = 0.0; debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dT()\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // Make sure the counter variables are setup counterIJ_setup(); @@ -2240,7 +2236,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // Store the ionic molality in the object for reference. m_IionicMolality = Is; sqrtIs = sqrt(Is); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelog(" Step 1: \n"); writelogf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); @@ -2252,11 +2248,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // Step 2: Find the coefficients E-theta and E-thetaprime for all // combinations of positive unlike charges up to 4 - debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 2: \n", m_debugCalc); for (int z1 = 1; z1 <=4; z1++) { for (int z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } @@ -2265,7 +2261,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const debuglog(" Step 3: \n" " Species Species g(x) hfunc(x) \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // calculate g(x) and hfunc(x) for each cation-anion pair MX // In the original literature, hfunc, was called gprime. However, @@ -2304,7 +2300,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_gfunc_IJ[counterIJ] = 0.0; m_hfunc_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { std::string sni = speciesName(i); std::string snj = speciesName(j); writelogf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), @@ -2318,7 +2314,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // quantities. debuglog(" Step 4: \n" " Species Species BMX BprimeMX BphiMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk - 1; i++) { for (size_t j = i+1; j < m_kk; j++) { @@ -2332,7 +2328,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_BMX_IJ_L[counterIJ] = m_Beta0MX_ij_L[counterIJ] + m_Beta1MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ] + m_Beta2MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf("%d %g: %g %g %g %g\n", counterIJ, m_BMX_IJ_L[counterIJ], m_Beta0MX_ij_L[counterIJ], m_Beta1MX_ij_L[counterIJ], m_Beta2MX_ij_L[counterIJ], m_gfunc_IJ[counterIJ]); @@ -2349,7 +2345,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_BprimeMX_IJ_L[counterIJ] = 0.0; m_BphiMX_IJ_L[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n", speciesName(i), speciesName(j), m_BMX_IJ_L[counterIJ], m_BprimeMX_IJ_L[counterIJ], m_BphiMX_IJ_L[counterIJ]); @@ -2359,8 +2355,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // --------- SUBSECTION TO CALCULATE CMX_L ---------- debuglog(" Step 5: \n" - " Species Species CMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + " Species Species CMX \n", m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -2375,7 +2370,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const } else { m_CMX_IJ_L[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f \n", speciesName(i), speciesName(j), m_CMX_IJ_L[counterIJ]); } @@ -2385,7 +2380,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- debuglog(" Step 6: \n" " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -2403,7 +2398,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_Phiprime_IJ[counterIJ] = 0.0; m_PhiPhi_IJ_L[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n", speciesName(i), speciesName(j), m_Phi_IJ_L[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_L[counterIJ]); @@ -2412,12 +2407,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const } // ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- - debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 7: \n", m_debugCalc); double dA_DebyedT = dA_DebyedT_TP(); double dAphidT = dA_DebyedT /3.0; double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" initial value of dFdT = %10.6f \n", dFdT); } for (size_t i = 1; i < m_kk-1; i++) { @@ -2437,12 +2432,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const if (charge(i)*charge(j) > 0) { dFdT += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ]; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" dFdT = %10.6f \n", dFdT); } } } - debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 8: \n", m_debugCalc); for (size_t i = 1; i < m_kk; i++) { // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS ----- @@ -2521,7 +2516,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_dlnActCoeffMolaldT_Unscaled[i] = zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5; d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -2598,7 +2593,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const m_dlnActCoeffMolaldT_Unscaled[i] = zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5; d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n", speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -2627,13 +2622,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i]; m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3; d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); } } } - debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 9: \n", m_debugCalc); // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT --------- double sum1 = 0.0; @@ -2755,7 +2750,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const d_osmotic_coef_dT = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", term1, sum1, sum2, sum3, sum4, sum5); @@ -2772,7 +2767,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const // We have just computed act_0. However, this routine returns // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { double d_wateract_dT = exp(d_lnwateract_dT); writelogf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n", d_lnwateract_dT, d_wateract_dT); @@ -2809,9 +2804,6 @@ void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const { -#ifdef DEBUG_MODE - m_debugCalc = 0; -#endif // HKM -> Assumption is made that the solvent is species 0. if (m_indexSolvent != 0) { throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2", @@ -2835,7 +2827,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const double molalitysum = 0.0; debuglog("\n Debugging information from s_Pitzer_d2lnMolalityActCoeff_dT2()\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // Make sure the counter variables are setup counterIJ_setup(); @@ -2853,7 +2845,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // Store the ionic molality in the object for reference. m_IionicMolality = Is; sqrtIs = sqrt(Is); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelog(" Step 1: \n"); writelogf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); @@ -2865,11 +2857,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // Step 2: Find the coefficients E-theta and E-thetaprime for all // combinations of positive unlike charges up to 4 - debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 2: \n", m_debugCalc); for (int z1 = 1; z1 <=4; z1++) { for (int z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } @@ -2878,7 +2870,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const debuglog(" Step 3: \n" " Species Species g(x) hfunc(x) \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // calculate gfunc(x) and hfunc(x) for each cation-anion pair MX. In the // original literature, hfunc, was called gprime. However, it's not the @@ -2917,7 +2909,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const m_gfunc_IJ[counterIJ] = 0.0; m_hfunc_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i), speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]); } @@ -2929,7 +2921,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // quantities. debuglog(" Step 4: \n" " Species Species BMX BprimeMX BphiMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk - 1; i++) { for (size_t j = i+1; j < m_kk; j++) { @@ -2943,7 +2935,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const m_BMX_IJ_LL[counterIJ] = m_Beta0MX_ij_LL[counterIJ] + m_Beta1MX_ij_LL[counterIJ] * m_gfunc_IJ[counterIJ] + m_Beta2MX_ij_LL[counterIJ] * m_g2func_IJ[counterIJ]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf("%d %g: %g %g %g %g\n", counterIJ, m_BMX_IJ_LL[counterIJ], m_Beta0MX_ij_LL[counterIJ], m_Beta1MX_ij_LL[counterIJ], m_Beta2MX_ij_LL[counterIJ], m_gfunc_IJ[counterIJ]); @@ -2960,7 +2952,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const m_BprimeMX_IJ_LL[counterIJ] = 0.0; m_BphiMX_IJ_LL[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n", speciesName(i), speciesName(j), m_BMX_IJ_LL[counterIJ], m_BprimeMX_IJ_LL[counterIJ], m_BphiMX_IJ_LL[counterIJ]); @@ -2970,8 +2962,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // --------- SUBSECTION TO CALCULATE CMX_LL ---------- debuglog(" Step 5: \n" - " Species Species CMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + " Species Species CMX \n", m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -2986,7 +2977,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } else { m_CMX_IJ_LL[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f \n", speciesName(i), speciesName(j), m_CMX_IJ_LL[counterIJ]); } @@ -2996,7 +2987,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- debuglog(" Step 6: \n" " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -3014,7 +3005,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const m_Phiprime_IJ[counterIJ] = 0.0; m_PhiPhi_IJ_LL[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n", speciesName(i), speciesName(j), m_Phi_IJ_LL[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_LL[counterIJ]); @@ -3023,11 +3014,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } // ----------- SUBSECTION FOR CALCULATION OF d2FdT2 --------------------- - debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 7: \n", m_debugCalc); double d2AphidT2 = d2A_DebyedT2_TP() / 3.0; double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" initial value of d2FdT2 = %10.6f \n", d2FdT2); } for (size_t i = 1; i < m_kk-1; i++) { @@ -3047,12 +3038,12 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const if (charge(i)*charge(j) > 0) { d2FdT2 += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ]; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" d2FdT2 = %10.6f \n", d2FdT2); } } } - debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 8: \n", m_debugCalc); for (size_t i = 1; i < m_kk; i++) { // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS ----- @@ -3128,7 +3119,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // solute activity coefficients (molality scale) m_d2lnActCoeffMolaldT2_Unscaled[i] = zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s d2lngammadT2[i]=%10.6f \n", speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -3205,7 +3196,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } m_d2lnActCoeffMolaldT2_Unscaled[i] = zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s d2lngammadT2[i]=%10.6f\n", speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -3233,14 +3224,14 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i]; m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelog(" %-16s d2lngammadT2[i]=%10.6f \n", speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]); } } } - debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 9: \n", m_debugCalc); // ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 --------- double sum1 = 0.0; @@ -3364,7 +3355,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } else { d2_osmotic_coef_dT2 = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", term1, sum1, sum2, sum3, sum4, sum5); @@ -3382,7 +3373,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { double d2_wateract_dT2 = exp(d2_lnwateract_dT2); writelogf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n", d2_lnwateract_dT2, d2_wateract_dT2); @@ -3415,9 +3406,6 @@ void HMWSoln::s_update_dlnMolalityActCoeff_dP() const void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const { -#ifdef DEBUG_MODE - m_debugCalc = 0; -#endif // HKM -> Assumption is made that the solvent is species 0. if (m_indexSolvent != 0) { throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP", @@ -3443,7 +3431,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const double currPres = pressure(); debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dP()\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // Make sure the counter variables are setup counterIJ_setup(); @@ -3461,7 +3449,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // Store the ionic molality in the object for reference. m_IionicMolality = Is; sqrtIs = sqrt(Is); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelog(" Step 1: \n"); writelogf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); @@ -3474,11 +3462,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // Step 2: Find the coefficients E-theta and E-thetaprime for all // combinations of positive unlike charges up to 4 - debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 2: \n", m_debugCalc); for (int z1 = 1; z1 <=4; z1++) { for (int z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } @@ -3487,7 +3475,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const debuglog(" Step 3: \n" " Species Species g(x) hfunc(x)\n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); // calculate g(x) and hfunc(x) for each cation-anion pair MX // In the original literature, hfunc, was called gprime. However, @@ -3526,7 +3514,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_gfunc_IJ[counterIJ] = 0.0; m_hfunc_IJ[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i), speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]); } @@ -3538,7 +3526,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // quantities. debuglog(" Step 4: \n" " Species Species BMX BprimeMX BphiMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk - 1; i++) { for (size_t j = i+1; j < m_kk; j++) { @@ -3552,7 +3540,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_BMX_IJ_P[counterIJ] = m_Beta0MX_ij_P[counterIJ] + m_Beta1MX_ij_P[counterIJ] * m_gfunc_IJ[counterIJ] + m_Beta2MX_ij_P[counterIJ] * m_g2func_IJ[counterIJ]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf("%d %g: %g %g %g %g\n", counterIJ, m_BMX_IJ_P[counterIJ], m_Beta0MX_ij_P[counterIJ], m_Beta1MX_ij_P[counterIJ], m_Beta2MX_ij_P[counterIJ], m_gfunc_IJ[counterIJ]); @@ -3569,7 +3557,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_BprimeMX_IJ_P[counterIJ] = 0.0; m_BphiMX_IJ_P[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n", speciesName(i), speciesName(j), m_BMX_IJ_P[counterIJ], m_BprimeMX_IJ_P[counterIJ], m_BphiMX_IJ_P[counterIJ]); @@ -3579,8 +3567,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // --------- SUBSECTION TO CALCULATE CMX_P ---------- debuglog(" Step 5: \n" - " Species Species CMX \n", - DEBUG_MODE_ENABLED && m_debugCalc); + " Species Species CMX \n", m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -3595,7 +3582,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } else { m_CMX_IJ_P[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %11.7f \n", speciesName(i), speciesName(j), m_CMX_IJ_P[counterIJ]); } @@ -3605,7 +3592,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- debuglog(" Step 6: \n" " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n", - DEBUG_MODE_ENABLED && m_debugCalc); + m_debugCalc); for (size_t i = 1; i < m_kk-1; i++) { for (size_t j = i+1; j < m_kk; j++) { // Find the counterIJ for the symmetric binary interaction @@ -3623,7 +3610,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_Phiprime_IJ[counterIJ] = 0.0; m_PhiPhi_IJ_P[counterIJ] = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n", speciesName(i), speciesName(j), m_Phi_IJ_P[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_P[counterIJ]); @@ -3632,12 +3619,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } // ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- - debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 7: \n", m_debugCalc); double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres); double dAphidP = dA_DebyedP /3.0; double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" initial value of dFdP = %10.6f \n", dFdP); } for (size_t i = 1; i < m_kk-1; i++) { @@ -3657,12 +3644,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const if (charge(i)*charge(j) > 0) { dFdP += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ]; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" dFdP = %10.6f \n", dFdP); } } } - debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 8: \n", m_debugCalc); for (size_t i = 1; i < m_kk; i++) { // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS ----- @@ -3740,7 +3727,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_dlnActCoeffMolaldP_Unscaled[i] = zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s lngamma[i]=%10.6f \n", speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -3817,7 +3804,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } m_dlnActCoeffMolaldP_Unscaled[i] = zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n", speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]); writelogf(" %12g %12g %12g %12g %12g %12g\n", @@ -3843,13 +3830,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i]; m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n", speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]); } } } - debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc); + debuglog(" Step 9: \n", m_debugCalc); // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP --------- double sum1 = 0.0; @@ -3973,7 +3960,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } else { d_osmotic_coef_dP = 0.0; } - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f " "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n", term1, sum1, sum2, sum3, sum4, sum5); @@ -3990,7 +3977,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const // We have just computed act_0. However, this routine returns // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n", d_lnwateract_dP, exp(d_lnwateract_dP)); } @@ -4007,7 +3994,7 @@ void HMWSoln::calc_lambdas(double is) const // aphi is the Debye-Huckel constant at 25 C double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528; double aphi = 0.392; /* Value at 25 C */ - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" Is = %g\n", is); } if (is < 1.0E-150) { @@ -4039,7 +4026,7 @@ void HMWSoln::calc_lambdas(double is) const elambda[ij] = zprod*jfunc / (4.0*is); // eqn 14 elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is)) - elambda[ij])/is; - if (DEBUG_MODE_ENABLED && m_debugCalc) { + if (m_debugCalc) { writelogf(" ij = %d, elambda = %g, elambda1 = %g\n", ij, elambda[ij], elambda1[ij]); } @@ -4324,11 +4311,7 @@ doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const int HMWSoln::debugPrinting() { -#ifdef DEBUG_MODE return m_debugCalc; -#else - return 0; -#endif } } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index b9944e4a2..4b87c7a22 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -453,15 +453,13 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) { NeutralMolecMoleFractions_[k] = 0.0; } - if (DEBUG_MODE_ENABLED) { - sum = -1.0; - for (size_t k = 0; k < m_kk; k++) { - sum += moleFractions_[k]; - } - if (fabs(sum) > 1.0E-11) { - throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", - "molefracts don't sum to one: {}", sum); - } + sum = -1.0; + for (size_t k = 0; k < m_kk; k++) { + sum += moleFractions_[k]; + } + if (fabs(sum) > 1.0E-11) { + throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", + "molefracts don't sum to one: {}", sum); } switch (ionSolnType_) { @@ -494,28 +492,26 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij; } - if (DEBUG_MODE_ENABLED) { + for (size_t k = 0; k < m_kk; k++) { + moleFractionsTmp_[k] = moleFractions_[k]; + } + for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { for (size_t k = 0; k < m_kk; k++) { - moleFractionsTmp_[k] = moleFractions_[k]; + fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; + moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut]; } - for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { - for (size_t k = 0; k < m_kk; k++) { - fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; - moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut]; - } - } - for (size_t k = 0; k < m_kk; k++) { - if (fabs(moleFractionsTmp_[k]) > 1.0E-13) { - // Check to see if we have in fact found the inverse. - if (anionList_[0] != k) { + } + for (size_t k = 0; k < m_kk; k++) { + if (fabs(moleFractionsTmp_[k]) > 1.0E-13) { + // Check to see if we have in fact found the inverse. + if (anionList_[0] != k) { + throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", + "neutral molecule calc error"); + } else { + // For the single anion case, we will allow some slippage + if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", - "neutral molecule calc error"); - } else { - // For the single anion case, we will allow some slippage - if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { - throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", - "neutral molecule calc error - anion"); - } + "neutral molecule calc error - anion"); } } } diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index fd2ae5063..1477852eb 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -230,13 +230,11 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const // At this point we can check against the mole fraction vector of the // underlying LatticePhase objects and get the same answer. - if (DEBUG_MODE_ENABLED) { - m_lattice[n]->getMoleFractions(&m_x[strt]); - for (size_t k = 0; k < nsp; k++) { - if (fabs((x + strt)[k] - m_x[strt+k]) > 1.0E-14) { - throw CanteraError("LatticeSolidPhase::getMoleFractions()", - "internal error"); - } + m_lattice[n]->getMoleFractions(&m_x[strt]); + for (size_t k = 0; k < nsp; k++) { + if (fabs((x + strt)[k] - m_x[strt+k]) > 1.0E-14) { + throw CanteraError("LatticeSolidPhase::getMoleFractions()", + "internal error"); } } strt += nsp; diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp index 3378668e8..96bf99cd7 100644 --- a/src/thermo/PDSS_HKFT.cpp +++ b/src/thermo/PDSS_HKFT.cpp @@ -205,13 +205,11 @@ doublereal PDSS_HKFT::enthalpy_mole() const return h; } -#ifdef DEBUG_MODE doublereal PDSS_HKFT::enthalpy_mole2() const { double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * 1.0E3 * 4.184; return deltaH() + enthTRPR; } -#endif doublereal PDSS_HKFT::intEnergy_mole() const { @@ -634,7 +632,6 @@ void PDSS_HKFT::constructPDSSFile(VPStandardStateTP* tp, size_t spindex, constructPDSSXML(tp, spindex, *s, *fxml_phase, true); } -#ifdef DEBUG_MODE doublereal PDSS_HKFT::deltaH() const { doublereal pbar = m_pres * 1.0E-5; @@ -683,7 +680,6 @@ doublereal PDSS_HKFT::deltaH() const // Convert to Joules / kmol return deltaH_calgmol * 1.0E3 * 4.184; } -#endif doublereal PDSS_HKFT::deltaG() const { diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp index d0ac1de98..783c7e366 100644 --- a/src/thermo/RedlichKisterVPSSTP.cpp +++ b/src/thermo/RedlichKisterVPSSTP.cpp @@ -571,7 +571,6 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) resizeNumInteractions(numBinaryInteractions_); } -#ifdef DEBUG_MODE void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut) { double XA; @@ -612,5 +611,5 @@ void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut) VintOut = Volts; voltsOut = Volts + termp; } -#endif + } diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index e402f313a..56b03bec3 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -169,7 +169,6 @@ doublereal RedlichKwongMFTP::cv_mole() const doublereal RedlichKwongMFTP::pressure() const { -#ifdef DEBUG_MODE _updateReferenceStateThermo(); // Get a copy of the private variables stored in the State object @@ -179,7 +178,7 @@ doublereal RedlichKwongMFTP::pressure() const if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) { throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe"); } -#endif + return m_Pcurrent; } diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp index 4da6ef4d9..c62d8518b 100644 --- a/src/transport/GasTransport.cpp +++ b/src/transport/GasTransport.cpp @@ -440,23 +440,15 @@ void GasTransport::setupMM() } // initialize the collision integral calculator for the desired T* range - if (DEBUG_MODE_ENABLED && m_log_level) { - writelog("*** collision_integrals ***\n"); - } + debuglog("*** collision_integrals ***\n", m_log_level); MMCollisionInt integrals; integrals.init(tstar_min, tstar_max, m_log_level); fitCollisionIntegrals(integrals); - if (DEBUG_MODE_ENABLED && m_log_level) { - writelog("*** end of collision_integrals ***\n"); - } + debuglog("*** end of collision_integrals ***\n", m_log_level); // make polynomial fits - if (DEBUG_MODE_ENABLED && m_log_level) { - writelog("*** property fits ***\n"); - } + debuglog("*** property fits ***\n", m_log_level); fitProperties(integrals); - if (DEBUG_MODE_ENABLED && m_log_level) { - writelog("*** end of property fits ***\n"); - } + debuglog("*** end of property fits ***\n", m_log_level); } void GasTransport::getTransportData() @@ -514,7 +506,7 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals) // Chemkin fits to sixth order polynomials int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE); - if (DEBUG_MODE_ENABLED && m_log_level) { + if (m_log_level) { writelog("tstar_fits\n" "fits to A*, B*, and C* vs. log(T*).\n" "These are done only for the required dstar(j,k) values.\n\n"); @@ -578,13 +570,13 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) vector_fp c(degree + 1), c2(degree + 1); // fit the pure-species viscosity and thermal conductivity for each species - if (DEBUG_MODE_ENABLED && m_log_level && m_log_level < 2) { + if (m_log_level && m_log_level < 2) { writelog("*** polynomial coefficients not printed (log_level < 2) ***\n"); } double sqrt_T, visc, err, relerr, mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0; - if (DEBUG_MODE_ENABLED && m_log_level) { + if (m_log_level) { writelog("Polynomial fits for viscosity:\n"); if (m_mode == CK_Mode) { writelog("log(viscosity) fit to cubic polynomial in log(T)\n"); @@ -695,11 +687,11 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) m_visccoeffs.push_back(c); m_condcoeffs.push_back(c2); - if (DEBUG_MODE_ENABLED && m_log_level >= 2) { + if (m_log_level >= 2) { writelog(m_thermo->speciesName(k) + ": [" + vec2str(c) + "]\n"); } } - if (DEBUG_MODE_ENABLED && m_log_level) { + if (m_log_level) { writelogf("Maximum viscosity absolute error: %12.6g\n", mxerr); writelogf("Maximum viscosity relative error: %12.6g\n", mxrelerr); writelog("\nPolynomial fits for conductivity:\n"); @@ -775,13 +767,13 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) mxrelerr = std::max(mxrelerr, fabs(relerr)); } m_diffcoeffs.push_back(c); - if (DEBUG_MODE_ENABLED && m_log_level >= 2) { + if (m_log_level >= 2) { writelog(m_thermo->speciesName(k) + "__" + m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n"); } } } - if (DEBUG_MODE_ENABLED && m_log_level) { + if (m_log_level) { writelogf("Maximum binary diffusion coefficient absolute error:" " %12.6g\n", mxerr); writelogf("Maximum binary diffusion coefficient relative error:" diff --git a/src/transport/MMCollisionInt.cpp b/src/transport/MMCollisionInt.cpp index 15ff98e51..b3f275c2d 100644 --- a/src/transport/MMCollisionInt.cpp +++ b/src/transport/MMCollisionInt.cpp @@ -216,9 +216,7 @@ double MMCollisionInt::cstar_table[39*8] = { void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) { m_loglevel = log_level; - if (DEBUG_MODE_ENABLED && m_loglevel > 0) { - writelog("Collision Integral Polynomial Fits\n"); - } + debuglog("Collision Integral Polynomial Fits\n", m_loglevel > 0); m_nmin = -1; m_nmax = -1; @@ -234,14 +232,14 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) m_nmin = 0; m_nmax = 36; } - if (DEBUG_MODE_ENABLED && m_loglevel > 0) { + if (m_loglevel > 0) { writelogf("T*_min = %g\n", tstar[m_nmin + 1]); writelogf("T*_max = %g\n", tstar[m_nmax + 1]); } m_logTemp.resize(37); doublereal rmserr, e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0; - if (DEBUG_MODE_ENABLED && m_loglevel > 0) { + if (m_loglevel > 0) { writelog("Collision integral fits at each tabulated T* vs. delta*.\n" "These polynomial fits are used to interpolate between " "columns (delta*)\n in the Monchick and Mason tables." @@ -256,7 +254,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) vector_fp c(DeltaDegree+1); rmserr = fitDelta(0, i, DeltaDegree, c.data()); - if (DEBUG_MODE_ENABLED && log_level > 3) { + if (log_level > 3) { writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]); writelog("omega22 = [" + vec2str(c) + "]\n"); } @@ -265,26 +263,26 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) rmserr = fitDelta(1, i, DeltaDegree, c.data()); m_apoly.push_back(c); - if (DEBUG_MODE_ENABLED && log_level > 3) { + if (log_level > 3) { writelog("A* = [" + vec2str(c) + "]\n"); } ea = std::max(ea, rmserr); rmserr = fitDelta(2, i, DeltaDegree, c.data()); m_bpoly.push_back(c); - if (DEBUG_MODE_ENABLED && log_level > 3) { + if (log_level > 3) { writelog("B* = [" + vec2str(c) + "]\n"); } eb = std::max(eb, rmserr); rmserr = fitDelta(3, i, DeltaDegree, c.data()); m_cpoly.push_back(c); - if (DEBUG_MODE_ENABLED && log_level > 3) { + if (log_level > 3) { writelog("C* = [" + vec2str(c) + "]\n"); } ec = std::max(ec, rmserr); - if (DEBUG_MODE_ENABLED && log_level > 0) { + if (log_level > 0) { writelogf("max RMS errors in fits vs. delta*:\n" " omega_22 = %12.6g \n" " A* = %12.6g \n" @@ -433,7 +431,7 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar, } w[0]= -1.0; rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, o22); - if (DEBUG_MODE_ENABLED && m_loglevel > 0 && rmserr > 0.01) { + if (m_loglevel > 0 && rmserr > 0.01) { writelogf("Warning: RMS error = %12.6g in omega_22 fit" "with delta* = %12.6g\n", rmserr, deltastar); } @@ -477,7 +475,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar, } w[0]= -1.0; rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, c); - if (DEBUG_MODE_ENABLED && m_loglevel > 2) { + if (m_loglevel > 2) { writelogf("\nT* fit at delta* = %.6g\n", deltastar); writelog("astar = [" + vec2str(vector_fp(a, a+degree+1))+ "]\n");