From a60217cfc6d54f887e05729aaf0e6781d94e38cd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 3 Apr 2016 17:40:38 -0400 Subject: [PATCH] [Equil] Make better use of local variables --- include/cantera/equil/vcs_internal.h | 2 +- src/equil/BasisOptimize.cpp | 148 ++++++++--------- src/equil/ChemEquil.cpp | 213 +++++++++++-------------- src/equil/MultiPhase.cpp | 200 +++++++++-------------- src/equil/MultiPhaseEquil.cpp | 199 ++++++++++------------- src/equil/vcs_MultiPhaseEquil.cpp | 48 +++--- src/equil/vcs_VolPhase.cpp | 4 +- src/equil/vcs_phaseStability.cpp | 36 ++--- src/equil/vcs_prob.cpp | 2 - src/equil/vcs_rxnadj.cpp | 7 +- src/equil/vcs_setMolesLinProg.cpp | 22 ++- src/equil/vcs_solve_TP.cpp | 96 +++++------ src/equil/vcs_solve_phaseStability.cpp | 10 +- src/equil/vcs_species_thermo.cpp | 18 +-- src/equil/vcs_util.cpp | 7 +- 15 files changed, 412 insertions(+), 600 deletions(-) diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index 9a3121513..b3f649336 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/include/cantera/equil/vcs_internal.h @@ -110,7 +110,7 @@ typedef double(*VCS_FUNC_PTR)(double xval, double Vtarget, * @param vec vector of doubles * @return Returns the l2 norm of the vector */ -double vcs_l2norm(const vector_fp vec); +double vcs_l2norm(const vector_fp& vec); //! Finds the location of the maximum component in a double vector /*! diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index a83da399a..fc5286b8c 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -10,6 +10,7 @@ using namespace std; namespace Cantera { int BasisOptimize_print_lvl = 0; +static const double USEDBEFORE = -1; //! Print a string within a given space limit. /*! @@ -28,40 +29,27 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, std::vector& orderVectorElements, vector_fp& formRxnMatrix) { - size_t j, jj, k=0, kk, i, jl, ml; - std::string ename; - std::string sname; - // Get the total number of elements defined in the multiphase object size_t ne = mphase->nElements(); // Get the total number of species in the multiphase object size_t nspecies = mphase->nSpecies(); - doublereal tmp; - doublereal const USEDBEFORE = -1; // Perhaps, initialize the element ordering if (orderVectorElements.size() < ne) { orderVectorElements.resize(ne); - for (j = 0; j < ne; j++) { - orderVectorElements[j] = j; - } + iota(orderVectorElements.begin(), orderVectorElements.end(), 0); } // Perhaps, initialize the species ordering if (orderVectorSpecies.size() != nspecies) { orderVectorSpecies.resize(nspecies); - for (k = 0; k < nspecies; k++) { - orderVectorSpecies[k] = k; - } + iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0); } if (BasisOptimize_print_lvl >= 1) { writelog(" "); - for (i=0; i<77; i++) { - writelog("-"); - } - writelog("\n"); + writeline('-', 77); writelog(" --- Subroutine BASOPT called to "); writelog("calculate the number of components and "); writelog("evaluate the formation matrix\n"); @@ -70,22 +58,22 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, writelog(" --- Formula Matrix used in BASOPT calculation\n"); writelog(" --- Species | Order | "); - for (j = 0; j < ne; j++) { - jj = orderVectorElements[j]; + for (size_t j = 0; j < ne; j++) { + size_t jj = orderVectorElements[j]; writelog(" "); - ename = mphase->elementName(jj); + std::string ename = mphase->elementName(jj); print_stringTrunc(ename.c_str(), 4, 1); writelogf("(%1d)", j); } writelog("\n"); - for (k = 0; k < nspecies; k++) { - kk = orderVectorSpecies[k]; + for (size_t k = 0; k < nspecies; k++) { + size_t kk = orderVectorSpecies[k]; writelog(" --- "); - sname = mphase->speciesName(kk); + std::string sname = mphase->speciesName(kk); print_stringTrunc(sname.c_str(), 11, 1); writelogf(" | %4d |", k); - for (j = 0; j < ne; j++) { - jj = orderVectorElements[j]; + for (size_t j = 0; j < ne; j++) { + size_t jj = orderVectorElements[j]; double num = mphase->nAtoms(kk,jj); writelogf("%6.1g ", num); } @@ -117,10 +105,8 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, } // For debugging purposes keep an unmodified copy of the array. - vector_fp molNumBase; - molNumBase = molNum; + vector_fp molNumBase = molNum; double molSave = 0.0; - size_t jr = 0; // Top of a loop of some sort based on the index JR. JR is the current @@ -128,11 +114,13 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, while (jr < nComponents) { // Top of another loop point based on finding a linearly independent // species + size_t k = npos; while (true) { // Search the remaining part of the mole number vector, molNum for // the largest remaining species. Return its identity. kk is the raw // number. k is the orderVectorSpecies index. - kk = max_element(molNum.begin(), molNum.end()) - molNum.begin(); + size_t kk = max_element(molNum.begin(), molNum.end()) - molNum.begin(); + size_t j; for (j = 0; j < nspecies; j++) { if (orderVectorSpecies[j] == kk) { k = j; @@ -162,9 +150,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Modified Gram-Schmidt Method, p. 202 Dalquist // QR factorization of a matrix without row pivoting. - jl = jr; + size_t jl = jr; for (j = 0; j < ne; ++j) { - jj = orderVectorElements[j]; + size_t jj = orderVectorElements[j]; sm[j + jr*ne] = mphase->nAtoms(kk,jj); } if (jl > 0) { @@ -173,7 +161,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // different than Dalquist) R_JA_JA = 1 for (j = 0; j < jl; ++j) { ss[j] = 0.0; - for (i = 0; i < ne; ++i) { + for (size_t i = 0; i < ne; ++i) { ss[j] += sm[i + jr*ne] * sm[i + j*ne]; } ss[j] /= sa[j]; @@ -191,8 +179,8 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Find the new length of the new column in Q. // It will be used in the denominator in future row calcs. sa[jr] = 0.0; - for (ml = 0; ml < ne; ++ml) { - tmp = sm[ml + jr*ne]; + for (size_t ml = 0; ml < ne; ++ml) { + double tmp = sm[ml + jr*ne]; sa[jr] += tmp * tmp; } @@ -205,9 +193,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // REARRANGE THE DATA if (jr != k) { if (BasisOptimize_print_lvl >= 1) { - kk = orderVectorSpecies[k]; + size_t kk = orderVectorSpecies[k]; writelogf(" --- %-12.12s", mphase->speciesName(kk)); - jj = orderVectorSpecies[jr]; + size_t jj = orderVectorSpecies[jr]; writelogf("(%9.2g) replaces %-12.12s", molSave, mphase->speciesName(jj)); writelogf("(%9.2g) as component %3d\n", molNum[jj], jr); @@ -253,19 +241,19 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Note the rearrangement of elements need only be done once in the problem. // It's actually very similar to the top of this program with ne being the // species and nc being the elements!! - for (k = 0; k < nComponents; ++k) { - kk = orderVectorSpecies[k]; - for (j = 0; j < nComponents; ++j) { - jj = orderVectorElements[j]; + for (size_t k = 0; k < nComponents; ++k) { + size_t kk = orderVectorSpecies[k]; + for (size_t j = 0; j < nComponents; ++j) { + size_t jj = orderVectorElements[j]; sm[j + k*ne] = mphase->nAtoms(kk, jj); } } - for (i = 0; i < nNonComponents; ++i) { - k = nComponents + i; - kk = orderVectorSpecies[k]; - for (j = 0; j < nComponents; ++j) { - jj = orderVectorElements[j]; + for (size_t i = 0; i < nNonComponents; ++i) { + size_t k = nComponents + i; + size_t kk = orderVectorSpecies[k]; + for (size_t j = 0; j < nComponents; ++j) { + size_t jj = orderVectorElements[j]; formRxnMatrix[j + i * ne] = - mphase->nAtoms(kk, jj); } } @@ -284,39 +272,36 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, writelogf(" --- Number of Components = %d\n", nComponents); writelog(" --- Formula Matrix:\n"); writelog(" --- Components: "); - for (k = 0; k < nComponents; k++) { - kk = orderVectorSpecies[k]; + for (size_t k = 0; k < nComponents; k++) { + size_t kk = orderVectorSpecies[k]; writelogf(" %3d (%3d) ", k, kk); } writelog("\n --- Components Moles: "); - for (k = 0; k < nComponents; k++) { - kk = orderVectorSpecies[k]; + for (size_t k = 0; k < nComponents; k++) { + size_t kk = orderVectorSpecies[k]; writelogf("%-11.3g", molNumBase[kk]); } writelog("\n --- NonComponent | Moles | "); - for (i = 0; i < nComponents; i++) { - kk = orderVectorSpecies[i]; + for (size_t i = 0; i < nComponents; i++) { + size_t kk = orderVectorSpecies[i]; writelogf("%-11.10s", mphase->speciesName(kk)); } writelog("\n"); - for (i = 0; i < nNonComponents; i++) { - k = i + nComponents; - kk = orderVectorSpecies[k]; + for (size_t i = 0; i < nNonComponents; i++) { + size_t k = i + nComponents; + size_t kk = orderVectorSpecies[k]; writelogf(" --- %3d (%3d) ", k, kk); writelogf("%-10.10s", mphase->speciesName(kk)); writelogf("|%10.3g|", molNumBase[kk]); // Print the negative of formRxnMatrix[]; it's easier to interpret. - for (j = 0; j < nComponents; j++) { + for (size_t j = 0; j < nComponents; j++) { writelogf(" %6.2f", - formRxnMatrix[j + i * ne]); } writelog("\n"); } writelog(" "); - for (i=0; i<77; i++) { - writelog("-"); - } - writelog("\n"); + writeline('-', 77); } return nComponents; @@ -368,19 +353,13 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, std::vector& orderVectorSpecies, std::vector& orderVectorElements) { - size_t j, k, i, jl, ml, jr, ielem, jj, kk=0; size_t nelements = mphase->nElements(); - std::string ename; // Get the total number of species in the multiphase object size_t nspecies = mphase->nSpecies(); - double test = -1.0E10; if (BasisOptimize_print_lvl > 0) { writelog(" "); - for (i=0; i<77; i++) { - writelog("-"); - } - writelog("\n"); + writeline('-', 77); writelog(" --- Subroutine ElemRearrange() called to "); writelog("check stoich. coefficient matrix\n"); writelog(" --- and to rearrange the element ordering once\n"); @@ -389,7 +368,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Perhaps, initialize the element ordering if (orderVectorElements.size() < nelements) { orderVectorElements.resize(nelements); - for (j = 0; j < nelements; j++) { + for (size_t j = 0; j < nelements; j++) { orderVectorElements[j] = j; } } @@ -398,7 +377,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // this ordering is assumed to yield the component species for the problem if (orderVectorSpecies.size() != nspecies) { orderVectorSpecies.resize(nspecies); - for (k = 0; k < nspecies; k++) { + for (size_t k = 0; k < nspecies; k++) { orderVectorSpecies[k] = k; } } @@ -408,9 +387,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // zero species to the end of the element ordering. vector_fp eAbund(nelements,0.0); if (elementAbundances.size() != nelements) { - for (j = 0; j < nelements; j++) { + for (size_t j = 0; j < nelements; j++) { eAbund[j] = 0.0; - for (k = 0; k < nspecies; k++) { + for (size_t k = 0; k < nspecies; k++) { eAbund[j] += fabs(mphase->nAtoms(k, j)); } } @@ -425,25 +404,26 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Top of a loop of some sort based on the index JR. JR is the current // number independent elements found. - jr = 0; + size_t jr = 0; while (jr < nComponents) { // Top of another loop point based on finding a linearly independent // element + size_t k = nelements; while (true) { // Search the element vector. We first locate elements that are // present in any amount. Then, we locate elements that are not // present in any amount. Return its identity in K. - k = nelements; - for (ielem = jr; ielem < nelements; ielem++) { + size_t kk; + for (size_t ielem = jr; ielem < nelements; ielem++) { kk = orderVectorElements[ielem]; - if (eAbund[kk] != test && eAbund[kk] > 0.0) { + if (eAbund[kk] != USEDBEFORE && eAbund[kk] > 0.0) { k = ielem; break; } } - for (ielem = jr; ielem < nelements; ielem++) { + for (size_t ielem = jr; ielem < nelements; ielem++) { kk = orderVectorElements[ielem]; - if (eAbund[kk] != test) { + if (eAbund[kk] != USEDBEFORE) { k = ielem; break; } @@ -460,21 +440,21 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Assign a large negative number to the element that we have // just found, in order to take it out of further consideration. - eAbund[kk] = test; + eAbund[kk] = USEDBEFORE; // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX // LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX // Modified Gram-Schmidt Method, p. 202 Dalquist // QR factorization of a matrix without row pivoting. - jl = jr; + size_t jl = jr; // Fill in the row for the current element, k, under consideration // The row will contain the Formula matrix value for that element // with respect to the vector of component species. (note j and k // indices are flipped compared to the previous routine) - for (j = 0; j < nComponents; ++j) { - jj = orderVectorSpecies[j]; + for (size_t j = 0; j < nComponents; ++j) { + size_t jj = orderVectorSpecies[j]; kk = orderVectorElements[k]; sm[j + jr*nComponents] = mphase->nAtoms(jj,kk); } @@ -482,9 +462,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Compute the coefficients of JA column of the the upper // triangular R matrix, SS(J) = R_J_JR (this is slightly // different than Dalquist) R_JA_JA = 1 - for (j = 0; j < jl; ++j) { + for (size_t j = 0; j < jl; ++j) { ss[j] = 0.0; - for (i = 0; i < nComponents; ++i) { + for (size_t i = 0; i < nComponents; ++i) { ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents]; } ss[j] /= sa[j]; @@ -492,7 +472,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Now make the new column, (*,JR), orthogonal to the // previous columns - for (j = 0; j < jl; ++j) { + for (size_t j = 0; j < jl; ++j) { for (size_t i = 0; i < nComponents; ++i) { sm[i + jr*nComponents] -= ss[j] * sm[i + j*nComponents]; } @@ -502,7 +482,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // Find the new length of the new column in Q. // It will be used in the denominator in future row calcs. sa[jr] = 0.0; - for (ml = 0; ml < nComponents; ++ml) { + for (size_t ml = 0; ml < nComponents; ++ml) { double tmp = sm[ml + jr*nComponents]; sa[jr] += tmp * tmp; } @@ -514,7 +494,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances, // REARRANGE THE DATA if (jr != k) { if (BasisOptimize_print_lvl > 0) { - kk = orderVectorElements[k]; + size_t kk = orderVectorElements[k]; writelog(" --- "); writelogf("%-2.2s", mphase->elementName(kk)); writelog("replaces "); diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index bd5110b4b..49655d0ea 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -96,15 +96,12 @@ void ChemEquil::initialize(thermo_t& s) // set up elemental composition matrix size_t mneg = npos; - doublereal na, ewt; for (size_t m = 0; m < m_mm; m++) { for (size_t k = 0; k < m_kk; k++) { - na = s.nAtoms(k,m); - // handle the case of negative atom numbers (used to // represent positive ions, where the 'element' is an // electron - if (na < 0.0) { + if (s.nAtoms(k,m) < 0.0) { // if negative atom numbers have already been specified // for some element other than this one, throw // an exception @@ -113,11 +110,10 @@ void ChemEquil::initialize(thermo_t& s) "negative atom numbers allowed for only one element"); } mneg = m; - ewt = s.atomicWeight(m); // the element should be an electron... if it isn't // print a warning. - if (ewt > 1.0e-3) { + if (s.atomicWeight(m) > 1.0e-3) { writelog("WARNING: species {} has {} atoms of element {}," " but this element is not an electron.\n", s.speciesName(k), s.nAtoms(k,m), s.elementName(m)); @@ -256,18 +252,16 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, m_orderVectorSpecies, m_orderVectorElements); s.getChemPotentials(mu_RT.data()); - doublereal rrt = 1.0/(GasConstant* s.temperature()); - scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt); + scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), + 1.0/(GasConstant* s.temperature())); 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)); } - double pres = s.pressure(); - double temp = s.temperature(); - writelogf("Pressure = %g\n", pres); - writelogf("Temperature = %g\n", temp); + writelogf("Pressure = %g\n", s.pressure()); + writelogf("Temperature = %g\n", s.temperature()); writelog(" id Name MF mu/RT \n"); for (size_t n = 0; n < s.nSpecies(); n++) { writelogf("%10d %15s %10.5g %10.5g\n", @@ -316,11 +310,9 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, int ChemEquil::equilibrate(thermo_t& s, const char* XY, bool useThermoPhaseElementPotentials, int loglevel) { - vector_fp elMolesGoal(s.nElements()); initialize(s); update(s); - copy(m_elementmolefracs.begin(), m_elementmolefracs.end(), - elMolesGoal.begin()); + vector_fp elMolesGoal = m_elementmolefracs; return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials, loglevel-1); } @@ -330,7 +322,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, bool useThermoPhaseElementPotentials, int loglevel) { - doublereal xval, yval, tmp; int fail = 0; bool tempFixed = true; int XY = _equilflag(XYstr); @@ -396,8 +387,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // Before we do anything to change the ThermoPhase object, we calculate and // store the two specified thermodynamic properties that we are after. - xval = m_p1->value(s); - yval = m_p2->value(s); + double xval = m_p1->value(s); + double yval = m_p2->value(s); size_t mm = m_mm; size_t nvar = mm + 1; @@ -409,10 +400,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // specified property calculation. // // We choose the equation of the element with the highest element abundance. - size_t m; - tmp = -1.0; + double tmp = -1.0; for (size_t im = 0; im < m_nComponents; im++) { - m = m_orderVectorElements[im]; + size_t m = m_orderVectorElements[im]; if (elMolesGoal[m] > tmp) { m_skip = m; tmp = elMolesGoal[m]; @@ -520,7 +510,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (s.temperature() < 100.) { writelog("we are here {:g}\n", s.temperature()); } - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { x[m] *= 1.0 / s.RT(); } } else { @@ -555,7 +545,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // containing that element. vector_fp above(nvar); vector_fp below(nvar); - for (m = 0; m < mm; m++) { + for (size_t m = 0; m < mm; m++) { above[m] = 200.0; below[m] = -2000.0; if (elMolesGoal[m] < m_elemFracCutoff && m != m_eloc) { @@ -571,20 +561,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp grad(nvar, 0.0); // gradient of f = F*F/2 vector_fp oldx(nvar, 0.0); // old solution vector_fp oldresid(nvar, 0.0); - doublereal f, oldf; - doublereal fctr = 1.0, newval; for (int iter = 0; iter < options.maxIterations; iter++) { // check for convergence. equilResidual(s, x, elMolesGoal, res_trial, xval, yval); - f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); - doublereal xx, yy, deltax, deltay; - xx = m_p1->value(s); - yy = m_p2->value(s); - deltax = (xx - xval)/xval; - deltay = (yy - yval)/yval; + double f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); + double xx = m_p1->value(s); + double yy = m_p2->value(s); + double deltax = (xx - xval)/xval; + double deltay = (yy - yval)/yval; bool passThis = true; - for (m = 0; m < nvar; m++) { + for (size_t m = 0; m < nvar; m++) { double tval = options.relTolerance; if (m < mm) { // Special case convergence requirements for electron element. @@ -612,7 +599,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (iter > 0 && passThis && fabs(deltax) < options.relTolerance && fabs(deltay) < options.relTolerance) { options.iterations = iter; - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { m_lambda[m] = x[m]* s.RT(); } @@ -641,7 +628,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (ChemEquil_print_lvl > 0) { writelogf("Jacobian matrix %d:\n", iter); - for (m = 0; m <= m_mm; m++) { + for (size_t m = 0; m <= m_mm; m++) { writelog(" [ "); for (size_t n = 0; n <= m_mm; n++) { writelog("{:10.5g} ", jac(m,n)); @@ -661,7 +648,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, } oldx = x; - oldf = f; + double oldf = f; scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0); // Solve the system @@ -677,9 +664,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, // find the factor by which the Newton step can be multiplied // to keep the solution within bounds. - fctr = 1.0; - for (m = 0; m < nvar; m++) { - newval = x[m] + res_trial[m]; + double fctr = 1.0; + for (size_t m = 0; m < nvar; m++) { + double newval = x[m] + res_trial[m]; if (newval > above[m]) { fctr = std::max(0.0, std::min(fctr,0.8*(above[m] - x[m])/(newval - x[m]))); @@ -729,11 +716,9 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, double& f, vector_fp& elmols, double xval, double yval) { - double damp; - // Carry out a delta damping approach on the dimensionless element // potentials. - damp = 1.0; + double damp = 1.0; for (size_t m = 0; m < m_mm; m++) { if (m == m_eloc) { if (step[m] > 1.25) { @@ -770,9 +755,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, const vector_fp& elmFracGoal, vector_fp& resid, doublereal xval, doublereal yval, int loglevel) { - doublereal xx, yy; - doublereal temp = exp(x[m_mm]); - setToEquilState(s, x, temp); + setToEquilState(s, x, exp(x[m_mm])); // residuals are the total element moles vector_fp& elmFrac = m_elementmolefracs; @@ -802,8 +785,8 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, } } - xx = m_p1->value(s); - yy = m_p2->value(s); + double xx = m_p1->value(s); + double yy = m_p2->value(s); resid[m_mm] = xx/xval - 1.0; resid[m_skip] = yy/yval - 1.0; @@ -823,25 +806,23 @@ void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, size_t len = x.size(); r0.resize(len); r1.resize(len); - size_t n, m; - doublereal rdx, dx, xsave; doublereal atol = 1.e-10; equilResidual(s, x, elmols, r0, xval, yval, loglevel-1); m_doResPerturb = false; - for (n = 0; n < len; n++) { - xsave = x[n]; - dx = std::max(atol, fabs(xsave) * 1.0E-7); + for (size_t n = 0; n < len; n++) { + double xsave = x[n]; + double dx = std::max(atol, fabs(xsave) * 1.0E-7); x[n] = xsave + dx; dx = x[n] - xsave; - rdx = 1.0/dx; + double rdx = 1.0/dx; // calculate perturbed residual equilResidual(s, x, elmols, r1, xval, yval, loglevel-1); // compute nth column of Jacobian - for (m = 0; m < x.size(); m++) { + for (size_t m = 0; m < x.size(); m++) { jac(m, n) = (r1[m] - r0[m])*rdx; } x[n] = xsave; @@ -855,7 +836,6 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, double pressureConst) { double n_t_calc = 0.0; - double tmp; // Calculate the activity coefficients of the solution, at the previous // solution state. @@ -865,7 +845,7 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, s.getActivityCoefficients(actCoeff.data()); for (size_t k = 0; k < m_kk; k++) { - tmp = - (m_muSS_RT[k] + log(actCoeff[k])); + double tmp = - (m_muSS_RT[k] + log(actCoeff[k])); for (size_t m = 0; m < m_mm; m++) { tmp += nAtoms(k,m) * x[m]; } @@ -893,11 +873,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // things go drastically wrong, we will restore the saved state. vector_fp state; s.saveState(state); - double tmp, sum; bool modifiedMatrix = false; size_t neq = m_mm+1; int retn = 1; - size_t m, n, k, im; DenseMatrix a1(neq, neq, 0.0); vector_fp b(neq, 0.0); vector_fp n_i(m_kk,0.0); @@ -920,14 +898,14 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, vector_fp eMolesCalc(m_mm, 0.0); vector_fp eMolesFix(m_mm, 0.0); double elMolesTotal = 0.0; - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { elMolesTotal += elMoles[m]; - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { eMolesFix[m] += nAtoms(k,m) * n_i[k]; } } - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { if (elMoles[m] > 1.0E-70) { x[m] = clip(x[m], -100.0, 50.0); } else { @@ -936,16 +914,15 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } double n_t = 0.0; - double sum2 = 0.0; double nAtomsMax = 1.0; s.setMoleFractions(Xmol_i_calc.data()); s.setPressure(pressureConst); s.getActivityCoefficients(actCoeff.data()); - for (k = 0; k < m_kk; k++) { - tmp = - (m_muSS_RT[k] + log(actCoeff[k])); - sum2 = 0.0; - for (m = 0; m < m_mm; m++) { - sum = nAtoms(k,m); + for (size_t k = 0; k < m_kk; k++) { + double tmp = - (m_muSS_RT[k] + log(actCoeff[k])); + double sum2 = 0.0; + for (size_t m = 0; m < m_mm; m++) { + double sum = nAtoms(k,m); tmp += sum * x[m]; sum2 += sum; nAtomsMax = std::max(nAtomsMax, sum2); @@ -959,25 +936,23 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, if (ChemEquil_print_lvl > 0) { writelog("estimateEP_Brinkley::\n\n"); - double temp = s.temperature(); - double pres = s.pressure(); - writelogf("temp = %g\n", temp); - writelogf("pres = %g\n", pres); + writelogf("temp = %g\n", s.temperature()); + writelogf("pres = %g\n", s.pressure()); writelog("Initial mole numbers and mu_SS:\n"); writelog(" Name MoleNum mu_SS actCoeff\n"); - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { writelogf("%15s %13.5g %13.5g %13.5g\n", s.speciesName(k), n_i[k], m_muSS_RT[k], actCoeff[k]); } writelogf("Initial n_t = %10.5g\n", n_t); writelog("Comparison of Goal Element Abundance with Initial Guess:\n"); writelog(" eName eCurrent eGoal\n"); - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { writelogf("%5s %13.5g %13.5g\n", s.elementName(m), eMolesFix[m], elMoles[m]); } } - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { if (m != m_eloc && elMoles[m] <= options.absElemTol) { x[m] = -200.; } @@ -986,7 +961,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // Main Loop. for (int iter = 0; iter < 20* options.maxIterations; iter++) { // Save the old solution - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { x_old[m] = x[m]; } x_old[m_mm] = n_t; @@ -998,29 +973,28 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc, pressureConst); - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; } if (ChemEquil_print_lvl > 0) { writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { writelogf("%15s: %10.5g %10.5g\n", s.speciesName(k), n_i_calc[k], Xmol_i_calc[k]); } writelogf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); writelogf("(iter %d) element moles bal: Goal Calculated\n", iter); - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { writelogf(" %8s: %10.5g %10.5g \n", s.elementName(m), elMoles[m], eMolesCalc[m]); } } - double nCutoff; bool normalStep = true; // Decide if we are to do a normal step or a modified step size_t iM = npos; - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { if (elMoles[m] > 0.001 * elMolesTotal) { if (eMolesCalc[m] > 1000. * elMoles[m]) { normalStep = false; @@ -1038,8 +1012,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, if (!normalStep) { beta = 1.0; resid[m_mm] = 0.0; - for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; + for (size_t im = 0; im < m_mm; im++) { + size_t m = m_orderVectorElements[im]; resid[m] = 0.0; if (im < m_nComponents && elMoles[m] > 0.001 * elMolesTotal) { if (eMolesCalc[m] > 1000. * elMoles[m]) { @@ -1082,25 +1056,25 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // Jordan factorization scheme in the future. For Example the scheme // below would fail for the set: HCl NH4Cl, NH3. Hopefully, it's // caught by the equal rows logic below. - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { lumpSum[m] = 1; } - nCutoff = 1.0E-9 * n_t_calc; + double nCutoff = 1.0E-9 * n_t_calc; if (ChemEquil_print_lvl > 0) { writelog(" Lump Sum Elements Calculation: \n"); } - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { size_t kMSp = npos; size_t kMSp2 = npos; int nSpeciesWithElem = 0; - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { if (n_i_calc[k] > nCutoff && fabs(nAtoms(k,m)) > 0.001) { nSpeciesWithElem++; if (kMSp != npos) { kMSp2 = k; double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m)); - for (n = 0; n < m_mm; n++) { + for (size_t n = 0; n < m_mm; n++) { if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { lumpSum[m] = 0; break; @@ -1118,19 +1092,19 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } // Formulate the matrix. - for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; + for (size_t im = 0; im < m_mm; im++) { + size_t m = m_orderVectorElements[im]; if (im < m_nComponents) { - for (n = 0; n < m_mm; n++) { + for (size_t n = 0; n < m_mm; n++) { a1(m,n) = 0.0; - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k]; } } a1(m,m_mm) = eMolesCalc[m]; a1(m_mm, m) = eMolesCalc[m]; } else { - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { a1(m,n) = 0.0; } a1(m,m) = 1.0; @@ -1140,9 +1114,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // Formulate the residual, resid, and the estimate for the // convergence criteria, sum - sum = 0.0; - for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; + double sum = 0.0; + for (size_t im = 0; im < m_mm; im++) { + size_t m = m_orderVectorElements[im]; if (im < m_nComponents) { resid[m] = elMoles[m] - eMolesCalc[m]; } else { @@ -1154,6 +1128,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // criteria by a condition limited by finite precision of // inverting a matrix. Other equations with just positive // coefficients aren't limited by this. + double tmp; if (m == m_eloc) { tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol); } else { @@ -1162,12 +1137,12 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, sum += tmp * tmp; } - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { if (a1(m,m) < 1.0E-50) { if (ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m); } - for (n = 0; n < m_mm; n++) { + for (size_t n = 0; n < m_mm; n++) { a1(m,n) = 0.0; } a1(m,m) = 1.0; @@ -1185,17 +1160,16 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, if (ChemEquil_print_lvl > 0) { writelog("Matrix:\n"); - for (m = 0; m <= m_mm; m++) { + for (size_t m = 0; m <= m_mm; m++) { writelog(" ["); - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { writelogf(" %10.5g", a1(m,n)); } writelogf("] = %10.5g\n", resid[m]); } } - tmp = resid[m_mm] /(n_t + 1.0E-15); - sum += tmp * tmp; + sum += pow(resid[m_mm] /(n_t + 1.0E-15), 2); if (ChemEquil_print_lvl > 0) { writelogf("(it %d) Convergence = %g\n", iter, sum); } @@ -1210,16 +1184,16 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } // Row Sum scaling - for (m = 0; m <= m_mm; m++) { - tmp = 0.0; - for (n = 0; n <= m_mm; n++) { + for (size_t m = 0; m <= m_mm; m++) { + double tmp = 0.0; + for (size_t n = 0; n <= m_mm; n++) { tmp += fabs(a1(m,n)); } if (m < m_mm && tmp < 1.0E-30) { if (ChemEquil_print_lvl > 0) { writelogf(" NOTE: Diagonalizing row %d\n", m); } - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { if (n != m) { a1(m,n) = 0.0; a1(n,m) = 0.0; @@ -1227,7 +1201,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } tmp = 1.0/tmp; - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { a1(m,n) *= tmp; } resid[m] *= tmp; @@ -1235,9 +1209,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, if (ChemEquil_print_lvl > 0) { writelog("Row Summed Matrix:\n"); - for (m = 0; m <= m_mm; m++) { + for (size_t m = 0; m <= m_mm; m++) { writelog(" ["); - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { writelogf(" %10.5g", a1(m,n)); } writelogf("] = %10.5g\n", resid[m]); @@ -1264,11 +1238,11 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // 1.0E-3. If two rows are anywhere close to being equivalent, the // algorithm can get stuck in an oscillatory mode. modifiedMatrix = false; - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { size_t sameAsRow = npos; for (size_t im = 0; im < m; im++) { bool theSame = true; - for (n = 0; n < m_mm; n++) { + for (size_t n = 0; n < m_mm; n++) { if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) { theSame = false; break; @@ -1287,7 +1261,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } modifiedMatrix = true; - for (n = 0; n < m_mm; n++) { + for (size_t n = 0; n < m_mm; n++) { if (n != m) { a1(m,m) += fabs(a1(m,n)); a1(m,n) = 0.0; @@ -1298,9 +1272,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, if (ChemEquil_print_lvl > 0 && modifiedMatrix) { writelog("Row Summed, MODIFIED Matrix:\n"); - for (m = 0; m <= m_mm; m++) { + for (size_t m = 0; m <= m_mm; m++) { writelog(" ["); - for (n = 0; n <= m_mm; n++) { + for (size_t n = 0; n <= m_mm; n++) { writelogf(" %10.5g", a1(m,n)); } writelogf("] = %10.5g\n", resid[m]); @@ -1320,7 +1294,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, // Figure out the damping coefficient: Use a delta damping // coefficient formulation: magnitude of change is capped to exp(1). beta = 1.0; - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { if (resid[m] > 1.0) { beta = std::min(beta, 1.0 / resid[m]); } @@ -1333,14 +1307,14 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } // Update the solution vector - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { x[m] += beta * resid[m]; } n_t *= exp(beta * resid[m_mm]); if (ChemEquil_print_lvl > 0) { writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); - for (m = 0; m < m_mm; m++) { + for (size_t m = 0; m < m_mm; m++) { writelogf(" %5s %10.5g %10.5g %10.5g\n", s.elementName(m), x_old[m], x[m], resid[m]); } @@ -1372,13 +1346,12 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) return; } s.getMoleFractions(m_molefractions.data()); - size_t k; size_t maxPosEloc = npos; size_t maxNegEloc = npos; double maxPosVal = -1.0; double maxNegVal = -1.0; if (ChemEquil_print_lvl > 0) { - for (k = 0; k < m_kk; k++) { + for (size_t 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]; maxPosEloc = k; @@ -1392,7 +1365,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) double sumPos = 0.0; double sumNeg = 0.0; - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0) { sumPos += nAtoms(k,m_eloc) * m_molefractions[k]; } @@ -1412,7 +1385,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) s.speciesName(maxPosEloc), m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); } - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0) { m_molefractions[k] *= factor; } @@ -1424,7 +1397,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) s.speciesName(maxNegEloc), m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); } - for (k = 0; k < m_kk; k++) { + for (size_t k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) < 0.0) { m_molefractions[k] *= factor; } diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index f5029750d..456616f65 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -68,8 +68,7 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right) void MultiPhase::addPhases(MultiPhase& mix) { - size_t n; - for (n = 0; n < mix.nPhases(); n++) { + for (size_t n = 0; n < mix.nPhases(); n++) { addPhase(mix.m_phase[n], mix.m_moles[n]); } } @@ -77,9 +76,7 @@ void MultiPhase::addPhases(MultiPhase& mix) void MultiPhase::addPhases(std::vector& phases, const vector_fp& phaseMoles) { - size_t np = phases.size(); - size_t n; - for (n = 0; n < np; n++) { + for (size_t n = 0; n < phases.size(); n++) { addPhase(phases[n], phaseMoles[n]); } init(); @@ -105,11 +102,9 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles) // determine if this phase has new elements for each new element, add an // entry in the map from names to index number + 1: - string ename; // iterate over the elements in this phase - size_t m, nel = p->nElements(); - for (m = 0; m < nel; m++) { - ename = p->elementName(m); + for (size_t m = 0; m < p->nElements(); m++) { + string ename = p->elementName(m); // if no entry is found for this element name, then it is a new element. // In this case, add the name to the list of names, increment the @@ -154,9 +149,6 @@ void MultiPhase::init() if (m_init) { return; } - size_t ip, kp, k = 0, nsp, m; - size_t mlocal; - string sym; // allocate space for the atomic composition matrix m_atoms.resize(m_nel, m_nsp, 0.0); @@ -165,15 +157,14 @@ void MultiPhase::init() // iterate over the elements // -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), m_spstart(ip) - for (m = 0; m < m_nel; m++) { - sym = m_enames[m]; - k = 0; + for (size_t m = 0; m < m_nel; m++) { + size_t k = 0; // iterate over the phases - for (ip = 0; ip < nPhases(); ip++) { + for (size_t ip = 0; ip < nPhases(); ip++) { ThermoPhase* p = m_phase[ip]; - nsp = p->nSpecies(); - mlocal = p->elementIndex(sym); - for (kp = 0; kp < nsp; kp++) { + size_t nsp = p->nSpecies(); + size_t mlocal = p->elementIndex(m_enames[m]); + for (size_t kp = 0; kp < nsp; kp++) { if (mlocal != npos) { m_atoms(m, k) = p->nAtoms(kp, mlocal); } @@ -189,18 +180,6 @@ void MultiPhase::init() } } - if (m_eloc != npos) { - doublereal esum; - for (k = 0; k < m_nsp; k++) { - esum = 0.0; - for (m = 0; m < m_nel; m++) { - if (m != m_eloc) { - esum += m_atoms(m,k) * m_atomicNumber[m]; - } - } - } - } - // set the initial composition within each phase to the // mole fractions stored in the phase objects m_init = true; @@ -241,13 +220,12 @@ doublereal MultiPhase::speciesMoles(size_t k) const doublereal MultiPhase::elementMoles(size_t m) const { - doublereal sum = 0.0, phasesum; - size_t i, k = 0, ik, nsp; - for (i = 0; i < nPhases(); i++) { - phasesum = 0.0; - nsp = m_phase[i]->nSpecies(); - for (ik = 0; ik < nsp; ik++) { - k = speciesIndex(ik, i); + doublereal sum = 0.0; + for (size_t i = 0; i < nPhases(); i++) { + double phasesum = 0.0; + size_t nsp = m_phase[i]->nSpecies(); + for (size_t ik = 0; ik < nsp; ik++) { + size_t k = speciesIndex(ik, i); phasesum += m_atoms(m,k)*m_moleFractions[k]; } sum += phasesum * m_moles[i]; @@ -258,8 +236,7 @@ doublereal MultiPhase::elementMoles(size_t m) const doublereal MultiPhase::charge() const { doublereal sum = 0.0; - size_t i; - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { sum += phaseCharge(i); } return sum; @@ -284,9 +261,9 @@ size_t MultiPhase::speciesIndex(const std::string& speciesName, const std::strin doublereal MultiPhase::phaseCharge(size_t p) const { doublereal phasesum = 0.0; - size_t ik, k, nsp = m_phase[p]->nSpecies(); - for (ik = 0; ik < nsp; ik++) { - k = speciesIndex(ik, p); + size_t nsp = m_phase[p]->nSpecies(); + for (size_t ik = 0; ik < nsp; ik++) { + size_t k = speciesIndex(ik, p); phasesum += m_phase[p]->charge(ik)*m_moleFractions[k]; } return Faraday*phasesum*m_moles[p]; @@ -294,9 +271,9 @@ doublereal MultiPhase::phaseCharge(size_t p) const void MultiPhase::getChemPotentials(doublereal* mu) const { - size_t i, loc = 0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + size_t loc = 0; + for (size_t i = 0; i < nPhases(); i++) { m_phase[i]->getChemPotentials(mu + loc); loc += m_phase[i]->nSpecies(); } @@ -305,10 +282,10 @@ void MultiPhase::getChemPotentials(doublereal* mu) const void MultiPhase::getValidChemPotentials(doublereal not_mu, doublereal* mu, bool standard) const { - size_t i, loc = 0; updatePhases(); // iterate over the phases - for (i = 0; i < nPhases(); i++) { + size_t loc = 0; + for (size_t i = 0; i < nPhases(); i++) { if (tempOK(i) || m_phase[i]->nSpecies() > 1) { if (!standard) { m_phase[i]->getChemPotentials(mu + loc); @@ -333,10 +310,9 @@ bool MultiPhase::solutionSpecies(size_t k) const doublereal MultiPhase::gibbs() const { - size_t i; doublereal sum = 0.0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { if (m_moles[i] > 0.0) { sum += m_phase[i]->gibbs_mole() * m_moles[i]; } @@ -346,10 +322,9 @@ doublereal MultiPhase::gibbs() const doublereal MultiPhase::enthalpy() const { - size_t i; doublereal sum = 0.0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { if (m_moles[i] > 0.0) { sum += m_phase[i]->enthalpy_mole() * m_moles[i]; } @@ -359,10 +334,9 @@ doublereal MultiPhase::enthalpy() const doublereal MultiPhase::IntEnergy() const { - size_t i; doublereal sum = 0.0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { if (m_moles[i] > 0.0) { sum += m_phase[i]->intEnergy_mole() * m_moles[i]; } @@ -372,10 +346,9 @@ doublereal MultiPhase::IntEnergy() const doublereal MultiPhase::entropy() const { - size_t i; doublereal sum = 0.0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { if (m_moles[i] > 0.0) { sum += m_phase[i]->entropy_mole() * m_moles[i]; } @@ -385,10 +358,9 @@ doublereal MultiPhase::entropy() const doublereal MultiPhase::cp() const { - size_t i; doublereal sum = 0.0; updatePhases(); - for (i = 0; i < nPhases(); i++) { + for (size_t i = 0; i < nPhases(); i++) { if (m_moles[i] > 0.0) { sum += m_phase[i]->cp_mole() * m_moles[i]; } @@ -430,13 +402,12 @@ void MultiPhase::getMoles(doublereal* molNum) const { // First copy in the mole fractions copy(m_moleFractions.begin(), m_moleFractions.end(), molNum); - size_t ik; doublereal* dtmp = molNum; for (size_t ip = 0; ip < nPhases(); ip++) { doublereal phasemoles = m_moles[ip]; ThermoPhase* p = m_phase[ip]; size_t nsp = p->nSpecies(); - for (ik = 0; ik < nsp; ik++) { + for (size_t ik = 0; ik < nsp; ik++) { *(dtmp++) *= phasemoles; } } @@ -447,14 +418,13 @@ void MultiPhase::setMoles(const doublereal* n) if (!m_init) { init(); } - size_t ip, loc = 0; - size_t ik, k = 0, nsp; - doublereal phasemoles; - for (ip = 0; ip < nPhases(); ip++) { + size_t loc = 0; + size_t k = 0; + for (size_t ip = 0; ip < nPhases(); ip++) { ThermoPhase* p = m_phase[ip]; - nsp = p->nSpecies(); - phasemoles = 0.0; - for (ik = 0; ik < nsp; ik++) { + size_t nsp = p->nSpecies(); + double phasemoles = 0.0; + for (size_t ik = 0; ik < nsp; ik++) { phasemoles += n[k]; k++; } @@ -502,9 +472,8 @@ void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres, void MultiPhase::getElemAbundances(doublereal* elemAbundances) const { - size_t eGlobal; calcElemAbundances(); - for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) { elemAbundances[eGlobal] = m_elemAbundances[eGlobal]; } } @@ -512,20 +481,18 @@ void MultiPhase::getElemAbundances(doublereal* elemAbundances) const void MultiPhase::calcElemAbundances() const { size_t loc = 0; - size_t eGlobal; - size_t ik, kGlobal; doublereal spMoles; - for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) { m_elemAbundances[eGlobal] = 0.0; } for (size_t ip = 0; ip < nPhases(); ip++) { ThermoPhase* p = m_phase[ip]; size_t nspPhase = p->nSpecies(); doublereal phasemoles = m_moles[ip]; - for (ik = 0; ik < nspPhase; ik++) { - kGlobal = loc + ik; + for (size_t ik = 0; ik < nspPhase; ik++) { + size_t kGlobal = loc + ik; spMoles = m_moleFractions[kGlobal] * phasemoles; - for (eGlobal = 0; eGlobal < m_nel; eGlobal++) { + for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) { m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles; } } @@ -535,9 +502,8 @@ void MultiPhase::calcElemAbundances() const doublereal MultiPhase::volume() const { - int i; doublereal sum = 0; - for (i = 0; i < int(nPhases()); i++) { + for (size_t i = 0; i < nPhases(); i++) { double vol = 1.0/m_phase[i]->molarDensity(); sum += m_moles[i] * vol; } @@ -549,14 +515,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, int loglevel) { bool strt = false; - doublereal dt; - doublereal h0; - int n; - doublereal hnow, herr = 1.0; - doublereal snow, s0; - doublereal Tlow = -1.0, Thigh = -1.0; - doublereal Hlow = Undef, Hhigh = Undef, tnew; - doublereal dta=0.0, dtmax, cpb; + doublereal dta = 0.0; if (!m_init) { init(); } @@ -566,10 +525,11 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, MultiPhaseEquil e(this); return e.equilibrate(XY, err, maxsteps, loglevel); } else if (XY == HP) { - h0 = enthalpy(); - Tlow = 0.5*m_Tmin; // lower bound on T - Thigh = 2.0*m_Tmax; // upper bound on T - for (n = 0; n < maxiter; n++) { + double h0 = enthalpy(); + double Tlow = 0.5*m_Tmin; // lower bound on T + double Thigh = 2.0*m_Tmax; // upper bound on T + doublereal Hlow = Undef, Hhigh = Undef; + for (int n = 0; n < maxiter; n++) { // if 'strt' is false, the current composition will be used as // the starting estimate; otherwise it will be estimated MultiPhaseEquil e(this, strt); @@ -578,7 +538,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, try { e.equilibrate(TP, err, maxsteps, loglevel); - hnow = enthalpy(); + double hnow = enthalpy(); // the equilibrium enthalpy monotonically increases with T; // if the current value is below the target, the we know the // current temperature is too low. Set @@ -595,25 +555,26 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, Hhigh = hnow; } } + double dt; if (Hlow != Undef && Hhigh != Undef) { - cpb = (Hhigh - Hlow)/(Thigh - Tlow); + double cpb = (Hhigh - Hlow)/(Thigh - Tlow); dt = (h0 - hnow)/cpb; dta = fabs(dt); - dtmax = 0.5*fabs(Thigh - Tlow); + double dtmax = 0.5*fabs(Thigh - Tlow); if (dta > dtmax) { dt *= dtmax/dta; } } else { - tnew = sqrt(Tlow*Thigh); + double tnew = sqrt(Tlow*Thigh); dt = tnew - m_temp; } - herr = fabs((h0 - hnow)/h0); + double herr = fabs((h0 - hnow)/h0); if (herr < err) { return err; } - tnew = m_temp + dt; + double tnew = m_temp + dt; if (tnew < 0.0) { tnew = 0.5*m_temp; } @@ -629,7 +590,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, if (!strt) { strt = true; } else { - tnew = 0.5*(m_temp + Thigh); + double tnew = 0.5*(m_temp + Thigh); if (fabs(tnew - m_temp) < 1.0) { tnew = m_temp + 1.0; } @@ -640,31 +601,31 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil", "No convergence for T"); } else if (XY == SP) { - s0 = entropy(); - Tlow = 1.0; // lower bound on T - Thigh = 1.0e6; // upper bound on T - for (n = 0; n < maxiter; n++) { + double s0 = entropy(); + double Tlow = 1.0; // lower bound on T + double Thigh = 1.0e6; // upper bound on T + for (int n = 0; n < maxiter; n++) { MultiPhaseEquil e(this, strt); try { e.equilibrate(TP, err, maxsteps, loglevel); - snow = entropy(); + double snow = entropy(); if (snow < s0) { Tlow = std::max(Tlow, m_temp); } else { Thigh = std::min(Thigh, m_temp); } - dt = (s0 - snow)*m_temp/cp(); - dtmax = 0.5*fabs(Thigh - Tlow); + double dt = (s0 - snow)*m_temp/cp(); + double dtmax = 0.5*fabs(Thigh - Tlow); dtmax = (dtmax > 500.0 ? 500.0 : dtmax); dta = fabs(dt); if (dta > dtmax) { dt *= dtmax/dta; } - if (herr < err || dta < 1.0e-4) { + if (dta < 1.0e-4) { return err; } - tnew = m_temp + dt; + double tnew = m_temp + dt; setTemperature(tnew); // if the size of Delta T is not too large, use @@ -676,7 +637,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, if (!strt) { strt = true; } else { - tnew = 0.5*(m_temp + Thigh); + double tnew = 0.5*(m_temp + Thigh); setTemperature(tnew); } } @@ -685,25 +646,22 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err, "No convergence for T"); } else if (XY == TV) { doublereal v0 = volume(); - doublereal dVdP; - int n; bool start = true; - doublereal vnow, pnow, verr; - for (n = 0; n < maxiter; n++) { - pnow = pressure(); + for (int n = 0; n < maxiter; n++) { + double pnow = pressure(); MultiPhaseEquil e(this, start); start = false; e.equilibrate(TP, err, maxsteps, loglevel); - vnow = volume(); - verr = fabs((v0 - vnow)/v0); + double vnow = volume(); + double verr = fabs((v0 - vnow)/v0); if (verr < err) { return err; } // find dV/dP setPressure(pnow*1.01); - dVdP = (volume() - vnow)/(0.01*pnow); + double dVdP = (volume() - vnow)/(0.01*pnow); setPressure(pnow + 0.5*(v0 - vnow)/dVdP); } } else { @@ -852,11 +810,8 @@ std::string MultiPhase::phaseName(const size_t iph) const int MultiPhase::phaseIndex(const std::string& pName) const { - std::string tmp; for (int iph = 0; iph < (int) nPhases(); iph++) { - const ThermoPhase* tptr = m_phase[iph]; - tmp = tptr->id(); - if (tmp == pName) { + if (m_phase[iph]->id() == pName) { return iph; } } @@ -890,8 +845,8 @@ bool MultiPhase::tempOK(const size_t p) const void MultiPhase::uploadMoleFractionsFromPhases() { - size_t ip, loc = 0; - for (ip = 0; ip < nPhases(); ip++) { + size_t loc = 0; + for (size_t ip = 0; ip < nPhases(); ip++) { ThermoPhase* p = m_phase[ip]; p->getMoleFractions(&m_moleFractions[loc]); loc += p->nSpecies(); @@ -901,11 +856,10 @@ void MultiPhase::uploadMoleFractionsFromPhases() void MultiPhase::updatePhases() const { - size_t p, nsp, loc = 0; - for (p = 0; p < nPhases(); p++) { - nsp = m_phase[p]->nSpecies(); + size_t loc = 0; + for (size_t p = 0; p < nPhases(); p++) { m_phase[p]->setState_TPX(m_temp, m_press, &m_moleFractions[loc]); - loc += nsp; + loc += m_phase[p]->nSpecies(); m_temp_OK[p] = true; if (m_temp < m_phase[p]->minTemp() || m_temp > m_phase[p]->maxTemp()) { m_temp_OK[p] = false; diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index f8c52efbb..edc7319e3 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -20,14 +20,13 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ m_press = mix->pressure(); m_temp = mix->temperature(); - size_t m, k; m_force = true; m_nel = 0; m_nsp = 0; m_eloc = 1000; m_incl_species.resize(m_nsp_mix,1); m_incl_element.resize(m_nel_mix,1); - for (m = 0; m < m_nel_mix; m++) { + for (size_t m = 0; m < m_nel_mix; m++) { string enm = mix->elementName(m); // element 'E' or 'e' represents an electron; this requires special // handling, so save its index for later use @@ -40,7 +39,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ // number of 'atoms' of electrons (positive ions). if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) { m_incl_element[m] = 0; - for (k = 0; k < m_nsp_mix; k++) { + for (size_t k = 0; k < m_nsp_mix; k++) { if (m_mix->nAtoms(k,m) != 0.0) { m_incl_species[k] = 0; } @@ -55,7 +54,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ m_nel++; } // add the included elements other than electrons - for (m = 0; m < m_nel_mix; m++) { + for (size_t m = 0; m < m_nel_mix; m++) { if (m_incl_element[m] == 1 && m != m_eloc) { m_nel++; m_element.push_back(m); @@ -71,9 +70,8 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ // only extend to 273.15 K, and give unphysical results above this // temperature, leading (incorrectly) to Gibbs free energies at high // temperature lower than for liquid water. - size_t ip; - for (k = 0; k < m_nsp_mix; k++) { - ip = m_mix->speciesPhaseIndex(k); + for (size_t k = 0; k < m_nsp_mix; k++) { + size_t ip = m_mix->speciesPhaseIndex(k); if (!m_mix->solutionSpecies(k) && !m_mix->tempOK(ip)) { m_incl_species[k] = 0; @@ -88,7 +86,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ } // Now build the list of all species to be included in the calculation. - for (k = 0; k < m_nsp_mix; k++) { + for (size_t k = 0; k < m_nsp_mix; k++) { if (m_incl_species[k] ==1) { m_nsp++; m_species.push_back(k); @@ -107,8 +105,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ m_dxi.resize(nFree()); // initialize the mole numbers to the mixture composition - size_t ik; - for (ik = 0; ik < m_nsp; ik++) { + for (size_t ik = 0; ik < m_nsp; ik++) { m_moles[ik] = m_mix->speciesMoles(m_species[ik]); } @@ -121,9 +118,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ m_A.resize(m_nel, m_nsp, 0.0); m_N.resize(m_nsp, nFree()); m_order.resize(std::max(m_nsp, m_nel), 0); - for (k = 0; k < m_nsp; k++) { - m_order[k] = k; - } + iota(m_order.begin(), m_order.begin() + m_nsp, 0); // if the 'start' flag is set, estimate the initial mole numbers by doing a // linear Gibbs minimization. In this case, only the elemental composition @@ -141,7 +136,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ unsort(m_work); } - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_moles[k] += m_work[k]; m_lastmoles[k] = m_moles[k]; if (m_mix->solutionSpecies(m_species[k])) { @@ -181,8 +176,7 @@ doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err, void MultiPhaseEquil::updateMixMoles() { fill(m_work3.begin(), m_work3.end(), 0.0); - size_t k; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_work3[m_species[k]] = m_moles[k]; } m_mix->setMoles(m_work3.data()); @@ -191,8 +185,7 @@ void MultiPhaseEquil::updateMixMoles() void MultiPhaseEquil::finish() { fill(m_work3.begin(), m_work3.end(), 0.0); - size_t k; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0); } m_mix->setMoles(m_work3.data()); @@ -200,13 +193,9 @@ void MultiPhaseEquil::finish() int MultiPhaseEquil::setInitialMoles(int loglevel) { - size_t ik, j; double not_mu = 1.0e12; m_mix->getValidChemPotentials(not_mu, m_mu.data(), true); - doublereal dg_rt; - int idir; - double nu; - double delta_xi, dxi_min = 1.0e10; + double dxi_min = 1.0e10; bool redo = true; int iter = 0; @@ -220,18 +209,18 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) } // loop over all reactions - for (j = 0; j < nFree(); j++) { - dg_rt = 0.0; + for (size_t j = 0; j < nFree(); j++) { + double dg_rt = 0.0; dxi_min = 1.0e10; - for (ik = 0; ik < m_nsp; ik++) { + for (size_t ik = 0; ik < m_nsp; ik++) { dg_rt += mu(ik) * m_N(ik,j); } // fwd or rev direction - idir = (dg_rt < 0.0 ? 1 : -1); + int idir = (dg_rt < 0.0 ? 1 : -1); - for (ik = 0; ik < m_nsp; ik++) { - nu = m_N(ik, j); + for (size_t ik = 0; ik < m_nsp; ik++) { + double nu = m_N(ik, j); // set max change in progress variable by // non-negativity requirement @@ -239,7 +228,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) // isn't zero. This causes differences between // optimized and debug versions of the code if (nu*idir < 0) { - delta_xi = fabs(0.99*moles(ik)/nu); + double delta_xi = fabs(0.99*moles(ik)/nu); // if a component has nearly zero moles, redo // with a new set of components if (!redo && delta_xi < 1.0e-10 && ik < m_nel) { @@ -249,7 +238,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) } } // step the composition by dxi_min - for (ik = 0; ik < m_nsp; ik++) { + for (size_t ik = 0; ik < m_nsp; ik++) { moles(ik) += m_N(ik, j) * idir*dxi_min; } } @@ -261,36 +250,33 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) void MultiPhaseEquil::getComponents(const std::vector& order) { - size_t m, k, j; - // if the input species array has the wrong size, ignore it // and consider the species for components in declaration order. if (order.size() != m_nsp) { - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_order[k] = k; } } else { - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_order[k] = order[k]; } } size_t nRows = m_nel; size_t nColumns = m_nsp; - doublereal fctr; // set up the atomic composition matrix - for (m = 0; m < nRows; m++) { - for (k = 0; k < nColumns; k++) { + for (size_t m = 0; m < nRows; m++) { + for (size_t k = 0; k < nColumns; k++) { m_A(m, k) = m_mix->nAtoms(m_species[m_order[k]], m_element[m]); } } // Do Gaussian elimination - for (m = 0; m < nRows; m++) { + for (size_t m = 0; m < nRows; m++) { // Check for rows that are zero bool isZeroRow = true; - for (k = m; k < nColumns; k++) { + for (size_t k = m; k < nColumns; k++) { if (fabs(m_A(m,k)) > sqrt(Tiny)) { isZeroRow = false; break; @@ -301,7 +287,7 @@ void MultiPhaseEquil::getComponents(const std::vector& order) size_t n = nRows - 1; bool foundSwapCandidate = false; for (; n > m; n--) { - for (k = m; k < nColumns; k++) { + for (size_t k = m; k < nColumns; k++) { if (fabs(m_A(n,k)) > sqrt(Tiny)) { foundSwapCandidate = true; break; @@ -313,7 +299,7 @@ void MultiPhaseEquil::getComponents(const std::vector& order) } if (m != n) { // Swap this row with the last non-zero row - for (k = 0; k < nColumns; k++) { + for (size_t k = 0; k < nColumns; k++) { std::swap(m_A(n,k), m_A(m,k)); } } else { @@ -334,7 +320,7 @@ void MultiPhaseEquil::getComponents(const std::vector& order) // satisfies these criteria. doublereal maxmoles = -999.0; size_t kmax = 0; - for (k = m+1; k < nColumns; k++) { + for (size_t k = m+1; k < nColumns; k++) { if (m_A(m,k) != 0.0 && fabs(m_moles[m_order[k]]) > maxmoles) { kmax = k; maxmoles = fabs(m_moles[m_order[k]]); @@ -352,8 +338,8 @@ void MultiPhaseEquil::getComponents(const std::vector& order) } // scale row m so that the diagonal element is unity - fctr = 1.0/m_A(m,m); - for (k = 0; k < nColumns; k++) { + double fctr = 1.0/m_A(m,m); + for (size_t k = 0; k < nColumns; k++) { m_A(m,k) *= fctr; } @@ -361,7 +347,7 @@ void MultiPhaseEquil::getComponents(const std::vector& order) // * (row m) from row n, so that A(n,m) = 0. for (size_t n = m+1; n < m_nel; n++) { fctr = m_A(n,m)/m_A(m,m); - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_A(n,k) -= m_A(m,k)*fctr; } } @@ -369,11 +355,11 @@ void MultiPhaseEquil::getComponents(const std::vector& order) // The left m_nel columns of A are now upper-diagonal. Now // reduce the m_nel columns to diagonal form by back-solving - for (m = std::min(nRows,nColumns)-1; m > 0; m--) { + for (size_t m = std::min(nRows,nColumns)-1; m > 0; m--) { for (size_t n = m-1; n != npos; n--) { if (m_A(n,m) != 0.0) { - fctr = m_A(n,m); - for (k = m; k < m_nsp; k++) { + double fctr = m_A(n,m); + for (size_t k = m; k < m_nsp; k++) { m_A(n,k) -= fctr*m_A(m,k); } } @@ -383,11 +369,11 @@ void MultiPhaseEquil::getComponents(const std::vector& order) // create stoichiometric coefficient matrix. for (size_t n = 0; n < m_nsp; n++) { if (n < m_nel) { - for (k = 0; k < nFree(); k++) { + for (size_t k = 0; k < nFree(); k++) { m_N(n, k) = -m_A(n, k + m_nel); } } else { - for (k = 0; k < nFree(); k++) { + for (size_t k = 0; k < nFree(); k++) { m_N(n, k) = 0.0; } m_N(n, n - m_nel) = 1.0; @@ -395,9 +381,9 @@ void MultiPhaseEquil::getComponents(const std::vector& order) } // find reactions involving solution phase species - for (j = 0; j < nFree(); j++) { + for (size_t j = 0; j < nFree(); j++) { m_solnrxn[j] = false; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { if (m_N(k, j) != 0 && m_mix->solutionSpecies(m_species[m_order[k]])) { m_solnrxn[j] = true; } @@ -408,8 +394,7 @@ void MultiPhaseEquil::getComponents(const std::vector& order) void MultiPhaseEquil::unsort(vector_fp& x) { m_work2 = x; - size_t k; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { x[m_order[k]] = m_work2[k]; } } @@ -417,19 +402,18 @@ void MultiPhaseEquil::unsort(vector_fp& x) void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN, int loglevel) { - size_t k, ik; if (omega < 0.0) { throw CanteraError("MultiPhaseEquil::step","negative omega"); } - for (ik = 0; ik < m_nel; ik++) { - k = m_order[ik]; + for (size_t ik = 0; ik < m_nel; ik++) { + size_t k = m_order[ik]; m_lastmoles[k] = m_moles[k]; m_moles[k] += omega * deltaN[k]; } - for (ik = m_nel; ik < m_nsp; ik++) { - k = m_order[ik]; + for (size_t ik = m_nel; ik < m_nsp; ik++) { + size_t k = m_order[ik]; m_lastmoles[k] = m_moles[k]; if (m_majorsp[k]) { m_moles[k] += omega * deltaN[k]; @@ -444,7 +428,6 @@ void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN, doublereal MultiPhaseEquil::stepComposition(int loglevel) { m_iter++; - size_t ik, k = 0; doublereal grad0 = computeReactionSteps(m_dxi); // compute the mole fraction changes. @@ -458,9 +441,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) // scale omega to keep the major species non-negative doublereal FCTR = 0.99; const doublereal MAJOR_THRESHOLD = 1.0e-12; - doublereal omega = 1.0, omax, omegamax = 1.0; - for (ik = 0; ik < m_nsp; ik++) { - k = m_order[ik]; + double omegamax = 1.0; + for (size_t ik = 0; ik < m_nsp; ik++) { + size_t k = m_order[ik]; if (ik < m_nel) { FCTR = 0.99; if (m_moles[k] < MAJOR_THRESHOLD) { @@ -477,7 +460,7 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) if (m_moles[k] < MAJOR_THRESHOLD) { m_force = true; } - omax = m_moles[k]*FCTR/(fabs(m_work[k]) + Tiny); + double omax = m_moles[k]*FCTR/(fabs(m_work[k]) + Tiny); if (m_work[k] < 0.0 && omax < omegamax) { omegamax = omax; if (omegamax < 1.0e-5) { @@ -490,7 +473,7 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) } } else { if (m_work[k] < 0.0 && m_moles[k] > 0.0) { - omax = -m_moles[k]/m_work[k]; + double omax = -m_moles[k]/m_work[k]; if (omax < omegamax) { omegamax = omax; if (omegamax < 1.0e-5) { @@ -510,14 +493,14 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) doublereal not_mu = 1.0e12; m_mix->getValidChemPotentials(not_mu, m_mu.data()); doublereal grad1 = 0.0; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { grad1 += m_work[k] * m_mu[m_species[k]]; } - omega = omegamax; + double omega = omegamax; if (grad1 > 0.0) { omega *= fabs(grad0) / (grad1 + fabs(grad0)); - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_moles[k] = m_lastmoles[k]; } step(omega, m_work); @@ -527,8 +510,6 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) { - size_t j, k, ik, kc, ip; - doublereal stoich, nmoles, csum, term1, fctr, rfctr; vector_fp nu; doublereal grad = 0.0; dxi.resize(nFree()); @@ -536,24 +517,23 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) doublereal not_mu = 1.0e12; m_mix->getValidChemPotentials(not_mu, m_mu.data()); - for (j = 0; j < nFree(); j++) { + for (size_t j = 0; j < nFree(); j++) { // get stoichiometric vector getStoichVector(j, nu); // compute Delta G doublereal dg_rt = 0.0; - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { dg_rt += m_mu[m_species[k]] * nu[k]; } dg_rt /= (m_temp * GasConstant); m_deltaG_RT[j] = dg_rt; - fctr = 1.0; + double fctr = 1.0; // if this is a formation reaction for a single-component phase, // check whether reaction should be included - ik = j + m_nel; - k = m_order[ik]; + size_t k = m_order[j + m_nel]; if (!m_dsoln[k]) { if (m_moles[k] <= 0.0 && dg_rt > 0.0) { fctr = 0.0; @@ -564,36 +544,34 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) fctr = 1.0; } else { // component sum - csum = 0.0; + double csum = 0.0; for (k = 0; k < m_nel; k++) { - kc = m_order[k]; - stoich = nu[kc]; - nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny; - csum += stoich*stoich*m_dsoln[kc]/nmoles; + size_t kc = m_order[k]; + double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny; + csum += pow(nu[kc], 2)*m_dsoln[kc]/nmoles; } // noncomponent term - kc = m_order[j + m_nel]; - nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny; - term1 = m_dsoln[kc]/nmoles; + size_t kc = m_order[j + m_nel]; + double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny; + double term1 = m_dsoln[kc]/nmoles; // sum over solution phases - doublereal sum = 0.0, psum; - for (ip = 0; ip < m_mix->nPhases(); ip++) { + doublereal sum = 0.0; + for (size_t ip = 0; ip < m_mix->nPhases(); ip++) { ThermoPhase& p = m_mix->phase(ip); if (p.nSpecies() > 1) { - psum = 0.0; + double psum = 0.0; for (k = 0; k < m_nsp; k++) { kc = m_species[k]; if (m_mix->speciesPhaseIndex(kc) == ip) { - stoich = nu[k]; - psum += stoich * stoich; + psum += pow(nu[k], 2); } } sum -= psum / (fabs(m_mix->phaseMoles(ip)) + Tiny); } } - rfctr = term1 + csum + sum; + double rfctr = term1 + csum + sum; if (fabs(rfctr) < Tiny) { fctr = 1.0; } else { @@ -602,8 +580,7 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) } dxi[j] = -fctr*dg_rt; - size_t m; - for (m = 0; m < m_nel; m++) { + for (size_t m = 0; m < m_nel; m++) { if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0)) { dxi[j] = 0.0; } @@ -628,7 +605,6 @@ void MultiPhaseEquil::computeN() m_sortindex[k] = moleFractions[k].second; } - bool ok; for (size_t m = 0; m < m_nel; m++) { size_t k = 0; for (size_t ik = 0; ik < m_nsp; ik++) { @@ -637,7 +613,7 @@ void MultiPhaseEquil::computeN() break; } } - ok = false; + bool ok = false; for (size_t ij = 0; ij < m_nel; ij++) { if (k == m_order[ij]) { ok = true; @@ -685,17 +661,10 @@ double MultiPhaseEquil::phaseMoles(size_t iph) const void MultiPhaseEquil::reportCSV(const std::string& reportFile) { - size_t k; - size_t istart; - size_t nSpecies; - double vol = 0.0; - string sName; FILE* FP = fopen(reportFile.c_str(), "w"); if (!FP) { throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file"); } - double Temp = m_mix->temperature(); - double pres = m_mix->pressure(); vector_fp mf(m_nsp_mix, 1.0); vector_fp fe(m_nsp_mix, 0.0); vector_fp VolPM; @@ -705,18 +674,18 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) vector_fp mu0; vector_fp molalities; - vol = 0.0; + double vol = 0.0; for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) { - istart = m_mix->speciesIndex(0, iphase); + size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); - nSpecies = tref.nSpecies(); + size_t nSpecies = tref.nSpecies(); VolPM.resize(nSpecies, 0.0); tref.getMoleFractions(&mf[istart]); tref.getPartialMolarVolumes(VolPM.data()); double TMolesPhase = phaseMoles(iphase); double VolPhaseVolumes = 0.0; - for (k = 0; k < nSpecies; k++) { + for (size_t k = 0; k < nSpecies; k++) { VolPhaseVolumes += VolPM[k] * mf[istart + k]; } VolPhaseVolumes *= TMolesPhase; @@ -724,18 +693,18 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) } fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT" " -----------------------------\n"); - fprintf(FP,"Temperature = %11.5g kelvin\n", Temp); - fprintf(FP,"Pressure = %11.5g Pascal\n", pres); + fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature()); + fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure()); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) { - istart = m_mix->speciesIndex(0, iphase); + size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); ThermoPhase* tp = &tref; tp->getMoleFractions(&mf[istart]); string phaseName = tref.name(); double TMolesPhase = phaseMoles(iphase); - nSpecies = tref.nSpecies(); + size_t nSpecies = tref.nSpecies(); activity.resize(nSpecies, 0.0); ac.resize(nSpecies, 0.0); mu0.resize(nSpecies, 0.0); @@ -749,7 +718,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) tp->getPartialMolarVolumes(VolPM.data()); tp->getChemPotentials(mu.data()); double VolPhaseVolumes = 0.0; - for (k = 0; k < nSpecies; k++) { + for (size_t k = 0; k < nSpecies; k++) { VolPhaseVolumes += VolPM[k] * mf[istart + k]; } VolPhaseVolumes *= TMolesPhase; @@ -768,8 +737,8 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) ", , ," " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n"); } - for (k = 0; k < nSpecies; k++) { - sName = tp->speciesName(k); + for (size_t k = 0; k < nSpecies; k++) { + string sName = tp->speciesName(k); fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e," "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n", sName.c_str(), @@ -789,11 +758,11 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) ", , ," " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n"); } - for (k = 0; k < nSpecies; k++) { + for (size_t k = 0; k < nSpecies; k++) { molalities[k] = 0.0; } - for (k = 0; k < nSpecies; k++) { - sName = tp->speciesName(k); + for (size_t k = 0; k < nSpecies; k++) { + string sName = tp->speciesName(k); fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, " "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n", sName.c_str(), diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index b6b3eddaa..9b5d0da89 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -68,7 +68,6 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, double P2 = 0.0; doublereal Tlow = 0.5 * m_mix->minTemp(); doublereal Thigh = 2.0 * m_mix->maxTemp(); - doublereal Vnow, Verr; int printLvlSub = std::max(0, printLvl - 1); for (int n = 0; n < maxiter; n++) { double Pnow = m_mix->pressure(); @@ -93,7 +92,7 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, break; } strt = false; - Vnow = m_mix->volume(); + double Vnow = m_mix->volume(); if (n == 0) { V2 = Vnow; P2 = Pnow; @@ -107,7 +106,7 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, V1 = Vnow; } - Verr = fabs((Vtarget - Vnow)/Vtarget); + double Verr = fabs((Vtarget - Vnow)/Vtarget); if (Verr < err) { goto done; } @@ -171,7 +170,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, Thigh = 2.0 * m_mix->maxTemp(); } - doublereal cpb = 1.0, Tnew; + doublereal cpb = 1.0; doublereal Hlow = Undef; doublereal Hhigh = Undef; doublereal Tnow = m_mix->temperature(); @@ -211,17 +210,17 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, Hhigh = Hnow; } } - double dT, dTa, dTmax, Tnew; + double dT; if (Hlow != Undef && Hhigh != Undef) { cpb = (Hhigh - Hlow)/(Thigh - Tlow); dT = (Htarget - Hnow)/cpb; - dTa = fabs(dT); - dTmax = 0.5*fabs(Thigh - Tlow); + double dTa = fabs(dT); + double dTmax = 0.5*fabs(Thigh - Tlow); if (dTa > dTmax) { dT *= dTmax/dTa; } } else { - Tnew = sqrt(Tlow*Thigh); + double Tnew = sqrt(Tlow*Thigh); dT = clip(Tnew - Tnow, -200.0, 200.0); } double acpb = std::max(fabs(cpb), 1.0E-6); @@ -244,7 +243,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, } goto done; } - Tnew = Tnow + dT; + double Tnew = Tnow + dT; if (Tnew < 0.0) { Tnew = 0.5*Tnow; } @@ -253,7 +252,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, if (!estimateEquil) { strt = -1; } else { - Tnew = 0.5*(Tnow + Thigh); + double Tnew = 0.5*(Tnow + Thigh); if (fabs(Tnew - Tnow) < 1.0) { Tnew = Tnow + 1.0; } @@ -286,7 +285,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, Thigh = 2.0 * m_mix->maxTemp(); } - doublereal cpb = 1.0, dT, dTa, dTmax, Tnew; + doublereal cpb = 1.0, dT; doublereal Slow = Undef; doublereal Shigh = Undef; doublereal Tnow = m_mix->temperature(); @@ -338,9 +337,9 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, if (Slow != Undef && Shigh != Undef) { cpb = (Shigh - Slow)/(Thigh - Tlow); dT = (Starget - Snow)/cpb; - Tnew = Tnow + dT; - dTa = fabs(dT); - dTmax = 0.5*fabs(Thigh - Tlow); + double Tnew = Tnow + dT; + double dTa = fabs(dT); + double dTmax = 0.5*fabs(Thigh - Tlow); if (Tnew > Thigh || Tnew < Tlow) { dTmax = 1.5*fabs(Thigh - Tlow); } @@ -349,7 +348,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, dT *= dTmax/dTa; } } else { - Tnew = sqrt(Tlow*Thigh); + double Tnew = sqrt(Tlow*Thigh); dT = Tnew - Tnow; } @@ -373,7 +372,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, } return iSuccess; } - Tnew = Tnow + dT; + double Tnew = Tnow + dT; if (Tnew < 0.0) { Tnew = 0.5*Tnow; } @@ -382,7 +381,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, if (!estimateEquil) { strt = -1; } else { - Tnew = 0.5*(Tnow + Thigh); + double Tnew = 0.5*(Tnow + Thigh); if (fabs(Tnew - Tnow) < 1.0) { Tnew = Tnow + 1.0; } @@ -464,13 +463,11 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, // Check obvious bounds on the temperature and pressure NOTE, we may want to // do more here with the real bounds given by the ThermoPhase objects. - double T = m_mix->temperature(); - if (T <= 0.0) { + if (m_mix->temperature() <= 0.0) { throw CanteraError("vcs_MultiPhaseEquil::equilibrate", "Temperature less than zero on input"); } - double pres = m_mix->pressure(); - if (pres <= 0.0) { + if (m_mix->pressure() <= 0.0) { throw CanteraError("vcs_MultiPhaseEquil::equilibrate", "Pressure less than zero on input"); } @@ -568,7 +565,6 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) { - double vol = 0.0; size_t nphase = m_vprob.NPhase; FILE* FP = fopen(reportFile.c_str(), "w"); @@ -576,8 +572,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "Failure to open file"); } - double Temp = m_mix->temperature(); - double pres = m_mix->pressure(); vector_fp& mf = m_vprob.mf; double* fe = &m_vprob.m_gibbsSpecies[0]; vector_fp VolPM; @@ -587,7 +581,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) vector_fp mu0; vector_fp molalities; - vol = 0.0; + double vol = 0.0; for (size_t iphase = 0; iphase < nphase; iphase++) { size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); @@ -607,8 +601,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT" " -----------------------------\n"); - fprintf(FP,"Temperature = %11.5g kelvin\n", Temp); - fprintf(FP,"Pressure = %11.5g Pascal\n", pres); + fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature()); + fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure()); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); fprintf(FP,"Number Basis optimizations = %d\n", m_vprob.m_NumBasisOptimizations); fprintf(FP,"Number VCS iterations = %d\n", m_vprob.m_Iterations); diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 2d0e9d96f..3e74f1a85 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -57,9 +57,7 @@ vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) : vcs_VolPhase::~vcs_VolPhase() { for (size_t k = 0; k < m_numSpecies; k++) { - vcs_SpeciesProperties* sp = ListSpeciesPtr[k]; - delete sp; - sp = 0; + delete ListSpeciesPtr[k]; } } diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index aa5c33b5f..142a7e7b8 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -142,7 +142,6 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() // out components which have a pos stoichiometric value with another species // in the phase. std::vector< std::vector > zeroedPhaseLinkedZeroComponents(m_numPhases); - vector_int linkedPhases; // The logic below calculates zeroedPhaseLinkedZeroComponents for (size_t iph = 0; iph < m_numPhases; iph++) { @@ -150,7 +149,6 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() iphList.clear(); vcs_VolPhase* Vphase = m_VolPhaseList[iph]; if (Vphase->exists() < 0) { - linkedPhases.clear(); size_t nsp = Vphase->nSpecies(); for (size_t k = 0; k < nsp; k++) { size_t kspec = Vphase->spGlobalIndexVCS(k); @@ -447,26 +445,18 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) vector_fp X_est(nsp, 0.0); vector_fp delFrac(nsp, 0.0); vector_fp E_phi(nsp, 0.0); - vector_fp fracDelta_new(nsp, 0.0); vector_fp fracDelta_old(nsp, 0.0); vector_fp fracDelta_raw(nsp, 0.0); vector creationGlobalRxnNumbers(nsp, npos); m_deltaGRxn_Deficient = m_deltaGRxn_old; - vector_fp m_feSpecies_Deficient(m_numComponents, 0.0); - doublereal damp = 1.0; - doublereal dampOld = 1.0; - doublereal normUpdate = 1.0; - doublereal normUpdateOld = 1.0; - doublereal sum = 0.0; - doublereal dirProd = 0.0; - doublereal dirProdOld = 0.0; + vector_fp feSpecies_Deficient = m_feSpecies_old; // get the activity coefficients Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]); // Get the stored estimate for the composition of the phase if // it gets created - fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); + vector_fp fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); std::vector componentList; for (size_t k = 0; k < nsp; k++) { @@ -476,11 +466,8 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } } - for (size_t k = 0; k < m_numComponents; k++) { - m_feSpecies_Deficient[k] = m_feSpecies_old[k]; - } - normUpdate = 0.1 * vcs_l2norm(fracDelta_new); - damp = 1.0E-2; + double normUpdate = 0.1 * vcs_l2norm(fracDelta_new); + double damp = 1.0E-2; if (doSuccessiveSubstitution) { int KP = 0; @@ -500,11 +487,12 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } } bool converged = false; + double dirProd = 0.0; for (int its = 0; its < 200 && (!converged); its++) { - dampOld = damp; - normUpdateOld = normUpdate; + double dampOld = damp; + double normUpdateOld = normUpdate; fracDelta_old = fracDelta_new; - dirProdOld = dirProd; + double dirProdOld = dirProd; // Given a set of fracDelta's, we calculate the fracDelta's // for the component species, if any @@ -552,10 +540,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) size_t kc = componentList[i]; size_t kc_spec = Vphase->spGlobalIndexVCS(kc); if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) { - m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]); } else { - m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF); } } @@ -571,14 +559,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) { m_deltaGRxn_Deficient[irxn] += - m_stoichCoeffRxnMatrix(kc_spec,irxn) * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); + m_stoichCoeffRxnMatrix(kc_spec,irxn) * (feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); } } } } // Calculate the E_phi's - sum = 0.0; + double sum = 0.0; funcPhaseStability = sum_Xcomp - 1.0; for (size_t k = 0; k < nsp; k++) { size_t kspec = Vphase->spGlobalIndexVCS(k); diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 20ca98750..815b38e66 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -99,11 +99,9 @@ VCS_PROB::~VCS_PROB() { for (size_t i = 0; i < nspecies; i++) { delete SpeciesThermo[i]; - SpeciesThermo[i] = 0; } for (size_t iph = 0; iph < NPhase; iph++) { delete VPhaseList[iph]; - VPhaseList[iph] = 0; } } diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 5a987aec6..8ad95d73d 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -358,12 +358,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } // Start of the regular processing - double s; - if (m_SSPhase[kspec]) { - s = 0.0; - } else { - s = 1.0 / m_molNumSpecies_old[kspec]; - } + double s = (m_SSPhase[kspec]) ? 0.0 : 1.0 / m_molNumSpecies_old[kspec]; for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 721faefc0..b6d8278b0 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -31,18 +31,13 @@ static void printProgress(const vector &spName, int VCS_SOLVE::vcs_setMolesLinProg() { - size_t ik, irxn; double test = -1.0E-10; if (m_debug_print_lvl >= 2) { plogf(" --- call setInitialMoles\n"); } - double dg_rt; - int idir; - double nu; - double delta_xi, dxi_min = 1.0e10; - bool redo = true; + double dxi_min = 1.0e10; int retn; int iter = 0; bool abundancesOK = true; @@ -53,7 +48,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() vector_fp wx(m_numElemConstraints, 0.0); vector_fp aw(m_numSpeciesTot, 0.0); - for (ik = 0; ik < m_numSpeciesTot; ik++) { + for (size_t ik = 0; ik < m_numSpeciesTot; ik++) { if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) { m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]); } @@ -63,6 +58,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); } + bool redo = true; while (redo) { if (!vcs_elabcheck(0)) { if (m_debug_print_lvl >= 2) { @@ -98,10 +94,10 @@ int VCS_SOLVE::vcs_setMolesLinProg() } // loop over all reactions - for (irxn = 0; irxn < m_numRxnTot; irxn++) { + for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) { // dg_rt is the Delta_G / RT value for the reaction - ik = m_numComponents + irxn; - dg_rt = m_SSfeSpecies[ik]; + size_t ik = m_numComponents + irxn; + double dg_rt = m_SSfeSpecies[ik]; dxi_min = 1.0e10; const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t jcomp = 0; jcomp < m_numElemConstraints; jcomp++) { @@ -110,17 +106,17 @@ int VCS_SOLVE::vcs_setMolesLinProg() // fwd or rev direction. // idir > 0 implies increasing the current species // idir < 0 implies decreasing the current species - idir = (dg_rt < 0.0 ? 1 : -1); + int idir = (dg_rt < 0.0 ? 1 : -1); if (idir < 0) { dxi_min = m_molNumSpecies_old[ik]; } for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) { - nu = sc_irxn[jcomp]; + double nu = sc_irxn[jcomp]; // set max change in progress variable by // non-negativity requirement if (nu*idir < 0) { - delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu); + double delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu); // 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)) { diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 6b7e3ca8f..503f4f3cf 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -48,10 +48,9 @@ void VCS_SOLVE::checkDelta1(double* const dsLocal, int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) { int stage = MAIN; - int solveFail; bool allMinorZeroedSpecies = false; size_t it1 = 0; - size_t npb, iti; + size_t iti; int rangeErrorFound = 0; bool giveUpOnElemAbund = false; int finalElemAbundAttempts = 0; @@ -75,7 +74,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) m_aw.assign(m_numSpeciesTot, 0.0); m_wx.assign(m_numElemConstraints, 0.0); - solveFail = false; + int solveFail = false; // Evaluate the elemental composition vcs_elab(); @@ -216,7 +215,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) // reason why we are here is because all of the non-component // species in the problem have been eliminated for one reason or // another. - npb = vcs_recheck_deleted(); + size_t npb = vcs_recheck_deleted(); // If we haven't found any species that needed adding we are done. if (npb <= 0) { @@ -232,7 +231,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) } } else if (stage == RETURN_A) { // CLEANUP AND RETURN BLOCK - npb = vcs_recheck_deleted(); + size_t npb = vcs_recheck_deleted(); // If we haven't found any species that needed adding we are done. if (npb > 0) { @@ -249,7 +248,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) } else if (stage == RETURN_B) { // Add back deleted species in non-zeroed phases. Estimate their // mole numbers. - npb = vcs_add_all_deleted(); + size_t npb = vcs_add_all_deleted(); if (npb > 0) { iti = 0; if (m_debug_print_lvl >= 1) { @@ -1427,12 +1426,8 @@ void VCS_SOLVE::solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec, double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, char* ANOTE) const { - double dx = 0.0, a; double w_kspec = m_molNumSpecies_old[kspec]; - double molNum_kspec_new; - double wTrial, tmp; double dg_irxn = m_deltaGRxn_old[irxn]; - doublereal s; size_t iph = m_phaseID[kspec]; *do_delete = false; @@ -1445,20 +1440,22 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, sprintf(ANOTE,"minor species alternative calc"); } if (dg_irxn >= 23.0) { - molNum_kspec_new = w_kspec * 1.0e-10; + double molNum_kspec_new = w_kspec * 1.0e-10; if (w_kspec < VCS_DELETE_MINORSPECIES_CUTOFF) { - goto L_ZERO_SPECIES; + // delete the species from the current list of active species, + // because its concentration has gotten too small. + *do_delete = true; + return - w_kspec; } return molNum_kspec_new - w_kspec; } else { if (fabs(dg_irxn) <= m_tolmin2) { - molNum_kspec_new = w_kspec; return 0.0; } } // get the diagonal of the activity coefficient Jacobian - s = m_np_dLnActCoeffdMolNum(kspec,kspec) / m_tPhaseMoles_old[iph]; + double s = m_np_dLnActCoeffdMolNum(kspec,kspec) / m_tPhaseMoles_old[iph]; // We fit it to a power law approximation of the activity coefficient // @@ -1469,10 +1466,10 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, // We then solve the resulting calculation: // // gamma * x = gamma_0 * x0 exp (-deltaG/RT); - a = clip(w_kspec * s, -1.0+1e-8, 100.0); - tmp = clip(-dg_irxn / (1.0 + a), -200.0, 200.0); - wTrial = w_kspec * exp(tmp); - molNum_kspec_new = wTrial; + double a = clip(w_kspec * s, -1.0+1e-8, 100.0); + double tmp = clip(-dg_irxn / (1.0 + a), -200.0, 200.0); + double wTrial = w_kspec * exp(tmp); + double molNum_kspec_new = wTrial; if (wTrial > 100. * w_kspec) { double molNumMax = 0.0001 * m_tPhaseMoles_old[iph]; @@ -1491,39 +1488,34 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, } if ((molNum_kspec_new) < VCS_DELETE_MINORSPECIES_CUTOFF) { - goto L_ZERO_SPECIES; + // delete the species from the current list of active species, + // because its concentration has gotten too small. + *do_delete = true; + return - w_kspec; } return molNum_kspec_new - w_kspec; - // Alternate return based for cases where we need to delete the species - // from the current list of active species, because its concentration - // has gotten too small. -L_ZERO_SPECIES: - ; - *do_delete = true; - return - w_kspec; } else { // Voltage calculation // Need to check the sign -> This is good for electrons - dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim; + double dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim; if (ANOTE) { sprintf(ANOTE,"voltage species alternative calc"); } + return dx; } - return dx; } int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) { size_t irxn = kspec - m_numComponents; int retn = 1; - double delta = *delta_ptr; AssertThrowMsg(kspec >= m_numComponents, "VCS_SOLVE::delta_species", "delete_species() ERROR: called for a component {}", kspec); if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { // Attempt the given dx. If it doesn't work, try to see if a smaller one // would work, - double dx = delta; + double dx = *delta_ptr; double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t j = 0; j < m_numComponents; ++j) { if (m_molNumSpecies_old[j] > 0.0) { @@ -1742,7 +1734,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) } } - double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0; + double dxPerm = 0.0, dxPerm2 = 0.0; for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { if (m_phaseID[kcomp] == iph) { if (m_debug_print_lvl >= 2) { @@ -1754,7 +1746,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) size_t irxn = kspec - m_numComponents; if (m_phaseID[kspec] != iph) { if (m_stoichCoeffRxnMatrix(kcomp,irxn) != 0.0) { - dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn); + double dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn); if (dxWant + m_molNumSpecies_old[kspec] < 0.0) { dxPerm = -m_molNumSpecies_old[kspec]; } @@ -1763,7 +1755,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) if (m_phaseID[jcomp] == iph) { dxPerm = 0.0; } else { - dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn); + double dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn); if (dj + m_molNumSpecies_old[kcomp] < 0.0) { dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(jcomp,irxn); } @@ -1938,7 +1930,6 @@ bool VCS_SOLVE::recheck_deleted_phase(const int iphase) size_t VCS_SOLVE::vcs_add_all_deleted() { - size_t retn; if (m_numSpeciesRdc == m_numSpeciesTot) { return 0; } @@ -1981,7 +1972,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() size_t iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { double dx = m_molNumSpecies_new[kspec]; - retn = delta_species(kspec, &dx); + size_t retn = delta_species(kspec, &dx); if (retn == 0) { if (m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", @@ -2012,7 +2003,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); vcs_deltag(0, true, VCS_STATECALC_OLD); - retn = 0; + size_t retn = 0; for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { size_t kspec = m_indexRxnToSpecies[irxn]; size_t iph = m_phaseID[kspec]; @@ -2149,7 +2140,6 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], size_t k; size_t juse = npos; size_t jlose = npos; - double* scrxn_ptr; clockWC tickTock; if (m_debug_print_lvl >= 2) { plogf(" "); @@ -2197,7 +2187,6 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], m_numComponents = ncTrial; *usedZeroedSpecies = false; vector_int ipiv(ncTrial); - int info; // Use a temporary work array for the mole numbers, aw[] std::copy(m_molNumSpecies_old.begin(), @@ -2457,6 +2446,7 @@ L_END_LOOP: } // Solve the linear system to calculate the reaction matrix, // m_stoichCoeffRxnMatrix. + int info; ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info); if (info) { plogf("vcs_solve_TP ERROR: Error factorizing stoichiometric coefficient matrix\n"); @@ -2618,7 +2608,7 @@ L_END_LOOP: // m_deltaMolNumPhase(iphase,irxn), and the phase participation array, // PhaseParticipation[irxn][iphase] for (size_t irxn = 0; irxn < m_numRxnTot; ++irxn) { - scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); + double* scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); size_t kspec = m_indexRxnToSpecies[irxn]; size_t iph = m_phaseID[kspec]; m_deltaMolNumPhase(iph,irxn) = 1.0; @@ -2702,8 +2692,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const size_t iph = m_phaseID[kspec]; int irxn = int(kspec) - int(m_numComponents); - vcs_VolPhase* VPhase = m_VolPhaseList[iph]; - int phaseExist = VPhase->exists(); + int phaseExist = m_VolPhaseList[iph]->exists(); // Treat zeroed out species first if (m_molNumSpecies_old[kspec] <= 0.0) { @@ -2714,8 +2703,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const // see if the species has an element which is so low that species will // always be zero for (size_t j = 0; j < m_numElemConstraints; ++j) { - int elType = m_elType[j]; - if (elType == VCS_ELEM_TYPE_ABSPOS) { + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { double atomComp = m_formulaMatrix(kspec,j); if (atomComp > 0.0) { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; @@ -2758,9 +2746,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const } } } else if (negChangeComp < 0.0) { - size_t jph = m_phaseID[j]; - vcs_VolPhase* jVPhase = m_VolPhaseList[jph]; - if (jVPhase->exists() <= 0) { + if (m_VolPhaseList[m_phaseID[j]]->exists() <= 0) { 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 " @@ -3395,7 +3381,6 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, const int vcsState, const bool alterZeroedPhases) { - int icase = 0; size_t irxnl = m_numRxnRdc; if (doDeleted) { irxnl = m_numRxnTot; @@ -3435,7 +3420,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { size_t kspec = irxn + m_numComponents; if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { - icase = 0; + int icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (kspec = 0; kspec < m_numComponents; ++kspec) { @@ -3452,7 +3437,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, } else if (L == 0) { // ALL REACTIONS for (size_t irxn = 0; irxn < irxnl; ++irxn) { - icase = 0; + int icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t kspec = 0; kspec < m_numComponents; ++kspec) { @@ -3471,7 +3456,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) { size_t kspec = irxn + m_numComponents; if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) { - icase = 0; + int icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (kspec = 0; kspec < m_numComponents; ++kspec) { @@ -3719,10 +3704,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc"); } - size_t irxnl = m_numRxnRdc; - if (doDeleted) { - irxnl = m_numRxnTot; - } + size_t irxnl = (doDeleted) ? m_numRxnTot : m_numRxnRdc; vcs_VolPhase* vPhase = m_VolPhaseList[iphase]; if (m_debug_print_lvl >= 2) { @@ -3900,19 +3882,17 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { return dx; } - double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF; // Check to make sure that species is zero in the solution vector // If it isn't, we don't know what's happening AssertThrowMsg(m_molNumSpecies_old[kspec] == 0.0, "VCS_SOLVE::vcs_birthGuess", "we shouldn't be here"); - int ss = m_SSPhase[kspec]; - if (!ss) { + if (!m_SSPhase[kspec]) { // Logic to handle species in multiple species phases. We cap the moles // here at 1.0E-15 kmol. bool soldel_ret; double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret); - dx = std::min(w_kspec + dxm, 1e-15); + dx = std::min(VCS_DELETE_MINORSPECIES_CUTOFF + dxm, 1e-15); } else { // Logic to handle single species phases. There is no real way to // estimate the moles. So we set it to a small number. diff --git a/src/equil/vcs_solve_phaseStability.cpp b/src/equil/vcs_solve_phaseStability.cpp index 4f47db49a..e1565e4f7 100644 --- a/src/equil/vcs_solve_phaseStability.cpp +++ b/src/equil/vcs_solve_phaseStability.cpp @@ -18,7 +18,6 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl { // ifunc determines the problem type int ifunc = 0; - int iStab = 0; // This function is called to create the private data using the public data. size_t nspecies0 = vprob->nspecies + 10; @@ -91,7 +90,7 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl // Solve the problem at a fixed Temperature and Pressure (all information // concerning Temperature and Pressure has already been derived. The free // energies are now in dimensionless form.) - iStab = vcs_solve_phaseStability(iphase, ifunc, feStable, printLvl); + int iStab = vcs_solve_phaseStability(iphase, ifunc, feStable, printLvl); // Redimensionalize the free energies using the reverse of vcs_nondim to add // back units. @@ -109,7 +108,6 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, { double test = -1.0E-10; bool usedZeroedSpecies; - int iStab = 0; vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0); vector_fp ss(m_numElemConstraints, 0.0); @@ -133,12 +131,10 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, m_deltaGRxn_Deficient = m_deltaGRxn_old; funcVal = vcs_phaseStabilityTest(iph); if (funcVal > 0.0) { - iStab = 1; + return 1; } else { - iStab = 0; + return 0; } - - return iStab; } } diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp index a0f82fe62..5e8c12ab6 100644 --- a/src/equil/vcs_species_thermo.cpp +++ b/src/equil/vcs_species_thermo.cpp @@ -94,26 +94,21 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin, throw CanteraError("VCS_SPECIES_THERMO::GStar_R_calc", "Possible inconsistency"); } - size_t kspec = IndexSpeciesPhase; OwningPhase->setState_TP(TKelvin, pres); - fe = OwningPhase->GStar_calc_one(kspec); + fe = OwningPhase->GStar_calc_one(IndexSpeciesPhase); double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - fe /= R; - return fe; + return fe / R; } double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin, double presPA) { - double vol; if (m_VCS_UnitsFormat != VCS_UNITS_MKS) { throw CanteraError("VCS_SPECIES_THERMO::VolStar_calc", "Possible inconsistency"); } - size_t kspec = IndexSpeciesPhase; OwningPhase->setState_TP(TKelvin, presPA); - vol = OwningPhase->VolStar_calc_one(kspec); - return vol; + return OwningPhase->VolStar_calc_one(IndexSpeciesPhase); } double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) @@ -128,9 +123,8 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin) throw CanteraError("VCS_SPECIES_THERMO::G0_R_calc", "Possible inconsistency"); } - size_t kspec = IndexSpeciesPhase; OwningPhase->setState_T(TKelvin); - double fe = OwningPhase->G0_calc_one(kspec); + double fe = OwningPhase->G0_calc_one(IndexSpeciesPhase); double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); fe /= R; SS0_feSave = fe; @@ -144,9 +138,7 @@ double VCS_SPECIES_THERMO::eval_ac(size_t kglob) // they are, then the currPhAC[] boolean may be used to reduce repeated // work. Just set currPhAC[iph], when the activity coefficients for all // species in the phase are reevaluated. - size_t kspec = IndexSpeciesPhase; - double ac = OwningPhase->AC_calc_one(kspec); - return ac; + return OwningPhase->AC_calc_one(IndexSpeciesPhase); } } diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index b88a579c9..719f94f01 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -20,17 +20,16 @@ using namespace std; namespace Cantera { -double vcs_l2norm(const vector_fp vec) +double vcs_l2norm(const vector_fp& vec) { - size_t len = vec.size(); - if (len == 0) { + if (vec.empty()) { return 0.0; } double sum = 0.0; for (const auto& val : vec) { sum += val * val; } - return std::sqrt(sum / len); + return std::sqrt(sum / vec.size()); } size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)