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

Eliminate debug_verbose option and DEBUG_MODE macro

Browse source 
Vebose printing can always be enabled at runtime by setting flags on various
classes.
This commit is contained in:
Ray Speth 2015-11-18 11:02:30 -05:00
parent ee95c60813
commit e58bd09859
31 changed files with 483 additions and 726 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
SConstruct
View file

@ -390,11 +390,6 @@ config_options = [
('FORTRANFLAGS',
'Compilation options for the Fortran (90) compiler.',
'-O3'),
BoolVariable(
'debug_verbose',
"""Enable extra printing to aid in debugging. This code is marked
by the preprocessor macros DEBUG_MODE and DEBUG_MODE_ENABLED.""",
False),
BoolVariable(
'coverage',
"""Enable collection of code coverage information with gcov.
@ -1201,8 +1196,6 @@ def cdefine(definevar, configvar, comp=True, value=1):
else:
configh[definevar] = None
cdefine('DEBUG_MODE', 'debug_verbose')
# Need to test all of these to see what platform.system() returns
configh['SOLARIS'] = 1 if env['OS'] == 'Solaris' else None
configh['DARWIN'] = 1 if env['OS'] == 'Darwin' else None

5
doc/sphinx/scons-options.txt
View file

@ -132,11 +132,6 @@ running 'scons build'. The format of this file is:
Compilation options for the Fortran (90) compiler.
- default: '-O3'
* debug_verbose: [ yes | no ]
Enable extra printing to aid in debugging. This code is marked by
the preprocessor macros DEBUG_MODE and DEBUG_MODE_ENABLED.
- default: 'no'
* coverage: [ yes | no ]
Enable collection of code coverage information with gcov. Available
only when compiling with gcc.

13
include/cantera/base/config.h.in
View file

@ -8,19 +8,6 @@
// Just the major + minor version (i.e. 2.2 instead of 2.2.0)
%(CANTERA_SHORT_VERSION)s
//------------------------ Development flags ------------------//
//
// Compile in additional debug printing where available.
// Note, the printing may need to be turned on via a switch.
// This just compiles in the code.
%(DEBUG_MODE)s
#ifdef DEBUG_MODE
#define DEBUG_MODE_ENABLED 1
#else
#define DEBUG_MODE_ENABLED 0
#endif
//------------------------ Fortran settings -------------------//
// define types doublereal, integer, and ftnlen to match the

8
include/cantera/equil/vcs_solve.h
View file

@ -509,9 +509,8 @@ public:
* flags with total moles in each phase.
*/
double vcs_tmoles();
#ifdef DEBUG_MODE
void check_tmoles() const;
#endif
//! This subroutine calculates reaction free energy changes for
//! all noncomponent formation reactions.
@ -1291,11 +1290,8 @@ private:
*/
double l2normdg(double dg[]) const;
#ifdef DEBUG_MODE
//! Print out and check the elemental abundance vector
void prneav() const;
#endif
void checkDelta1(double* const ds, double* const delTPhMoles, size_t kspec);
@ -1895,8 +1891,6 @@ public:
* * 5 each iteration in solve_TP gets a report with one line per species
* * 6 Each decision in solve_TP gets a line per species in addition to 4
* * 10 Additionally Hessian matrix is printed out
*
* Levels of printing above 4 are only accessible when DEBUG_MODE is turned on
*/
int m_debug_print_lvl;

8
include/cantera/thermo/HMWSoln.h
View file

@ -2789,16 +2789,10 @@ private:
static int interp_est(const std::string& estString);
public:
/*!
* Turn on copious debug printing when this
* is true and DEBUG_MODE is turned on.
*/
//! Turn on copious debug printing when this is true
mutable int m_debugCalc;
//! Return int specifying the amount of debug printing
/*!
* This will return 0 if DEBUG_MODE is not turned on
*/
int debugPrinting();
};

4
include/cantera/thermo/PDSS_HKFT.h
View file

@ -86,7 +86,6 @@ public:
virtual doublereal enthalpy_mole() const;
#ifdef DEBUG_MODE
//! Return the molar enthalpy in units of J kmol-1
/*!
* Returns the species standard state enthalpy in J kmol-1 at the
@ -97,7 +96,6 @@ public:
* @return returns the species standard state enthalpy in J kmol-1
*/
doublereal enthalpy_mole2() const;
#endif
virtual doublereal intEnergy_mole() const;
virtual doublereal entropy_mole() const;
@ -221,14 +219,12 @@ private:
*/
doublereal deltaS() const;
#ifdef DEBUG_MODE
//! Routine that actually calculates the enthalpy difference
//! between the reference state at Tr, Pr and T,P
/*!
* This is an extra routine that was added to check the arithmetic
*/
doublereal deltaH() const;
#endif
//! Internal formula for the calculation of a_g()
/*!

2
include/cantera/thermo/RedlichKisterVPSSTP.h
View file

@ -432,7 +432,6 @@ private:
*/
void s_update_dlnActCoeff_dX_() const;
#ifdef DEBUG_MODE
public:
//! Utility routine that calculates a literature expression
/*!
@ -442,7 +441,6 @@ public:
* nonideal term and mf term
*/
void Vint(double& VintOut, double& voltsOut);
#endif
protected:
//! number of binary interaction expressions

10
src/base/ctml.cpp
View file

@ -214,18 +214,8 @@ doublereal getFloatCurrent(const XML_Node& node, const std::string& type)
fctr = toSI(units);
} else if (type != "" && units != "") {
fctr = toSI(units);
#ifdef DEBUG_MODE
writelog("\nWarning: conversion toSI() was done on node value " + node.name() +
"but wasn't explicitly requested. Type was \"" + type + "\"\n");
#endif
#ifdef DEBUG_MODE_MORE
} else if (type == "" && units != "") {
// Note, below currently produces a lot of output due to transport blocks.
// This needs to be addressed.
writelog("\nWarning: XML node " + node.name() +
"has a units attribute, \"" + units + "\","
"but no conversion was done because the getFloat() command didn't have a type\n");
#endif
}
return fctr*x;
}

18
src/equil/BasisOptimize.cpp
View file

@ -56,7 +56,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
}
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
if (BasisOptimize_print_lvl >= 1) {
writelog(" ");
for (i=0; i<77; i++) {
writelog("-");
@ -118,9 +118,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// For debugging purposes keep an unmodified copy of the array.
vector_fp molNumBase;
if (DEBUG_MODE_ENABLED) {
molNumBase = molNum;
}
molNumBase = molNum;
double molSave = 0.0;
size_t jr = 0;
@ -157,9 +155,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// Assign a small negative number to the component that we have
// just found, in order to take it out of further consideration.
#ifdef DEBUG_MODE
molSave = molNum[kk];
#endif
molNum[kk] = USEDBEFORE;
// CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES
@ -208,7 +204,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
if (BasisOptimize_print_lvl >= 1) {
kk = orderVectorSpecies[k];
writelogf(" --- %-12.12s", mphase->speciesName(kk));
jj = orderVectorSpecies[jr];
@ -283,7 +279,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne,
&ipiv[0], &formRxnMatrix[0], ne, info);
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
if (BasisOptimize_print_lvl >= 1) {
writelog(" ---\n");
writelogf(" --- Number of Components = %d\n", nComponents);
writelog(" --- Formula Matrix:\n");
@ -379,7 +375,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
size_t nspecies = mphase->nSpecies();
double test = -1.0E10;
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
if (BasisOptimize_print_lvl > 0) {
writelog(" ");
for (i=0; i<77; i++) {
writelog("-");
@ -456,7 +452,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
if (k == nelements) {
// When we are here, there is an error usually.
// We haven't found the number of elements necessary.
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
if (BasisOptimize_print_lvl > 0) {
writelogf("Error exit: returning with nComponents = %d\n", jr);
}
throw CanteraError("ElemRearrange", "Required number of elements not found.");
@ -517,7 +513,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
}
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
if (BasisOptimize_print_lvl > 0) {
kk = orderVectorElements[k];
writelog(" --- ");
writelogf("%-2.2s", mphase->elementName(kk));

54
src/equil/ChemEquil.cpp
View file

@ -204,7 +204,7 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal,
// and element abundance vectors kept within the ChemEquil object.
update(s);
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog("setInitialMoles: Estimated Mole Fractions\n");
writelogf(" Temperature = %g\n", s.temperature());
writelogf(" Pressure = %g\n", s.pressure());
@ -259,7 +259,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT,
doublereal rrt = 1.0/(GasConstant* s.temperature());
scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt);
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
for (size_t m = 0; m < m_nComponents; m++) {
size_t isp = m_component[m];
writelogf("isp = %d, %s\n", isp, s.speciesName(isp));
@ -293,7 +293,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT,
lambda_RT[m_orderVectorElements[m]] = 0.0;
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog(" id CompSpecies ChemPot EstChemPot Diff\n");
for (size_t m = 0; m < m_nComponents; m++) {
size_t isp = m_component[m];
@ -639,7 +639,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// Compute the Jacobian matrix
equilJacobian(s, x, elMolesGoal, jac, xval, yval);
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf("Jacobian matrix %d:\n", iter);
for (m = 0; m <= m_mm; m++) {
writelog(" [ ");
@ -698,7 +698,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
fctr = std::min(fctr, 0.2/fabs(res_trial[mm]));
}
}
if (fctr != 1.0 && DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (fctr != 1.0 && ChemEquil_print_lvl > 0) {
writelogf("WARNING Soln Damping because of bounds: %g\n", fctr);
}
@ -756,7 +756,7 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
for (size_t m = 0; m < x.size(); m++) {
x[m] = oldx[m] + damp * step[m];
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf("Solution Unknowns: damp = %g\n", damp);
writelog(" X_new X_old Step\n");
for (size_t m = 0; m < m_mm; m++) {
@ -794,7 +794,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !m_doResPerturb) {
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
writelog("Residual: ElFracGoal ElFracCurrent Resid\n");
for (size_t n = 0; n < m_mm; n++) {
writelogf(" % -14.7E % -14.7E % -10.5E\n",
@ -807,7 +807,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
resid[m_mm] = xx/xval - 1.0;
resid[m_skip] = yy/yval - 1.0;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !m_doResPerturb) {
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
writelog(" Goal Xvalue Resid\n");
writelogf(" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]);
writelogf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
@ -957,7 +957,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog("estimateEP_Brinkley::\n\n");
double temp = s.temperature();
double pres = s.pressure();
@ -991,7 +991,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
x_old[m_mm] = n_t;
// Calculate the mole numbers of species
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf("START ITERATION %d:\n", iter);
}
// Calculate the mole numbers of species and elements.
@ -1002,7 +1002,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n");
for (k = 0; k < m_kk; k++) {
writelogf("%15s: %10.5g %10.5g\n",
@ -1032,7 +1032,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !normalStep) {
if (ChemEquil_print_lvl > 0 && !normalStep) {
writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
}
if (!normalStep) {
@ -1087,7 +1087,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
nCutoff = 1.0E-9 * n_t_calc;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog(" Lump Sum Elements Calculation: \n");
}
for (m = 0; m < m_mm; m++) {
@ -1111,7 +1111,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf(" %5s %3d : %5d %5d\n",
s.elementName(m), lumpSum[m], kMSp, kMSp2);
}
@ -1164,7 +1164,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
for (m = 0; m < m_mm; m++) {
if (a1(m,m) < 1.0E-50) {
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m);
}
for (n = 0; n < m_mm; n++) {
@ -1183,7 +1183,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
resid[m_mm] = n_t - n_t_calc;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog("Matrix:\n");
for (m = 0; m <= m_mm; m++) {
writelog(" [");
@ -1196,7 +1196,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
tmp = resid[m_mm] /(n_t + 1.0E-15);
sum += tmp * tmp;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf("(it %d) Convergence = %g\n", iter, sum);
}
@ -1216,7 +1216,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
tmp += fabs(a1(m,n));
}
if (m < m_mm && tmp < 1.0E-30) {
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf(" NOTE: Diagonalizing row %d\n", m);
}
for (n = 0; n <= m_mm; n++) {
@ -1233,7 +1233,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
resid[m] *= tmp;
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelog("Row Summed Matrix:\n");
for (m = 0; m <= m_mm; m++) {
writelog(" [");
@ -1279,7 +1279,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
if (sameAsRow != npos || lumpSum[m]) {
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
if (lumpSum[m]) {
writelogf("Lump summing row %d, due to rank deficiency analysis\n", m);
} else if (sameAsRow != npos) {
@ -1296,7 +1296,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && modifiedMatrix) {
if (ChemEquil_print_lvl > 0 && modifiedMatrix) {
writelog("Row Summed, MODIFIED Matrix:\n");
for (m = 0; m <= m_mm; m++) {
writelog(" [");
@ -1328,7 +1328,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
beta = std::min(beta, -1.0 / resid[m]);
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && beta != 1.0) {
if (ChemEquil_print_lvl > 0 && beta != 1.0) {
writelogf("(it %d) Beta = %g\n", iter, beta);
}
}
@ -1338,7 +1338,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
n_t *= exp(beta * resid[m_mm]);
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
for (m = 0; m < m_mm; m++) {
writelogf(" %5s %10.5g %10.5g %10.5g\n",
@ -1347,7 +1347,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
writelogf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm]));
}
}
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
double temp = s.temperature();
double pres = s.pressure();
@ -1377,7 +1377,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
size_t maxNegEloc = npos;
double maxPosVal = -1.0;
double maxNegVal = -1.0;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
if (ChemEquil_print_lvl > 0) {
for (k = 0; k < m_kk; k++) {
if (nAtoms(k,m_eloc) > 0.0 && m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) {
maxPosVal = m_molefractions[k];
@ -1407,7 +1407,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
return;
}
double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) {
if (ChemEquil_print_lvl > 0 && factor < 0.9999999999) {
writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
s.speciesName(maxPosEloc),
m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor);
@ -1419,7 +1419,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
}
} else {
double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg;
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && factor < 0.9999999999) {
if (ChemEquil_print_lvl > 0 && factor < 0.9999999999) {
writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
s.speciesName(maxNegEloc),
m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor);

4
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -579,9 +579,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
double Temp = m_mix->temperature();
double pres = m_mix->pressure();
vector_fp& mf = m_vprob.mf;
#ifdef DEBUG_MODE
double* fe = &m_vprob.m_gibbsSpecies[0];
#endif
vector_fp VolPM;
vector_fp activity;
vector_fp ac;
@ -692,7 +690,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
}
}
#ifdef DEBUG_MODE
// Check consistency: These should be equal
tref.getChemPotentials(fe+istart);
for (size_t k = 0; k < nSpecies; k++) {
@ -702,7 +699,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "incompatibility!");
}
}
#endif
}
fclose(FP);

19
src/equil/vcs_VolPhase.cpp
View file

@ -430,9 +430,10 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_old[0];
} else if (stateCalc == VCS_STATECALC_NEW) {
molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_new[0];
} else if (DEBUG_MODE_ENABLED) {
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here"); }
} else if (DEBUG_MODE_ENABLED && m_owningSolverObject) {
} else {
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
}
} else if (m_owningSolverObject) {
if (stateCalc == VCS_STATECALC_OLD) {
if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_old[0]) {
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
@ -890,20 +891,16 @@ void vcs_VolPhase::setExistence(const int existence)
{
if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
if (v_totalMoles != 0.0) {
if (DEBUG_MODE_ENABLED) {
throw CanteraError("vcs_VolPhase::setExistence",
"setting false existence for phase with moles");
} else {
v_totalMoles = 0.0;
}
throw CanteraError("vcs_VolPhase::setExistence",
"setting false existence for phase with moles");
}
} else if (DEBUG_MODE_ENABLED && m_totalMolesInert == 0.0) {
} else if (m_totalMolesInert == 0.0) {
if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) {
throw CanteraError("vcs_VolPhase::setExistence",
"setting true existence for phase with no moles");
}
}
if (DEBUG_MODE_ENABLED && m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) {
if (m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) {
throw CanteraError("vcs_VolPhase::setExistence",
"Trying to set existence of an electron phase to false");
}

8
src/equil/vcs_elem.cpp
View file

@ -94,7 +94,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
{
int retn = 0;
#ifdef DEBUG_MODE
vector_fp ga_save(m_elemAbundances);
if (m_debug_print_lvl >= 2) {
plogf(" --- vcsc_elcorr: Element abundances correction routine");
@ -112,7 +111,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
l2before += x[i] * x[i];
}
l2before = sqrt(l2before/m_numElemConstraints);
#endif
// Special section to take out single species, single component,
// moles. These are species which have non-zero entries in the
@ -190,7 +188,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[i] / atomComp;
if (m_molNumSpecies_old[kspec] > maxPermissible) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) {
if (m_debug_print_lvl >= 3) {
plogf(" --- vcs_elcorr: Reduced species %s from %g to %g "
"due to %s max bounds constraint\n",
m_speciesName[kspec], m_molNumSpecies_old[kspec],
@ -205,7 +203,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
} else {
m_speciesStatus[kspec] = VCS_SPECIES_ACTIVEBUTZERO;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_elcorr: Zeroed species %s and changed "
"status to %d due to max bounds constraint\n",
m_speciesName[kspec], m_speciesStatus[kspec]);
@ -422,7 +420,6 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
L_CLEANUP:
;
vcs_tmoles();
#ifdef DEBUG_MODE
double l2after = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) {
l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2);
@ -439,7 +436,6 @@ L_CLEANUP:
plogf(" --- Diff_Norm: %20.12E %20.12E\n",
l2before, l2after);
}
#endif
return retn;
}

4
src/equil/vcs_elem_rearrange.cpp
View file

@ -21,7 +21,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
double* const sm, double* const ss)
{
size_t ncomponents = m_numComponents;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ");
for (size_t i=0; i<77; i++) {
plogf("-");
@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
}
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-2.2s", m_elementName[k]);
plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]);

22
src/equil/vcs_inest.cpp
View file

@ -28,7 +28,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
// molNum(I) .GE. 0.0. Note, both of these programs do this.
vcs_setMolesLinProg();
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
pprefix);
plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
@ -116,7 +116,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
}
}
vcs_deltag(0, true, VCS_STATECALC_NEW);
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s", pprefix);
plogf("%-12.12s", m_speciesName[kspec]);
@ -171,7 +171,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
}
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf("%sdirection (", pprefix);
@ -269,7 +269,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
lt = 1;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%s Final Mole Numbers produced by inest:\n",
pprefix);
plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
@ -290,14 +290,14 @@ int VCS_SOLVE::vcs_inest_TP()
// Calculate the elemental abundances
vcs_elab();
if (vcs_elabcheck(0)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%s Initial guess passed element abundances on input\n", pprefix);
plogf("%s m_doEstimateEquil = 1 so will use the input mole "
"numbers as estimates", pprefix);
plogendl();
}
return retn;
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf("%s Initial guess failed element abundances on input\n", pprefix);
plogf("%s m_doEstimateEquil = 1 so will discard input "
"mole numbers and find our own estimate", pprefix);
@ -312,7 +312,7 @@ int VCS_SOLVE::vcs_inest_TP()
vector_fp aw(m_numSpeciesTot+ m_numElemConstraints, 0.0);
// Go get the estimate of the solution
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%sGo find an initial estimate for the equilibrium problem",
pprefix);
plogendl();
@ -332,7 +332,7 @@ int VCS_SOLVE::vcs_inest_TP()
// matrix.
bool rangeCheck = vcs_elabcheck(1);
if (!vcs_elabcheck(0)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%sInitial guess failed element abundances\n", pprefix);
plogf("%sCall vcs_elcorr to attempt fix", pprefix);
plogendl();
@ -347,7 +347,7 @@ int VCS_SOLVE::vcs_inest_TP()
plogendl();
retn = -1;
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (rangeCheck) {
plogf("%sInitial guess now satisfies element abundances", pprefix);
plogendl();
@ -361,7 +361,7 @@ int VCS_SOLVE::vcs_inest_TP()
}
}
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (rangeCheck) {
plogf("%sInitial guess satisfies element abundances", pprefix);
plogendl();
@ -375,7 +375,7 @@ int VCS_SOLVE::vcs_inest_TP()
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_new[0],
&m_tPhaseMoles_old[0]));

4
src/equil/vcs_nondim.cpp
View file

@ -105,7 +105,7 @@ void VCS_SOLVE::vcs_nondim_TP()
if (m_totalMoleScale != 1.0) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale);
plogendl();
}
@ -149,7 +149,7 @@ void VCS_SOLVE::vcs_redim_TP()
}
if (m_totalMoleScale != 1.0) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale);
plogendl();
}

40
src/equil/vcs_phaseStability.cpp
View file

@ -205,24 +205,18 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
doublereal FephaseMax = -1.0E30;
doublereal Fephase = -1.0E30;
#ifdef DEBUG_MODE
char anote[128];
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_popPhaseID() called\n");
plogf(" --- Phase Status F_e MoleNum\n");
plogf(" --------------------------------------------------------------------------\n");
}
#else
char* anote;
#endif
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (DEBUG_MODE_ENABLED) {
strcpy(anote, "");
}
strcpy(anote, "");
if (existence > 0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d NA %11.3e\n",
Vphase->PhaseName, existence, m_tPhaseMoles_old[iph]);
}
@ -234,24 +228,20 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
doublereal deltaGRxn = m_deltaGRxn_old[irxn];
Fephase = exp(-deltaGRxn) - 1.0;
if (Fephase > 0.0) {
if (DEBUG_MODE_ENABLED) {
strcpy(anote," (ready to be birthed)");
}
strcpy(anote," (ready to be birthed)");
if (Fephase > FephaseMax) {
iphasePop = iph;
FephaseMax = Fephase;
if (DEBUG_MODE_ENABLED) {
strcpy(anote," (chosen to be birthed)");
}
strcpy(anote," (chosen to be birthed)");
}
}
if (DEBUG_MODE_ENABLED && Fephase < 0.0) {
if (Fephase < 0.0) {
strcpy(anote," (not stable)");
AssertThrowMsg(m_tPhaseMoles_old[iph] <= 0.0,
"VCS_SOLVE::vcs_popPhaseID", "shouldn't be here");
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d %10.3g %10.3g %s\n",
Vphase->PhaseName, existence, Fephase,
m_tPhaseMoles_old[iph], anote);
@ -268,13 +258,13 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
} else {
FephaseMax = std::max(FephaseMax, Fephase);
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d %11.3g %11.3g\n",
Vphase->PhaseName, existence, Fephase,
m_tPhaseMoles_old[iph]);
}
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %18s %5d blocked %11.3g\n",
Vphase->PhaseName,
existence, m_tPhaseMoles_old[iph]);
@ -290,7 +280,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
// Insert logic here to figure out if phase pops are linked together. Only
// do one linked pop at a time.
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ---------------------------------------------------------------------\n");
}
return iphasePop;
@ -312,7 +302,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes",
"called for a phase that exists!");
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
Vphase->PhaseName, iphasePop);
}
@ -494,13 +484,13 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
if (doSuccessiveSubstitution) {
int KP = 0;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_phaseStabilityTest() called\n");
plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
" normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
plogf(" --------------------------------------------------------------"
"--------------------------------------------------------\n");
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) {
} else if (m_debug_print_lvl == 1) {
plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
}
@ -679,7 +669,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
fracDelta_new[k] = fracDelta_old[k] + damp * delFrac[k];
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
}
@ -703,10 +693,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
} else {
throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet");
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ------------------------------------------------------------"
"-------------------------------------------------------------\n");
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) {
} else if (m_debug_print_lvl == 1) {
if (funcPhaseStability > 0.0) {
plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
} else {

10
src/equil/vcs_prob.cpp
View file

@ -460,16 +460,6 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
}
}
if (DEBUG_MODE_ENABLED) {
// Check consistency: These should be equal
tp->getChemPotentials(&m_gibbsSpecies[0]+istart);
for (size_t k = 0; k < nSpeciesPhase; k++) {
if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) {
fclose(FP);
throw CanteraError("VCS_PROB::reportCSV", "incompatibility");
}
}
}
iK += nSpeciesPhase;
}
fclose(FP);

225
src/equil/vcs_rxnadj.cpp
View file

@ -21,7 +21,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
{
size_t iphDel = npos;
size_t k = 0;
#ifdef DEBUG_MODE
char ANOTE[128];
if (m_debug_print_lvl >= 2) {
plogf(" ");
@ -38,9 +37,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
" | Comment\n");
}
#else
char* ANOTE = 0;
#endif
// We update the matrix dlnActCoeffdmolNumber[][] at the top of the loop,
// when necessary
@ -50,16 +46,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
// LOOP OVER THE FORMATION REACTIONS
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Normal Calc");
}
sprintf(ANOTE, "Normal Calc");
size_t kspec = m_indexRxnToSpecies[irxn];
if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) {
m_deltaMolNumSpecies[kspec] = 0.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed");
}
sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed");
} else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
// MULTISPECIES PHASE WITH total moles equal to zero
@ -78,37 +70,27 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
m_deltaMolNumSpecies[kspec] = 0.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
} else {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
}
}
} else {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E",
sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
}
} else {
sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
m_deltaMolNumSpecies[kspec] = 0.0;
if (Vphase->exists() > 0 && trphmoles > 0.0) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE,
"MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE,
"MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
}
}
} else {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
m_deltaMolNumSpecies[kspec] = 0.0;
}
} else {
@ -118,14 +100,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
// bother if superconvergence has already been achieved in this
// mode.
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
continue;
}
@ -133,14 +113,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
// Don't calculate for minor or nonexistent species if their
// values are to be decreasing anyway.
if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
continue;
}
@ -170,10 +148,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
if (m_useActCoeffJac) {
double s_old = s;
s = vcs_Hessian_diag_adj(irxn, s_old);
if (DEBUG_MODE_ENABLED && s_old != s) {
sprintf(ANOTE, "Normal calc: diag adjusted from %g "
"to %g due to act coeff", s_old, s);
}
sprintf(ANOTE, "Normal calc: diag adjusted from %g "
"to %g due to act coeff", s_old, s);
}
m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
@ -184,18 +160,14 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
double negChangeComp = -stoicC * m_deltaMolNumSpecies[kspec];
if (negChangeComp > m_molNumSpecies_old[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
}
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[j] / stoicC;
} else {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
}
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
m_deltaMolNumSpecies[kspec] = 0.0;
}
}
@ -204,11 +176,9 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
// Implement a damping term that limits m_deltaMolNumSpecies
// to the size of the mole number
if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Delta damped from %g "
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec],
m_speciesName[kspec].c_str());
}
sprintf(ANOTE, "Delta damped from %g "
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec],
m_speciesName[kspec].c_str());
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
}
} else {
@ -261,12 +231,10 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
// having all of the mole numbers of that phases. it
// seems that we can suggest a zero of the species
// and the code will recover.
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
}
sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
@ -287,31 +255,29 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
iphDel = m_phaseID[k];
kSpecial = k;
if (DEBUG_MODE_ENABLED) {
if (k != kspec) {
sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel,
m_speciesName[k].c_str());
} else {
sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel);
}
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s",
m_speciesName[k]);
plogendl();
}
if (k != kspec) {
sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel,
m_speciesName[k].c_str());
} else {
sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel);
}
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s",
m_speciesName[k]);
plogendl();
}
if (k != kspec) {
forceComponentCalc = 1;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Force a component recalculation \n");
plogendl();
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ");
writeline('-', 82);
}
@ -319,7 +285,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
}
} // End of regular processing
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
@ -327,7 +293,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
} // End of loop over m_speciesUnknownType
} // End of loop over non-component stoichiometric formation reactions
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ");
writeline('-', 82);
}
@ -337,7 +303,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
int VCS_SOLVE::vcs_rxn_adj_cg()
{
int soldel = 0;
#ifdef DEBUG_MODE
char ANOTE[128];
plogf(" ");
for (size_t j = 0; j < 77; j++) {
@ -345,17 +310,12 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
}
plogf("\n --- Subroutine rxn_adj_cg() called\n");
plogf(" --- Species Moles Rxn_Adjustment | Comment\n");
#else
char* ANOTE = 0;
#endif
// Precalculation loop -> we calculate quantities based on loops over the
// number of species. We also evaluate whether the matrix is appropriate for
// this algorithm. If not, we bail out.
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Normal Calc");
}
sprintf(ANOTE, "Normal Calc");
size_t kspec = m_indexRxnToSpecies[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
@ -364,16 +324,12 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
// HKM -> the statment below presupposes units in m_deltaGRxn_new[].
// It probably should be replaced with something more relativistic
if (m_deltaGRxn_new[irxn] < -1.0e-4) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]);
m_deltaMolNumSpecies[kspec] = 1.0e-10;
m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
--m_numRxnMinorZeroed;
} else {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]);
m_deltaMolNumSpecies[kspec] = 0.0;
}
} else {
@ -382,26 +338,22 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
// First take care of cases where we want to bail out. Don't bother
// if superconvergence has already been achieved in this mode.
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
}
sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
continue;
}
// Don't calculate for minor or nonexistent species if their values
// are to be decreasing anyway.
if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
}
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
continue;
}
@ -474,11 +426,9 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
}
m_molNumSpecies_old[k] = 0.0;
m_tPhaseMoles_old[m_phaseID[k]] = 0.0;
if (DEBUG_MODE_ENABLED) {
plogf(" --- vcs_st2 Special section to delete ");
plogf("%-12.12s", m_speciesName[k]);
plogf("\n --- Immediate return - Restart iteration\n");
}
plogf(" --- vcs_st2 Special section to delete ");
plogf("%-12.12s", m_speciesName[k]);
plogf("\n --- Immediate return - Restart iteration\n");
// We need to immediately recompute the component basis,
// because we just zeroed it out.
@ -491,25 +441,21 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
}
}
} // End of regular processing
if (DEBUG_MODE_ENABLED) {
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
}
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
} // End of loop over non-component stoichiometric formation reactions
// When we form the Hessian we must be careful to ensure that it is a
// symmetric positive definite matrix, still. This means zeroing out columns
// when we zero out rows as well. I suggest writing a small program to make
// sure of this property.
if (DEBUG_MODE_ENABLED) {
plogf(" ");
for (size_t j = 0; j < 77; j++) {
plogf("-");
}
plogf("\n");
plogf(" ");
for (size_t j = 0; j < 77; j++) {
plogf("-");
}
plogf("\n");
return soldel;
}
@ -613,7 +559,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char*
if (deltaGOrig > 0.0) {
if (dx_orig > 0.0) {
dx = 0.0;
if (DEBUG_MODE_ENABLED && ANOTE) {
if (ANOTE) {
sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx>0 dg > 0");
}
return dx;
@ -621,7 +567,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char*
} else if (deltaGOrig < 0.0) {
if (dx_orig < 0.0) {
dx = 0.0;
if (DEBUG_MODE_ENABLED && ANOTE) {
if (ANOTE) {
sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx<0 dg < 0");
}
return dx;
@ -699,17 +645,12 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char*
finalize:
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
if (its >= MAXITS) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Rxn reduced to zero step size from %g to %g (MAXITS)", dx_orig, dx);
return dx;
#endif
}
#ifdef DEBUG_MODE
if (dx != dx_orig) {
sprintf(ANOTE, "Line Search reduced step size from %g to %g", dx_orig, dx);
}
#endif
return dx;
}

14
src/equil/vcs_setMolesLinProg.cpp
View file

@ -34,7 +34,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
size_t ik, irxn;
double test = -1.0E-10;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- call setInitialMoles\n");
}
@ -59,13 +59,13 @@ int VCS_SOLVE::vcs_setMolesLinProg()
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
}
while (redo) {
if (!vcs_elabcheck(0)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- seMolesLinProg Mole numbers failing element abundances\n");
plogf(" --- seMolesLinProg Call vcs_elcorr to attempt fix\n");
}
@ -88,7 +88,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
return retn;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("iteration %d\n", iter);
}
redo = false;
@ -124,7 +124,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Component too small: %s\n", m_speciesName[jcomp]);
}
redo = true;
@ -152,12 +152,12 @@ int VCS_SOLVE::vcs_setMolesLinProg()
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) {
if (m_debug_print_lvl == 1) {
printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
plogf(" --- setInitialMoles end\n");
}

4
src/equil/vcs_solve.cpp
View file

@ -362,11 +362,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
// NPhase = number of phases
m_numPhases = nph;
#ifdef DEBUG_MODE
m_debug_print_lvl = pub->vcs_debug_print_lvl;
#else
m_debug_print_lvl = std::min(2, pub->vcs_debug_print_lvl);
#endif
// FormulaMatrix[] -> Copy the formula matrix over
for (size_t i = 0; i < nspecies; i++) {

301
src/equil/vcs_solve_TP.cpp
View file

@ -58,13 +58,9 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
bool uptodate_minors = true;
int forceComponentCalc = 1;
#ifdef DEBUG_MODE
char ANOTE[128];
// Set the debug print lvl to the same as the print lvl.
m_debug_print_lvl = printDetails;
#else
char* ANOTE = 0;
#endif
if (printDetails > 0 && print_lvl == 0) {
print_lvl = 1;
}
@ -256,7 +252,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
npb = vcs_add_all_deleted();
if (npb > 0) {
iti = 0;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) {
if (m_debug_print_lvl >= 1) {
plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!");
plogendl();
}
@ -328,7 +324,7 @@ int VCS_SOLVE::solve_tp_component_calc(bool& allMinorZeroedSpecies)
// EVALUATE THE ELELEMT ABUNDANCE CHECK
if (! vcs_elabcheck(0)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Element Abundance check failed");
plogendl();
}
@ -338,7 +334,7 @@ int VCS_SOLVE::solve_tp_component_calc(bool& allMinorZeroedSpecies)
// Update the phase objects with the contents of the soln vector
vcs_updateVP(VCS_STATECALC_OLD);
vcs_deltag(0, false, VCS_STATECALC_OLD);
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf(" --- Element Abundance check passed");
plogendl();
}
@ -384,9 +380,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// Calculate the total moles in each phase -> old solution
// -> Needed for numerical stability when phases disappear.
// -> the phase moles tend to drift off without this step.
#ifdef DEBUG_MODE
check_tmoles();
#endif
vcs_tmoles();
// COPY OLD into NEW and ZERO VECTORS
// Copy the old solution into the new solution as an initial guess
@ -411,7 +405,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
int soldel = vcs_popPhaseRxnStepSizes(iphasePop);
if (soldel == 3) {
iphasePop = npos;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich "
"prevented phase %d popping\n");
}
@ -428,7 +422,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// minor species in the future too)
kspec = npos;
iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec);
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_RxnStepSizes not called because alternative"
"phase creation delta was used instead\n");
}
@ -453,7 +447,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// m_molNumSpecies_old[kspec] m_molNumSpecies_new[kspec] Species Mole Numbers
// m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers
if (iphaseDelete != npos) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Main Loop Treatment -> Circumvented due to Phase Deletion ");
plogendl();
}
@ -490,7 +484,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// only if ITI NE 0
vcs_deltag(0, false, VCS_STATECALC_NEW);
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Main Loop Treatment of each non-component species ");
if (iti == 0) {
plogf("- Full Calculation:\n");
@ -505,18 +499,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (DEBUG_MODE_ENABLED) {
ANOTE[0] = '\0';
}
ANOTE[0] = '\0';
double dx;
if (iphasePop != npos) {
if (iph == iphasePop) {
dx = m_deltaMolNumSpecies[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Phase pop");
}
sprintf(ANOTE, "Phase pop");
} else {
dx = 0.0;
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
@ -529,7 +519,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
} else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) {
// ZEROED OUT SPECIES
bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0);
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) {
if (m_debug_print_lvl >= 3) {
plogf(" --- %s currently zeroed (SpStatus=%-2d):",
m_speciesName[kspec], m_speciesStatus[kspec]);
plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
@ -540,17 +530,15 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
m_deltaMolNumSpecies[kspec] = 0.0;
resurrect = false;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
if (m_deltaGRxn_new[irxn] < 0.0) {
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
sprintf(ANOTE, "Species stays zeroed even though dg neg due to "
"STOICH/PHASEPOP constraint: DG = %11.4E",
m_deltaGRxn_new[irxn]);
} else {
sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
if (m_deltaGRxn_new[irxn] < 0.0) {
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
sprintf(ANOTE, "Species stays zeroed even though dg neg due to "
"STOICH/PHASEPOP constraint: DG = %11.4E",
m_deltaGRxn_new[irxn]);
} else {
sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
m_deltaGRxn_new[irxn]);
}
}
} else {
@ -561,11 +549,9 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Species stays zeroed even though dG "
"neg, because of %s elemAbund",
m_elementName[j].c_str());
}
sprintf(ANOTE, "Species stays zeroed even though dG "
"neg, because of %s elemAbund",
m_elementName[j].c_str());
resurrect = false;
break;
}
@ -577,14 +563,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// Resurrect the species
if (resurrect) {
if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Zeroed species changed to major: ");
plogf("%-12s\n", m_speciesName[kspec]);
}
m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
allMinorZeroedSpecies = false;
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Zeroed species changed to minor: ");
plogf("%-12s\n", m_speciesName[kspec]);
}
@ -598,9 +584,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
}
m_deltaMolNumSpecies[kspec] = dx;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
}
sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
} else {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
m_deltaMolNumSpecies[kspec] = 0.0;
@ -615,16 +599,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
m_deltaMolNumSpecies[kspec] = 0.0;
dx = 0.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE,"minor species not considered");
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-12s", m_speciesName[kspec]);
plogf("%3d%11.4E%11.4E%11.4E | %s",
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
plogendl();
}
sprintf(ANOTE,"minor species not considered");
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-12s", m_speciesName[kspec]);
plogf("%3d%11.4E%11.4E%11.4E | %s",
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
plogendl();
}
continue;
}
@ -648,7 +630,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
if (soldel_ret) {
// DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF
// MOLE NUMBER
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Delete minor species in multispec phase: %-12s",
m_speciesName[kspec]);
plogendl();
@ -675,9 +657,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
} else {
// MAJOR SPECIES
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "Normal Major Calc");
}
sprintf(ANOTE, "Normal Major Calc");
// Check for superconvergence of the formation reaction. Do
// nothing if it is superconverged. Skip to the end of the irxn
@ -686,16 +666,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
m_deltaMolNumSpecies[kspec] = 0.0;
dx = 0.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "major species is converged");
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-12s", m_speciesName[kspec]);
plogf("%3d%11.4E%11.4E%11.4E | %s",
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
plogendl();
}
sprintf(ANOTE, "major species is converged");
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-12s", m_speciesName[kspec]);
plogf("%3d%11.4E%11.4E%11.4E | %s",
m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
plogendl();
}
continue;
}
@ -712,10 +690,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
} else {
dx = 0.0;
m_deltaMolNumSpecies[kspec] = 0.0;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
"changed initial point");
}
sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
"changed initial point");
}
//Form a tentative value of the new species moles
@ -726,10 +702,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// the outcome, we branch to sections below, or we restart the
// entire iteration.
if (m_molNumSpecies_new[kspec] <= 0.0) {
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
m_molNumSpecies_new[kspec]);
}
sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
m_molNumSpecies_new[kspec]);
// NON-POSITIVE MOLES OF MAJOR SPECIES
//
// We are here when a tentative value of a mole fraction
@ -768,19 +742,15 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
if (m_molNumSpecies_new[kspec] > 0.0) {
m_deltaMolNumSpecies[kspec] = dx;
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE,
"zeroing SS phase created a neg component species "
"-> reducing step size instead");
}
sprintf(ANOTE,
"zeroing SS phase created a neg component species "
"-> reducing step size instead");
} else {
// We are going to zero the single species phase.
// Set the existence flag
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
if (DEBUG_MODE_ENABLED) {
sprintf(ANOTE, "zeroing out SS phase: ");
}
sprintf(ANOTE, "zeroing out SS phase: ");
// Change the base mole numbers for the iteration.
// We need to do this here, because we have decided
@ -789,13 +759,11 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[kspec] = 0.0;
doPhaseDeleteIph = iph;
if (DEBUG_MODE_ENABLED) {
doPhaseDeleteKspec = kspec;
if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) {
plogf(" --- SS species changed to zeroedss: ");
plogf("%-12s", m_speciesName[kspec]);
plogendl();
}
doPhaseDeleteKspec = kspec;
if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) {
plogf(" --- SS species changed to zeroedss: ");
plogf("%-12s", m_speciesName[kspec]);
plogendl();
}
m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
++m_numRxnMinorZeroed;
@ -842,10 +810,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
}
if (DEBUG_MODE_ENABLED) {
checkDelta1(&m_deltaMolNumSpecies[0],
&m_deltaPhaseMoles[0], kspec+1);
}
checkDelta1(&m_deltaMolNumSpecies[0],
&m_deltaPhaseMoles[0], kspec+1);
// Branch point for returning
if (m_debug_print_lvl >= 2) {
@ -860,7 +826,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
if (doPhaseDeleteIph != npos) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- ");
plogf("%-12.12s Main Loop Special Case deleting phase with species: ",
m_speciesName[doPhaseDeleteKspec]);
@ -869,7 +835,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
break;
}
} // END OF MAIN LOOP OVER FORMATION REACTIONS
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
for (size_t k = 0; k < m_numComponents; k++) {
plogf(" --- ");
plogf("%-12.12s", m_speciesName[k]);
@ -888,7 +854,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria
// to limit its magnitude.
double par = 0.5;
size_t ll; // only used in DEBUG_MODE
size_t ll;
for (size_t k = 0; k < m_numComponents; ++k) {
if (m_molNumSpecies_old[k] > 0.0) {
double xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k];
@ -908,7 +874,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
if (par <= 1.01 && par > 0.0) {
// Reduce the size of the step by the multiplicative factor, par
par *= 0.99;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Reduction in step size due to component ");
plogf("%s", m_speciesName[ll]);
plogf(" going negative = %11.3E", par);
@ -923,10 +889,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
} else {
par = 1.0;
}
if (DEBUG_MODE_ENABLED) {
checkDelta1(&m_deltaMolNumSpecies[0],
&m_deltaPhaseMoles[0], m_numSpeciesTot);
}
checkDelta1(&m_deltaMolNumSpecies[0],
&m_deltaPhaseMoles[0], m_numSpeciesTot);
// Now adjust the wt[kspec]'s so that the reflect the decrease in the
// overall length of m_deltaMolNumSpecies[kspec] just calculated. At the
@ -1074,7 +1038,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0],
&m_tPhaseMoles_new[0]));
plogendl();
if (DEBUG_MODE_ENABLED && m_VCount->Its > 550) {
if (m_VCount->Its > 550) {
plogf(" --- Troublesome solve");
plogendl();
}
@ -1101,7 +1065,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// Increment the iteration counters
++m_VCount->Its;
++it1;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Increment counter increased, step is accepted: %4d",
m_VCount->Its);
plogendl();
@ -1117,7 +1081,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
for (size_t iph = 0; iph < m_numPhases; iph++) {
if (!m_VolPhaseList[iph]->m_singleSpecies && m_tPhaseMoles_old[iph] != 0.0 &&
m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) {
if (m_debug_print_lvl >= 1) {
plogf(" --- Setting microscopic phase %d to zero", iph);
plogendl();
}
@ -1146,12 +1110,12 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
// CHECK FOR ELEMENT ABUNDANCE
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Normal element abundance check");
}
vcs_elab();
if (! vcs_elabcheck(0)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" - failed -> redoing element abundances.");
plogendl();
}
@ -1160,7 +1124,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true, VCS_STATECALC_OLD);
uptodate_minors = true;
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf(" - passed");
plogendl();
}
@ -1192,7 +1156,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
}
if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Get a new basis because ");
plogf("%s", m_speciesName[k]);
plogf(" is better than comp ");
@ -1206,7 +1170,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
}
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Check for an optimum basis passed");
plogendl();
}
@ -1219,7 +1183,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
if (iti != 0) {
return;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Reevaluate major-minor status of noncomponents:\n");
}
m_numRxnMinorZeroed = 0;
@ -1227,13 +1191,13 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
size_t kspec = m_indexRxnToSpecies[irxn];
int speciesType = vcs_species_type(kspec);
if (speciesType < VCS_SPECIES_MINOR) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) {
if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) {
plogf(" --- major/minor species is now zeroed out: %s\n",
m_speciesName[kspec]);
}
++m_numRxnMinorZeroed;
} else if (speciesType == VCS_SPECIES_MINOR) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) {
plogf(" --- Noncomponent turned from major to minor: ");
} else if (kspec < m_numComponents) {
@ -1247,7 +1211,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
++m_numRxnMinorZeroed;
} else if (speciesType == VCS_SPECIES_MAJOR) {
if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
plogf(" --- Noncomponent turned from minor to major: ");
} else if (kspec < m_numComponents) {
@ -1275,7 +1239,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
int& stage, bool& lec)
{
if (! allMinorZeroedSpecies) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Equilibrium check for major species: ");
}
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
@ -1288,7 +1252,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
// -> we have run out of iterations!
stage = RETURN_A;
} else {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]);
}
// Convergence amongst major species has not been achieved.
@ -1299,11 +1263,11 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
return;
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" MAJOR SPECIES CONVERGENCE achieved");
plogendl();
}
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf(" MAJOR SPECIES CONVERGENCE achieved "
"(because there are no major species)");
plogendl();
@ -1320,7 +1284,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
vcs_deltag(1, false, VCS_STATECALC_OLD);
uptodate_minors = true;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Equilibrium check for minor species: ");
}
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
@ -1333,7 +1297,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
stage = RETURN_A;
return;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]);
}
@ -1344,7 +1308,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
return;
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" CONVERGENCE achieved\n");
}
}
@ -1357,14 +1321,14 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
// LEC is only true when we are near the end game
if (lec) {
if (!giveUpOnElemAbund) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Check the Full Element Abundances: ");
}
// Final element abundance check: If we fail then we need to go back
// and correct the element abundances, and then go do a major step
if (! vcs_elabcheck(1)) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (! vcs_elabcheck(0)) {
plogf(" failed\n");
} else {
@ -1375,7 +1339,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
stage = ELEM_ABUND_CHECK;
return;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" passed\n");
}
}
@ -1437,7 +1401,7 @@ void VCS_SOLVE::solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec,
} else if (ncAfter) {
if (!neAfter) {
// Probably an unrecoverable range error
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n");
plogf(" --- vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n");
plogf(" --- vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n");
@ -1477,7 +1441,7 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
}
dg_irxn = std::max(dg_irxn, -200.0);
if (DEBUG_MODE_ENABLED && ANOTE) {
if (ANOTE) {
sprintf(ANOTE,"minor species alternative calc");
}
if (dg_irxn >= 23.0) {
@ -1542,7 +1506,7 @@ L_ZERO_SPECIES:
// Voltage calculation
// Need to check the sign -> This is good for electrons
dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
if (DEBUG_MODE_ENABLED && ANOTE) {
if (ANOTE) {
sprintf(ANOTE,"voltage species alternative calc");
}
}
@ -1608,7 +1572,7 @@ int VCS_SOLVE::vcs_zero_species(const size_t kspec)
double dx = -m_molNumSpecies_old[kspec];
if (dx != 0.0) {
retn = delta_species(kspec, &dx);
if (DEBUG_MODE_ENABLED && !retn && m_debug_print_lvl >= 1) {
if (!retn && m_debug_print_lvl >= 1) {
plogf("vcs_zero_species: Couldn't zero the species %d, "
"did delta of %g. orig conc of %g",
kspec, dx, m_molNumSpecies_old[kspec] + dx);
@ -1683,7 +1647,7 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec)
void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
{
size_t iph = m_phaseID[kspec];
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Add back a deleted species: %-12s\n", m_speciesName[kspec]);
}
@ -1738,7 +1702,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
// set the phase existence flag to dead
Vphase->setTotalMoles(0.0);
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName);
}
@ -1753,7 +1717,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
int retn = delta_species(kspec, &dxTent);
if (retn != 1) {
successful = false;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
iph, Vphase->PhaseName, m_speciesName[kspec]);
plogf(" --- delta attempted: %g achieved: %g "
@ -1781,7 +1745,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0;
for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
if (m_phaseID[kcomp] == iph) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n",
kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]);
}
@ -1818,7 +1782,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
}
if (m_molNumSpecies_old[kcomp] != 0.0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n",
kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]);
plogf(" --- zeroing it \n");
@ -1843,7 +1807,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
++m_numRxnRdc;
++m_numSpeciesRdc;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Make %s", m_speciesName[kspec]);
plogf(" an active but zeroed species because its phase "
"was zeroed\n");
@ -1868,7 +1832,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
int VCS_SOLVE::vcs_recheck_deleted()
{
vector_fp& xtcutoff = m_TmpPhase;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Start rechecking deleted species in multispec phases\n");
}
if (m_numSpeciesRdc == m_numSpeciesTot) {
@ -2019,20 +1983,20 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
double dx = m_molNumSpecies_new[kspec];
retn = delta_species(kspec, &dx);
if (retn == 0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl) {
if (m_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec], kspec, dx);
}
if (dx > 1.0E-50) {
dx = 1.0E-50;
retn = delta_species(kspec, &dx);
if (DEBUG_MODE_ENABLED && retn == 0 && m_debug_print_lvl) {
if (retn == 0 && m_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec], kspec, dx);
}
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (retn != 0) {
plogf(" --- add_deleted(): species %s added back in with mol number %g",
m_speciesName[kspec], dx);
@ -2057,7 +2021,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
VCS_DELETE_MINORSPECIES_CUTOFF) ||
(m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) {
retn++;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
m_speciesName[kspec], m_molNumSpecies_old[kspec],
m_deltaGRxn_old[irxn]);
@ -2092,13 +2056,13 @@ bool VCS_SOLVE::vcs_globStepDamp()
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1);
plogf(" --- subroutine FORCE: End Slope = %g\n", s2);
}
if (s1 > 0.0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE produced no adjustments,");
if (s1 < 1.0E-40) {
plogf(" s1 positive but really small");
@ -2111,7 +2075,7 @@ bool VCS_SOLVE::vcs_globStepDamp()
}
if (s2 <= 0.0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE produced no adjustments, s2 < 0");
plogendl();
}
@ -2124,17 +2088,17 @@ bool VCS_SOLVE::vcs_globStepDamp()
al = s1 / (s1 - s2);
}
if (al >= 0.95 || al < 0.0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE produced no adjustments (al = %g)\n", al);
}
return false;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE produced a damping factor = %g\n", al);
}
// ADJUST MOLE NUMBERS, CHEM. POT
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
m_deltaGRxn_tmp = m_deltaGRxn_new;
}
@ -2148,7 +2112,7 @@ bool VCS_SOLVE::vcs_globStepDamp()
}
vcs_updateVP(VCS_STATECALC_NEW);
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE adjusted the mole "
"numbers, AL = %10.3f\n", al);
}
@ -2172,7 +2136,7 @@ bool VCS_SOLVE::vcs_globStepDamp()
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- subroutine FORCE: Adj End Slope = %g", s2);
plogendl();
}
@ -2187,7 +2151,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
size_t jlose = npos;
double* scrxn_ptr;
clockWC tickTock;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ");
for (size_t i=0; i<77; i++) {
plogf("-");
@ -2363,7 +2327,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
m_indexRxnToSpecies[i] = ncTrial + i;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Total number of components found = %3d (ne = %d)\n ",
ncTrial, m_numElemConstraints);
}
@ -2417,7 +2381,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[k]);
if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf("(Volts = %9.2g)", m_molNumSpecies_old[k]);
@ -2434,7 +2398,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
}
vcs_switch_pos(false, jr, k);
std::swap(aw[jr], aw[k]);
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[k]);
if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf("(Volts = %9.2g) remains ", m_molNumSpecies_old[k]);
@ -2560,7 +2524,7 @@ L_END_LOOP:
m_scSize[i] = szTmp;
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Components:");
for (size_t j = 0; j < ncTrial; j++) {
plogf(" %3d", j);
@ -2756,7 +2720,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %s can not be nonzero because"
" needed element %s is zero\n",
m_speciesName[kspec], m_elementName[j]);
@ -2782,7 +2746,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
double negChangeComp = - stoicC;
if (negChangeComp > 0.0) {
if (m_molNumSpecies_old[j] < 1.0E-60) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %s is prevented from popping into existence because"
" a needed component to be consumed, %s, has a zero mole number\n",
m_speciesName[kspec], m_speciesName[j]);
@ -2797,7 +2761,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
size_t jph = m_phaseID[j];
vcs_VolPhase* jVPhase = m_VolPhaseList[jph];
if (jVPhase->exists() <= 0) {
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %s is prevented from popping into existence because"
" a needed component %s is in a zeroed-phase that would be "
"popped into existence at the same time\n",
@ -2954,7 +2918,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
tPhMoles_ptr = &m_tPhaseMoles_new[0];
actCoeff_ptr = &m_actCoeffSpecies_new[0];
molNum = &m_molNumSpecies_new[0];
} else if (DEBUG_MODE_ENABLED) {
} else {
throw CanteraError("VCS_SOLVE::vcs_dfe",
"Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc);
}
@ -2962,7 +2926,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
AssertThrowMsg(m_unitsState != VCS_DIMENSIONAL_G, "VCS_SOLVE::vcs_dfe",
"called with wrong units state");
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
if (ll == 0) {
if (lbot != 0) {
plogf(" --- Subroutine vcs_dfe called for one species: ");
@ -3239,7 +3203,6 @@ void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const
writeline('-', 132);
}
#ifdef DEBUG_MODE
void VCS_SOLVE::prneav() const
{
vector_fp eav(m_numElemConstraints, 0.0);
@ -3276,7 +3239,6 @@ void VCS_SOLVE::prneav() const
plogf("--------------------------------------------------");
plogendl();
}
#endif
double VCS_SOLVE::l2normdg(double dgLocal[]) const
{
@ -3320,7 +3282,6 @@ double VCS_SOLVE::vcs_tmoles()
return m_totalMolNum;
}
#ifdef DEBUG_MODE
void VCS_SOLVE::check_tmoles() const
{
double sum = 0.0;
@ -3341,7 +3302,6 @@ void VCS_SOLVE::check_tmoles() const
}
}
}
#endif
void VCS_SOLVE::vcs_updateVP(const int vcsState)
{
@ -3355,7 +3315,7 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState)
Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW,
&m_molNumSpecies_new[0],
&m_tPhaseMoles_new[0]);
} else if (DEBUG_MODE_ENABLED) {
} else {
throw CanteraError("VCS_SOLVE::vcs_updateVP",
"wrong stateCalc value: {}", vcsState);
}
@ -3365,15 +3325,15 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState)
bool VCS_SOLVE::vcs_evaluate_speciesType()
{
m_numRxnMinorZeroed = 0;
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Species Status decision is reevaluated: All species are minor except for:\n");
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) {
} else if (m_debug_print_lvl >= 5) {
plogf(" --- Species Status decision is reevaluated");
plogendl();
}
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_speciesStatus[kspec] = vcs_species_type(kspec);
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) {
if (m_debug_print_lvl >= 5) {
plogf(" --- %-16s: ", m_speciesName[kspec]);
if (kspec < m_numComponents) {
plogf("(COMP) ");
@ -3383,7 +3343,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType()
plogf(" %10.3g ", m_molNumSpecies_old[kspec]);
const char* sString = vcs_speciesType_string(m_speciesStatus[kspec], 100);
plogf("%s\n", sString);
} else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
} else if (m_debug_print_lvl >= 2) {
if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
switch (m_speciesStatus[kspec]) {
case VCS_SPECIES_COMPONENT:
@ -3425,7 +3385,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType()
++m_numRxnMinorZeroed;
}
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" ---");
plogendl();
}
@ -3459,7 +3419,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
throw CanteraError("VCS_SOLVE::vcs_deltag", "bad vcsState");
}
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Subroutine vcs_deltag called for ");
if (L < 0) {
plogf("major noncomponents\n");
@ -3755,7 +3715,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
feSpecies = &m_feSpecies_old[0];
deltaGRxn = &m_deltaGRxn_old[0];
actCoeffSpecies = &m_actCoeffSpecies_old[0];
} else if (DEBUG_MODE_ENABLED) {
} else {
throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc");
}
@ -3765,7 +3725,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
}
vcs_VolPhase* vPhase = m_VolPhaseList[iphase];
if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
if (m_debug_print_lvl >= 2) {
plogf(" --- Subroutine vcs_deltag_Phase called for phase %d\n",
iphase);
}
@ -3858,8 +3818,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
if (k1 == k2) {
return;
}
if (DEBUG_MODE_ENABLED && (k1 >= m_numSpeciesTot ||
k2 >= m_numSpeciesTot)) {
if (k1 >= m_numSpeciesTot || k2 >= m_numSpeciesTot) {
plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
k1, k2);
}
@ -3915,7 +3874,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
// Find the Rxn indices corresponding to the two species
size_t i1 = k1 - m_numComponents;
size_t i2 = k2 - m_numComponents;
if (DEBUG_MODE_ENABLED && (i1 > m_numRxnTot || i2 >= m_numRxnTot)) {
if (i1 > m_numRxnTot || i2 >= m_numRxnTot) {
plogf("switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n",
i1 , i2);
}

50
src/numerics/RootFind.cpp
View file

@ -8,14 +8,6 @@
*/
#include "cantera/numerics/RootFind.h"
// turn on debugging for now
#ifndef DEBUG_MODE
#define DEBUG_MODE
#undef DEBUG_MODE_ENABLED
#define DEBUG_MODE_ENABLED 1
#endif
#include "cantera/base/utilities.h"
#include "cantera/base/stringUtils.h"
@ -201,7 +193,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
rfT.reasoning = "First Point: ";
callNum++;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fp = fopen(fmt::format("RootFind_%d.log", callNum).c_str(), "w");
fprintf(fp, " Iter TP_its xval Func_val | Reasoning\n");
fprintf(fp, "-----------------------------------------------------"
@ -256,7 +248,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
int its = 1;
f1 = func(x1);
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
print_funcEval(fp, x1, f1, its);
fprintf(fp, "%-5d %-5d %-15.5E %-15.5E\n", -2, 0, x1, f1);
}
@ -305,7 +297,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
deltaX2 = x2 - x1;
its++;
f2 = func(x2);
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
print_funcEval(fp, x2, f2, its);
fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2);
}
@ -359,7 +351,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
while (!converged && its < itmax) {
// Find an estimate of the next point, xnew, to try based on a linear
// approximation from the last two points.
if (DEBUG_MODE_ENABLED && fabs(x2 - x1) < 1.0E-14) {
if (fabs(x2 - x1) < 1.0E-14) {
writelogf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1);
}
doublereal delXtmp = deltaXControlled(x2, x1);
@ -387,7 +379,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
}
rfT.reasoning += "Slope is good. ";
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | xlin = %-11.5E", xnew);
}
deltaXnew = xnew - x2;
@ -478,13 +470,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
xDelMin = fabs(x2 - x1) / 10.;
if (fabs(xnew - x1) < xDelMin) {
xnew = x1 + sign(xnew-x1) * xDelMin;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | x10%% = %-11.5E", xnew);
}
}
if (fabs(xnew - x2) < 0.1 * xDelMin) {
xnew = x2 + sign(xnew-x2) * 0.1 * xDelMin;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | x10%% = %-11.5E", xnew);
}
}
@ -498,7 +490,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
}
if (fabs(xDelMax) < fabs(xnew - x2)) {
xnew = x2 + sign(xnew-x2) * xDelMax;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | xlimitsize = %-11.5E", xnew);
}
}
@ -510,13 +502,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
xDelMin = 0.1 * fabs(x2 - x1);
if (fabs(xnew - x2) < xDelMin) {
xnew = x2 + sign(xnew - x2) * xDelMin;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | x10%% = %-11.5E", xnew);
}
}
if (fabs(xnew - x1) < xDelMin) {
xnew = x1 + sign(xnew - x1) * xDelMin;
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | x10%% = %-11.5E", xnew);
}
}
@ -558,7 +550,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
}
}
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_ && xorig != xnew) {
if (printLvl >= 3 && writeLogAllowed_ && xorig != xnew) {
fprintf(fp, " | xstraddle = %-11.5E", xnew);
}
}
@ -596,7 +588,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
xnew = xmax;
}
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | xlimitmax = %-11.5E", xnew);
}
}
@ -620,7 +612,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
xnew = xmin;
}
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | xlimitmin = %-11.5E", xnew);
}
}
@ -628,7 +620,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun
its++;
fnew = func(xnew);
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp,"\n");
print_funcEval(fp, xnew, fnew, its);
fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew);
@ -947,7 +939,7 @@ done:
if (printLvl >= 1) {
writelogf("RootFind success: convergence achieved\n");
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, " | RootFind success in %d its, fnorm = %g\n", its, fnorm);
}
rfHistory_.push_back(rfT);
@ -972,7 +964,7 @@ done:
}
rfT.reasoning += "Maximum iterations exceeded without convergence, cause unknown";
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fprintf(fp, "\nRootFind failure in %d its\n", its);
}
@ -982,7 +974,7 @@ done:
rfT.fval = f2;
rfHistory_.push_back(rfT);
}
if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) {
if (printLvl >= 3 && writeLogAllowed_) {
fclose(fp);
}
@ -996,13 +988,9 @@ done:
doublereal RootFind::func(doublereal x)
{
doublereal r;
if (DEBUG_MODE_ENABLED) {
checkFinite(x);
}
checkFinite(x);
m_residFunc->evalSS(0.0, &x, &r);
if (DEBUG_MODE_ENABLED) {
checkFinite(r);
}
checkFinite(r);
doublereal ff = r - m_funcTargetValue;
if (x >= x_maxTried_) {
x_maxTried_ = x;

249
src/thermo/HMWSoln.cpp
View file

@ -1492,7 +1492,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// molalitysum is the sum of the molalities over all solutes, even those
// with zero charge.
double molalitysumUncropped = 0.0;
debuglog("\n Debugging information from hmw_act \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog("\n Debugging information from hmw_act \n", m_debugCalc);
// Make sure the counter variables are setup
counterIJ_setup();
@ -1510,7 +1510,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// Store the ionic molality in the object for reference.
m_IionicMolality = Is;
sqrtIs = sqrt(Is);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelog(" Step 1: \n");
writelogf(" ionic strenth = %14.7le \n total molar "
"charge = %14.7le \n", Is, molarcharge);
@ -1522,11 +1522,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// Step 2: Find the coefficients E-theta and E-thetaprime for all
// combinations of positive unlike charges up to 4
debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 2: \n", m_debugCalc);
for (int z1 = 1; z1 <=4; z1++) {
for (int z2 =1; z2 <=4; z2++) {
calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
}
@ -1535,7 +1535,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
debuglog(" Step 3: \n"
" Species Species g(x) hfunc(x)\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// calculate g(x) and hfunc(x) for each cation-anion pair MX. In the
// original literature, hfunc, was called gprime. However, it's not the
@ -1574,7 +1574,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
m_gfunc_IJ[counterIJ] = 0.0;
m_hfunc_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i),
speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
}
@ -1585,7 +1585,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// Agrees with Pitzer, Eq. (49), (51), (55)
debuglog(" Step 4: \n"
" Species Species BMX BprimeMX BphiMX\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk - 1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
@ -1600,7 +1600,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
+ m_Beta1MX_ij[counterIJ] * m_gfunc_IJ[counterIJ]
+ m_Beta2MX_ij[counterIJ] * m_g2func_IJ[counterIJ];
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf("%d %g: %g %g %g %g\n",
counterIJ, m_BMX_IJ[counterIJ], m_Beta0MX_ij[counterIJ],
m_Beta1MX_ij[counterIJ], m_Beta2MX_ij[counterIJ], m_gfunc_IJ[counterIJ]);
@ -1617,7 +1617,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
m_BprimeMX_IJ[counterIJ] = 0.0;
m_BphiMX_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
speciesName(i), speciesName(j),
m_BMX_IJ[counterIJ], m_BprimeMX_IJ[counterIJ], m_BphiMX_IJ[counterIJ]);
@ -1628,8 +1628,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// SUBSECTION TO CALCULATE CMX
// Agrees with Pitzer, Eq. (53).
debuglog(" Step 5: \n"
" Species Species CMX\n",
DEBUG_MODE_ENABLED && m_debugCalc);
" Species Species CMX\n", m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -1644,7 +1643,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
} else {
m_CMX_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f \n",
speciesName(i), speciesName(j), m_CMX_IJ[counterIJ]);
}
@ -1655,7 +1654,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// Agrees with Pitzer, Eq. 72, 73, 74
debuglog(" Step 6: \n"
" Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -1675,7 +1674,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
m_Phiprime_IJ[counterIJ] = 0.0;
m_PhiPhi_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
speciesName(i), speciesName(j),
m_Phi_IJ[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ[counterIJ]);
@ -1685,11 +1684,11 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// SUBSECTION FOR CALCULATION OF F
// Agrees with Pitzer Eqn. (65)
debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 7: \n", m_debugCalc);
double Aphi = A_Debye_TP() / 3.0;
double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
+ (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" initial value of F = %10.6f \n", F);
}
for (size_t i = 1; i < m_kk-1; i++) {
@ -1709,13 +1708,13 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
if (charge(i)*charge(j) > 0) {
F += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" F = %10.6f \n", F);
}
}
}
debuglog(" Step 8: Summing in All Contributions to Activity Coefficients \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk; i++) {
@ -1723,12 +1722,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// equations agree with my notes, Eqn. (118).
// Equations agree with Pitzer, eqn.(63)
if (charge(i) > 0.0) {
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
}
// species i is the cation (positive) to calc the actcoeff
double zsqF = charge(i)*charge(i)*F;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Unary term: z*z*F = %10.5f\n", zsqF);
}
double sum1 = 0.0;
@ -1745,7 +1744,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// sum over all anions
sum1 += molality[j] *
(2.0*m_BMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
std::string snj = speciesName(j) + ":";
writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj,
molality[j]*2.0*m_BMX_IJ[counterIJ]);
@ -1761,7 +1760,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
if (charge(k) < 0.0) {
n = k + j * m_kk + i * m_kk * m_kk;
sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) {
if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1775,7 +1774,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// sum over all cations
if (j != i) {
sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
std::string snj = speciesName(j) + ":";
writelogf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj,
molality[j]*(2.0*m_Phi_IJ[counterIJ]));
@ -1786,7 +1785,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// two inner sums over anions
n = k + j * m_kk + i * m_kk * m_kk;
sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) {
if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1797,7 +1796,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
size_t counterIJ2 = m_CounterIJ[n];
sum4 += (fabs(charge(i))*
molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj,
fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
@ -1809,7 +1808,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// Handle neutral j species
if (charge(j) == 0) {
sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
std::string snj = speciesName(j) + ":";
writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj,
molality[j]*2.0*m_Lambda_nj(j,i));
@ -1824,7 +1823,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
double zeta = m_Psi_ijk[n];
if (zeta != 0.0) {
sum5 += molality[j]*molality[k]*zeta;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1839,7 +1838,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// activity coefficients (molality scale)
m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
}
@ -1849,12 +1848,12 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// equations agree with my notes, Eqn. (119).
// Equations agree with Pitzer, eqn.(64)
if (charge(i) < 0) {
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
}
// species i is an anion (negative)
double zsqF = charge(i)*charge(i)*F;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Unary term: z*z*F = %10.5f\n", zsqF);
}
double sum1 = 0.0;
@ -1871,7 +1870,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
if (charge(j) > 0) {
sum1 += molality[j]*
(2.0*m_BMX_IJ[counterIJ]+molarcharge*m_CMX_IJ[counterIJ]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
std::string snj = speciesName(j) + ":";
writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj,
molality[j]*2.0*m_BMX_IJ[counterIJ]);
@ -1884,7 +1883,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
if (charge(k) > 0) {
n = k + j * m_kk + i * m_kk * m_kk;
sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) {
if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1899,7 +1898,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// sum over all anions
if (j != i) {
sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
std::string snj = speciesName(j) + ":";
writelogf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj,
molality[j]*(2.0*m_Phi_IJ[counterIJ]));
@ -1910,7 +1909,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// two inner sums over cations
n = k + j * m_kk + i * m_kk * m_kk;
sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) {
if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1920,7 +1919,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
size_t counterIJ2 = m_CounterIJ[n];
sum4 += fabs(charge(i))*
molality[j]*molality[k]*m_CMX_IJ[counterIJ2];
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj,
fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
@ -1932,7 +1931,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// for Anions, do the neutral species interaction
if (charge(j) == 0.0) {
sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
std::string snj = speciesName(j) + ":";
writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj,
molality[j]*2.0*m_Lambda_nj(j,i));
@ -1947,7 +1946,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
double zeta = m_Psi_ijk[n];
if (zeta != 0.0) {
sum5 += molality[j]*molality[k]*zeta;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1959,7 +1958,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
}
m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
}
@ -1969,14 +1968,14 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
// equations agree with my notes,
// Equations agree with Pitzer,
if (charge(i) == 0.0) {
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
}
double sum1 = 0.0;
double sum3 = 0.0;
for (size_t j = 1; j < m_kk; j++) {
sum1 += molality[j]*2.0*m_Lambda_nj(i,j);
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Lambda_nj(i,j) != 0.0) {
if (m_debugCalc && m_Lambda_nj(i,j) != 0.0) {
std::string snj = speciesName(j) + ":";
writelogf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj,
molality[j]*2.0*m_Lambda_nj(i,j));
@ -1987,7 +1986,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
if (charge(k) < 0.0) {
size_t n = k + j * m_kk + i * m_kk * m_kk;
sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Psi_ijk[n] != 0.0) {
if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
std::string snj = speciesName(j) + "," + speciesName(k) + ":";
writelogf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
molality[j]*molality[k]*m_Psi_ijk[n]);
@ -1997,19 +1996,19 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
}
}
double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
if (DEBUG_MODE_ENABLED && m_debugCalc && m_Mu_nnn[i] != 0.0) {
if (m_debugCalc && m_Mu_nnn[i] != 0.0) {
writelogf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
}
m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
}
}
}
debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 9: \n", m_debugCalc);
// SUBSECTION FOR CALCULATING THE OSMOTIC COEFF
// equations agree with my notes, Eqn. (117).
@ -2133,7 +2132,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
} else {
osmotic_coef = 1.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
"sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
term1, sum1, sum2, sum3, sum4, sum5);
@ -2152,7 +2151,7 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
double xmolSolvent = moleFraction(m_indexSolvent);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
double wateract = exp(lnwateract);
writelogf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
writelogf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
@ -2195,9 +2194,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// immediately preceding the calling of this routine. Therefore, some
// quantities do not need to be recalculated in this routine.
#ifdef DEBUG_MODE
m_debugCalc = 0;
#endif
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
@ -2222,7 +2218,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
double molalitysum = 0.0;
debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dT()\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// Make sure the counter variables are setup
counterIJ_setup();
@ -2240,7 +2236,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// Store the ionic molality in the object for reference.
m_IionicMolality = Is;
sqrtIs = sqrt(Is);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelog(" Step 1: \n");
writelogf(" ionic strenth = %14.7le \n total molar "
"charge = %14.7le \n", Is, molarcharge);
@ -2252,11 +2248,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// Step 2: Find the coefficients E-theta and E-thetaprime for all
// combinations of positive unlike charges up to 4
debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 2: \n", m_debugCalc);
for (int z1 = 1; z1 <=4; z1++) {
for (int z2 =1; z2 <=4; z2++) {
calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
}
@ -2265,7 +2261,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
debuglog(" Step 3: \n"
" Species Species g(x) hfunc(x) \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// calculate g(x) and hfunc(x) for each cation-anion pair MX
// In the original literature, hfunc, was called gprime. However,
@ -2304,7 +2300,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_gfunc_IJ[counterIJ] = 0.0;
m_hfunc_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
std::string sni = speciesName(i);
std::string snj = speciesName(j);
writelogf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
@ -2318,7 +2314,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// quantities.
debuglog(" Step 4: \n"
" Species Species BMX BprimeMX BphiMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk - 1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
@ -2332,7 +2328,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_BMX_IJ_L[counterIJ] = m_Beta0MX_ij_L[counterIJ]
+ m_Beta1MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ]
+ m_Beta2MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ];
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf("%d %g: %g %g %g %g\n",
counterIJ, m_BMX_IJ_L[counterIJ], m_Beta0MX_ij_L[counterIJ],
m_Beta1MX_ij_L[counterIJ], m_Beta2MX_ij_L[counterIJ], m_gfunc_IJ[counterIJ]);
@ -2349,7 +2345,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_BprimeMX_IJ_L[counterIJ] = 0.0;
m_BphiMX_IJ_L[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
speciesName(i), speciesName(j),
m_BMX_IJ_L[counterIJ], m_BprimeMX_IJ_L[counterIJ], m_BphiMX_IJ_L[counterIJ]);
@ -2359,8 +2355,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// --------- SUBSECTION TO CALCULATE CMX_L ----------
debuglog(" Step 5: \n"
" Species Species CMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
" Species Species CMX \n", m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -2375,7 +2370,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
} else {
m_CMX_IJ_L[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f \n",
speciesName(i), speciesName(j), m_CMX_IJ_L[counterIJ]);
}
@ -2385,7 +2380,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
debuglog(" Step 6: \n"
" Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -2403,7 +2398,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_Phiprime_IJ[counterIJ] = 0.0;
m_PhiPhi_IJ_L[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
speciesName(i), speciesName(j),
m_Phi_IJ_L[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_L[counterIJ]);
@ -2412,12 +2407,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
}
// ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 7: \n", m_debugCalc);
double dA_DebyedT = dA_DebyedT_TP();
double dAphidT = dA_DebyedT /3.0;
double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
+ (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" initial value of dFdT = %10.6f \n", dFdT);
}
for (size_t i = 1; i < m_kk-1; i++) {
@ -2437,12 +2432,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
if (charge(i)*charge(j) > 0) {
dFdT += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" dFdT = %10.6f \n", dFdT);
}
}
}
debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 8: \n", m_debugCalc);
for (size_t i = 1; i < m_kk; i++) {
// -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
@ -2521,7 +2516,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_dlnActCoeffMolaldT_Unscaled[i] =
zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -2598,7 +2593,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
m_dlnActCoeffMolaldT_Unscaled[i] =
zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -2627,13 +2622,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
}
}
}
debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 9: \n", m_debugCalc);
// ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
double sum1 = 0.0;
@ -2755,7 +2750,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
d_osmotic_coef_dT = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
"sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
term1, sum1, sum2, sum3, sum4, sum5);
@ -2772,7 +2767,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
// We have just computed act_0. However, this routine returns
// ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
double d_wateract_dT = exp(d_lnwateract_dT);
writelogf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
d_lnwateract_dT, d_wateract_dT);
@ -2809,9 +2804,6 @@ void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
{
#ifdef DEBUG_MODE
m_debugCalc = 0;
#endif
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
@ -2835,7 +2827,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
double molalitysum = 0.0;
debuglog("\n Debugging information from s_Pitzer_d2lnMolalityActCoeff_dT2()\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// Make sure the counter variables are setup
counterIJ_setup();
@ -2853,7 +2845,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// Store the ionic molality in the object for reference.
m_IionicMolality = Is;
sqrtIs = sqrt(Is);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelog(" Step 1: \n");
writelogf(" ionic strenth = %14.7le \n total molar "
"charge = %14.7le \n", Is, molarcharge);
@ -2865,11 +2857,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// Step 2: Find the coefficients E-theta and E-thetaprime for all
// combinations of positive unlike charges up to 4
debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 2: \n", m_debugCalc);
for (int z1 = 1; z1 <=4; z1++) {
for (int z2 =1; z2 <=4; z2++) {
calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
}
@ -2878,7 +2870,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
debuglog(" Step 3: \n"
" Species Species g(x) hfunc(x) \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// calculate gfunc(x) and hfunc(x) for each cation-anion pair MX. In the
// original literature, hfunc, was called gprime. However, it's not the
@ -2917,7 +2909,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
m_gfunc_IJ[counterIJ] = 0.0;
m_hfunc_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i), speciesName(j),
m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
}
@ -2929,7 +2921,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// quantities.
debuglog(" Step 4: \n"
" Species Species BMX BprimeMX BphiMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk - 1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
@ -2943,7 +2935,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
m_BMX_IJ_LL[counterIJ] = m_Beta0MX_ij_LL[counterIJ]
+ m_Beta1MX_ij_LL[counterIJ] * m_gfunc_IJ[counterIJ]
+ m_Beta2MX_ij_LL[counterIJ] * m_g2func_IJ[counterIJ];
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf("%d %g: %g %g %g %g\n",
counterIJ, m_BMX_IJ_LL[counterIJ], m_Beta0MX_ij_LL[counterIJ],
m_Beta1MX_ij_LL[counterIJ], m_Beta2MX_ij_LL[counterIJ], m_gfunc_IJ[counterIJ]);
@ -2960,7 +2952,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
m_BprimeMX_IJ_LL[counterIJ] = 0.0;
m_BphiMX_IJ_LL[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
speciesName(i), speciesName(j),
m_BMX_IJ_LL[counterIJ], m_BprimeMX_IJ_LL[counterIJ], m_BphiMX_IJ_LL[counterIJ]);
@ -2970,8 +2962,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// --------- SUBSECTION TO CALCULATE CMX_LL ----------
debuglog(" Step 5: \n"
" Species Species CMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
" Species Species CMX \n", m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -2986,7 +2977,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
} else {
m_CMX_IJ_LL[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f \n",
speciesName(i), speciesName(j), m_CMX_IJ_LL[counterIJ]);
}
@ -2996,7 +2987,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
debuglog(" Step 6: \n"
" Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -3014,7 +3005,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
m_Phiprime_IJ[counterIJ] = 0.0;
m_PhiPhi_IJ_LL[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
speciesName(i), speciesName(j),
m_Phi_IJ_LL[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_LL[counterIJ]);
@ -3023,11 +3014,11 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
}
// ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 7: \n", m_debugCalc);
double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
+ (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" initial value of d2FdT2 = %10.6f \n", d2FdT2);
}
for (size_t i = 1; i < m_kk-1; i++) {
@ -3047,12 +3038,12 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
if (charge(i)*charge(j) > 0) {
d2FdT2 += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" d2FdT2 = %10.6f \n", d2FdT2);
}
}
}
debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 8: \n", m_debugCalc);
for (size_t i = 1; i < m_kk; i++) {
// -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
@ -3128,7 +3119,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// solute activity coefficients (molality scale)
m_d2lnActCoeffMolaldT2_Unscaled[i] =
zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s d2lngammadT2[i]=%10.6f \n",
speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -3205,7 +3196,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
}
m_d2lnActCoeffMolaldT2_Unscaled[i] =
zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s d2lngammadT2[i]=%10.6f\n",
speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -3233,14 +3224,14 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
}
double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelog(" %-16s d2lngammadT2[i]=%10.6f \n",
speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
}
}
}
debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 9: \n", m_debugCalc);
// ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
double sum1 = 0.0;
@ -3364,7 +3355,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
} else {
d2_osmotic_coef_dT2 = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
"sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
term1, sum1, sum2, sum3, sum4, sum5);
@ -3382,7 +3373,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
// ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
writelogf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
d2_lnwateract_dT2, d2_wateract_dT2);
@ -3415,9 +3406,6 @@ void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
{
#ifdef DEBUG_MODE
m_debugCalc = 0;
#endif
// HKM -> Assumption is made that the solvent is species 0.
if (m_indexSolvent != 0) {
throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
@ -3443,7 +3431,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
double currPres = pressure();
debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dP()\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// Make sure the counter variables are setup
counterIJ_setup();
@ -3461,7 +3449,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// Store the ionic molality in the object for reference.
m_IionicMolality = Is;
sqrtIs = sqrt(Is);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelog(" Step 1: \n");
writelogf(" ionic strenth = %14.7le \n total molar "
"charge = %14.7le \n", Is, molarcharge);
@ -3474,11 +3462,11 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// Step 2: Find the coefficients E-theta and E-thetaprime for all
// combinations of positive unlike charges up to 4
debuglog(" Step 2: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 2: \n", m_debugCalc);
for (int z1 = 1; z1 <=4; z1++) {
for (int z2 =1; z2 <=4; z2++) {
calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
}
@ -3487,7 +3475,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
debuglog(" Step 3: \n"
" Species Species g(x) hfunc(x)\n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
// calculate g(x) and hfunc(x) for each cation-anion pair MX
// In the original literature, hfunc, was called gprime. However,
@ -3526,7 +3514,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
m_gfunc_IJ[counterIJ] = 0.0;
m_hfunc_IJ[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i),
speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
}
@ -3538,7 +3526,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// quantities.
debuglog(" Step 4: \n"
" Species Species BMX BprimeMX BphiMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk - 1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
@ -3552,7 +3540,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
m_BMX_IJ_P[counterIJ] = m_Beta0MX_ij_P[counterIJ]
+ m_Beta1MX_ij_P[counterIJ] * m_gfunc_IJ[counterIJ]
+ m_Beta2MX_ij_P[counterIJ] * m_g2func_IJ[counterIJ];
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf("%d %g: %g %g %g %g\n",
counterIJ, m_BMX_IJ_P[counterIJ], m_Beta0MX_ij_P[counterIJ],
m_Beta1MX_ij_P[counterIJ], m_Beta2MX_ij_P[counterIJ], m_gfunc_IJ[counterIJ]);
@ -3569,7 +3557,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
m_BprimeMX_IJ_P[counterIJ] = 0.0;
m_BphiMX_IJ_P[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
speciesName(i), speciesName(j),
m_BMX_IJ_P[counterIJ], m_BprimeMX_IJ_P[counterIJ], m_BphiMX_IJ_P[counterIJ]);
@ -3579,8 +3567,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// --------- SUBSECTION TO CALCULATE CMX_P ----------
debuglog(" Step 5: \n"
" Species Species CMX \n",
DEBUG_MODE_ENABLED && m_debugCalc);
" Species Species CMX \n", m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -3595,7 +3582,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
} else {
m_CMX_IJ_P[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %11.7f \n",
speciesName(i), speciesName(j), m_CMX_IJ_P[counterIJ]);
}
@ -3605,7 +3592,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
debuglog(" Step 6: \n"
" Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
DEBUG_MODE_ENABLED && m_debugCalc);
m_debugCalc);
for (size_t i = 1; i < m_kk-1; i++) {
for (size_t j = i+1; j < m_kk; j++) {
// Find the counterIJ for the symmetric binary interaction
@ -3623,7 +3610,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
m_Phiprime_IJ[counterIJ] = 0.0;
m_PhiPhi_IJ_P[counterIJ] = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
speciesName(i), speciesName(j),
m_Phi_IJ_P[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_P[counterIJ]);
@ -3632,12 +3619,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
}
// ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
debuglog(" Step 7: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 7: \n", m_debugCalc);
double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
double dAphidP = dA_DebyedP /3.0;
double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
+ (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" initial value of dFdP = %10.6f \n", dFdP);
}
for (size_t i = 1; i < m_kk-1; i++) {
@ -3657,12 +3644,12 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
if (charge(i)*charge(j) > 0) {
dFdP += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" dFdP = %10.6f \n", dFdP);
}
}
}
debuglog(" Step 8: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 8: \n", m_debugCalc);
for (size_t i = 1; i < m_kk; i++) {
// -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
@ -3740,7 +3727,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
m_dlnActCoeffMolaldP_Unscaled[i] =
zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s lngamma[i]=%10.6f \n",
speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -3817,7 +3804,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
}
m_dlnActCoeffMolaldP_Unscaled[i] =
zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
writelogf(" %12g %12g %12g %12g %12g %12g\n",
@ -3843,13 +3830,13 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
}
double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
}
}
}
debuglog(" Step 9: \n", DEBUG_MODE_ENABLED && m_debugCalc);
debuglog(" Step 9: \n", m_debugCalc);
// ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
double sum1 = 0.0;
@ -3973,7 +3960,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
} else {
d_osmotic_coef_dP = 0.0;
}
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
"sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
term1, sum1, sum2, sum3, sum4, sum5);
@ -3990,7 +3977,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
// We have just computed act_0. However, this routine returns
// ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
d_lnwateract_dP, exp(d_lnwateract_dP));
}
@ -4007,7 +3994,7 @@ void HMWSoln::calc_lambdas(double is) const
// aphi is the Debye-Huckel constant at 25 C
double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
double aphi = 0.392; /* Value at 25 C */
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" Is = %g\n", is);
}
if (is < 1.0E-150) {
@ -4039,7 +4026,7 @@ void HMWSoln::calc_lambdas(double is) const
elambda[ij] = zprod*jfunc / (4.0*is); // eqn 14
elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
- elambda[ij])/is;
if (DEBUG_MODE_ENABLED && m_debugCalc) {
if (m_debugCalc) {
writelogf(" ij = %d, elambda = %g, elambda1 = %g\n",
ij, elambda[ij], elambda1[ij]);
}
@ -4324,11 +4311,7 @@ doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const
int HMWSoln::debugPrinting()
{
#ifdef DEBUG_MODE
return m_debugCalc;
#else
return 0;
#endif
}
}

52
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -453,15 +453,13 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
NeutralMolecMoleFractions_[k] = 0.0;
}
if (DEBUG_MODE_ENABLED) {
sum = -1.0;
for (size_t k = 0; k < m_kk; k++) {
sum += moleFractions_[k];
}
if (fabs(sum) > 1.0E-11) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"molefracts don't sum to one: {}", sum);
}
sum = -1.0;
for (size_t k = 0; k < m_kk; k++) {
sum += moleFractions_[k];
}
if (fabs(sum) > 1.0E-11) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"molefracts don't sum to one: {}", sum);
}
switch (ionSolnType_) {
@ -494,28 +492,26 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
}
if (DEBUG_MODE_ENABLED) {
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
}
for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
}
for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
}
}
for (size_t k = 0; k < m_kk; k++) {
if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
// Check to see if we have in fact found the inverse.
if (anionList_[0] != k) {
}
for (size_t k = 0; k < m_kk; k++) {
if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
// Check to see if we have in fact found the inverse.
if (anionList_[0] != k) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"neutral molecule calc error");
} else {
// For the single anion case, we will allow some slippage
if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"neutral molecule calc error");
} else {
// For the single anion case, we will allow some slippage
if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"neutral molecule calc error - anion");
}
"neutral molecule calc error - anion");
}
}
}

12
src/thermo/LatticeSolidPhase.cpp
View file

@ -230,13 +230,11 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
// At this point we can check against the mole fraction vector of the
// underlying LatticePhase objects and get the same answer.
if (DEBUG_MODE_ENABLED) {
m_lattice[n]->getMoleFractions(&m_x[strt]);
for (size_t k = 0; k < nsp; k++) {
if (fabs((x + strt)[k] - m_x[strt+k]) > 1.0E-14) {
throw CanteraError("LatticeSolidPhase::getMoleFractions()",
"internal error");
}
m_lattice[n]->getMoleFractions(&m_x[strt]);
for (size_t k = 0; k < nsp; k++) {
if (fabs((x + strt)[k] - m_x[strt+k]) > 1.0E-14) {
throw CanteraError("LatticeSolidPhase::getMoleFractions()",
"internal error");
}
}
strt += nsp;

4
src/thermo/PDSS_HKFT.cpp
View file

@ -205,13 +205,11 @@ doublereal PDSS_HKFT::enthalpy_mole() const
return h;
}
#ifdef DEBUG_MODE
doublereal PDSS_HKFT::enthalpy_mole2() const
{
double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * 1.0E3 * 4.184;
return deltaH() + enthTRPR;
}
#endif
doublereal PDSS_HKFT::intEnergy_mole() const
{
@ -634,7 +632,6 @@ void PDSS_HKFT::constructPDSSFile(VPStandardStateTP* tp, size_t spindex,
constructPDSSXML(tp, spindex, *s, *fxml_phase, true);
}
#ifdef DEBUG_MODE
doublereal PDSS_HKFT::deltaH() const
{
doublereal pbar = m_pres * 1.0E-5;
@ -683,7 +680,6 @@ doublereal PDSS_HKFT::deltaH() const
// Convert to Joules / kmol
return deltaH_calgmol * 1.0E3 * 4.184;
}
#endif
doublereal PDSS_HKFT::deltaG() const
{

3
src/thermo/RedlichKisterVPSSTP.cpp
View file

@ -571,7 +571,6 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
resizeNumInteractions(numBinaryInteractions_);
}
#ifdef DEBUG_MODE
void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut)
{
double XA;
@ -612,5 +611,5 @@ void RedlichKisterVPSSTP::Vint(double& VintOut, double& voltsOut)
VintOut = Volts;
voltsOut = Volts + termp;
}
#endif
}

3
src/thermo/RedlichKwongMFTP.cpp
View file

@ -169,7 +169,6 @@ doublereal RedlichKwongMFTP::cv_mole() const
doublereal RedlichKwongMFTP::pressure() const
{
#ifdef DEBUG_MODE
_updateReferenceStateThermo();
// Get a copy of the private variables stored in the State object
@ -179,7 +178,7 @@ doublereal RedlichKwongMFTP::pressure() const
if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) {
throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe");
}
#endif
return m_Pcurrent;
}

30
src/transport/GasTransport.cpp
View file

@ -440,23 +440,15 @@ void GasTransport::setupMM()
}
// initialize the collision integral calculator for the desired T* range
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** collision_integrals ***\n");
}
debuglog("*** collision_integrals ***\n", m_log_level);
MMCollisionInt integrals;
integrals.init(tstar_min, tstar_max, m_log_level);
fitCollisionIntegrals(integrals);
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** end of collision_integrals ***\n");
}
debuglog("*** end of collision_integrals ***\n", m_log_level);
// make polynomial fits
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** property fits ***\n");
}
debuglog("*** property fits ***\n", m_log_level);
fitProperties(integrals);
if (DEBUG_MODE_ENABLED && m_log_level) {
writelog("*** end of property fits ***\n");
}
debuglog("*** end of property fits ***\n", m_log_level);
}
void GasTransport::getTransportData()
@ -514,7 +506,7 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
// Chemkin fits to sixth order polynomials
int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
if (DEBUG_MODE_ENABLED && m_log_level) {
if (m_log_level) {
writelog("tstar_fits\n"
"fits to A*, B*, and C* vs. log(T*).\n"
"These are done only for the required dstar(j,k) values.\n\n");
@ -578,13 +570,13 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
vector_fp c(degree + 1), c2(degree + 1);
// fit the pure-species viscosity and thermal conductivity for each species
if (DEBUG_MODE_ENABLED && m_log_level && m_log_level < 2) {
if (m_log_level && m_log_level < 2) {
writelog("*** polynomial coefficients not printed (log_level < 2) ***\n");
}
double sqrt_T, visc, err, relerr,
mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
if (DEBUG_MODE_ENABLED && m_log_level) {
if (m_log_level) {
writelog("Polynomial fits for viscosity:\n");
if (m_mode == CK_Mode) {
writelog("log(viscosity) fit to cubic polynomial in log(T)\n");
@ -695,11 +687,11 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
m_visccoeffs.push_back(c);
m_condcoeffs.push_back(c2);
if (DEBUG_MODE_ENABLED && m_log_level >= 2) {
if (m_log_level >= 2) {
writelog(m_thermo->speciesName(k) + ": [" + vec2str(c) + "]\n");
}
}
if (DEBUG_MODE_ENABLED && m_log_level) {
if (m_log_level) {
writelogf("Maximum viscosity absolute error: %12.6g\n", mxerr);
writelogf("Maximum viscosity relative error: %12.6g\n", mxrelerr);
writelog("\nPolynomial fits for conductivity:\n");
@ -775,13 +767,13 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
mxrelerr = std::max(mxrelerr, fabs(relerr));
}
m_diffcoeffs.push_back(c);
if (DEBUG_MODE_ENABLED && m_log_level >= 2) {
if (m_log_level >= 2) {
writelog(m_thermo->speciesName(k) + "__" +
m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
}
}
}
if (DEBUG_MODE_ENABLED && m_log_level) {
if (m_log_level) {
writelogf("Maximum binary diffusion coefficient absolute error:"
" %12.6g\n", mxerr);
writelogf("Maximum binary diffusion coefficient relative error:"

22
src/transport/MMCollisionInt.cpp
View file

@ -216,9 +216,7 @@ double MMCollisionInt::cstar_table[39*8] = {
void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
{
m_loglevel = log_level;
if (DEBUG_MODE_ENABLED && m_loglevel > 0) {
writelog("Collision Integral Polynomial Fits\n");
}
debuglog("Collision Integral Polynomial Fits\n", m_loglevel > 0);
m_nmin = -1;
m_nmax = -1;
@ -234,14 +232,14 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
m_nmin = 0;
m_nmax = 36;
}
if (DEBUG_MODE_ENABLED && m_loglevel > 0) {
if (m_loglevel > 0) {
writelogf("T*_min = %g\n", tstar[m_nmin + 1]);
writelogf("T*_max = %g\n", tstar[m_nmax + 1]);
}
m_logTemp.resize(37);
doublereal rmserr, e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
if (DEBUG_MODE_ENABLED && m_loglevel > 0) {
if (m_loglevel > 0) {
writelog("Collision integral fits at each tabulated T* vs. delta*.\n"
"These polynomial fits are used to interpolate between "
"columns (delta*)\n in the Monchick and Mason tables."
@ -256,7 +254,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
vector_fp c(DeltaDegree+1);
rmserr = fitDelta(0, i, DeltaDegree, c.data());
if (DEBUG_MODE_ENABLED && log_level > 3) {
if (log_level > 3) {
writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]);
writelog("omega22 = [" + vec2str(c) + "]\n");
}
@ -265,26 +263,26 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
rmserr = fitDelta(1, i, DeltaDegree, c.data());
m_apoly.push_back(c);
if (DEBUG_MODE_ENABLED && log_level > 3) {
if (log_level > 3) {
writelog("A* = [" + vec2str(c) + "]\n");
}
ea = std::max(ea, rmserr);
rmserr = fitDelta(2, i, DeltaDegree, c.data());
m_bpoly.push_back(c);
if (DEBUG_MODE_ENABLED && log_level > 3) {
if (log_level > 3) {
writelog("B* = [" + vec2str(c) + "]\n");
}
eb = std::max(eb, rmserr);
rmserr = fitDelta(3, i, DeltaDegree, c.data());
m_cpoly.push_back(c);
if (DEBUG_MODE_ENABLED && log_level > 3) {
if (log_level > 3) {
writelog("C* = [" + vec2str(c) + "]\n");
}
ec = std::max(ec, rmserr);
if (DEBUG_MODE_ENABLED && log_level > 0) {
if (log_level > 0) {
writelogf("max RMS errors in fits vs. delta*:\n"
" omega_22 = %12.6g \n"
" A* = %12.6g \n"
@ -433,7 +431,7 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
}
w[0]= -1.0;
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, o22);
if (DEBUG_MODE_ENABLED && m_loglevel > 0 && rmserr > 0.01) {
if (m_loglevel > 0 && rmserr > 0.01) {
writelogf("Warning: RMS error = %12.6g in omega_22 fit"
"with delta* = %12.6g\n", rmserr, deltastar);
}
@ -477,7 +475,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
}
w[0]= -1.0;
rmserr = polyfit(n, logT, values.data(), w.data(), degree, ndeg, 0.0, c);
if (DEBUG_MODE_ENABLED && m_loglevel > 2) {
if (m_loglevel > 2) {
writelogf("\nT* fit at delta* = %.6g\n", deltastar);
writelog("astar = [" + vec2str(vector_fp(a, a+degree+1))+ "]\n");