diff --git a/include/cantera/base/global.h b/include/cantera/base/global.h index 584e72380..beb2901e3 100644 --- a/include/cantera/base/global.h +++ b/include/cantera/base/global.h @@ -271,6 +271,11 @@ inline T clip(const T& value, const T& lower, const T& upper) return std::max(lower, std::min(upper, value)); } +//! Sign of a number. Returns -1 if x < 0, 1 if x > 0 and 0 if x == 0. +template int sign(T x) { + return (T(0) < x) - (x < T(0)); +} + } #endif diff --git a/include/cantera/equil/vcs_defs.h b/include/cantera/equil/vcs_defs.h index 433b30867..b256ef4c0 100644 --- a/include/cantera/equil/vcs_defs.h +++ b/include/cantera/equil/vcs_defs.h @@ -14,21 +14,6 @@ namespace VCSnonideal { -/* - * COMMON DEFINITIONS -> Protect them against redefinitions - */ -//@{ - -#ifndef SQUARE -# define SQUARE(x) ((x) * (x)) -#endif - -#ifndef DSIGN -# define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 )) -#endif - -//@} - /*! * ERROR CODES * diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 022f0fef5..73e660f8b 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -473,7 +473,7 @@ L_CLEANUP: #ifdef DEBUG_MODE double l2after = 0.0; for (size_t i = 0; i < m_numElemConstraints; ++i) { - l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]); + l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2); } l2after = sqrt(l2after/m_numElemConstraints); if (m_debug_print_lvl >= 2) { diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 38747aaab..f78716d81 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -137,7 +137,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, */ sa[jr] = 0.0; for (size_t ml = 0; ml < ncomponents; ++ml) { - sa[jr] += SQUARE(sm[ml + jr*ncomponents]); + sa[jr] += pow(sm[ml + jr*ncomponents], 2); } /* **************************************************** */ /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */ diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 96291d982..38e79fc88 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -384,7 +384,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { if (m_molNumSpecies_old[j] > 0.0) { - s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; + s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; } } } @@ -392,7 +392,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) Vphase = m_VolPhaseList[j]; if (! Vphase->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } } diff --git a/src/equil/vcs_rank.cpp b/src/equil/vcs_rank.cpp index f7693d46a..cfc2cc88f 100644 --- a/src/equil/vcs_rank.cpp +++ b/src/equil/vcs_rank.cpp @@ -178,7 +178,7 @@ namespace VCSnonideal { */ sa[jr] = 0.0; for (ml = 0; ml < numElemConstraints; ++ml) { - sa[jr] += SQUARE(sm[ml + jr * numElemConstraints]); + sa[jr] += pow(sm[ml + jr * numElemConstraints], 2); } /* **************************************************** */ /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */ @@ -253,7 +253,7 @@ namespace VCSnonideal { sa[jr] = 0.0; for (ml = 0; ml < numSpecies; ++ml) { - sa[jr] += SQUARE(sm[ml + jr * numSpecies]); + sa[jr] += pow(sm[ml + jr * numSpecies], 2); } if (sa[jr] < 1.0e-6) lindep = true; diff --git a/src/equil/vcs_root1d.cpp b/src/equil/vcs_root1d.cpp index ceed270b8..4caeb098a 100644 --- a/src/equil/vcs_root1d.cpp +++ b/src/equil/vcs_root1d.cpp @@ -165,9 +165,9 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, c[3] = x0; c[4] = x1; c[5] = x2; - c[6] = SQUARE(x0); - c[7] = SQUARE(x1); - c[8] = SQUARE(x2); + c[6] = pow(x0, 2); + c[7] = pow(x1, 2); + c[8] = pow(x2, 2); f[0] = f0; f[1] = f1; f[2] = f2; @@ -199,8 +199,8 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax, * Pick out situations where the convergence may be * accelerated. */ - if ((DSIGN(xnew - x2) == DSIGN(x2 - x1)) && - (DSIGN(x2 - x1) == DSIGN(x1 - x0))) { + if ((sign(xnew - x2) == sign(x2 - x1)) && + (sign(x2 - x1) == sign(x1 - x0))) { xnew += xnew - x2; if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xquada = %-9.4g", xnew); @@ -225,13 +225,13 @@ QUAD_BAIL: */ slope = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < slope) { - xnew = x1 + DSIGN(xnew-x1) * slope; + xnew = x1 + sign(xnew-x1) * slope; if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | x10%% = %-9.4g", xnew); } } if (fabs(xnew - x2) < slope) { - xnew = x2 + DSIGN(xnew-x2) * slope; + xnew = x2 + sign(xnew-x2) * slope; if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | x10%% = %-9.4g", xnew); } @@ -243,7 +243,7 @@ QUAD_BAIL: */ slope = 2.0 * fabs(x2 - x1); if (fabs(slope) < fabs(xnew - x2)) { - xnew = x2 + DSIGN(xnew-x2) * slope; + xnew = x2 + sign(xnew-x2) * slope; if (DEBUG_MODE_ENABLED && printLvl >= 3) { fprintf(fp, " | xlimitsize = %-9.4g", xnew); } diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 01fec624b..f91c65931 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -176,7 +176,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { if (m_molNumSpecies_old[j] > 0.0) { - s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; + s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; } } } @@ -184,7 +184,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) vcs_VolPhase* Vphase = m_VolPhaseList[j]; if (!Vphase->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } } @@ -462,13 +462,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } for (size_t j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { - s += SQUARE(m_stoichCoeffRxnMatrix(j,irxn)) / m_molNumSpecies_old[j]; + s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; } } for (size_t j = 0; j < m_numPhases; j++) { if (!(m_VolPhaseList[j])->m_singleSpecies) { if (m_tPhaseMoles_old[j] > 0.0) { - s -= SQUARE(m_deltaMolNumPhase(j,irxn)) / m_tPhaseMoles_old[j]; + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } } diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 36f647286..2b70a379d 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -2767,7 +2767,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], */ sa[jr] = 0.0; for (size_t ml = 0; ml < m_numElemConstraints; ++ml) { - sa[jr] += SQUARE(sm[ml + jr*m_numElemConstraints]); + sa[jr] += pow(sm[ml + jr*m_numElemConstraints], 2); } /* **************************************************** */ /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */ diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 3713035a1..f7dc32d64 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -26,14 +26,6 @@ using namespace std; namespace Cantera { -#ifndef SQUARE -# define SQUARE(x) ( (x) * (x) ) -#endif - -#ifndef DSIGN -#define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 )) -#endif - //! Print out a form for the current function evaluation /*! * @param fp Pointer to the FILE object @@ -524,13 +516,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun */ xDelMin = fabs(x2 - x1) / 10.; if (fabs(xnew - x1) < xDelMin) { - xnew = x1 + DSIGN(xnew-x1) * xDelMin; + xnew = x1 + sign(xnew-x1) * xDelMin; if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } if (fabs(xnew - x2) < 0.1 * xDelMin) { - xnew = x2 + DSIGN(xnew-x2) * 0.1 * xDelMin; + xnew = x2 + sign(xnew-x2) * 0.1 * xDelMin; if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } @@ -548,7 +540,7 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } if (fabs(xDelMax) < fabs(xnew - x2)) { - xnew = x2 + DSIGN(xnew-x2) * xDelMax; + xnew = x2 + sign(xnew-x2) * xDelMax; if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | xlimitsize = %-11.5E", xnew); } @@ -560,13 +552,13 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun */ xDelMin = 0.1 * fabs(x2 - x1); if (fabs(xnew - x2) < xDelMin) { - xnew = x2 + DSIGN(xnew - x2) * xDelMin; + xnew = x2 + sign(xnew - x2) * xDelMin; if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); } } if (fabs(xnew - x1) < xDelMin) { - xnew = x1 + DSIGN(xnew - x1) * xDelMin; + xnew = x1 + sign(xnew - x1) * xDelMin; if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { fprintf(fp, " | x10%% = %-11.5E", xnew); }