diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h index 7569202fb..0459e1906 100644 --- a/include/cantera/equil/MultiPhase.h +++ b/include/cantera/equil/MultiPhase.h @@ -757,14 +757,12 @@ size_t ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, std::vector& orderVectorSpecies, std::vector& orderVectorElements); -#ifdef DEBUG_MODE //! External int that is used to turn on debug printing for the //! BasisOptimze program. /*! * Set this to 1 if you want debug printing from BasisOptimize. */ extern int BasisOptimize_print_lvl; -#endif } #endif diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 1f59513eb..1a57848cb 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -1364,9 +1364,9 @@ private: //! Print out and check the elemental abundance vector void prneav() const; +#endif void checkDelta1(double* const ds, double* const delTPhMoles, int kspec); -#endif //! Estimate equilibrium compositions /*! diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index bd162d879..694417994 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -9,7 +9,6 @@ using namespace Cantera; using namespace std; -#ifdef DEBUG_MODE namespace Cantera { int BasisOptimize_print_lvl = 0; @@ -26,7 +25,6 @@ int BasisOptimize_print_lvl = 0; * 2 left aligned */ static void print_stringTrunc(const char* str, int space, int alignment); -#endif //! Finds the location of the maximum component in a vector *x* /*! @@ -78,9 +76,7 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, } } -#ifdef DEBUG_MODE - double molSave = 0.0; - if (BasisOptimize_print_lvl >= 1) { + if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) { writelog(" "); for (i=0; i<77; i++) { writelog("-"); @@ -118,7 +114,6 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, writelog(" --- \n"); } } -#endif /* * Calculate the maximum value of the number of components possible @@ -148,13 +143,14 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, formRxnMatrix.resize(nspecies*ne, 0.0); } -#ifdef DEBUG_MODE /* * For debugging purposes keep an unmodified copy of the array. */ - vector_fp molNumBase(molNum); -#endif - + vector_fp molNumBase; + if (DEBUG_MODE_ENABLED) { + molNumBase = molNum; + } + double molSave = 0.0; size_t jr = npos; /* @@ -260,8 +256,7 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, /* **** REARRANGE THE DATA ****************** */ /* ****************************************** */ if (jr != k) { -#ifdef DEBUG_MODE - if (BasisOptimize_print_lvl >= 1) { + if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) { kk = orderVectorSpecies[k]; sname = mphase->speciesName(kk); writelogf(" --- %-12.12s", sname.c_str()); @@ -270,7 +265,6 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str()); writelogf("(%9.2g) as component %3d\n", molNum[jj], jr); } -#endif std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]); } @@ -346,8 +340,7 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne, &ipiv[0], &formRxnMatrix[0], ne, info); -#ifdef DEBUG_MODE - if (Cantera::BasisOptimize_print_lvl >= 1) { + if (DEBUG_MODE_ENABLED && Cantera::BasisOptimize_print_lvl >= 1) { writelog(" ---\n"); writelogf(" --- Number of Components = %d\n", nComponents); writelog(" --- Formula Matrix:\n"); @@ -390,14 +383,10 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, } writelog("\n"); } -#endif return nComponents; } /* basopt() ************************************************************/ - - -#ifdef DEBUG_MODE static void print_stringTrunc(const char* str, int space, int alignment) /*********************************************************************** @@ -442,7 +431,6 @@ static void print_stringTrunc(const char* str, int space, int alignment) } } } -#endif /* * Finds the location of the maximum component in a double vector @@ -482,8 +470,7 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda size_t nspecies = mphase->nSpecies(); double test = -1.0E10; -#ifdef DEBUG_MODE - if (BasisOptimize_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { writelog(" "); for (i=0; i<77; i++) { writelog("-"); @@ -493,7 +480,6 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda writelog("check stoich. coefficient matrix\n"); writelog(" --- and to rearrange the element ordering once\n"); } -#endif /* * Perhaps, initialize the element ordering @@ -578,11 +564,9 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda // When we are here, there is an error usually. // We haven't found the number of elements necessary. // This is signalled by returning jr != nComponents. -#ifdef DEBUG_MODE - if (BasisOptimize_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { writelogf("Error exit: returning with nComponents = %d\n", jr); } -#endif return jr; } @@ -659,8 +643,7 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda /* **** REARRANGE THE DATA ****************** */ /* ****************************************** */ if (jr != k) { -#ifdef DEBUG_MODE - if (BasisOptimize_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) { kk = orderVectorElements[k]; ename = mphase->elementName(kk); writelog(" --- "); @@ -671,7 +654,6 @@ size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbunda writelogf("%-2.2s", ename.c_str()); writelogf(" as element %3d\n", jr); } -#endif std::swap(orderVectorElements[jr], orderVectorElements[k]); } diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index c9c61e43f..bf3717dfc 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -219,8 +219,7 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, */ update(s); -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog("setInitialMoles: Estimated Mole Fractions\n"); writelogf(" Temperature = %g\n", s.temperature()); writelogf(" Pressure = %g\n", s.pressure()); @@ -236,7 +235,6 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]); } } -#endif iok = 0; } catch (CanteraError& err) { @@ -292,8 +290,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); -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { for (size_t m = 0; m < m_nComponents; m++) { int isp = m_component[m]; string nnn = s.speciesName(isp); @@ -310,7 +307,6 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, n, nnn.c_str(), xMF_est[n], mu_RT[n]); } } -#endif DenseMatrix aa(m_nComponents, m_nComponents, 0.0); for (size_t m = 0; m < m_nComponents; m++) { for (size_t n = 0; n < m_nComponents; n++) { @@ -330,8 +326,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, lambda_RT[m_orderVectorElements[m]] = 0.0; } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog(" id CompSpecies ChemPot EstChemPot Diff\n"); for (size_t m = 0; m < m_nComponents; m++) { int isp = m_component[m]; @@ -350,7 +345,6 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, writelogf(" %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]); } } -#endif return info; } @@ -388,10 +382,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, "Input ThermoPhase is incompatible with initialization"); } -#ifdef DEBUG_MODE - int n; - const vector& eNames = s.elementNames(); -#endif initialize(s); update(s); switch (XY) { @@ -745,18 +735,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // Compute the Jacobian matrix equilJacobian(s, x, elMolesGoal, jac, xval, yval); -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf("Jacobian matrix %d:\n", iter); for (m = 0; m <= m_mm; m++) { writelog(" [ "); - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { writelogf("%10.5g ", jac(m,n)); } writelog(" ]"); char xName[32]; if (m < m_mm) { - string nnn = eNames[m]; + string nnn = s.elementName(m); sprintf(xName, "x_%-10s", nnn.c_str()); } else { sprintf(xName, "x_XX"); @@ -771,7 +760,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, writelogf(" = - (%10.5g)\n", res_trial[m]); } } -#endif copy(x.begin(), x.end(), oldx.begin()); oldf = f; @@ -817,12 +805,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, } } } - if (fctr != 1.0) { -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { - writelogf("WARNING Soln Damping because of bounds: %g\n", fctr); - } -#endif + if (fctr != 1.0 && DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + writelogf("WARNING Soln Damping because of bounds: %g\n", fctr); } // multiply the step by the scaling factor @@ -883,15 +867,13 @@ 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]; } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && 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++) { writelogf(" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]); } } -#endif return 1; } @@ -926,28 +908,24 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, } } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !m_doResPerturb) { writelog("Residual: ElFracGoal ElFracCurrent Resid\n"); for (int n = 0; n < m_mm; n++) { writelogf(" % -14.7E % -14.7E % -10.5E\n", elmFracGoal[n], elmFrac[n], resid[n]); } } -#endif xx = m_p1->value(s); yy = m_p2->value(s); resid[m_mm] = xx/xval - 1.0; resid[m_skip] = yy/yval - 1.0; -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { + if (DEBUG_MODE_ENABLED && 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]); } -#endif } void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, @@ -1121,10 +1099,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } - -#ifdef DEBUG_MODE - const vector& eNames = s.elementNames(); - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog("estimateEP_Brinkley::\n\n"); double temp = s.temperature(); double pres = s.pressure(); @@ -1145,7 +1120,6 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, writelogf("%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]); } } -#endif for (m = 0; m < m_mm; m++) { if (m != m_eloc) { if (elMoles[m] <= options.absElemTol) { @@ -1169,11 +1143,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, /* * Calculate the mole numbers of species */ -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf("START ITERATION %d:\n", iter); } -#endif /* * Calculate the mole numbers of species and elements. */ @@ -1184,9 +1156,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; } - -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); for (k = 0; k < m_kk; k++) { string nnn = s.speciesName(k); @@ -1195,11 +1165,10 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, writelogf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); writelogf("(iter %d) element moles bal: Goal Calculated\n", iter); for (m = 0; m < m_mm; m++) { - string nnn = eNames[m]; + string nnn = s.elementName(m); writelogf(" %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]); } } -#endif double nCutoff; @@ -1283,9 +1252,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } nCutoff = 1.0E-9 * n_t_calc; -#ifdef DEBUG_MODE - writelog(" Lump Sum Elements Calculation: \n", ChemEquil_print_lvl); -#endif + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + writelog(" Lump Sum Elements Calculation: \n"); + } for (m = 0; m < m_mm; m++) { size_t kMSp = npos; size_t kMSp2 = npos; @@ -1309,12 +1278,10 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { - string nnn = eNames[m]; + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { + string nnn = s.elementName(m); writelogf(" %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2); } -#endif } /* @@ -1367,11 +1334,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, for (m = 0; m < m_mm; m++) { if (a1(m,m) < 1.0E-50) { -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m); } -#endif for (n = 0; n < m_mm; n++) { a1(m,n) = 0.0; } @@ -1389,8 +1354,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, resid[m_mm] = n_t - n_t_calc; -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog("Matrix:\n"); for (m = 0; m <= m_mm; m++) { writelog(" ["); @@ -1400,15 +1364,12 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, writelogf("] = %10.5g\n", resid[m]); } } -#endif tmp = resid[m_mm] /(n_t + 1.0E-15); sum += tmp * tmp; -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf("(it %d) Convergence = %g\n", iter, sum); } -#endif /* * Insist on 20x accuracy compared to the top routine. * There are instances, for ill-conditioned or @@ -1430,11 +1391,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, tmp += fabs(a1(m,n)); } if (m < m_mm && tmp < 1.0E-30) { -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing row %d\n", m); } -#endif for (n = 0; n <= m_mm; n++) { if (n != m) { a1(m,n) = 0.0; @@ -1449,8 +1408,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, resid[m] *= tmp; } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelog("Row Summed Matrix:\n"); for (m = 0; m <= m_mm; m++) { writelog(" ["); @@ -1460,7 +1418,6 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, writelogf("] = %10.5g\n", resid[m]); } } -#endif /* * Next Step: We have row-summed the equations. * However, there are some degenerate cases where two @@ -1498,15 +1455,13 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } if (sameAsRow != npos || lumpSum[m]) { -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { if (lumpSum[m]) { writelogf("Lump summing row %d, due to rank deficiency analysis\n", m); } else if (sameAsRow != npos) { writelogf("Identified that rows %d and %d are the same\n", m, sameAsRow); } } -#endif modifiedMatrix = true; for (n = 0; n < m_mm; n++) { if (n != m) { @@ -1532,9 +1487,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, solve(a1, DATA_PTR(resid)); } catch (CanteraError& err) { err.save(); -#ifdef DEBUG_MODE - writelog("Matrix is SINGULAR.ERROR\n", ChemEquil_print_lvl); -#endif + if (DEBUG_MODE_ENABLED) { + writelog("Matrix is SINGULAR.ERROR\n", ChemEquil_print_lvl); + } s.restoreState(state); throw CanteraError("equilibrate:estimateEP_Brinkley()", "Jacobian is singular. \nTry adding more species, " @@ -1562,13 +1517,11 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { if (beta != 1.0) { writelogf("(it %d) Beta = %g\n", iter, beta); } } -#endif } /* * Update the solution vector @@ -1578,19 +1531,16 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } n_t *= exp(beta * resid[m_mm]); -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); for (m = 0; m < m_mm; m++) { - string eee = eNames[m]; + string eee = s.elementName(m); writelogf(" %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]); } writelogf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); } -#endif } -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { double temp = s.temperature(); double pres = s.pressure(); @@ -1602,7 +1552,6 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, temp, pres); } } -#endif return retn; } @@ -1618,12 +1567,11 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) s.getMoleFractions(DATA_PTR(m_molefractions)); size_t k; -#ifdef DEBUG_MODE int maxPosEloc = -1; int maxNegEloc = -1; double maxPosVal = -1.0; double maxNegVal = -1.0; - if (ChemEquil_print_lvl > 0) { + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { for (k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0) { if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { @@ -1639,7 +1587,6 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) } } } -#endif double sumPos = 0.0; double sumNeg = 0.0; @@ -1658,16 +1605,12 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) return; } double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos; -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { - if (factor < 0.9999999999) { - string nnn = s.speciesName(maxPosEloc); - writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", - nnn.c_str(), - m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); - } + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) { + string nnn = s.speciesName(maxPosEloc); + writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); } -#endif for (k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0) { m_molefractions[k] *= factor; @@ -1675,16 +1618,12 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) } } else { double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg; -#ifdef DEBUG_MODE - if (ChemEquil_print_lvl > 0) { - if (factor < 0.9999999999) { - string nnn = s.speciesName(maxNegEloc); - writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", - nnn.c_str(), - m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); - } + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) { + string nnn = s.speciesName(maxNegEloc); + writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); } -#endif for (k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) < 0.0) { m_molefractions[k] *= factor; diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 5b2da4c38..312608906 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -202,16 +202,10 @@ void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies, const size_t numElem, const char* const phaseName, const double molesInert) { -#ifdef DEBUG_MODE - if (nspecies <= 0) { + if (DEBUG_MODE_ENABLED && nspecies <= 0) { plogf("nspecies Error\n"); exit(EXIT_FAILURE); } - if (phaseNum < 0) { - plogf("phaseNum should be greater than 0\n"); - exit(EXIT_FAILURE); - } -#endif setTotalMolesInert(molesInert); m_phi = 0.0; m_phiVarIndex = npos; @@ -483,41 +477,31 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc, v_totalMoles = m_totalMolesInert; if (molesSpeciesVCS == 0) { -#ifdef DEBUG_MODE - if (m_owningSolverObject == 0) { + if (DEBUG_MODE_ENABLED && m_owningSolverObject == 0) { printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); exit(EXIT_FAILURE); } -#endif if (stateCalc == VCS_STATECALC_OLD) { molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old); } else if (stateCalc == VCS_STATECALC_NEW) { molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new); - } -#ifdef DEBUG_MODE - else { + } else if (DEBUG_MODE_ENABLED) { printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); exit(EXIT_FAILURE); } -#endif - } -#ifdef DEBUG_MODE - else { - if (m_owningSolverObject) { - if (stateCalc == VCS_STATECALC_OLD) { - if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old)) { - printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); - exit(EXIT_FAILURE); - } - } else if (stateCalc == VCS_STATECALC_NEW) { - if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new)) { - printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); - exit(EXIT_FAILURE); - } + } else if (DEBUG_MODE_ENABLED && m_owningSolverObject) { + if (stateCalc == VCS_STATECALC_OLD) { + if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old)) { + printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); + exit(EXIT_FAILURE); + } + } else if (stateCalc == VCS_STATECALC_NEW) { + if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new)) { + printf("vcs_VolPhase::setMolesFromVCS shouldn't be here\n"); + exit(EXIT_FAILURE); } } } -#endif for (size_t k = 0; k < m_numSpecies; k++) { if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { @@ -937,14 +921,12 @@ void vcs_VolPhase::setTotalMoles(const double totalMols) v_totalMoles = totalMols; if (m_totalMolesInert > 0.0) { m_existence = VCS_PHASE_EXIST_ALWAYS; -#ifdef DEBUG_MODE - if (totalMols < m_totalMolesInert) { + if (DEBUG_MODE_ENABLED && totalMols < m_totalMolesInert) { printf(" vcs_VolPhase::setTotalMoles:: ERROR totalMoles " "less than inert moles: %g %g\n", totalMols, m_totalMolesInert); exit(EXIT_FAILURE); } -#endif } else { if (m_singleSpecies && (m_phiVarIndex == 0)) { m_existence = VCS_PHASE_EXIST_ALWAYS; @@ -1046,30 +1028,24 @@ void vcs_VolPhase::setExistence(const int existence) { if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) { if (v_totalMoles != 0.0) { -#ifdef DEBUG_MODE - plogf("vcs_VolPhase::setExistence setting false existence for phase with moles"); - plogendl(); - exit(EXIT_FAILURE); -#else - v_totalMoles = 0.0; -#endif + if (DEBUG_MODE_ENABLED) { + plogf("vcs_VolPhase::setExistence setting false existence for phase with moles"); + plogendl(); + exit(EXIT_FAILURE); + } else { + v_totalMoles = 0.0; + } } - } -#ifdef DEBUG_MODE - else { - if (m_totalMolesInert == 0.0) { - if (v_totalMoles == 0.0) { - if (!m_singleSpecies || m_phiVarIndex != 0) { - plogf("vcs_VolPhase::setExistence setting true existence for phase with no moles"); - plogendl(); - exit(EXIT_FAILURE); - } + } else if (DEBUG_MODE_ENABLED && m_totalMolesInert == 0.0) { + if (v_totalMoles == 0.0) { + if (!m_singleSpecies || m_phiVarIndex != 0) { + plogf("vcs_VolPhase::setExistence setting true existence for phase with no moles"); + plogendl(); + exit(EXIT_FAILURE); } } } -#endif -#ifdef DEBUG_MODE - if (m_singleSpecies) { + if (DEBUG_MODE_ENABLED && m_singleSpecies) { if (m_phiVarIndex == 0) { if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) { plogf("vcs_VolPhase::Trying to set existence of an electron phase to false"); @@ -1078,7 +1054,6 @@ void vcs_VolPhase::setExistence(const int existence) } } } -#endif m_existence = existence; } diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 2269fcf30..4146eb298 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -211,15 +211,12 @@ 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) { - -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 3) { + if (DEBUG_MODE_ENABLED && 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].c_str(), m_molNumSpecies_old[kspec], maxPermissible, m_elementName[i].c_str()); } -#endif m_molNumSpecies_old[kspec] = maxPermissible; changed = true; if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) { @@ -229,13 +226,11 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) } else { m_speciesStatus[kspec] = VCS_SPECIES_ACTIVEBUTZERO; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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].c_str(), m_speciesStatus[kspec]); } -#endif } } } diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 3ae62e4f6..88b31c30b 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -20,8 +20,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, double* const sm, double* const ss) { size_t ncomponents = m_numComponents; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); for (size_t i=0; i<77; i++) { plogf("-"); @@ -32,7 +31,6 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, plogf(" --- and to rearrange the element ordering once"); plogendl(); } -#endif /* * Use a temporary work array for the element numbers @@ -152,8 +150,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, /* **** REARRANGE THE DATA ****************** */ /* ****************************************** */ if (jr != k) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-2.2s", (m_elementName[k]).c_str()); plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]); @@ -161,7 +158,6 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, plogf("(%9.2g) as element %3d", m_elemAbundancesGoal[jr], jr); plogendl(); } -#endif vcs_switch_elem_pos(jr, k); std::swap(aw[jr], aw[k]); } @@ -180,15 +176,13 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) if (ipos == jpos) { return; } -#ifdef DEBUG_MODE - if (ipos > (m_numElemConstraints - 1) || + if (DEBUG_MODE_ENABLED && ipos > (m_numElemConstraints - 1) || jpos > (m_numElemConstraints - 1)) { plogf("vcs_switch_elem_pos: ifunc = 0: inappropriate args: %d %d\n", ipos, jpos); plogendl(); exit(EXIT_FAILURE); } -#endif /* * Change the element Global Index list in each vcs_VolPhase object * to reflect the switch in the element positions. diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 53db9c5f2..67df2f3ab 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -33,8 +33,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, */ vcs_setMolesLinProg(); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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); @@ -61,7 +60,6 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } plogendl(); } -#endif /* * Make sure all species have positive definite mole numbers @@ -135,8 +133,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } vcs_deltag(0, true, VCS_STATECALC_NEW); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { for (size_t kspec = 0; kspec < nspecies; ++kspec) { plogf("%s", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str()); @@ -148,7 +145,6 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, m_feSpecies_new[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]); } } -#endif /* ********************************************************** */ /* **** ESTIMATE REACTION ADJUSTMENTS *********************** */ /* ********************************************************** */ @@ -196,8 +192,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { for (size_t kspec = 0; kspec < nspecies; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("%sdirection (", pprefix); @@ -214,7 +209,6 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } } -#endif /* *********************************************************** */ /* **** KEEP COMPONENT SPECIES POSITIVE ********************** */ /* *********************************************************** */ @@ -316,8 +310,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } while (!finished); finished: ; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("%s Final Mole Numbers produced by inest:\n", pprefix); plogf("%s SPECIES MOLE_NUMBER\n", pprefix); @@ -328,7 +321,6 @@ finished: plogendl(); } } -#endif } int VCS_SOLVE::vcs_inest_TP() @@ -342,24 +334,18 @@ int VCS_SOLVE::vcs_inest_TP() */ vcs_elab(); if (vcs_elabcheck(0)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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(); } -#endif return retn; -#ifdef DEBUG_MODE - } 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); - plogendl(); - } -#endif + } else if (DEBUG_MODE_ENABLED && 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); + plogendl(); } } @@ -378,13 +364,11 @@ int VCS_SOLVE::vcs_inest_TP() /* * Go get the estimate of the solution */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("%sGo find an initial estimate for the equilibrium problem", pprefix); plogendl(); } -#endif double test = -1.0E20; vcs_inest(VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test); @@ -448,14 +432,12 @@ int VCS_SOLVE::vcs_inest_TP() } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix, vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_tPhaseMoles_old))); plogendl(); } -#endif /* * Record time diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index 7761a9ff9..e17e8675a 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -123,12 +123,10 @@ void VCS_SOLVE::vcs_nondim_TP() if (m_totalMoleScale != 1.0) { if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale); plogendl(); } -#endif for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale); @@ -170,12 +168,10 @@ void VCS_SOLVE::vcs_redim_TP(void) } if (m_totalMoleScale != 1.0) { if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale); plogendl(); } -#endif for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { m_molNumSpecies_old[i] *= m_totalMoleScale; diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 3e6f1ccd7..ffee1388b 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -51,13 +51,11 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const */ for (size_t k = 0; k < Vphase->nSpecies(); k++) { size_t kspec = Vphase->spGlobalIndexVCS(k); -#ifdef DEBUG_MODE - if (m_molNumSpecies_old[kspec] > 0.0) { + if (DEBUG_MODE_ENABLED && m_molNumSpecies_old[kspec] > 0.0) { printf("ERROR vcs_popPhasePossible we shouldn't be here %lu %g > 0.0", kspec, m_molNumSpecies_old[kspec]); exit(-1); } -#endif size_t irxn = kspec - m_numComponents; if (kspec >= m_numComponents) { bool iPopPossible = true; @@ -264,23 +262,22 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) 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(); -#ifdef DEBUG_MODE - strcpy(anote, ""); -#endif + if (DEBUG_MODE_ENABLED) { + strcpy(anote, ""); + } if (existence > 0) { - -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %18s %5d NA %11.3e\n", Vphase->PhaseName.c_str(), existence, m_tPhaseMoles_old[iph]); } -#endif } else { if (Vphase->m_singleSpecies) { /*********************************************************************** @@ -293,35 +290,31 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) doublereal deltaGRxn = m_deltaGRxn_old[irxn]; Fephase = exp(-deltaGRxn) - 1.0; if (Fephase > 0.0) { -#ifdef DEBUG_MODE - strcpy(anote," (ready to be birthed)"); -#endif + if (DEBUG_MODE_ENABLED) { + strcpy(anote," (ready to be birthed)"); + } if (Fephase > FephaseMax) { iphasePop = iph; FephaseMax = Fephase; -#ifdef DEBUG_MODE - strcpy(anote," (chosen to be birthed)"); -#endif + if (DEBUG_MODE_ENABLED) { + strcpy(anote," (chosen to be birthed)"); + } } } -#ifdef DEBUG_MODE - if (Fephase < 0.0) { + if (DEBUG_MODE_ENABLED && Fephase < 0.0) { strcpy(anote," (not stable)"); if (m_tPhaseMoles_old[iph] > 0.0) { printf("shouldn't be here\n"); exit(-1); } } -#endif -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %18s %5d %10.3g %10.3g %s\n", Vphase->PhaseName.c_str(), existence, Fephase, m_tPhaseMoles_old[iph], anote); } -#endif } else { /*********************************************************************** @@ -341,22 +334,18 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) FephaseMax = Fephase; } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %18s %5d %11.3g %11.3g\n", Vphase->PhaseName.c_str(), existence, Fephase, m_tPhaseMoles_old[iph]); } -#endif } else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %18s %5d blocked %11.3g\n", Vphase->PhaseName.c_str(), existence, m_tPhaseMoles_old[iph]); } -#endif } } } @@ -371,11 +360,9 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) * pop at a time. */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" ---------------------------------------------------------------------\n"); } -#endif return iphasePop; } @@ -607,19 +594,16 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) damp = 1.0E-2; if (doSuccessiveSubstitution) { - -#ifdef DEBUG_MODE int KP = 0; - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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 (m_debug_print_lvl == 1) { + } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph); } -#endif for (size_t k = 0; k < nsp; k++) { if (fracDelta_new[k] < 1.0E-13) { @@ -828,12 +812,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]); } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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); } -#endif if (normUpdate < 1.0E-5 * damp) { converged = true; @@ -861,18 +843,16 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) printf("not done yet\n"); exit(-1); } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" ------------------------------------------------------------" "-------------------------------------------------------------\n"); - } else if (m_debug_print_lvl == 1) { + } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { if (funcPhaseStability > 0.0) { plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability); } else { plogf(" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability); } } -#endif return funcPhaseStability; } //==================================================================================================================== diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 44918d186..30194c229 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -348,11 +348,9 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) double const* const* const fm = volPhase->getFormulaMatrix(); for (size_t eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) { size_t e = volPhase->elemGlobalIndex(eVP); -#ifdef DEBUG_MODE - if (e == npos) { + if (DEBUG_MODE_ENABLED && e == npos) { exit(EXIT_FAILURE); } -#endif FormulaMatrix[e][kT] = fm[eVP][k]; } /* @@ -490,20 +488,20 @@ void VCS_PROB::reportCSV(const std::string& reportFile) } } -#ifdef DEBUG_MODE - /* - * Check consistency: These should be equal - */ - tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart); - for (size_t k = 0; k < nSpeciesPhase; k++) { - if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) { - fprintf(FP,"ERROR: incompatibility!\n"); - fclose(FP); - plogf("ERROR: incompatibility!\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + /* + * Check consistency: These should be equal + */ + tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart); + for (size_t k = 0; k < nSpeciesPhase; k++) { + if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) { + fprintf(FP,"ERROR: incompatibility!\n"); + fclose(FP); + plogf("ERROR: incompatibility!\n"); + exit(EXIT_FAILURE); + } } } -#endif iK += nSpeciesPhase; } fclose(FP); diff --git a/src/equil/vcs_rank.cpp b/src/equil/vcs_rank.cpp index 4406fed6a..5a669dd8d 100644 --- a/src/equil/vcs_rank.cpp +++ b/src/equil/vcs_rank.cpp @@ -58,8 +58,7 @@ namespace VCSnonideal { vector ss(numSpecies); double test = -0.2512345E298; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); for(i=0; i<77; i++) plogf("-"); plogf("\n"); plogf(" --- Subroutine vcs_rank called to "); plogf("calculate the rank and independent rows /colums of the following matrix\n"); @@ -88,8 +87,6 @@ namespace VCSnonideal { plogendl(); } } -#endif - /* * Calculate the maximum value of the number of components possible * It's equal to the minimum of the number of elements and the @@ -269,8 +266,7 @@ namespace VCSnonideal { numComponents = jr; LE_CLEANUP: ; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- vcs_rank found rank %d\n", numComponents); if (m_debug_print_lvl >= 5) { if (compRes.size() == elemComp.size()) { @@ -288,7 +284,6 @@ namespace VCSnonideal { } } } -#endif if (numComponentsR != numComponents) { printf("vcs_rank ERROR: number of components are different: %d %d\n", numComponentsR, numComponents); diff --git a/src/equil/vcs_root1d.cpp b/src/equil/vcs_root1d.cpp index 6fe2025f8..334604399 100644 --- a/src/equil/vcs_root1d.cpp +++ b/src/equil/vcs_root1d.cpp @@ -20,7 +20,6 @@ namespace VCSnonideal #define TOL_CONV 1.0E-5 -#ifdef DEBUG_MODE static void print_funcEval(FILE* fp, double xval, double fval, int its) { fprintf(fp,"\n"); @@ -32,7 +31,6 @@ static void print_funcEval(FILE* fp, double xval, double fval, int its) fprintf(fp,"...............................................................\n"); fprintf(fp,"\n"); } -#endif int vcsUtil_root1d(double xmin, double xmax, size_t itmax, VCS_FUNC_PTR func, void* fptrPassthrough, @@ -44,10 +42,14 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, const char* strw = "vcs_root1d WARNING: "; bool converged = false; int err = 0; + #ifdef DEBUG_MODE char fileName[80]; - FILE* fp = 0; +#else + char* fileName; #endif + FILE* fp = 0; + double x1, x2, xnew, f1, f2, fnew, slope; size_t its = 0; int posStraddle = 0; @@ -62,19 +64,15 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, double c[9], f[3], xn1, xn2, x0 = 0.0, f0 = 0.0, root, theta, xquad; callNum++; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { sprintf(fileName, "rootfd_%d.log", callNum); fp = fopen(fileName, "w"); fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n"); fprintf(fp, "-----------------------------------------------------" "-------------------------------\n"); - } -#else - if (printLvl >= 3) { + } else if (printLvl >= 3) { plogf("WARNING: vcsUtil_root1d: printlvl >= 3, but debug mode not turned on\n"); } -#endif if (xmax <= xmin) { plogf("%sxmin and xmax are bad: %g %g\n", stre, xmin, xmax); return VCS_PUB_BAD; @@ -84,12 +82,10 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, x1 = (xmin + xmax) / 2.0; } f1 = func(x1, FuncTargVal, varID, fptrPassthrough, &err); -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { print_funcEval(fp, x1, f1, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1); } -#endif if (f1 == 0.0) { *xbest = x1; return VCS_SUCCESS; @@ -105,12 +101,10 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, x2 = x1 - (xmax - xmin) / 100.; } f2 = func(x2, FuncTargVal, varID, fptrPassthrough, &err); -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { print_funcEval(fp, x2, f2, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2); } -#endif if (FuncTargVal != 0.0) { fnorm = fabs(FuncTargVal) + 1.0E-13; @@ -155,11 +149,9 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, } else { xnew = x2 - f2 / slope; } -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xlin = %-9.4g", xnew); } -#endif /* * Do a quadratic fit -> Note this algorithm seems @@ -197,13 +189,11 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, theta = fabs(xquad - xnew) / fabs(xnew - x2); theta = std::min(1.0, theta); xnew = theta * xnew + (1.0 - theta) * xquad; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { if (theta != 1.0) { fprintf(fp, " | xquad = %-9.4g", xnew); } } -#endif } else { /* * Pick out situations where the convergence may be @@ -212,11 +202,9 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, if ((DSIGN(xnew - x2) == DSIGN(x2 - x1)) && (DSIGN(x2 - x1) == DSIGN(x1 - x0))) { xnew += xnew - x2; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xquada = %-9.4g", xnew); } -#endif } } } @@ -238,19 +226,15 @@ QUAD_BAIL: slope = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < slope) { xnew = x1 + DSIGN(xnew-x1) * slope; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | x10%% = %-9.4g", xnew); } -#endif } if (fabs(xnew - x2) < slope) { xnew = x2 + DSIGN(xnew-x2) * slope; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | x10%% = %-9.4g", xnew); } -#endif } } else { /* @@ -260,30 +244,24 @@ QUAD_BAIL: slope = 2.0 * fabs(x2 - x1); if (fabs(slope) < fabs(xnew - x2)) { xnew = x2 + DSIGN(xnew-x2) * slope; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xlimitsize = %-9.4g", xnew); } -#endif } } if (xnew > xmax) { xnew = x2 + (xmax - x2) / 2.0; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xlimitmax = %-9.4g", xnew); } -#endif } if (xnew < xmin) { xnew = x2 + (x2 - xmin) / 2.0; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xlimitmin = %-9.4g", xnew); } -#endif } if (foundStraddle) { #ifdef DEBUG_MODE @@ -322,23 +300,19 @@ QUAD_BAIL: } } } -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { if (slope != xnew) { fprintf(fp, " | xstraddle = %-9.4g", xnew); } } -#endif } fnew = func(xnew, FuncTargVal, varID, fptrPassthrough, &err); -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp,"\n"); print_funcEval(fp, xnew, fnew, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew); } -#endif if (foundStraddle) { if (posStraddle) { @@ -397,28 +371,22 @@ QUAD_BAIL: if (printLvl >= 1) { plogf("vcs_root1d success: convergence achieved\n"); } -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | vcs_root1d success in %d its, fnorm = %g\n", its, fnorm); } -#endif } else { retn = VCS_FAILED_CONVERGENCE; if (printLvl >= 1) { plogf("vcs_root1d ERROR: maximum iterations exceeded without convergence\n"); } -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, "\nvcs_root1d failure in %lu its\n", its); } -#endif } *xbest = x2; -#ifdef DEBUG_MODE - if (printLvl >= 3) { + if (DEBUG_MODE_ENABLED && printLvl >= 3) { fclose(fp); } -#endif return retn; } diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 9852603de..bff2b25cb 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -41,6 +41,8 @@ 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 @@ -54,17 +56,17 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) ************************************************************************/ for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Normal Calc"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Normal Calc"); + } size_t kspec = m_indexRxnToSpecies[irxn]; if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) { m_deltaMolNumSpecies[kspec] = 0.0; -#ifdef DEBUG_MODE - sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); + } } else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; @@ -91,38 +93,38 @@ 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; -#ifdef DEBUG_MODE - 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]); -#endif + 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]); + } } else { m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.0; -#ifdef DEBUG_MODE - sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E", - vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]); -#endif + 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 { -#ifdef DEBUG_MODE - 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]); -#endif + if (DEBUG_MODE_ENABLED) { + 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.; -#ifdef DEBUG_MODE - 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]); -#endif + 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]); + } } } } else { -#ifdef DEBUG_MODE - sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15), - m_deltaGRxn_new[irxn]); -#endif + if (DEBUG_MODE_ENABLED) { + 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; } @@ -138,15 +140,15 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) * in this mode. */ if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]); - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); - plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + 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].c_str()); + plogf(" %12.4E %12.4E %12.4E | %s\n", + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); + } } -#endif continue; } /* @@ -154,15 +156,15 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) * their values are to be decreasing anyway. */ if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) { -#ifdef DEBUG_MODE - 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].c_str()); - plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], - m_deltaGRxn_new[irxn], ANOTE); + 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].c_str()); + plogf(" %12.4E %12.4E %12.4E | %s\n", + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); + } } -#endif continue; } /* @@ -198,12 +200,10 @@ 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); -#ifdef DEBUG_MODE - if (s_old != s) { + if (DEBUG_MODE_ENABLED && s_old != s) { sprintf(ANOTE, "Normal calc: diag adjusted from %g " "to %g due to act coeff", s_old, s); } -#endif } m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s; @@ -214,18 +214,18 @@ 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) { -#ifdef DEBUG_MODE - 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()); -#endif + 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()); + } m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[j] / stoicC; } else { -#ifdef DEBUG_MODE - 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()); -#endif + 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()); + } m_deltaMolNumSpecies[kspec] = 0.0; } } @@ -233,11 +233,11 @@ 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]) { -#ifdef DEBUG_MODE - 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()); -#endif + 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()); + } m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; } @@ -297,19 +297,17 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) * mole numbers of that phases. it seems that we can suggest a zero of the species * and the code will recover. */ -#ifdef DEBUG_MODE - sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); -#endif + 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()); + } m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], m_deltaGRxn_new[irxn], ANOTE); } -#endif continue; } /* @@ -326,58 +324,50 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) iphDel = m_phaseID[k]; kSpecial = k; -#ifdef DEBUG_MODE - 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 (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].c_str()); + 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].c_str()); + plogendl(); + } } - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); - 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].c_str()); - plogendl(); - } -#endif if (k != kspec) { forceComponentCalc = 1; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Force a component recalculation \n"); plogendl(); } -#endif } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); writeline('-', 82); } -#endif return iphDel; } } } /* End of regular processing */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], m_deltaGRxn_new[irxn], ANOTE); } -#endif } /* End of loop over m_speciesUnknownType */ } /* End of loop over non-component stoichiometric formation reactions */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); writeline('-', 82); } -#endif return iphDel; } @@ -392,6 +382,8 @@ 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 /* @@ -401,9 +393,9 @@ int VCS_SOLVE::vcs_rxn_adj_cg() * this algorithm. If not, we bail out. */ for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Normal Calc"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Normal Calc"); + } size_t kspec = m_indexRxnToSpecies[irxn]; double* dnPhase_irxn = m_deltaMolNumPhase[irxn]; @@ -417,16 +409,16 @@ int VCS_SOLVE::vcs_rxn_adj_cg() * should be replaced with something more relativistic */ if (m_deltaGRxn_new[irxn] < -1.0e-4) { -#ifdef DEBUG_MODE - (void) sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]); -#endif + if (DEBUG_MODE_ENABLED) { + 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 { -#ifdef DEBUG_MODE - (void) sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]); + } m_deltaMolNumSpecies[kspec] = 0.0; } } else { @@ -441,13 +433,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() * in this mode. */ if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); - plogf(" --- "); - plogf("%-12.12s", m_speciesName[kspec].c_str()); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]); + plogf(" --- "); + plogf("%-12.12s", m_speciesName[kspec].c_str()); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + } continue; } /* @@ -455,13 +447,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() * their values are to be decreasing anyway. */ if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) { -#ifdef DEBUG_MODE - 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].c_str()); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); -#endif + 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].c_str()); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + } continue; } /* @@ -543,11 +535,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } m_molNumSpecies_old[k] = 0.0; m_tPhaseMoles_old[m_phaseID[k]] = 0.0; -#ifdef DEBUG_MODE - plogf(" --- vcs_st2 Special section to delete "); - plogf("%-12.12s", m_speciesName[k].c_str()); - plogf("\n --- Immediate return - Restart iteration\n"); -#endif + if (DEBUG_MODE_ENABLED) { + plogf(" --- vcs_st2 Special section to delete "); + plogf("%-12.12s", m_speciesName[k].c_str()); + plogf("\n --- Immediate return - Restart iteration\n"); + } /* * We need to immediately recompute the * component basis, because we just zeroed @@ -562,12 +554,12 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } } } /* End of regular processing */ -#ifdef DEBUG_MODE - plogf(" --- "); - plogf("%-12.12s", m_speciesName[kspec].c_str()); - plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); -#endif + if (DEBUG_MODE_ENABLED) { + plogf(" --- "); + plogf("%-12.12s", m_speciesName[kspec].c_str()); + plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + } } /* End of loop over non-component stoichiometric formation reactions */ /* @@ -579,13 +571,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() * property. */ -#ifdef DEBUG_MODE - plogf(" "); - for (size_t j = 0; j < 77; j++) { - plogf("-"); + if (DEBUG_MODE_ENABLED) { + plogf(" "); + for (size_t j = 0; j < 77; j++) { + plogf("-"); + } + plogf("\n"); } - plogf("\n"); -#endif return soldel; } diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 2d067d6a7..39cbb99a0 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -17,7 +17,6 @@ using namespace std; namespace VCSnonideal { -#ifdef DEBUG_MODE static void printProgress(const vector &spName, const vector &soln, const vector &ff) @@ -33,20 +32,15 @@ static void printProgress(const vector &spName, } plogf(" --- Total sum to be minimized = %g\n", sum); } -#endif int VCS_SOLVE::vcs_setMolesLinProg() { size_t ik, irxn; double test = -1.0E-10; -#ifdef DEBUG_MODE - std::string pprefix(" --- seMolesLinProg "); - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- call setInitialMoles\n"); } -#endif - // m_mu are standard state chemical potentials // Boolean on the end specifies standard chem potentials @@ -75,21 +69,17 @@ int VCS_SOLVE::vcs_setMolesLinProg() } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); } -#endif while (redo) { if (!vcs_elabcheck(0)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf("%s Mole numbers failing element abundances\n", pprefix.c_str()); - plogf("%sCall vcs_elcorr to attempt fix\n", pprefix.c_str()); + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" --- seMolesLinProg Mole numbers failing element abundances\n"); + plogf(" --- seMolesLinProg Call vcs_elcorr to attempt fix\n"); } -#endif retn = vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); if (retn >= 2) { abundancesOK = false; @@ -111,11 +101,9 @@ int VCS_SOLVE::vcs_setMolesLinProg() return retn; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("iteration %d\n", iter); } -#endif redo = false; iter++; if (iter > 15) { @@ -152,11 +140,9 @@ int VCS_SOLVE::vcs_setMolesLinProg() // with a new set of components if (!redo) { if (delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Component too small: %s\n", m_speciesName[jcomp].c_str()); } -#endif redo = true; } } @@ -187,19 +173,15 @@ int VCS_SOLVE::vcs_setMolesLinProg() } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); } -#endif } -#ifdef DEBUG_MODE - if (m_debug_print_lvl == 1) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); plogf(" --- setInitialMoles end\n"); } -#endif retn = 0; if (!abundancesOK) { retn = -1; diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index c0d600871..9f36f08ed 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -35,7 +35,6 @@ static int prnfm(void); #endif /*****************************************************************************/ -#ifdef DEBUG_MODE void VCS_SOLVE::checkDelta1(double* const dsLocal, double* const delTPhMoles, int kspec) { @@ -54,7 +53,6 @@ void VCS_SOLVE::checkDelta1(double* const dsLocal, } } } -#endif int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) { @@ -77,13 +75,15 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) std::vector phasePopPhaseIDs(0); size_t doPhaseDeleteIph = npos; -#ifdef DEBUG_MODE size_t doPhaseDeleteKspec = npos; +#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; @@ -115,9 +115,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) solveFail = false; -#if DEBUG_MODE - int ll; -#endif + int ll; // only used in DEBUG_MODE /* ****************************************************** */ /* **** Evaluate the elemental composition ****** */ /* ****************************************************** */ @@ -266,27 +264,20 @@ L_COMPONENT_CALC: /************** EVALUATE THE ELELEMT ABUNDANCE CHECK ******************/ /*************************************************************************/ if (! vcs_elabcheck(0)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Element Abundance check failed"); plogendl(); } -#endif vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); // 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) { + plogf(" --- Element Abundance check passed"); + plogendl(); } -#ifdef DEBUG_MODE - else { - if (m_debug_print_lvl >= 2) { - plogf(" --- Element Abundance check passed"); - plogendl(); - } - } -#endif iti = 0; goto L_MAINLOOP_ALL_SPECIES; @@ -380,12 +371,10 @@ L_MAINLOOP_ALL_SPECIES: soldel = vcs_popPhaseRxnStepSizes(iphasePop); if (soldel == 3) { iphasePop = npos; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich " "prevented phase %d popping\n"); } -#endif } } @@ -403,16 +392,11 @@ L_MAINLOOP_ALL_SPECIES: */ kspec = npos; iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec); - } -#ifdef DEBUG_MODE - else { - if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_RxnStepSizes not called because alternative" - "phase creation delta was used instead\n"); - } + } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" --- vcs_RxnStepSizes not called because alternative" + "phase creation delta was used instead\n"); } doPhaseDeleteKspec = npos; -#endif lec = false; doPhaseDeleteIph = npos; /* @@ -452,13 +436,10 @@ L_MAINLOOP_ALL_SPECIES: * */ if (iphaseDelete != npos) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Main Loop Treatment -> Circumvented due to Phase Deletion "); plogendl(); } -#endif - for (size_t k = 0; k < m_numSpeciesTot; k++) { m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k]; size_t iph = m_phaseID[k]; @@ -497,8 +478,7 @@ L_MAINLOOP_ALL_SPECIES: */ vcs_deltag(0, false, VCS_STATECALC_NEW); } else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Main Loop Treatment of each non-component species "); if (iti == 0) { plogf("- Full Calculation:\n"); @@ -508,23 +488,22 @@ L_MAINLOOP_ALL_SPECIES: plogf(" --- Species IC "); plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n"); } -#endif for (size_t irxn = 0; irxn < m_numRxnRdc; irxn++) { kspec = m_indexRxnToSpecies[irxn]; double* sc_irxn = m_stoichCoeffRxnMatrix[irxn]; size_t iph = m_phaseID[kspec]; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; -#ifdef DEBUG_MODE - ANOTE[0] = '\0'; -#endif + if (DEBUG_MODE_ENABLED) { + 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]; -#ifdef DEBUG_MODE - sprintf(ANOTE, "Phase pop"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Phase pop"); + } } else { dx = 0.0; m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; @@ -547,33 +526,31 @@ L_MAINLOOP_ALL_SPECIES: /********************** ZEROED OUT SPECIES **************************/ /********************************************************************/ bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 3) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) { plogf(" --- %s currently zeroed (SpStatus=%-2d):", m_speciesName[kspec].c_str(), m_speciesStatus[kspec]); plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n", irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec], m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]); } -#endif if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; resurrect = false; -#ifdef DEBUG_MODE - 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]); + 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]); + } } } -#endif } else { for (size_t j = 0; j < m_numElemConstraints; ++j) { int elType = m_elType[j]; @@ -582,11 +559,11 @@ L_MAINLOOP_ALL_SPECIES: if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Species stays zeroed even though dG " - "neg, because of %s elemAbund", - m_elementName[j].c_str()); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Species stays zeroed even though dG " + "neg, because of %s elemAbund", + m_elementName[j].c_str()); + } resurrect = false; break; } @@ -604,21 +581,17 @@ L_MAINLOOP_ALL_SPECIES: } if (phaseResurrected) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Zeroed species changed to major: "); plogf("%-12s\n", m_speciesName[kspec].c_str()); } -#endif m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; allMinorZeroedSpecies = false; } else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Zeroed species changed to minor: "); plogf("%-12s\n", m_speciesName[kspec].c_str()); } -#endif m_speciesStatus[kspec] = VCS_SPECIES_MINOR; } if (m_deltaMolNumSpecies[kspec] > 0.0) { @@ -629,9 +602,9 @@ L_MAINLOOP_ALL_SPECIES: dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec]; } m_deltaMolNumSpecies[kspec] = dx; -#ifdef DEBUG_MODE - sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]); -#endif + if (DEBUG_MODE_ENABLED) { + 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; @@ -649,17 +622,17 @@ L_MAINLOOP_ALL_SPECIES: m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; -#ifdef DEBUG_MODE - sprintf(ANOTE,"minor species not considered"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); - plogf("%-12s", m_speciesName[kspec].c_str()); - 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(); + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE,"minor species not considered"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); + plogf("%-12s", m_speciesName[kspec].c_str()); + 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(); + } } -#endif continue; } /* @@ -691,13 +664,11 @@ L_MAINLOOP_ALL_SPECIES: /***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */ /***** MOLE NUMBER */ /*******************************************************************/ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Delete minor species in multispec phase: %-12s", m_speciesName[kspec].c_str()); plogendl(); } -#endif m_deltaMolNumSpecies[kspec] = 0.0; /* * Delete species, kspec. The alternate return is for the case @@ -728,9 +699,9 @@ L_MAINLOOP_ALL_SPECIES: /********************************************************************/ /*********************** MAJOR SPECIES ******************************/ /********************************************************************/ -#ifdef DEBUG_MODE - sprintf(ANOTE, "Normal Major Calc"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "Normal Major Calc"); + } /* * Check for superconvergence of the formation reaction. Do * nothing if it is superconverged. Skip to the end of the @@ -740,17 +711,17 @@ L_MAINLOOP_ALL_SPECIES: m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; -#ifdef DEBUG_MODE - sprintf(ANOTE, "major species is converged"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); - plogf("%-12s", m_speciesName[kspec].c_str()); - 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(); + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "major species is converged"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); + plogf("%-12s", m_speciesName[kspec].c_str()); + 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(); + } } -#endif continue; } /* @@ -767,10 +738,10 @@ L_MAINLOOP_ALL_SPECIES: } else { dx = 0.0; m_deltaMolNumSpecies[kspec] = 0.0; -#ifdef DEBUG_MODE - sprintf(ANOTE, "dx set to 0, DG flipped sign due to " - "changed initial point"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "dx set to 0, DG flipped sign due to " + "changed initial point"); + } } /* * Form a tentative value of the new species moles @@ -784,10 +755,10 @@ L_MAINLOOP_ALL_SPECIES: * or we restart the entire iteration. */ if (m_molNumSpecies_new[kspec] <= 0.0) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "initial nonpos kmoles= %11.3E", - m_molNumSpecies_new[kspec]); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "initial nonpos kmoles= %11.3E", + m_molNumSpecies_new[kspec]); + } /* ************************************************* */ /* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */ /* ************************************************* */ @@ -834,11 +805,11 @@ L_MAINLOOP_ALL_SPECIES: m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; if (m_molNumSpecies_new[kspec] > 0.0) { m_deltaMolNumSpecies[kspec] = dx; -#ifdef DEBUG_MODE - sprintf(ANOTE, - "zeroing SS phase created a neg component species " - "-> reducing step size instead"); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, + "zeroing SS phase created a neg component species " + "-> reducing step size instead"); + } } else { /* * We are going to zero the single species phase. @@ -847,9 +818,9 @@ L_MAINLOOP_ALL_SPECIES: iph = m_phaseID[kspec]; Vphase = m_VolPhaseList[iph]; //Vphase->setExistence(0); -#ifdef DEBUG_MODE - sprintf(ANOTE, "zeroing out SS phase: "); -#endif + if (DEBUG_MODE_ENABLED) { + sprintf(ANOTE, "zeroing out SS phase: "); + } /* * Change the base mole numbers for the iteration. * We need to do this here, because we have decided @@ -859,16 +830,16 @@ L_MAINLOOP_ALL_SPECIES: m_molNumSpecies_new[kspec] = 0.0; doPhaseDeleteIph = iph; -#ifdef DEBUG_MODE - doPhaseDeleteKspec = kspec; - if (m_debug_print_lvl >= 2) { - if (m_speciesStatus[kspec] >= 0) { - plogf(" --- SS species changed to zeroedss: "); - plogf("%-12s", m_speciesName[kspec].c_str()); - plogendl(); + if (DEBUG_MODE_ENABLED) { + doPhaseDeleteKspec = kspec; + if (m_debug_print_lvl >= 2) { + if (m_speciesStatus[kspec] >= 0) { + plogf(" --- SS species changed to zeroedss: "); + plogf("%-12s", m_speciesName[kspec].c_str()); + plogendl(); + } } } -#endif m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; ++m_numRxnMinorZeroed; allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); @@ -926,8 +897,7 @@ L_MAINLOOP_ALL_SPECIES: * to the reaction delta that we just computed. * This should keep the amount of material constant. */ -#ifdef DEBUG_MODE - if (fabs(m_deltaMolNumSpecies[kspec] -dx) > + if (DEBUG_MODE_ENABLED && fabs(m_deltaMolNumSpecies[kspec] -dx) > 1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32)) { plogf(" ds[kspec] = %20.16g dx = %20.16g , kspec = %d\n", m_deltaMolNumSpecies[kspec], dx, kspec); @@ -935,7 +905,6 @@ L_MAINLOOP_ALL_SPECIES: plogendl(); exit(EXIT_FAILURE); } -#endif for (size_t k = 0; k < m_numComponents; ++k) { m_deltaMolNumSpecies[k] += sc_irxn[k] * dx; } @@ -949,10 +918,10 @@ L_MAINLOOP_ALL_SPECIES: } } -#ifdef DEBUG_MODE - checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), - VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1); -#endif + if (DEBUG_MODE_ENABLED) { + checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), + VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1); + } /* * Branch point for returning - */ @@ -975,20 +944,17 @@ L_MAIN_LOOP_END_NO_PRINT: ; #endif if (doPhaseDeleteIph != npos) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-12.12s Main Loop Special Case deleting phase with species: ", m_speciesName[doPhaseDeleteKspec].c_str()); plogendl(); } -#endif break; } } /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { for (size_t k = 0; k < m_numComponents; k++) { plogf(" --- "); plogf("%-12.12s", m_speciesName[k].c_str()); @@ -1001,7 +967,6 @@ L_MAIN_LOOP_END_NO_PRINT: plogf(" --- Finished Main Loop"); plogendl(); } -#endif /*************************************************************************/ /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/ @@ -1018,9 +983,7 @@ L_MAIN_LOOP_END_NO_PRINT: double xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k]; if (par < xx) { par = xx; -#ifdef DEBUG_MODE ll = k; -#endif } } else { if (m_deltaMolNumSpecies[k] < 0.0) { @@ -1038,14 +1001,12 @@ L_MAIN_LOOP_END_NO_PRINT: if (par <= 1.01 && par > 0.0) { /* Reduce the size of the step by the multiplicative factor, par */ par *= 0.99; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Reduction in step size due to component "); plogf("%s", m_speciesName[ll].c_str()); plogf(" going negative = %11.3E", par); plogendl(); } -#endif for (size_t i = 0; i < m_numSpeciesTot; ++i) { m_deltaMolNumSpecies[i] *= par; } @@ -1055,10 +1016,10 @@ L_MAIN_LOOP_END_NO_PRINT: } else { par = 1.0; } -#ifdef DEBUG_MODE - checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), - VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot); -#endif + if (DEBUG_MODE_ENABLED) { + checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), + VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot); + } /* * Now adjust the wt[kspec]'s so that the reflect the decrease in @@ -1222,12 +1183,10 @@ L_MAIN_LOOP_END_NO_PRINT: vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_tPhaseMoles_new))); plogendl(); -#ifdef DEBUG_MODE - if (m_VCount->Its > 550) { + if (DEBUG_MODE_ENABLED && m_VCount->Its > 550) { plogf(" --- Troublesome solve"); plogendl(); } -#endif } /*************************************************************************/ @@ -1260,13 +1219,11 @@ L_MAIN_LOOP_END_NO_PRINT: */ ++(m_VCount->Its); ++it1; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Increment counter increased, step is accepted: %4d", m_VCount->Its); plogendl(); } -#endif /*************************************************************************/ /******************* HANDLE DELETION OF MULTISPECIES PHASES **************/ /*************************************************************************/ @@ -1284,12 +1241,10 @@ L_MAIN_LOOP_END_NO_PRINT: m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) { soldel = 1; if (soldel) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 1) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) { plogf(" --- Setting microscopic phase %d to zero", iph); plogendl(); } -#endif justDeletedMultiPhase = true; vcs_delete_multiphase(iph); } @@ -1308,9 +1263,9 @@ L_MAIN_LOOP_END_NO_PRINT: VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); if (retn != VCS_SUCCESS) { -#ifdef DEBUG_MODE + if (DEBUG_MODE_ENABLED) { plogf(" --- BASOPT returned with an error condition\n"); -#endif + } exit(EXIT_FAILURE); } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); @@ -1324,33 +1279,24 @@ L_MAIN_LOOP_END_NO_PRINT: /*************************************************************************/ /***************** CHECK FOR ELEMENT ABUNDANCE****************************/ /*************************************************************************/ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Normal element abundance check"); } -#endif vcs_elab(); if (! vcs_elabcheck(0)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" - failed -> redoing element abundances."); plogendl(); } -#endif vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); 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) { + plogf(" - passed"); + plogendl(); } -#ifdef DEBUG_MODE - else { - if (m_debug_print_lvl >= 2) { - plogf(" - passed"); - plogendl(); - } - } -#endif /*************************************************************************/ /***************** CHECK FOR OPTIMUM BASIS *******************************/ /*************************************************************************/ @@ -1401,15 +1347,13 @@ L_MAIN_LOOP_END_NO_PRINT: } if (doSwap) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because %s", m_speciesName[l].c_str()); plogf(" is better than comp %s", m_speciesName[j].c_str()); plogf(" and share nonzero stoic: %-9.1f", m_stoichCoeffRxnMatrix[i][j]); plogendl(); } -#endif goto L_COMPONENT_CALC; } } else { @@ -1420,15 +1364,13 @@ L_MAIN_LOOP_END_NO_PRINT: if (m_molNumSpecies_old[j] == 0.0) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { if (dg[i] < 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because %s", m_speciesName[l].c_str()); plogf(" has dg < 0.0 and comp %s has zero mole num", m_speciesName[j].c_str()); plogf(" and share nonzero stoic: %-9.1f", m_stoichCoeffRxnMatrix[i][j]); plogendl(); } -#endif goto L_COMPONENT_CALC; } } @@ -1466,8 +1408,7 @@ L_MAIN_LOOP_END_NO_PRINT: } if (doSwap) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because "); plogf("%s", m_speciesName[l].c_str()); plogf(" is better than comp "); @@ -1476,7 +1417,6 @@ L_MAIN_LOOP_END_NO_PRINT: m_stoichCoeffRxnMatrix[i][j]); plogendl(); } -#endif goto L_COMPONENT_CALC; } } @@ -1485,8 +1425,7 @@ L_MAIN_LOOP_END_NO_PRINT: if (m_molNumSpecies_old[j] == 0.0) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { if (dg[i] < 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because %s", m_speciesName[l].c_str()); plogf(" has dg < 0.0 and comp %s has zero mole num", m_speciesName[j].c_str()); @@ -1494,7 +1433,6 @@ L_MAIN_LOOP_END_NO_PRINT: m_stoichCoeffRxnMatrix[i][j]); plogendl(); } -#endif goto L_COMPONENT_CALC; } } @@ -1505,12 +1443,10 @@ L_MAIN_LOOP_END_NO_PRINT: } } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Check for an optimum basis passed"); plogendl(); } -#endif /*************************************************************************/ /********************** RE-EVALUATE MAJOR-MINOR VECTOR IF NECESSARY ******/ /*************************************************************************/ @@ -1519,29 +1455,24 @@ L_MAIN_LOOP_END_NO_PRINT: * Go right to the check equilibrium section */ if (iti == 0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Reevaluate major-minor status of noncomponents:\n"); } -#endif m_numRxnMinorZeroed = 0; for (size_t irxn = 0; irxn < m_numRxnRdc; irxn++) { kspec = m_indexRxnToSpecies[irxn]; int speciesType = vcs_species_type(kspec); if (speciesType < VCS_SPECIES_MINOR) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { plogf(" --- major/minor species is now zeroed out: %s\n", m_speciesName[kspec].c_str()); } } -#endif ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MINOR) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) { plogf(" --- Noncomponent turned from major to minor: "); @@ -1554,12 +1485,10 @@ L_MAIN_LOOP_END_NO_PRINT: plogf("%s\n", m_speciesName[kspec].c_str()); } } -#endif ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MAJOR) { if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { plogf(" --- Noncomponent turned from minor to major: "); } else if (kspec < m_numComponents) { @@ -1569,7 +1498,6 @@ L_MAIN_LOOP_END_NO_PRINT: } plogf("%s\n", m_speciesName[kspec].c_str()); } -#endif m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; /* * For this special case, we must reevaluate thermo functions @@ -1595,11 +1523,9 @@ L_MAIN_LOOP_END_NO_PRINT: L_EQUILIB_CHECK: ; if (! allMinorZeroedSpecies) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Equilibrium check for major species: "); } -#endif for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { kspec = irxn + m_numComponents; if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) { @@ -1611,11 +1537,9 @@ L_EQUILIB_CHECK: */ goto L_RETURN_BLOCK; } else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str()); } -#endif // Convergence amongst major species has not been achieved /* * Go back and do another iteration with variable ITI @@ -1627,22 +1551,16 @@ L_EQUILIB_CHECK: } } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" MAJOR SPECIES CONVERGENCE achieved"); plogendl(); } -#endif } -#ifdef DEBUG_MODE - else { - if (m_debug_print_lvl >= 2) { - plogf(" MAJOR SPECIES CONVERGENCE achieved " - "(because there are no major species)"); - plogendl(); - } + else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" MAJOR SPECIES CONVERGENCE achieved " + "(because there are no major species)"); + plogendl(); } -#endif // Convergence amongst major species has been achieved /*************************************************************************/ @@ -1659,11 +1577,9 @@ L_EQUILIB_CHECK: vcs_deltag(1, false, VCS_STATECALC_OLD); uptodate_minors = true; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Equilibrium check for minor species: "); } -#endif for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { kspec = irxn + m_numComponents; if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) { @@ -1675,11 +1591,9 @@ L_EQUILIB_CHECK: */ goto L_RETURN_BLOCK; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str()); } -#endif /* * Set iti to zero to force a full calculation, and go back * to the main loop to do another iteration. @@ -1691,11 +1605,9 @@ L_EQUILIB_CHECK: goto L_MAINLOOP_ALL_SPECIES; } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" CONVERGENCE achieved\n"); } -#endif } /*************************************************************************/ /*********************** FINAL ELEMENTAL ABUNDANCE CHECK *****************/ @@ -1709,34 +1621,28 @@ L_EQUILIB_CHECK: /* LEC is only true when we are near the end game */ if (lec) { if (!giveUpOnElemAbund) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Check the Full Element Abundances: "); } -#endif /* * 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)) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (! vcs_elabcheck(0)) { plogf(" failed\n"); } else { plogf(" passed for NC but failed for NE: RANGE ERROR\n"); } } -#endif // delete? goto L_ELEM_ABUND_CHECK; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" passed\n"); } -#endif } /* * If we have deleted a species then we need to recheck the @@ -1811,15 +1717,13 @@ L_ELEM_ABUND_CHECK: /* * Probably an unrecoverable range error */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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"); plogf(" --- vcs_solve_tp: - Returning the calculated equilibrium condition "); plogendl(); } -#endif rangeErrorFound = 1; giveUpOnElemAbund = true; goto L_EQUILIB_CHECK; @@ -1895,12 +1799,10 @@ L_RETURN_BLOCK_B: npb = vcs_add_all_deleted(); if (npb > 0) { iti = 0; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 1) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) { plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!"); plogendl(); } -#endif goto L_MAINLOOP_ALL_SPECIES; } @@ -2080,13 +1982,11 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) size_t irxn = kspec - m_numComponents; int retn = 1; double delta = *delta_ptr; -#ifdef DEBUG_MODE - if (kspec < m_numComponents) { + if (DEBUG_MODE_ENABLED && kspec < m_numComponents) { plogf(" --- delete_species() ERROR: called for a component %d", kspec); plogendl(); exit(EXIT_FAILURE); } -#endif if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { /* * Attempt the given dx. If it doesn't work, try to see if a smaller @@ -2148,16 +2048,12 @@ 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); -#ifdef DEBUG_MODE - if (!retn) { - if (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); - plogendl(); - } + if (DEBUG_MODE_ENABLED && !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); + plogendl(); } -#endif } } return retn; @@ -2244,11 +2140,9 @@ 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]; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Add back a deleted species: %-12s\n", m_speciesName[kspec].c_str()); } -#endif /* * Set the species back to minor species status * this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[] @@ -2311,11 +2205,9 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) * set the phase existence flag to dead */ Vphase->setTotalMoles(0.0); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName.c_str()); } -#endif /* * Loop over all of the species in the phase. @@ -2332,14 +2224,12 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) int retn = delta_species(kspec, &dxTent); if (retn != 1) { successful = false; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n", iph, Vphase->PhaseName.c_str(), m_speciesName[kspec].c_str()); plogf(" --- delta attempted: %g achieved: %g " " Zeroing it manually\n", dx, dxTent); } -#endif m_molNumSpecies_old[kspec] = 0.0; m_molNumSpecies_new[kspec] = 0.0; m_deltaMolNumSpecies[kspec] = 0.0; @@ -2365,12 +2255,10 @@ 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) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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].c_str(), m_molNumSpecies_old[kcomp]); } -#endif if (m_molNumSpecies_old[kcomp] != 0.0) { for (size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { size_t irxn = kspec - m_numComponents; @@ -2404,13 +2292,11 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) } if (m_molNumSpecies_old[kcomp] != 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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].c_str(), m_molNumSpecies_old[kcomp]); plogf(" --- zeroing it \n"); } -#endif m_molNumSpecies_old[kcomp] = 0.0; } m_speciesStatus[kcomp] = VCS_SPECIES_ZEROEDMS; @@ -2436,13 +2322,11 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) ++(m_numRxnRdc); ++(m_numSpeciesRdc); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Make %s", m_speciesName[kspec].c_str()); plogf(" an active but zeroed species because its phase " "was zeroed\n"); } -#endif if (kspec != (m_numSpeciesRdc - 1)) { /* * Rearrange both the species and the non-component global data @@ -2469,11 +2353,9 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) int VCS_SOLVE::vcs_recheck_deleted() { double* xtcutoff = VCS_DATA_PTR(m_TmpPhase); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Start rechecking deleted species in multispec phases\n"); } -#endif if (m_numSpeciesRdc == m_numSpeciesTot) { return 0; } @@ -2638,27 +2520,20 @@ size_t VCS_SOLVE::vcs_add_all_deleted() double dx = m_molNumSpecies_new[kspec]; retn = delta_species(kspec, &dx); if (retn == 0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", m_speciesName[kspec].c_str(), kspec, dx); } -#endif if (dx > 1.0E-50) { dx = 1.0E-50; retn = delta_species(kspec, &dx); -#ifdef DEBUG_MODE - if (retn == 0) { - if (m_debug_print_lvl) { - plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", - m_speciesName[kspec].c_str(), kspec, dx); - } + if (DEBUG_MODE_ENABLED && retn == 0 && m_debug_print_lvl) { + plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", + m_speciesName[kspec].c_str(), kspec, dx); } -#endif } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (retn != 0) { plogf(" --- add_deleted(): species %s added back in with mol number %g", m_speciesName[kspec].c_str(), dx); @@ -2668,7 +2543,6 @@ size_t VCS_SOLVE::vcs_add_all_deleted() plogendl(); } } -#endif } } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); @@ -2685,14 +2559,12 @@ size_t VCS_SOLVE::vcs_add_all_deleted() VCS_DELETE_MINORSPECIES_CUTOFF) || (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { retn++; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], m_deltaGRxn_old[irxn]); plogendl(); } -#endif } } } @@ -2728,16 +2600,13 @@ bool VCS_SOLVE::vcs_globStepDamp() } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1); plogf(" --- subroutine FORCE: End Slope = %g\n", s2); } -#endif if (s1 > 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments,"); if (s1 < 1.0E-40) { plogf(" s1 positive but really small"); @@ -2746,17 +2615,14 @@ bool VCS_SOLVE::vcs_globStepDamp() } plogendl(); } -#endif return false; } if (s2 <= 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments, s2 < 0"); plogendl(); } -#endif return false; } @@ -2768,28 +2634,22 @@ bool VCS_SOLVE::vcs_globStepDamp() al = s1 / (s1 - s2); } if (al >= 0.95 || al < 0.0) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments (al = %g)\n", al); } -#endif return false; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced a damping factor = %g\n", al); } -#endif /* *************************************************** */ /* **** ADJUST MOLE NUMBERS, CHEM. POT *************** */ /* *************************************************** */ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_tmp), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc); } -#endif dptr = VCS_DATA_PTR(m_molNumSpecies_new); for (size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { @@ -2801,12 +2661,10 @@ bool VCS_SOLVE::vcs_globStepDamp() } vcs_updateVP(VCS_STATECALC_NEW); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE adjusted the mole " "numbers, AL = %10.3f\n", al); } -#endif /* * Because we changed the mole numbers, we need to * calculate the chemical potentials again. If a major- @@ -2831,13 +2689,10 @@ bool VCS_SOLVE::vcs_globStepDamp() } } - -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE: Adj End Slope = %g", s2); plogendl(); } -#endif return true; } @@ -2850,8 +2705,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], size_t jlose = npos; double* dptr, *scrxn_ptr; Cantera::clockWC tickTock; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" "); for (size_t i=0; i<77; i++) { plogf("-"); @@ -2889,7 +2743,6 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], plogendl(); } } -#endif /* * Calculate the maximum value of the number of components possible @@ -3048,12 +2901,10 @@ 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; } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Total number of components found = %3d (ne = %d)\n ", ncTrial, m_numElemConstraints); } -#endif goto L_END_LOOP; } /* @@ -3113,8 +2964,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], /* **** REARRANGE THE DATA ****************** */ /* ****************************************** */ if (jr != k) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", (m_speciesName[k]).c_str()); if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf("(Volts = %9.2g)", m_molNumSpecies_old[k]); @@ -3129,23 +2979,18 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], } plogf(" as component %3d\n", jr); } -#endif vcs_switch_pos(false, jr, k); std::swap(aw[jr], aw[k]); } -#ifdef DEBUG_MODE - else { - if (m_debug_print_lvl >= 2) { - plogf(" --- %-12.12s", m_speciesName[k].c_str()); - if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf("(Volts = %9.2g) remains ", m_molNumSpecies_old[k]); - } else { - plogf("(%9.2g) remains ", m_molNumSpecies_old[k]); - } - plogf(" as component %3d\n", jr); + else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" --- %-12.12s", m_speciesName[k].c_str()); + if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf("(Volts = %9.2g) remains ", m_molNumSpecies_old[k]); + } else { + plogf("(%9.2g) remains ", m_molNumSpecies_old[k]); } + plogf(" as component %3d\n", jr); } -#endif /* - entry point from up above */ L_END_LOOP: ; @@ -3284,9 +3129,7 @@ L_END_LOOP: m_scSize[i] = szTmp; } - -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Components:"); for (size_t j = 0; j < ncTrial; j++) { plogf(" %3d", j); @@ -3370,7 +3213,6 @@ L_END_LOOP: } plogf("\n"); } -#endif /* **************************************************** */ /* **** EVALUATE DELTA N VALUES *********************** */ /* **************************************************** */ @@ -3508,13 +3350,11 @@ 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) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %s can not be nonzero because" " needed element %s is zero\n", m_speciesName[kspec].c_str(), (m_elementName[j]).c_str()); } -#endif if (m_SSPhase[kspec]) { return VCS_SPECIES_ZEROEDSS; } else { @@ -3541,13 +3381,11 @@ 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) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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].c_str(), m_speciesName[j].c_str()); } -#endif if (m_SSPhase[kspec]) { return VCS_SPECIES_ZEROEDSS; } else { @@ -3558,14 +3396,12 @@ 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) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && 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", m_speciesName[kspec].c_str(), m_speciesName[j].c_str()); } -#endif if (m_SSPhase[kspec]) { return VCS_SPECIES_ZEROEDSS; } else { @@ -3702,16 +3538,16 @@ void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc, } } if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { -#ifdef DEBUG_MODE - if (molNum[kspec] != phi) { - plogf("We have an inconsistency!\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + if (molNum[kspec] != phi) { + plogf("We have an inconsistency!\n"); + exit(EXIT_FAILURE); + } + if (m_chargeSpecies[kspec] != -1.0) { + plogf("We have an unexpected situation!\n"); + exit(EXIT_FAILURE); + } } - if (m_chargeSpecies[kspec] != -1.0) { - plogf("We have an unexpected situation!\n"); - exit(EXIT_FAILURE); - } -#endif mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi; } else { if (m_SSPhase[kspec]) { @@ -3744,25 +3580,19 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new); actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new); molNum = VCS_DATA_PTR(m_molNumSpecies_new); - } -#ifdef DEBUG_MODE - else { + } else if (DEBUG_MODE_ENABLED) { plogf("vcs_dfe: wrong stateCalc value"); plogf(" --- Subroutine vcs_dfe called with bad stateCalc value: %d", stateCalc); plogendl(); exit(EXIT_FAILURE); } -#endif -#ifdef DEBUG_MODE - if (m_unitsState == VCS_DIMENSIONAL_G) { + if (DEBUG_MODE_ENABLED && m_unitsState == VCS_DIMENSIONAL_G) { printf("vcs_dfe: called with wrong units state\n"); exit(EXIT_FAILURE); } -#endif -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (ll == 0) { if (lbot != 0) { plogf(" --- Subroutine vcs_dfe called for one species: "); @@ -3780,7 +3610,6 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, } plogendl(); } -#endif double* tlogMoles = VCS_DATA_PTR(m_TmpPhase); /* @@ -3798,15 +3627,15 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, tlogMoles[iph] += molNum[kspec]; } } -#ifdef DEBUG_MODE - for (size_t iph = 0; iph < m_numPhases; iph++) { - if (! vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph])) { - plogf("phase Moles may be off, iph = %d, %20.14g %20.14g \n", - iph, tlogMoles[iph], tPhMoles_ptr[iph]); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + for (size_t iph = 0; iph < m_numPhases; iph++) { + if (! vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph])) { + plogf("phase Moles may be off, iph = %d, %20.14g %20.14g \n", + iph, tlogMoles[iph], tPhMoles_ptr[iph]); + exit(EXIT_FAILURE); + } } } -#endif vcs_dzero(tlogMoles, m_numPhases); for (size_t iph = 0; iph < m_numPhases; iph++) { if (tPhMoles_ptr[iph] > 0.0) { @@ -3848,16 +3677,15 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, for (size_t kspec = l1; kspec < l2; ++kspec) { size_t iphase = m_phaseID[kspec]; if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { -#ifdef DEBUG_MODE - if (molNum[kspec] != m_phasePhi[iphase]) { - plogf("We have an inconsistency!\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + if (molNum[kspec] != m_phasePhi[iphase]) { + plogf("We have an inconsistency!\n"); + exit(EXIT_FAILURE); + } else if (m_chargeSpecies[kspec] != -1.0) { + plogf("We have an unexpected situation!\n"); + exit(EXIT_FAILURE); + } } - if (m_chargeSpecies[kspec] != -1.0) { - plogf("We have an unexpected situation!\n"); - exit(EXIT_FAILURE); - } -#endif feSpecies[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; } else { @@ -3898,16 +3726,15 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { size_t iphase = m_phaseID[kspec]; if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { -#ifdef DEBUG_MODE - if (molNum[kspec] != m_phasePhi[iphase]) { - plogf("We have an inconsistency!\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + if (molNum[kspec] != m_phasePhi[iphase]) { + plogf("We have an inconsistency!\n"); + exit(EXIT_FAILURE); + } else if (m_chargeSpecies[kspec] != -1.0) { + plogf("We have an unexpected situation!\n"); + exit(EXIT_FAILURE); + } } - if (m_chargeSpecies[kspec] != -1.0) { - plogf("We have an unexpected situation!\n"); - exit(EXIT_FAILURE); - } -#endif feSpecies[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; } else { @@ -3949,16 +3776,15 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { size_t iphase = m_phaseID[kspec]; if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { -#ifdef DEBUG_MODE - if (molNum[kspec] != m_phasePhi[iphase]) { - plogf("We have an inconsistency!\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + if (molNum[kspec] != m_phasePhi[iphase]) { + plogf("We have an inconsistency!\n"); + exit(EXIT_FAILURE); + } else if (m_chargeSpecies[kspec] != -1.0) { + plogf("We have an unexpected situation!\n"); + exit(EXIT_FAILURE); + } } - if (m_chargeSpecies[kspec] != -1.0) { - plogf("We have an unexpected situation!\n"); - exit(EXIT_FAILURE); - } -#endif feSpecies[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ; } else { @@ -4189,32 +4015,26 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState) Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_tPhaseMoles_new)); - } -#ifdef DEBUG_MODE - else { + } else if (DEBUG_MODE_ENABLED) { plogf("vcs_updateVP ERROR: wrong stateCalc value: %d", vcsState); plogendl(); exit(EXIT_FAILURE); } -#endif } } bool VCS_SOLVE::vcs_evaluate_speciesType() { m_numRxnMinorZeroed = 0; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Species Status decision is reevaluated: All species are minor except for:\n"); - } else if (m_debug_print_lvl >= 5) { + } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) { plogf(" --- Species Status decision is reevaluated"); plogendl(); } -#endif for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { m_speciesStatus[kspec] = vcs_species_type(kspec); -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 5) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) { plogf(" --- %-16s: ", m_speciesName[kspec].c_str()); if (kspec < m_numComponents) { plogf("(COMP) "); @@ -4225,7 +4045,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() const char* sString = vcs_speciesType_string(m_speciesStatus[kspec], 100); plogf("%s\n", sString); - } else if (m_debug_print_lvl >= 2) { + } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { switch (m_speciesStatus[kspec]) { case VCS_SPECIES_COMPONENT: @@ -4264,20 +4084,16 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() } } } -#endif if (kspec >= m_numComponents) { if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { ++m_numRxnMinorZeroed; } } } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" ---"); plogendl(); } -#endif - return (m_numRxnMinorZeroed >= m_numRxnRdc); } @@ -4323,8 +4139,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, exit(EXIT_FAILURE); } -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Subroutine vcs_deltag called for "); if (l < 0) { plogf("major noncomponents\n"); @@ -4334,7 +4149,6 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, plogf("minor noncomponents\n"); } } -#endif /* ************************************************* */ /* **** MAJORS and ZEROED SPECIES ONLY ************* */ /* ************************************************* */ @@ -4646,14 +4460,11 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, feSpecies = VCS_DATA_PTR(m_feSpecies_old); deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old); actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old); - } -#ifdef DEBUG_MODE - else { + } else if (DEBUG_MODE_ENABLED) { plogf("vcs_deltag_Phase: we shouldn't be here\n"); plogendl(); exit(EXIT_FAILURE); } -#endif size_t irxnl = m_numRxnRdc; if (doDeleted) { @@ -4661,24 +4472,20 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, } vcs_VolPhase* vPhase = m_VolPhaseList[iphase]; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- Subroutine vcs_deltag_Phase called for phase %d\n", iphase); } -#endif /* * Single species Phase */ if (vPhase->m_singleSpecies) { size_t kspec = vPhase->spGlobalIndexVCS(0); -#ifdef DEBUG_MODE - if (iphase != m_phaseID[kspec]) { + if (DEBUG_MODE_ENABLED && iphase != m_phaseID[kspec]) { plogf("vcs_deltag_Phase index error\n"); exit(EXIT_FAILURE); } -#endif if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; deltaGRxn[irxn] = feSpecies[kspec]; @@ -4784,13 +4591,11 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k if (k1 == k2) { return; } -#ifdef DEBUG_MODE - if (k1 > (m_numSpeciesTot - 1) || + if (DEBUG_MODE_ENABLED && k1 > (m_numSpeciesTot - 1) || k2 > (m_numSpeciesTot - 1)) { plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n", k1, k2); } -#endif /* * Handle the index pointer in the phase structures first */ @@ -4799,16 +4604,15 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k size_t kp1 = m_speciesLocalPhaseIndex[k1]; size_t kp2 = m_speciesLocalPhaseIndex[k2]; -#ifdef DEBUG_MODE - if (pv1->spGlobalIndexVCS(kp1) != k1) { - plogf("Indexing error in program\n"); - exit(EXIT_FAILURE); + if (DEBUG_MODE_ENABLED) { + if (pv1->spGlobalIndexVCS(kp1) != k1) { + plogf("Indexing error in program\n"); + exit(EXIT_FAILURE); + } else if (pv2->spGlobalIndexVCS(kp2) != k2) { + plogf("Indexing error in program\n"); + exit(EXIT_FAILURE); + } } - if (pv2->spGlobalIndexVCS(kp2) != k2) { - plogf("Indexing error in program\n"); - exit(EXIT_FAILURE); - } -#endif pv1->setSpGlobalIndexVCS(kp1, k2); pv2->setSpGlobalIndexVCS(kp2, k1); std::swap(m_speciesName[k1], m_speciesName[k2]); @@ -4851,13 +4655,11 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k */ size_t i1 = k1 - m_numComponents; size_t i2 = k2 - m_numComponents; -#ifdef DEBUG_MODE - if (i1 > (m_numRxnTot - 1) || + if (DEBUG_MODE_ENABLED && i1 > (m_numRxnTot - 1) || i2 > (m_numRxnTot - 1)) { plogf("switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n", i1 , i2); } -#endif for (size_t j = 0; j < m_numComponents; ++j) { std::swap(m_stoichCoeffRxnMatrix[i1][j], m_stoichCoeffRxnMatrix[i2][j]); } @@ -4881,15 +4683,14 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) return dx; } double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF; -#ifdef DEBUG_MODE + // Check to make sure that species is zero in the solution vector // If it isn't, we don't know what's happening - if (m_molNumSpecies_old[kspec] != 0.0) { + if (DEBUG_MODE_ENABLED && m_molNumSpecies_old[kspec] != 0.0) { w_kspec = 0.0; plogf("vcs_birthGuess:: we shouldn't be here\n"); exit(EXIT_FAILURE); } -#endif int ss = m_SSPhase[kspec]; if (!ss) { /* diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp index eb8b2858b..621516a71 100644 --- a/src/equil/vcs_species_thermo.cpp +++ b/src/equil/vcs_species_thermo.cpp @@ -185,9 +185,6 @@ double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin, double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) { -#ifdef DEBUG_MODE - char yo[] = "VS_SPECIES_THERMO::G0_R_calc "; -#endif double fe, H, S; if (SS0_Model == VCS_SS0_CONSTANT) { return SS0_feSave; @@ -216,9 +213,7 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) fe = H - TKelvin * S; break; default: -#ifdef DEBUG_MODE - plogf("%sERROR: unknown model\n", yo); -#endif + plogf("VS_SPECIES_THERMO::G0_R_calc ERROR: unknown model\n"); exit(EXIT_FAILURE); } } @@ -229,9 +224,6 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) double VCS_SPECIES_THERMO::eval_ac(size_t kglob) { -#ifdef DEBUG_MODE - char yo[] = "VCS_SPECIES_THERMO::eval_ac "; -#endif double ac; /* * Activity coefficients are frequently evaluated on a per phase @@ -248,9 +240,7 @@ double VCS_SPECIES_THERMO::eval_ac(size_t kglob) ac = 1.0; break; default: -#ifdef DEBUG_MODE - plogf("%sERROR: unknown model\n", yo); -#endif + plogf("VCS_SPECIES_THERMO::eval_ac ERROR: unknown model\n"); exit(EXIT_FAILURE); } } diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index 28aa108fe..136cf4d28 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -381,11 +381,9 @@ void NonlinearSolver::createSolnWeights(const doublereal* const y) { for (size_t i = 0; i < neq_; i++) { m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i]; -#ifdef DEBUG_MODE - if (m_ewt[i] <= 0.0) { + if (DEBUG_MODE_ENABLED && m_ewt[i] <= 0.0) { throw CanteraError(" NonlinearSolver::createSolnWeights()", "ewts <= 0.0"); } -#endif } } @@ -497,19 +495,19 @@ doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const doublereal sum_norm = 0.0, error; for (size_t i = 0; i < neq_; i++) { -#ifdef DEBUG_MODE - checkFinite(resid[i]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(resid[i]); + } error = resid[i] / m_residWts[i]; -#ifdef DEBUG_MODE - checkFinite(error); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(error); + } sum_norm += (error * error); } sum_norm = sqrt(sum_norm / neq_); -#ifdef DEBUG_MODE - checkFinite(sum_norm); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(sum_norm); + } if (printLargest) { const int num_entries = printLargest; doublereal dmax1, normContrib; @@ -694,9 +692,9 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, } else { m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol]; } -#ifdef DEBUG_MODE - checkFinite(m_rowWtScales[irow]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(m_rowWtScales[irow]); + } jptr++; } } @@ -718,9 +716,9 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, } else { m_rowWtScales[irow] += vv * m_ewt[jcol]; } -#ifdef DEBUG_MODE - checkFinite(m_rowWtScales[irow]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(m_rowWtScales[irow]); + } } } } @@ -1280,9 +1278,9 @@ doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac) } deltaX_CP_[j] -= m_resid[i] * jac(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j] / (m_residWts[i] * m_residWts[i]); -#ifdef DEBUG_MODE - checkFinite(deltaX_CP_[j]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(deltaX_CP_[j]); + } } } @@ -2664,10 +2662,8 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c doublereal stepNorm_1; doublereal stepNorm_2; -#ifdef DEBUG_MODE int legBest; doublereal alphaBest; -#endif bool trInit = false; copy(y_comm, y_comm + neq_, m_y_n_curr.begin()); @@ -2733,31 +2729,11 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c /* * If we are far enough away from the solution, redo the solution weights and the trust vectors. */ - if (m_normDeltaSoln_Newton > 1.0E2) { + if (m_normDeltaSoln_Newton > 1.0E2 || (num_newt_its % 5) == 1) { createSolnWeights(DATA_PTR(m_y_n_curr)); -#ifdef DEBUG_MODE - if (trInit) { + if (trInit && (DEBUG_MODE_ENABLED || doDogLeg_)) { readjustTrustVector(); } -#else - if (doDogLeg_ && trInit) { - readjustTrustVector(); - } -#endif - } else { - // Do this stuff every 5 iterations - if ((num_newt_its % 5) == 1) { - createSolnWeights(DATA_PTR(m_y_n_curr)); -#ifdef DEBUG_MODE - if (trInit) { - readjustTrustVector(); - } -#else - if (doDogLeg_ && trInit) { - readjustTrustVector(); - } -#endif - } } /* @@ -2829,21 +2805,12 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c } } - -#ifdef DEBUG_MODE - if (m_print_flag > 3) { - printf("\t solve_nonlinear_problem(): Calculate the steepest descent direction and Cauchy Point\n"); - } - m_normDeltaSoln_CP = doCauchyPointSolve(jac); - -#else - if (doDogLeg_) { + if (DEBUG_MODE_ENABLED || doDogLeg_) { if (m_print_flag > 3) { printf("\t solve_nonlinear_problem(): Calculate the steepest descent direction and Cauchy Point\n"); } m_normDeltaSoln_CP = doCauchyPointSolve(jac); } -#endif // compute the undamped Newton step if (doAffineSolve_) { @@ -2889,25 +2856,17 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c } } - - if (doDogLeg_) { - - - -#ifdef DEBUG_MODE + if (DEBUG_MODE_ENABLED && doDogLeg_ && m_print_flag >= 4) { doublereal trustD = calcTrustDistance(m_step_1); - if (m_print_flag >= 4) { - if (trustD > trustDelta_) { - printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n", - trustD, trustDelta_); - printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n", - trustD, trustDelta_); - } else { - printf("\t\t Newton's method step size, %g trustVectorUnits, smaller than trust region, %g trustVectorUnits\n", - trustD, trustDelta_); - } + if (trustD > trustDelta_) { + printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n", + trustD, trustDelta_); + printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n", + trustD, trustDelta_); + } else { + printf("\t\t Newton's method step size, %g trustVectorUnits, smaller than trust region, %g trustVectorUnits\n", + trustD, trustDelta_); } -#endif } /* @@ -2930,26 +2889,19 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c if (doDogLeg_) { setupDoubleDogleg(); -#ifdef DEBUG_MODE - if (s_print_DogLeg && m_print_flag >= 5) { + if (DEBUG_MODE_ENABLED && s_print_DogLeg && m_print_flag >= 5) { printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n"); residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest); } -#endif if (m_print_flag >= 4) { printf("\t solve_nonlinear_problem(): Calculate damping along dog-leg path to ensure residual decrease\n"); } retnDamp = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), m_step_1, DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial), stepNorm_1, stepNorm_2, jac, i_numTrials); + } else if (DEBUG_MODE_ENABLED && s_print_DogLeg && m_print_flag >= 5) { + printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n"); + residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest); } -#ifdef DEBUG_MODE - else { - if (s_print_DogLeg && m_print_flag >= 5) { - printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n"); - residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest); - } - } -#endif // Damp the Newton step /* @@ -3336,11 +3288,11 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, return info; } m_nJacEval++; -#ifdef DEBUG_MODE - for (int ii = 0; ii < neq_; ii++) { - checkFinite(f[ii]); + if (DEBUG_MODE_ENABLED) { + for (int ii = 0; ii < neq_; ii++) { + checkFinite(f[ii]); + } } -#endif /* * Malloc a vector and call the function object to return a set of @@ -3403,14 +3355,14 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, JacDelta_ResidEval, j, dy); m_nfe++; -#ifdef DEBUG_MODE - if (fabs(dy) < 1.0E-300) { - throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); + if (DEBUG_MODE_ENABLED) { + if (fabs(dy) < 1.0E-300) { + throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); + } + for (int ii = 0; ii < neq_; ii++) { + checkFinite(m_wksp[ii]); + } } - for (int ii = 0; ii < neq_; ii++) { - checkFinite(m_wksp[ii]); - } -#endif if (info != 1) { return info; @@ -3483,14 +3435,14 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp), JacDelta_ResidEval, static_cast(j), dy); m_nfe++; -#ifdef DEBUG_MODE - if (fabs(dy) < 1.0E-300) { - throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); + if (DEBUG_MODE_ENABLED) { + if (fabs(dy) < 1.0E-300) { + throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero"); + } + for (int ii = 0; ii < neq_; ii++) { + checkFinite(m_wksp[ii]); + } } - for (int ii = 0; ii < neq_; ii++) { - checkFinite(m_wksp[ii]); - } -#endif if (info != 1) { return info; } @@ -3613,17 +3565,17 @@ void NonlinearSolver::computeResidWts() if (checkUserResidualTols_ == 1) { for (size_t i = 0; i < neq_; i++) { m_residWts[i] = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / rtol_; -#ifdef DEBUG_MODE - checkFinite(m_residWts[i]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(m_residWts[i]); + } } } else { doublereal sum = 0.0; for (size_t i = 0; i < neq_; i++) { m_residWts[i] = m_rowWtScales[i]; -#ifdef DEBUG_MODE - checkFinite(m_residWts[i]); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(m_residWts[i]); + } sum += m_residWts[i]; } sum /= neq_; diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 582511ca4..3713035a1 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -14,6 +14,8 @@ // turn on debugging for now #ifndef DEBUG_MODE #define DEBUG_MODE +#undef DEBUG_MODE_ENABLED +#define DEBUG_MODE_ENABLED 1 #endif #include "cantera/base/global.h" @@ -32,7 +34,6 @@ namespace Cantera #define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 )) #endif -#ifdef DEBUG_MODE //! Print out a form for the current function evaluation /*! * @param fp Pointer to the FILE object @@ -51,7 +52,6 @@ static void print_funcEval(FILE* fp, doublereal xval, doublereal fval, int its) fprintf(fp,"...............................................................\n"); fprintf(fp,"\n"); } -#endif RootFind::RootFind(ResidEval* resid) : m_residFunc(resid), @@ -196,10 +196,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun int converged = 0; int bottomBump = 0; int topBump = 0; -#ifdef DEBUG_MODE - char fileName[80]; FILE* fp = 0; -#endif int doFinalFuncCall = 0; doublereal x1, x2, xnew, f1, f2, fnew, slope; doublereal deltaX2 = 0.0, deltaXnew = 0.0; @@ -223,19 +220,16 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun rfT.reasoning = "First Point: "; callNum++; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { + char fileName[80]; sprintf(fileName, "RootFind_%d.log", callNum); fp = fopen(fileName, "w"); fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n"); fprintf(fp, "-----------------------------------------------------" "-------------------------------\n"); - } -#else - if (printLvl >= 3) { + } else if (printLvl >= 3) { writelog("WARNING: RootFind: printlvl >= 3, but debug mode not turned on\n"); } -#endif if (xmax <= xmin) { writelogf("%sxmin and xmax are bad: %g %g\n", stre, xmin, xmax); funcTargetValue = func(*xbest); @@ -289,13 +283,10 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun int its = 1; f1 = func(x1); -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { print_funcEval(fp, x1, f1, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1); } -#endif - if (f1 == 0.0) { *xbest = x1; @@ -345,12 +336,10 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun deltaX2 = x2 - x1; its++; f2 = func(x2); -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { print_funcEval(fp, x2, f2, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2); } -#endif /* * Calculate the norm of the function, this is the nominal value of f. We try @@ -410,12 +399,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun * Find an estimate of the next point, xnew, to try based on * a linear approximation from the last two points. */ -#ifdef DEBUG_MODE - if (fabs(x2 - x1) < 1.0E-14) { + if (DEBUG_MODE_ENABLED && fabs(x2 - x1) < 1.0E-14) { printf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1); } -#endif - doublereal delXtmp = deltaXControlled(x2, x1); slope = (f2 - f1) / delXtmp; rfT.slope = slope; @@ -441,11 +427,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } rfT.reasoning += "Slope is good. "; } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlin = %-11.5E", xnew); } -#endif deltaXnew = xnew - x2; /* * If the suggested step size is too big, throw out step @@ -541,19 +525,15 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xDelMin = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < xDelMin) { xnew = x1 + DSIGN(xnew-x1) * xDelMin; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } -#endif } if (fabs(xnew - x2) < 0.1 * xDelMin) { xnew = x2 + DSIGN(xnew-x2) * 0.1 * xDelMin; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } -#endif } } else { /* @@ -569,11 +549,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } if (fabs(xDelMax) < fabs(xnew - x2)) { xnew = x2 + DSIGN(xnew-x2) * xDelMax; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitsize = %-11.5E", xnew); } -#endif } /* * If we are doing a jump outside the two previous points, make sure @@ -583,28 +561,22 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xDelMin = 0.1 * fabs(x2 - x1); if (fabs(xnew - x2) < xDelMin) { xnew = x2 + DSIGN(xnew - x2) * xDelMin; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } -#endif } if (fabs(xnew - x1) < xDelMin) { xnew = x1 + DSIGN(xnew - x1) * xDelMin; -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } -#endif } } /* * HKM -> Not sure this section is needed */ if (foundStraddle) { -#ifdef DEBUG_MODE double xorig = xnew; -#endif if (posStraddle) { if (f2 > 0.0) { if (xnew > x2) { @@ -638,13 +610,11 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { if (xorig != xnew) { fprintf(fp, " | xstraddle = %-11.5E", xnew); } } -#endif } /* @@ -686,11 +656,9 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = xmax; } } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitmax = %-11.5E", xnew); } -#endif } if (xnew < xmin) { bottomBump++; @@ -712,23 +680,19 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun xnew = xmin; } } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitmin = %-11.5E", xnew); } -#endif } its++; fnew = func(xnew); -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp,"\n"); print_funcEval(fp, xnew, fnew, its); fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew); } -#endif rfT.xval = xnew; rfT.fval = fnew; rfT.its = its; @@ -1063,11 +1027,9 @@ done: if (printLvl >= 1) { writelogf("RootFind success: convergence achieved\n"); } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | RootFind success in %d its, fnorm = %g\n", its, fnorm); } -#endif rfHistory_.push_back(rfT); } else { rfT.reasoning = "FAILED CONVERGENCE "; @@ -1090,11 +1052,9 @@ done: } rfT.reasoning += "Maximum iterations exceeded without convergence, cause unknown"; } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, "\nRootFind failure in %d its\n", its); } -#endif *xbest = x2; funcTargetValue = f2 + m_funcTargetValue; @@ -1102,11 +1062,9 @@ done: rfT.fval = f2; rfHistory_.push_back(rfT); } -#ifdef DEBUG_MODE - if (printLvl >= 3 && writeLogAllowed_) { + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fclose(fp); } -#endif if (printLvl >= 2) { printTable(); @@ -1118,13 +1076,13 @@ done: doublereal RootFind::func(doublereal x) { doublereal r; -#ifdef DEBUG_MODE - checkFinite(x); -#endif + if (DEBUG_MODE_ENABLED) { + checkFinite(x); + } m_residFunc->evalSS(0.0, &x, &r); -#ifdef DEBUG_MODE - checkFinite(r); -#endif + if (DEBUG_MODE_ENABLED) { + 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 60619e8cc..4482c58d1 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -1912,11 +1912,9 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const int z1, z2; size_t n, i, j, m, counterIJ, counterIJ2; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("\n Debugging information from hmw_act \n"); } -#endif /* * Make sure the counter variables are setup */ @@ -1939,13 +1937,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ m_IionicMolality = Is; sqrtIs = sqrt(Is); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 1: \n"); printf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); } -#endif /* * The following call to calc_lambdas() calculates all 16 elements @@ -1959,30 +1955,24 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * E-thetaprime for all combinations of positive * unlike charges up to 4 */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 2: \n"); } -#endif for (z1 = 1; z1 <=4; z1++) { for (z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 3: \n"); printf(" Species Species g(x) " " hfunc(x) \n"); } -#endif /* * calculate g(x) and hfunc(x) for each cation-anion pair MX @@ -2029,14 +2019,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const gfunc[counterIJ] = 0.0; hfunc[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), gfunc[counterIJ], hfunc[counterIJ]); } -#endif } } @@ -2044,13 +2032,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * --------- SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX ---------- * --------- Agrees with Pitzer, Eq. (49), (51), (55) */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 4: \n"); printf(" Species Species BMX " "BprimeMX BphiMX \n"); } -#endif for (i = 1; i < m_kk - 1; i++) { for (j = i+1; j < m_kk; j++) { @@ -2068,13 +2054,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const BMX[counterIJ] = beta0MX[counterIJ] + beta1MX[counterIJ] * gfunc[counterIJ] + beta2MX[counterIJ] * g2func[counterIJ]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("%d %g: %g %g %g %g\n", (int) counterIJ, BMX[counterIJ], beta0MX[counterIJ], beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]); } -#endif if (Is > 1.0E-150) { BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is + beta2MX[counterIJ] * h2func[counterIJ]/Is); @@ -2087,15 +2072,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const BprimeMX[counterIJ] = 0.0; BphiMX[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", sni.c_str(), snj.c_str(), BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ]); } -#endif } } @@ -2103,12 +2086,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * --------- SUBSECTION TO CALCULATE CMX ---------- * --------- Agrees with Pitzer, Eq. (53). */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 5: \n"); printf(" Species Species CMX \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -2126,14 +2107,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } else { CMX[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), CMX[counterIJ]); } -#endif } } @@ -2141,13 +2120,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- * --------- Agrees with Pitzer, Eq. 72, 73, 74 */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 6: \n"); printf(" Species Species Phi_ij " " Phiprime_ij Phi^phi_ij \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -2170,15 +2147,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const Phiprime[counterIJ] = 0.0; Phiphi[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", sni.c_str(), snj.c_str(), Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ]); } -#endif } } @@ -2186,19 +2161,15 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * ------------- SUBSECTION FOR CALCULATION OF F ---------------------- * ------------ Agrees with Pitzer Eqn. (65) -------------------------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 7: \n"); } -#endif Aphi = A_Debye_TP() / 3.0; F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" initial value of F = %10.6f \n", F); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -2220,18 +2191,14 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(i)*charge(j) > 0) { F = F + molality[i]*molality[j] * Phiprime[counterIJ]; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" F = %10.6f \n", F); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 8: Summing in All Contributions to Activity Coefficients \n"); } -#endif for (i = 1; i < m_kk; i++) { @@ -2241,20 +2208,15 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * -> Equations agree with Pitzer, eqn.(63) */ if (charge(i) > 0.0) { - -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str()); } -#endif // species i is the cation (positive) to calc the actcoeff zsqF = charge(i)*charge(i)*F; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Unary term: z*z*F = %10.5f\n", zsqF); } -#endif sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; @@ -2271,15 +2233,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all anions sum1 = sum1 + molality[j]* (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { snj = speciesName(j) + ":"; printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(), molality[j]*2.0*BMX[counterIJ]); printf(" m_j Z CMX = %10.5f\n", molality[j]* molarcharge*CMX[counterIJ]); } -#endif if (j < m_kk-1) { /* * This term is the ternary interaction involving the @@ -2291,15 +2251,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (psi_ijk[n] != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } } -#endif } } } @@ -2310,15 +2268,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all cations if (j != i) { sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j] * Phi[counterIJ])!= 0.0) { snj = speciesName(j) + ":"; printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(), molality[j]*(2.0*Phi[counterIJ])); } } -#endif } for (size_t k = 1; k < m_kk; k++) { if (charge(k) < 0.0) { @@ -2326,15 +2282,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const n = k + j * m_kk + i * m_kk * m_kk; sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (psi_ijk[n] != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } } -#endif /* * Find the counterIJ for the j,k interaction */ @@ -2342,15 +2296,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const counterIJ2 = m_CounterIJ[n]; sum4 = sum4 + (fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); } } -#endif } } } @@ -2360,15 +2312,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ if (charge(j) == 0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { snj = speciesName(j) + ":"; printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), molality[j]*2.0*m_Lambda_nj(j,i)); } } -#endif /* * Zeta interaction term */ @@ -2380,13 +2330,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double zeta = psi_ijk[n]; if (zeta != 0.0) { sum5 = sum5 + molality[j]*molality[k]*zeta; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } -#endif } } } @@ -2398,13 +2346,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n", sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } -#endif } /* @@ -2413,21 +2359,15 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * -> Equations agree with Pitzer, eqn.(64) */ if (charge(i) < 0) { - -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str()); } -#endif - // species i is an anion (negative) zsqF = charge(i)*charge(i)*F; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Unary term: z*z*F = %10.5f\n", zsqF); } -#endif sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; @@ -2446,30 +2386,26 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(j) > 0) { sum1 = sum1 + molality[j]* (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { snj = speciesName(j) + ":"; printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(), molality[j]*2.0*BMX[counterIJ]); printf(" m_j Z CMX = %10.5f\n", molality[j]* molarcharge*CMX[counterIJ]); } -#endif if (j < m_kk-1) { for (size_t k = j+1; k < m_kk; k++) { // an inner sum over all cations if (charge(k) > 0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (psi_ijk[n] != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } } -#endif } } } @@ -2482,30 +2418,26 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all anions if (j != i) { sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j] * Phi[counterIJ])!= 0.0) { snj = speciesName(j) + ":"; printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(), molality[j]*(2.0*Phi[counterIJ])); } } -#endif } for (size_t k = 1; k < m_kk; k++) { if (charge(k) > 0.0) { // two inner sums over cations n = k + j * m_kk + i * m_kk * m_kk; sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (psi_ijk[n] != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } } -#endif /* * Find the counterIJ for the symmetric binary interaction */ @@ -2514,15 +2446,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const sum4 = sum4 + (fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); } } -#endif } } } @@ -2532,15 +2462,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ if (charge(j) == 0.0) { sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { snj = speciesName(j) + ":"; printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), molality[j]*2.0*m_Lambda_nj(j,i)); } } -#endif /* * Zeta interaction term */ @@ -2553,13 +2481,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double zeta = psi_ijk[n]; if (zeta != 0.0) { sum5 = sum5 + molality[j]*molality[k]*zeta; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } -#endif } } } @@ -2567,13 +2493,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n", sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } -#endif } /* * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- @@ -2581,26 +2505,21 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const * -> Equations agree with Pitzer, */ if (charge(i) == 0.0) { -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str()); } -#endif sum1 = 0.0; sum3 = 0.0; for (j = 1; j < m_kk; j++) { sum1 = sum1 + molality[j]*2.0*m_Lambda_nj(i,j); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (m_Lambda_nj(i,j) != 0.0) { snj = speciesName(j) + ":"; printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(), molality[j]*2.0*m_Lambda_nj(i,j)); } } -#endif - /* * Zeta term -> we piggyback on the psi term */ @@ -2609,46 +2528,37 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (psi_ijk[n] != 0.0) { snj = speciesName(j) + "," + speciesName(k) + ":"; printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), molality[j]*molality[k]*psi_ijk[n]); } } -#endif } } } } sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { if (m_Mu_nnn[i] != 0.0) { printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n", 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]); } } -#endif - m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n", sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 9: \n"); } -#endif /* * -------- SUBSECTION FOR CALCULATING THE OSMOTIC COEFF --------- * -------- -> equations agree with my notes, Eqn. (117). @@ -2790,19 +2700,14 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const } else { osmotic_coef = 1.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" 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); printf(" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n", sum_m_phi_minus_1, osmotic_coef); - } - - if (m_debugCalc) { printf(" Step 10: \n"); } -#endif lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef; /* @@ -2816,15 +2721,13 @@ 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); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { double wateract = exp(lnwateract); printf(" Weight of Solvent = %16.7g\n", m_weightSolvent); printf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped); printf(" ln_a_water=%10.6f a_water=%10.6f\n\n", lnwateract, wateract); } -#endif } void HMWSoln::s_update_dlnMolalityActCoeff_dT() const @@ -2933,12 +2836,10 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const int z1, z2; size_t n, i, j, m, counterIJ, counterIJ2; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("\n Debugging information from " "s_Pitzer_dlnMolalityActCoeff_dT()\n"); } -#endif /* * Make sure the counter variables are setup */ @@ -2961,13 +2862,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const */ m_IionicMolality = Is; sqrtIs = sqrt(Is); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 1: \n"); printf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); } -#endif /* * The following call to calc_lambdas() calculates all 16 elements @@ -2981,30 +2880,24 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * E-thetaprime for all combinations of positive * unlike charges up to 4 */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 2: \n"); } -#endif for (z1 = 1; z1 <=4; z1++) { for (z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 3: \n"); printf(" Species Species g(x) " " hfunc(x) \n"); } -#endif /* * calculate g(x) and hfunc(x) for each cation-anion pair MX @@ -3050,14 +2943,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const gfunc[counterIJ] = 0.0; hfunc[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), gfunc[counterIJ], hfunc[counterIJ]); } -#endif } } @@ -3066,13 +2957,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * ------- These are now temperature derivatives of the * previously calculated quantities. */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 4: \n"); printf(" Species Species BMX " "BprimeMX BphiMX \n"); } -#endif for (i = 1; i < m_kk - 1; i++) { for (j = i+1; j < m_kk; j++) { @@ -3089,13 +2978,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const BMX_L[counterIJ] = beta0MX_L[counterIJ] + beta1MX_L[counterIJ] * gfunc[counterIJ] + beta2MX_L[counterIJ] * gfunc[counterIJ]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("%d %g: %g %g %g %g\n", (int) counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ], beta1MX_L[counterIJ], beta2MX_L[counterIJ], gfunc[counterIJ]); } -#endif if (Is > 1.0E-150) { BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is + beta2MX_L[counterIJ] * h2func[counterIJ]/Is); @@ -3108,15 +2995,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const BprimeMX_L[counterIJ] = 0.0; BphiMX_L[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", sni.c_str(), snj.c_str(), BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]); } -#endif } } @@ -3124,12 +3009,10 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * --------- SUBSECTION TO CALCULATE CMX_L ---------- * --------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 5: \n"); printf(" Species Species CMX \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -3147,14 +3030,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const } else { CMX_L[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), CMX_L[counterIJ]); } -#endif } } @@ -3162,13 +3043,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- * -------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 6: \n"); printf(" Species Species Phi_ij " " Phiprime_ij Phi^phi_ij \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -3193,35 +3072,29 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const Phiprime[counterIJ] = 0.0; Phiphi_L[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", sni.c_str(), snj.c_str(), Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ]); } -#endif } } /* * ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 7: \n"); } -#endif double dA_DebyedT = dA_DebyedT_TP(); double dAphidT = dA_DebyedT /3.0; dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" initial value of dFdT = %10.6f \n", dFdT); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -3243,18 +3116,14 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const if (charge(i)*charge(j) > 0) { dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ]; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" dFdT = %10.6f \n", dFdT); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 8: \n"); } -#endif for (i = 1; i < m_kk; i++) { @@ -3348,15 +3217,13 @@ 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]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqdFdT, sum1, sum2, sum3, sum4, sum5); } -#endif } /* @@ -3442,15 +3309,13 @@ 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]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n", sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqdFdT, sum1, sum2, sum3, sum4, sum5); } -#endif } /* * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- @@ -3477,21 +3342,17 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const 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]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n", sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 9: \n"); } -#endif /* * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT --------- * @@ -3633,19 +3494,14 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const } else { d_osmotic_coef_dT = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" 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); printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n", sum_m_phi_minus_1, d_osmotic_coef_dT); - } - - if (m_debugCalc) { printf(" Step 10: \n"); } -#endif d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT; /* @@ -3657,13 +3513,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). */ m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { double d_wateract_dT = exp(d_lnwateract_dT); printf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n", d_lnwateract_dT, d_wateract_dT); } -#endif } void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const @@ -3764,12 +3618,10 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const int z1, z2; size_t n, i, j, m, counterIJ, counterIJ2; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("\n Debugging information from " "s_Pitzer_d2lnMolalityActCoeff_dT2()\n"); } -#endif /* * Make sure the counter variables are setup */ @@ -3793,13 +3645,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const */ m_IionicMolality = Is; sqrtIs = sqrt(Is); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 1: \n"); printf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); } -#endif /* * The following call to calc_lambdas() calculates all 16 elements @@ -3813,30 +3663,24 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const * E-thetaprime for all combinations of positive * unlike charges up to 4 */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 2: \n"); } -#endif for (z1 = 1; z1 <=4; z1++) { for (z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 3: \n"); printf(" Species Species g(x) " " hfunc(x) \n"); } -#endif /* * @@ -3883,14 +3727,12 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const gfunc[counterIJ] = 0.0; hfunc[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), gfunc[counterIJ], hfunc[counterIJ]); } -#endif } } /* @@ -3898,13 +3740,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const * ------- These are now temperature derivatives of the * previously calculated quantities. */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 4: \n"); printf(" Species Species BMX " "BprimeMX BphiMX \n"); } -#endif for (i = 1; i < m_kk - 1; i++) { for (j = i+1; j < m_kk; j++) { @@ -3921,13 +3761,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const BMX_LL[counterIJ] = beta0MX_LL[counterIJ] + beta1MX_LL[counterIJ] * gfunc[counterIJ] + beta2MX_LL[counterIJ] * g2func[counterIJ]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("%d %g: %g %g %g %g\n", (int) counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ], beta1MX_LL[counterIJ], beta2MX_LL[counterIJ], gfunc[counterIJ]); } -#endif if (Is > 1.0E-150) { BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is + beta2MX_LL[counterIJ] * h2func[counterIJ]/Is); @@ -3940,15 +3778,13 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const BprimeMX_LL[counterIJ] = 0.0; BphiMX_LL[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", sni.c_str(), snj.c_str(), BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]); } -#endif } } @@ -3956,12 +3792,10 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const * --------- SUBSECTION TO CALCULATE CMX_LL ---------- * --------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 5: \n"); printf(" Species Species CMX \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -3979,14 +3813,12 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } else { CMX_LL[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), CMX_LL[counterIJ]); } -#endif } } @@ -3994,13 +3826,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- * -------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 6: \n"); printf(" Species Species Phi_ij " " Phiprime_ij Phi^phi_ij \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -4023,34 +3853,28 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const Phiprime[counterIJ] = 0.0; Phiphi_LL[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", sni.c_str(), snj.c_str(), Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]); } -#endif } } /* * ----------- SUBSECTION FOR CALCULATION OF d2FdT2 --------------------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 7: \n"); } -#endif double d2AphidT2 = d2A_DebyedT2_TP() / 3.0; d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" initial value of d2FdT2 = %10.6f \n", d2FdT2); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -4072,18 +3896,14 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const if (charge(i)*charge(j) > 0) { d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ]; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" d2FdT2 = %10.6f \n", d2FdT2); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 8: \n"); } -#endif for (i = 1; i < m_kk; i++) { @@ -4176,15 +3996,13 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const */ m_d2lnActCoeffMolaldT2_Unscaled[i] = zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s d2lngammadT2[i]=%10.6f \n", sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqd2FdT2, sum1, sum2, sum3, sum4, sum5); } -#endif } @@ -4273,15 +4091,13 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } m_d2lnActCoeffMolaldT2_Unscaled[i] = zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s d2lngammadT2[i]=%10.6f\n", sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqd2FdT2, sum1, sum2, sum3, sum4, sum5); } -#endif } /* * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF ------- @@ -4307,21 +4123,17 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i]; m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s d2lngammadT2[i]=%10.6f \n", sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 9: \n"); } -#endif /* * ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 --------- @@ -4464,19 +4276,14 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const } else { d2_osmotic_coef_dT2 = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" 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); printf(" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n", sum_m_phi_minus_1, d2_osmotic_coef_dT2); - } - - if (m_debugCalc) { printf(" Step 10: \n"); } -#endif d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2; /* @@ -4489,13 +4296,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const */ m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { double d2_wateract_dT2 = exp(d2_lnwateract_dT2); printf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n", d2_lnwateract_dT2, d2_wateract_dT2); } -#endif } void HMWSoln::s_update_dlnMolalityActCoeff_dP() const @@ -4589,12 +4394,10 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const double currTemp = temperature(); double currPres = pressure(); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("\n Debugging information from " "s_Pitzer_dlnMolalityActCoeff_dP()\n"); } -#endif /* * Make sure the counter variables are setup */ @@ -4617,13 +4420,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const */ m_IionicMolality = Is; sqrtIs = sqrt(Is); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 1: \n"); printf(" ionic strenth = %14.7le \n total molar " "charge = %14.7le \n", Is, molarcharge); } -#endif /* * The following call to calc_lambdas() calculates all 16 elements @@ -4638,30 +4439,24 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const * E-thetaprime for all combinations of positive * unlike charges up to 4 */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 2: \n"); } -#endif for (z1 = 1; z1 <=4; z1++) { for (z2 =1; z2 <=4; z2++) { calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n", z1, z2, etheta[z1][z2], etheta_prime[z1][z2]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 3: \n"); printf(" Species Species g(x) " " hfunc(x) \n"); } -#endif /* * @@ -4708,14 +4503,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const gfunc[counterIJ] = 0.0; hfunc[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(), gfunc[counterIJ], hfunc[counterIJ]); } -#endif } } @@ -4724,13 +4517,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const * ------- These are now temperature derivatives of the * previously calculated quantities. */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 4: \n"); printf(" Species Species BMX " "BprimeMX BphiMX \n"); } -#endif for (i = 1; i < m_kk - 1; i++) { for (j = i+1; j < m_kk; j++) { @@ -4747,13 +4538,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const BMX_P[counterIJ] = beta0MX_P[counterIJ] + beta1MX_P[counterIJ] * gfunc[counterIJ] + beta2MX_P[counterIJ] * g2func[counterIJ]; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf("%d %g: %g %g %g %g\n", (int) counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ], beta1MX_P[counterIJ], beta2MX_P[counterIJ], gfunc[counterIJ]); } -#endif if (Is > 1.0E-150) { BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is + beta2MX_P[counterIJ] * h2func[counterIJ]/Is); @@ -4766,27 +4555,23 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const BprimeMX_P[counterIJ] = 0.0; BphiMX_P[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f %11.7f %11.7f \n", sni.c_str(), snj.c_str(), BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]); } -#endif } } /* * --------- SUBSECTION TO CALCULATE CMX_P ---------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 5: \n"); printf(" Species Species CMX \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -4804,14 +4589,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } else { CMX_P[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(), CMX_P[counterIJ]); } -#endif } } @@ -4819,13 +4602,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ---------- * -------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 6: \n"); printf(" Species Species Phi_ij " " Phiprime_ij Phi^phi_ij \n"); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -4848,35 +4629,29 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const Phiprime[counterIJ] = 0.0; Phiphi_P[counterIJ] = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); snj = speciesName(j); printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", sni.c_str(), snj.c_str(), Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ]); } -#endif } } /* * ----------- SUBSECTION FOR CALCULATION OF dFdT --------------------- */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 7: \n"); } -#endif double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres); double dAphidP = dA_DebyedP /3.0; dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" initial value of dFdP = %10.6f \n", dFdP); } -#endif for (i = 1; i < m_kk-1; i++) { for (j = i+1; j < m_kk; j++) { /* @@ -4898,19 +4673,14 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const if (charge(i)*charge(j) > 0) { dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ]; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" dFdP = %10.6f \n", dFdP); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 8: \n"); } -#endif - for (i = 1; i < m_kk; i++) { @@ -5004,15 +4774,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const m_dlnActCoeffMolaldP_Unscaled[i] = zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s lngamma[i]=%10.6f \n", sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqdFdP, sum1, sum2, sum3, sum4, sum5); } -#endif } @@ -5100,15 +4868,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } m_dlnActCoeffMolaldP_Unscaled[i] = zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n", sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]); printf(" %12g %12g %12g %12g %12g %12g\n", zsqdFdP, sum1, sum2, sum3, sum4, sum5); } -#endif } /* @@ -5133,21 +4899,17 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i]; m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { sni = speciesName(i); printf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n", sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]); } -#endif } } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Step 9: \n"); } -#endif /* * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP --------- @@ -5292,19 +5054,14 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const } else { d_osmotic_coef_dP = 0.0; } -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" 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); printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n", sum_m_phi_minus_1, d_osmotic_coef_dP); - } - - if (m_debugCalc) { printf(" Step 10: \n"); } -#endif d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP; @@ -5317,14 +5074,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). */ m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { double d_wateract_dP = exp(d_lnwateract_dP); printf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n", d_lnwateract_dP, d_wateract_dP); } -#endif - } void HMWSoln::calc_lambdas(double is) const @@ -5345,11 +5099,9 @@ void HMWSoln::calc_lambdas(double is) const double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528; aphi = 0.392; /* Value at 25 C */ -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" Is = %g\n", is); } -#endif if (is < 1.0E-150) { for (i = 0; i < 17; i++) { elambda[i] = 0.0; @@ -5383,12 +5135,10 @@ 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; -#ifdef DEBUG_MODE - if (m_debugCalc) { + if (DEBUG_MODE_ENABLED && m_debugCalc) { printf(" ij = %d, elambda = %g, elambda1 = %g\n", ij, elambda[ij], elambda1[ij]); } -#endif } } } @@ -5406,12 +5156,10 @@ void HMWSoln::calc_thetas(int z1, int z2, i = abs(z1); j = abs(z2); -#ifdef DEBUG_MODE - if (i > 4 || j > 4) { + if (DEBUG_MODE_ENABLED && i > 4 || j > 4) { printf("we shouldn't be here\n"); exit(EXIT_FAILURE); } -#endif if ((i == 0) || (j == 0)) { printf("ERROR calc_thetas called with one species being neutral\n"); diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index 1938d7dc3..81c86d826 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -582,16 +582,16 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) { NeutralMolecMoleFractions_[k] = 0.0; } -#ifdef DEBUG_MODE - sum = -1.0; - for (size_t k = 0; k < m_kk; k++) { - sum += moleFractions_[k]; + 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: " + fp2str(sum)); + } } - if (fabs(sum) > 1.0E-11) { - throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", - "molefracts don't sum to one: " + fp2str(sum)); - } -#endif switch (ionSolnType_) { @@ -624,30 +624,30 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij; } -#ifdef DEBUG_MODE - for (size_t k = 0; k < m_kk; k++) { - moleFractionsTmp_[k] = moleFractions_[k]; - } - for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { + if (DEBUG_MODE_ENABLED) { for (size_t k = 0; k < m_kk; k++) { - fmij = fm_neutralMolec_ions_[k + jNeut * m_kk]; - moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut]; + moleFractionsTmp_[k] = moleFractions_[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("", "neutral molecule calc error"); - } else { - //! For the single anion case, we will allow some slippage - if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { - throw CanteraError("", "neutral molecule calc error - anion"); + 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) { + throw CanteraError("", "neutral molecule calc error"); + } else { + //! For the single anion case, we will allow some slippage + if (fabs(moleFractionsTmp_[k]) > 1.0E-5) { + throw CanteraError("", "neutral molecule calc error - anion"); + } } } } } -#endif // Normalize the Neutral Molecule mole fractions sum = 0.0; diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index be71d8e32..9977ccd5e 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -242,15 +242,15 @@ 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. */ -#ifdef DEBUG_MODE - 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"); + 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"); + } } } -#endif strt += nsp; } } diff --git a/src/transport/MMCollisionInt.cpp b/src/transport/MMCollisionInt.cpp index 7a8e296b0..404b36b6d 100644 --- a/src/transport/MMCollisionInt.cpp +++ b/src/transport/MMCollisionInt.cpp @@ -226,23 +226,22 @@ MMCollisionInt::~MMCollisionInt() void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, int log_level) { -#ifdef DEBUG_MODE - if (!xml) { - throw CanteraError("MMCollisionInt::init", "pointer to xml file is zero"); + ostream* logfile = 0; + if (DEBUG_MODE_ENABLED) { + if (!xml) { + throw CanteraError("MMCollisionInt::init", "pointer to xml file is zero"); + } + logfile = &xml->output(); + m_xml = xml; + } else { + m_xml = 0; + log_level = 0; } - ostream& logfile = xml->output(); - m_xml = xml; -#else - m_xml = 0; - log_level = 0; -#endif m_loglevel = log_level; -#ifdef DEBUG_MODE - if (m_loglevel > 0) { - m_xml->XML_comment(logfile, "Collision Integral Polynomial Fits"); + if (DEBUG_MODE_ENABLED && m_loglevel > 0) { + m_xml->XML_comment(*logfile, "Collision Integral Polynomial Fits"); } char p[200]; -#endif m_nmin = -1; m_nmax = -1; @@ -258,29 +257,25 @@ void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, i m_nmin = 0; m_nmax = 36; } -#ifdef DEBUG_MODE - if (m_loglevel > 0) { - m_xml->XML_item(logfile, "Tstar_min", tstar[m_nmin + 1]); - m_xml->XML_item(logfile, "Tstar_max", tstar[m_nmax + 1]); + if (DEBUG_MODE_ENABLED && m_loglevel > 0) { + m_xml->XML_item(*logfile, "Tstar_min", tstar[m_nmin + 1]); + m_xml->XML_item(*logfile, "Tstar_max", tstar[m_nmax + 1]); } -#endif m_logTemp.resize(37); doublereal rmserr, e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0; -#ifdef DEBUG_MODE - if (m_loglevel > 0) { - m_xml->XML_open(logfile, "dstar_fits"); - m_xml->XML_comment(logfile, "Collision integral fits at each " + if (DEBUG_MODE_ENABLED && m_loglevel > 0) { + m_xml->XML_open(*logfile, "dstar_fits"); + m_xml->XML_comment(*logfile, "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." " They are only used for nonzero delta*."); if (log_level < 4) { - m_xml->XML_comment(logfile, + m_xml->XML_comment(*logfile, "polynomial coefficients not printed (log_level < 4)"); } } -#endif string indent = " "; for (int i = 0; i < 37; i++) { @@ -288,15 +283,13 @@ void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, i vector_fp c(DeltaDegree+1); rmserr = fitDelta(0, i, DeltaDegree, DATA_PTR(c)); -#ifdef DEBUG_MODE - if (log_level > 3) { + if (DEBUG_MODE_ENABLED && log_level > 3) { sprintf(p, " Tstar=\"%12.6g\"", tstar[i+1]); - m_xml->XML_open(logfile, "dstar_fit", p); - m_xml->XML_item(logfile, "Tstar", tstar[i+1]); - m_xml->XML_writeVector(logfile, indent, "omega22", + m_xml->XML_open(*logfile, "dstar_fit", p); + m_xml->XML_item(*logfile, "Tstar", tstar[i+1]); + m_xml->XML_writeVector(*logfile, indent, "omega22", c.size(), DATA_PTR(c)); } -#endif m_o22poly.push_back(c); if (rmserr > e22) { e22 = rmserr; @@ -304,54 +297,46 @@ void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, i rmserr = fitDelta(1, i, DeltaDegree, DATA_PTR(c)); m_apoly.push_back(c); -#ifdef DEBUG_MODE - if (log_level > 3) - m_xml->XML_writeVector(logfile, indent, "astar", + if (DEBUG_MODE_ENABLED && log_level > 3) + m_xml->XML_writeVector(*logfile, indent, "astar", c.size(), DATA_PTR(c)); -#endif if (rmserr > ea) { ea = rmserr; } rmserr = fitDelta(2, i, DeltaDegree, DATA_PTR(c)); m_bpoly.push_back(c); -#ifdef DEBUG_MODE - if (log_level > 3) - m_xml->XML_writeVector(logfile, indent, "bstar", + if (DEBUG_MODE_ENABLED && log_level > 3) + m_xml->XML_writeVector(*logfile, indent, "bstar", c.size(), DATA_PTR(c)); -#endif if (rmserr > eb) { eb = rmserr; } rmserr = fitDelta(3, i, DeltaDegree, DATA_PTR(c)); m_cpoly.push_back(c); -#ifdef DEBUG_MODE - if (log_level > 3) { - m_xml->XML_writeVector(logfile, indent, "cstar", + if (DEBUG_MODE_ENABLED && log_level > 3) { + m_xml->XML_writeVector(*logfile, indent, "cstar", c.size(), DATA_PTR(c)); } -#endif if (rmserr > ec) { ec = rmserr; } -#ifdef DEBUG_MODE - if (log_level > 3) { - m_xml->XML_close(logfile, "dstar_fit"); + if (DEBUG_MODE_ENABLED && log_level > 3) { + m_xml->XML_close(*logfile, "dstar_fit"); } - if (log_level > 0) { + if (DEBUG_MODE_ENABLED && log_level > 0) { sprintf(p, "max RMS errors in fits vs. delta*:\n" " omega_22 = %12.6g \n" " A* = %12.6g \n" " B* = %12.6g \n" " C* = %12.6g \n", e22, ea, eb, ec); - m_xml->XML_comment(logfile, p); - m_xml->XML_close(logfile, "dstar_fits"); + m_xml->XML_comment(*logfile, p); + m_xml->XML_close(*logfile, "dstar_fits"); } -#endif } } diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index b133d96ba..380d237be 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -507,31 +507,23 @@ void TransportFactory::setupMM(std::ostream& flog, const std::vectorXML_open(flog, "collision_integrals"); } -#endif MMCollisionInt integrals; integrals.init(tr.xml, tstar_min, tstar_max, log_level); fitCollisionIntegrals(flog, tr, integrals); -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_close(flog, "collision_integrals"); } -#endif // make polynomial fits -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_open(flog, "property fits"); } -#endif fitProperties(tr, integrals, flog); -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_close(flog, "property fits"); } -#endif } void TransportFactory::setupLiquidTransport(std::ostream& flog, thermo_t* thermo, int log_level, @@ -717,8 +709,7 @@ void TransportFactory::fitCollisionIntegrals(ostream& logfile, // Chemkin fits to sixth order polynomials int degree = (mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE); -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_open(logfile, "tstar_fits"); tr.xml->XML_comment(logfile, "fits to A*, B*, and C* vs. log(T*).\n" "These are done only for the required dstar(j,k) values."); @@ -726,7 +717,6 @@ void TransportFactory::fitCollisionIntegrals(ostream& logfile, tr.xml->XML_comment(logfile, "*** polynomial coefficients not printed (log_level < 3) ***"); } } -#endif for (i = 0; i < nsp; i++) { for (j = i; j < nsp; j++) { // Chemkin fits only delta* = 0 @@ -766,11 +756,9 @@ void TransportFactory::fitCollisionIntegrals(ostream& logfile, tr.poly[j][i] = tr.poly[i][j]; } } -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_close(logfile, "tstar_fits"); } -#endif } void TransportFactory::getTransportData(const std::vector &xspecies, @@ -1178,9 +1166,6 @@ void TransportFactory::fitProperties(GasTransportParams& tr, { doublereal tstar; int ndeg = 0; -#ifdef DEBUG_MODE - char s[100]; -#endif // number of points to use in generating fit data const size_t np = 50; @@ -1206,29 +1191,26 @@ void TransportFactory::fitProperties(GasTransportParams& tr, // fit the pure-species viscosity and thermal conductivity for // each species -#ifdef DEBUG_MODE - if (tr.log_level < 2 && m_verbose) { + if (DEBUG_MODE_ENABLED && tr.log_level < 2 && m_verbose) { tr.xml->XML_comment(logfile, "*** polynomial coefficients not printed (log_level < 2) ***"); } -#endif doublereal sqrt_T, visc, err, relerr, mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0; -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { tr.xml->XML_open(logfile, "viscosity"); tr.xml->XML_comment(logfile,"Polynomial fits for viscosity"); if (mode == CK_Mode) { tr.xml->XML_comment(logfile,"log(viscosity) fit to cubic " "polynomial in log(T)"); } else { + char s[100]; sprintf(s, "viscosity/sqrt(T) fit to " "polynomial of degree %d in log(T)",degree); tr.xml->XML_comment(logfile,s); } } -#endif doublereal cp_R, cond, w_RT, f_int, A_factor, B_factor, c1, cv_rot, cv_int, f_rot, f_trans, om11; @@ -1348,15 +1330,13 @@ void TransportFactory::fitProperties(GasTransportParams& tr, tr.visccoeffs.push_back(c); tr.condcoeffs.push_back(c2); -#ifdef DEBUG_MODE - if (tr.log_level >= 2 && m_verbose) { + if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) { tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k), c.size(), DATA_PTR(c)); } -#endif } -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { + char s[100]; sprintf(s, "Maximum viscosity absolute error: %12.6g", mxerr); tr.xml->XML_comment(logfile,s); sprintf(s, "Maximum viscosity relative error: %12.6g", mxrelerr); @@ -1398,7 +1378,6 @@ void TransportFactory::fitProperties(GasTransportParams& tr, tr.xml->XML_comment(logfile,s); } } -#endif mxerr = 0.0, mxrelerr = 0.0; vector_fp diff(np + 1); @@ -1457,16 +1436,14 @@ void TransportFactory::fitProperties(GasTransportParams& tr, } } tr.diffcoeffs.push_back(c); -#ifdef DEBUG_MODE - if (tr.log_level >= 2 && m_verbose) { + if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) { tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k) + "__"+tr.thermo->speciesName(j), c.size(), DATA_PTR(c)); } -#endif } } -#ifdef DEBUG_MODE - if (m_verbose) { + if (DEBUG_MODE_ENABLED && m_verbose) { + char s[100]; sprintf(s,"Maximum binary diffusion coefficient absolute error:" " %12.6g", mxerr); tr.xml->XML_comment(logfile,s); @@ -1475,7 +1452,6 @@ void TransportFactory::fitProperties(GasTransportParams& tr, tr.xml->XML_comment(logfile,s); tr.xml->XML_close(logfile, "binary_diffusion_coefficients"); } -#endif } Transport* newTransportMgr(const std::string& transportModel, thermo_t* thermo, int loglevel, TransportFactory* f, int ndim)