From b5396837b208253af884a1ddf1a060f79a4512c8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 3 Jun 2014 16:52:43 +0000 Subject: [PATCH] Make use of std::min, std::max, and clip --- include/cantera/base/utilities.h | 7 +- include/cantera/base/vec_functions.h | 9 +-- include/cantera/thermo/SimpleThermo.h | 9 +-- include/cantera/zeroD/flowControllers.h | 12 +-- src/equil/ChemEquil.cpp | 73 ++++--------------- src/equil/MultiPhase.cpp | 22 ++---- src/equil/MultiPhaseEquil.cpp | 8 +- src/equil/vcs_MultiPhaseEquil.cpp | 16 +--- src/equil/vcs_VolPhase.cpp | 4 +- src/equil/vcs_elem.cpp | 25 ++----- src/equil/vcs_phaseStability.cpp | 27 ++----- src/equil/vcs_prob.cpp | 4 +- src/equil/vcs_rxnadj.cpp | 5 +- src/equil/vcs_setMolesLinProg.cpp | 4 +- src/equil/vcs_solve.cpp | 4 +- src/equil/vcs_solve_TP.cpp | 60 +++------------ src/kinetics/AqueousKinetics.cpp | 8 +- src/kinetics/GasKinetics.cpp | 8 +- src/kinetics/ImplicitSurfChem.cpp | 4 +- src/kinetics/InterfaceKinetics.cpp | 12 +-- src/kinetics/ReactionPath.cpp | 25 ++----- src/kinetics/solveSP.cpp | 10 +-- src/numerics/BandMatrix.cpp | 16 +--- src/numerics/NonlinearSolver.cpp | 4 +- src/numerics/SquareMatrix.cpp | 8 +- src/oneD/MultiNewton.cpp | 4 +- src/oneD/OneDim.cpp | 4 +- src/oneD/Sim1D.cpp | 4 +- src/oneD/boundaries1D.cpp | 4 +- src/thermo/DebyeHuckel.cpp | 10 +-- src/thermo/HMWSoln.cpp | 5 +- src/thermo/MolalityVPSSTP.cpp | 10 +-- src/thermo/NasaThermo.cpp | 8 +- src/thermo/PDSS_HKFT.cpp | 4 +- src/thermo/RedlichKisterVPSSTP.cpp | 11 +-- src/thermo/RedlichKwongMFTP.cpp | 6 +- src/thermo/ShomateThermo.h | 9 +-- src/thermo/SingleSpeciesTP.cpp | 28 +------ src/thermo/VPSSMgr.cpp | 10 +-- src/thermo/VPSSMgr_General.cpp | 16 +--- src/thermo/WaterPropsIAPWSphi.cpp | 5 +- src/transport/MMCollisionInt.cpp | 36 ++------- src/transport/TransportFactory.cpp | 36 +++------ src/zeroD/ConstPressureReactor.cpp | 4 +- src/zeroD/IdealGasConstPressureReactor.cpp | 4 +- src/zeroD/IdealGasReactor.cpp | 4 +- src/zeroD/Reactor.cpp | 4 +- .../cathermo/HMW_dupl_test/HMW_dupl_test.cpp | 9 +-- .../HMW_graph_CpvT/HMW_graph_CpvT.cpp | 9 +-- .../cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp | 6 +- .../cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp | 9 +-- .../cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp | 9 +-- .../surfSolverTest/surfaceSolver.cpp | 6 +- .../surfSolverTest/surfaceSolver2.cpp | 6 +- 54 files changed, 140 insertions(+), 524 deletions(-) diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h index 8506e1c27..ed039d11f 100644 --- a/include/cantera/base/utilities.h +++ b/include/cantera/base/utilities.h @@ -295,10 +295,9 @@ template inline doublereal absmax(InputIter begin, InputIter end) { doublereal amax = 0.0; - for (; begin != end; ++begin) - if (fabs(*begin) > amax) { - amax = fabs(*begin); - } + for (; begin != end; ++begin) { + amax = std::max(fabs(*begin), amax); + } return amax; } diff --git a/include/cantera/base/vec_functions.h b/include/cantera/base/vec_functions.h index 9f3f7ccd5..e544d0578 100644 --- a/include/cantera/base/vec_functions.h +++ b/include/cantera/base/vec_functions.h @@ -163,16 +163,9 @@ template inline T absmax(const std::vector& v) { int n = v.size(); - T val; T maxval = 0.0; for (int i = 0; i < n; i++) { - val = v[i]; - if (val < 0) { - val = -val; - } - if (val > maxval) { - maxval = val; - } + maxval = std::max(std::abs(v[i]), maxval); } return maxval; } diff --git a/include/cantera/thermo/SimpleThermo.h b/include/cantera/thermo/SimpleThermo.h index d490f4756..6f6419cee 100644 --- a/include/cantera/thermo/SimpleThermo.h +++ b/include/cantera/thermo/SimpleThermo.h @@ -144,13 +144,8 @@ public: m_nspData++; doublereal tlow = minTemp_; doublereal thigh = maxTemp_; - - if (tlow > m_tlow_max) { - m_tlow_max = tlow; - } - if (thigh < m_thigh_min) { - m_thigh_min = thigh; - } + m_tlow_max = std::max(tlow, m_tlow_max); + m_thigh_min = std::min(thigh, m_thigh_min); if (m_tlow.size() < index + 1) { m_tlow.resize(index + 1, tlow); diff --git a/include/cantera/zeroD/flowControllers.h b/include/cantera/zeroD/flowControllers.h index 8d005975f..73b8807cd 100644 --- a/include/cantera/zeroD/flowControllers.h +++ b/include/cantera/zeroD/flowControllers.h @@ -37,9 +37,7 @@ public: if (m_func) { m_mdot = m_func->eval(time); } - if (m_mdot < 0.0) { - m_mdot = 0.0; - } + m_mdot = std::max(m_mdot, 0.0); } }; @@ -67,9 +65,7 @@ public: doublereal master_mdot = m_master->massFlowRate(time); m_mdot = master_mdot + m_coeffs[0]*(in().pressure() - out().pressure()); - if (m_mdot < 0.0) { - m_mdot = 0.0; - } + m_mdot = std::max(m_mdot, 0.0); } protected: @@ -101,9 +97,7 @@ public: } else { m_mdot = m_coeffs[0]*delta_P; } - if (m_mdot < 0.0) { - m_mdot = 0.0; - } + m_mdot = std::max(m_mdot, 0.0); } }; diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index b917d0555..0a645a829 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -206,9 +206,7 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, e.setInitialMixMoles(loglevel-1); // store component indices - if (m_nComponents > m_kk) { - m_nComponents = m_kk; - } + m_nComponents = std::min(m_nComponents, m_kk); for (size_t m = 0; m < m_nComponents; m++) { m_component[m] = e.componentIndex(m); } @@ -253,9 +251,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, s.getMoleFractions(DATA_PTR(xMF_est)); for (size_t n = 0; n < s.nSpecies(); n++) { - if (xMF_est[n] < 1.0E-20) { - xMF_est[n] = 1.0E-20; - } + xMF_est[n] = std::max(xMF_est[n], 1e-20); } s.setMoleFractions(DATA_PTR(xMF_est)); s.getMoleFractions(DATA_PTR(xMF_est)); @@ -272,9 +268,7 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, for (size_t m = 0; m < m_nComponents; m++) { size_t k = m_orderVectorSpecies[m]; m_component[m] = k; - if (xMF_est[k] < 1.0E-8) { - xMF_est[k] = 1.0E-8; - } + xMF_est[k] = std::max(xMF_est[k], 1e-8); } s.setMoleFractions(DATA_PTR(xMF_est)); s.getMoleFractions(DATA_PTR(xMF_est)); @@ -489,20 +483,11 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, doublereal tminPhase = s.minTemp(); // loop to estimate T if (!tempFixed) { - doublereal tmin; - doublereal tmax; - - tmin = s.temperature(); - if (tmin < tminPhase) { - tmin = tminPhase; - } + doublereal tmin = std::max(s.temperature(), tminPhase); if (tmin > tmaxPhase) { tmin = tmaxPhase - 20; } - tmax = tmin + 10.; - if (tmax > tmaxPhase) { - tmax = tmaxPhase; - } + doublereal tmax = std::min(tmin + 10., tmaxPhase); if (tmax < tminPhase) { tmax = tminPhase + 20; } @@ -554,12 +539,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, if (fabs(dt) < 50.0) { break; } - if (dt > 200.) { - dt = 200.; - } - if (dt < -200.) { - dt = -200.; - } + dt = clip(dt, -200.0, 200.0); if ((t0 + dt) < tminPhase) { dt = 0.5*((t0) + tminPhase) - t0; } @@ -938,7 +918,7 @@ void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, r0.resize(len); r1.resize(len); size_t n, m; - doublereal rdx, dx, xsave, dx2; + doublereal rdx, dx, xsave; doublereal atol = 1.e-10; equilResidual(s, x, elmols, r0, xval, yval, loglevel-1); @@ -946,11 +926,7 @@ void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, m_doResPerturb = false; for (n = 0; n < len; n++) { xsave = x[n]; - dx = atol; - dx2 = fabs(xsave) * 1.0E-7; - if (dx2 > dx) { - dx = dx2; - } + dx = std::max(atol, fabs(xsave) * 1.0E-7); x[n] = xsave + dx; dx = x[n] - xsave; rdx = 1.0/dx; @@ -990,9 +966,7 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, for (size_t m = 0; m < m_mm; m++) { tmp += nAtoms(k,m) * x[m]; } - if (tmp > 100.) { - tmp = 100.; - } + tmp = std::min(tmp, 100.0); if (tmp < -300.) { n_i_calc[k] = 0.0; } else { @@ -1060,17 +1034,10 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } for (m = 0; m < m_mm; m++) { - if (x[m] > 50.0) { - x[m] = 50.; - } if (elMoles[m] > 1.0E-70) { - if (x[m] < -100) { - x[m] = -100.; - } + x[m] = clip(x[m], -100.0, 50.0); } else { - if (x[m] < -1000.) { - x[m] = -1000.; - } + x[m] = clip(x[m], -1000.0, 50.0); } } @@ -1088,9 +1055,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, sum = nAtoms(k,m); tmp += sum * x[m]; sum2 += sum; - if (sum2 > nAtomsMax) { - nAtomsMax = sum2; - } + nAtomsMax = std::max(nAtomsMax, sum2); } if (tmp > 100.) { n_t += 2.8E43; @@ -1218,9 +1183,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, resid[m_mm] = 0.1; } } else if (n_t > elMolesTotal) { - if (resid[m_mm] > 0.0) { - resid[m_mm] = 0.0; - } + resid[m_mm] = std::min(resid[m_mm], 0.0); } } else { /* @@ -1505,16 +1468,10 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, beta = 1.0; for (m = 0; m < m_mm; m++) { if (resid[m] > 1.0) { - double betat = 1.0 / resid[m]; - if (betat < beta) { - beta = betat; - } + beta = std::min(beta, 1.0 / resid[m]); } if (resid[m] < -1.0) { - double betat = -1.0 / resid[m]; - if (betat < beta) { - beta = betat; - } + beta = std::min(beta, -1.0 / resid[m]); } } if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index 56daeede1..7c3a9cb41 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -157,14 +157,8 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles) // water, only one of which should be present if the mixture // represents an equilibrium state. if (p->nSpecies() > 1) { - double t = p->minTemp(); - if (t > m_Tmin) { - m_Tmin = t; - } - t = p->maxTemp(); - if (t < m_Tmax) { - m_Tmax = t; - } + m_Tmin = std::max(p->minTemp(), m_Tmin); + m_Tmax = std::min(p->maxTemp(), m_Tmax); } } @@ -506,9 +500,7 @@ void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) vector_fp tmpMoles(m_nsp, 0.0); getMoles(DATA_PTR(tmpMoles)); tmpMoles[indexS] += addedMoles; - if (tmpMoles[indexS] < 0.0) { - tmpMoles[indexS] = 0.0; - } + tmpMoles[indexS] = std::max(tmpMoles[indexS], 0.0); setMoles(DATA_PTR(tmpMoles)); } @@ -692,13 +684,9 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, e.equilibrate(TP, err, maxsteps, loglevel); snow = entropy(); if (snow < s0) { - if (m_temp > Tlow) { - Tlow = m_temp; - } + Tlow = std::max(Tlow, m_temp); } else { - if (m_temp < Thigh) { - Thigh = m_temp; - } + Thigh = std::min(Thigh, m_temp); } dt = (s0 - snow)*m_temp/cp(); dtmax = 0.5*fabs(Thigh - Tlow); diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index 277d7e2fc..e6145c515 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -264,9 +264,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) if (!redo && delta_xi < 1.0e-10 && ik < m_nel) { redo = true; } - if (delta_xi < dxi_min) { - dxi_min = delta_xi; - } + dxi_min = std::min(dxi_min, delta_xi); } } // step the composition by dxi_min @@ -706,9 +704,7 @@ doublereal MultiPhaseEquil::error() } else { err = fabs(m_deltaG_RT[j]); } - if (err > maxerr) { - maxerr = err; - } + maxerr = std::max(maxerr, err); } return maxerr; } diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 152dec054..139753042 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -235,13 +235,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, } } else { Tnew = sqrt(Tlow*Thigh); - dT = Tnew - Tnow; - if (dT < -200.) { - dT = -200; - } - if (dT > 200.) { - dT = 200.; - } + dT = clip(Tnew - Tnow, -200.0, 200.0); } double acpb = std::max(fabs(cpb), 1.0E-6); double denom = std::max(fabs(Htarget), acpb); @@ -311,12 +305,8 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, doublereal Slow = Undef; doublereal Shigh = Undef; doublereal Tnow = m_mix->temperature(); - if (Tnow < Tlow) { - Tlow = Tnow; - } - if (Tnow > Thigh) { - Thigh = Tnow; - } + Tlow = std::min(Tnow, Tlow); + Thigh = std::max(Tnow, Thigh); int printLvlSub = std::max(printLvl - 1, 0); for (int n = 0; n < maxiter; n++) { diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index d68a2949f..692d716a5 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -434,9 +434,7 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles, } else { m_UpToDate = true; m_vcsStateStatus = vcsStateStatus; - if (m_existence > VCS_PHASE_EXIST_NO) { - m_existence = VCS_PHASE_EXIST_NO; - } + m_existence = std::min(m_existence, VCS_PHASE_EXIST_NO); } double fractotal = 1.0; v_totalMoles = totalMoles; diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 911895be0..022f0fef5 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -273,15 +273,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) double par = 0.5; for (size_t i = 0; i < m_numComponents; ++i) { if (m_molNumSpecies_old[i] > 0.0) { - double xx = -x[i] / m_molNumSpecies_old[i]; - if (par < xx) { - par = xx; - } + par = std::max(par, -x[i] / m_molNumSpecies_old[i]); } } - if (par > 100.0) { - par = 100.0; - } + par = std::min(par, 100.0); par = 1.0 / par; if (par < 1.0 && par > 0.0) { retn = 2; @@ -398,9 +393,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (m_elemAbundances[i] > 0.0) { if (m_formulaMatrix(kspec,i) < 0.0) { m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i) ; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); vcs_elab(); break; } @@ -408,9 +401,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (m_elemAbundances[i] < 0.0) { if (m_formulaMatrix(kspec,i) > 0.0) { m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i); - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); vcs_elab(); break; } @@ -453,9 +444,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (m_formulaMatrix(kspec,i) < 0.0) { double delta = dev / m_formulaMatrix(kspec,i) ; m_molNumSpecies_old[kspec] += delta; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); vcs_elab(); break; } @@ -464,9 +453,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (m_formulaMatrix(kspec,i) > 0.0) { double delta = dev / m_formulaMatrix(kspec,i) ; m_molNumSpecies_old[kspec] += delta; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); vcs_elab(); break; } diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 045f470e9..98a29285b 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -12,6 +12,7 @@ #include using namespace std; +using namespace Cantera; namespace VCSnonideal { @@ -330,9 +331,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) FephaseMax = Fephase; } } else { - if (Fephase > FephaseMax) { - FephaseMax = Fephase; - } + FephaseMax = std::max(FephaseMax, Fephase); } if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { plogf(" --- %18s %5d %11.3g %11.3g\n", @@ -484,10 +483,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) ratioComp = 1.0; if (deltaJ > 0.0) { double delta0 = m_molNumSpecies_old[j]; - double dampj = delta0 / deltaJ * 0.9; - if (dampj < damp) { - damp = dampj; - } + damp = std::min(damp, delta0 / deltaJ * 0.9); } } else { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { @@ -605,9 +601,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } for (size_t k = 0; k < nsp; k++) { - if (fracDelta_new[k] < 1.0E-13) { - fracDelta_new[k] = 1.0E-13; - } + fracDelta_new[k] = std::min(fracDelta_new[k], 1e-13); } bool converged = false; for (int its = 0; its < 200 && (!converged); its++) { @@ -704,13 +698,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) size_t kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; - double deltaGRxn = m_deltaGRxn_Deficient[irxn]; - if (deltaGRxn > 50.0) { - deltaGRxn = 50.0; - } - if (deltaGRxn < -50.0) { - deltaGRxn = -50.0; - } + double deltaGRxn = clip(m_deltaGRxn_Deficient[irxn], -50.0, 50.0); E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec]; sum += E_phi[k]; funcPhaseStability += E_phi[k]; @@ -802,10 +790,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } } } - if (damp < 0.000001) { - damp = 0.000001; - } - + damp = std::max(damp, 0.000001); for (size_t k = 0; k < nsp; k++) { fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]); } diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 027650679..93520fb42 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -148,9 +148,7 @@ void VCS_PROB::resizeElements(size_t nel, int force) m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS); ElActive.resize(nel, 1); NE0 = nel; - if (ne > NE0) { - ne = NE0; - } + ne = std::min(ne, NE0); } } diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 038505e1f..505b84115 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -598,10 +598,7 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) { size_t kspec = m_indexRxnToSpecies[irxn]; size_t kph = m_phaseID[kspec]; - double np_kspec = m_tPhaseMoles_old[kph]; - if (np_kspec < 1.0E-13) { - np_kspec = 1.0E-13; - } + double np_kspec = std::max(m_tPhaseMoles_old[kph], 1e-13); double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); /* * First the diagonal term of the Jacobian diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 161a0bfd7..79050e348 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -145,9 +145,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() redo = true; } } - if (delta_xi < dxi_min) { - dxi_min = delta_xi; - } + dxi_min = std::min(dxi_min, delta_xi); } } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 3dc1a6243..2184929bc 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -248,9 +248,7 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) Cantera::clockWC tickTock; int iprintTime = std::max(ipr, ip1); - if (m_timing_print_lvl < iprintTime) { - iprintTime = m_timing_print_lvl ; - } + iprintTime = std::min(iprintTime, m_timing_print_lvl); if (ifunc > 2) { plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n", diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index e596cf3f7..e99f5f706 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -192,9 +192,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) plogendl(); size_t iph = m_phaseID[i]; double tmp = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10; - if (VCS_DELETE_MINORSPECIES_CUTOFF*10. > tmp) { - tmp = VCS_DELETE_MINORSPECIES_CUTOFF*10.; - } + tmp = std::max(tmp, VCS_DELETE_MINORSPECIES_CUTOFF*10.); m_molNumSpecies_old[i] = tmp; } } @@ -1729,9 +1727,7 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, if (w_kspec <= 0.0) { w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF; } - if (dg_irxn < -200.) { - dg_irxn = -200.; - } + dg_irxn = std::max(dg_irxn, -200.0); if (DEBUG_MODE_ENABLED && ANOTE) { sprintf(ANOTE,"minor species alternative calc"); } @@ -1766,18 +1762,8 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete, * * */ - a = w_kspec * s; - if (a < (-1.0 + 1.0E-8)) { - a = -1.0 + 1.0E-8; - } else if (a > 100.0) { - a = 100.0; - } - tmp = -dg_irxn / (1.0 + a); - if (tmp < -200.) { - tmp = -200.; - } else if (tmp > 200.) { - tmp = 200.; - } + 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; @@ -1879,9 +1865,7 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) m_molNumSpecies_old[j] += tmp; m_tPhaseMoles_old[iph] += tmp; vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD); - if (m_molNumSpecies_old[j] < 0.0) { - m_molNumSpecies_old[j] = 0.0; - } + m_molNumSpecies_old[j] = std::max(m_molNumSpecies_old[j], 0.0); } } } @@ -2307,12 +2291,7 @@ bool VCS_SOLVE::recheck_deleted_phase(const int iphase) for (size_t kk = 0; kk < Vphase->nSpecies(); kk++) { size_t kspec = Vphase->spGlobalIndexVCS(kk); size_t irxn = kspec + m_numComponents; - if (m_deltaGRxn_old[irxn] > 50.0) { - m_deltaGRxn_old[irxn] = 50.0; - } - if (m_deltaGRxn_old[irxn] < -50.0) { - m_deltaGRxn_old[irxn] = -50.0; - } + m_deltaGRxn_old[irxn] = clip(m_deltaGRxn_old[irxn], -50.0, 50.0); phaseDG -= exp(-m_deltaGRxn_old[irxn]); } @@ -2710,12 +2689,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], } } } - if (nonZeroesKspec < minNonZeroes) { - minNonZeroes = nonZeroesKspec; - } - if (maxConcPossKspec > maxConcPoss) { - maxConcPoss = maxConcPossKspec; - } + minNonZeroes = std::min(minNonZeroes, nonZeroesKspec); + maxConcPoss = std::max(maxConcPoss, maxConcPossKspec); } } } @@ -4115,12 +4090,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, // We may need to look at deltaGRxn for components! if (kspec >= m_numComponents) { size_t irxn = kspec - m_numComponents; - if (deltaGRxn[irxn] > 50.0) { - deltaGRxn[irxn] = 50.0; - } - if (deltaGRxn[irxn] < -50.0) { - deltaGRxn[irxn] = -50.0; - } + deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0); poly += exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec]; } } @@ -4395,12 +4365,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, for (size_t irxn = 0; irxn < irxnl; ++irxn) { size_t kspec = m_indexRxnToSpecies[irxn]; if (m_phaseID[kspec] == iphase) { - if (deltaGRxn[irxn] > 50.0) { - deltaGRxn[irxn] = 50.0; - } - if (deltaGRxn[irxn] < -50.0) { - deltaGRxn[irxn] = -50.0; - } + deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0); phaseDG -= exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec]; } } @@ -4536,10 +4501,7 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) */ bool soldel_ret; double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret); - dx = w_kspec + dxm; - if (dx > 1.0E-15) { - dx = 1.0E-15; - } + dx = std::min(w_kspec + dxm, 1e-15); } else { /* * Logic to handle single species phases diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index a870c6e35..144f34ac7 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -400,9 +400,7 @@ void AqueousKinetics::installReagents(const ReactionData& r) reactantGlobalOrder += nsFlt; ns = (size_t) nsFlt; if ((doublereal) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } if (r.rstoich[n] != 0.0) { m_rrxn[r.reactants[n]][rnum] += r.rstoich[n]; @@ -420,9 +418,7 @@ void AqueousKinetics::installReagents(const ReactionData& r) productGlobalOrder += nsFlt; ns = (size_t) nsFlt; if ((double) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } if (r.pstoich[n] != 0.0) { m_prxn[r.products[n]][rnum] += r.pstoich[n]; diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 94bab6322..b41088045 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -616,9 +616,7 @@ void GasKinetics::installReagents(const ReactionData& r) reactantGlobalOrder += nsFlt; ns = (size_t) nsFlt; if ((doublereal) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } if (r.rstoich[n] != 0.0) { m_rrxn[r.reactants[n]][rnum] += r.rstoich[n]; @@ -636,9 +634,7 @@ void GasKinetics::installReagents(const ReactionData& r) productGlobalOrder += nsFlt; ns = (size_t) nsFlt; if ((double) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } if (r.pstoich[n] != 0.0) { m_prxn[r.products[n]][rnum] += r.pstoich[n]; diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp index 52e58a2f7..97b332518 100644 --- a/src/kinetics/ImplicitSurfChem.cpp +++ b/src/kinetics/ImplicitSurfChem.cpp @@ -51,9 +51,7 @@ ImplicitSurfChem::ImplicitSurfChem(vector k) : m_nsp.push_back(nsp); m_nv += m_nsp.back(); nt = k[n]->nTotalSpecies(); - if (nt > ntmax) { - ntmax = nt; - } + ntmax = std::max(nt, ntmax); m_specStartIndex.push_back(kinSpIndex); kinSpIndex += nsp; diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 6b9b908f7..9f44d6bc1 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -862,9 +862,7 @@ void InterfaceKinetics::installReagents(const ReactionData& r) nsFlt = r.rstoich[n]; ns = (size_t) nsFlt; if ((doublereal) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } /* * Add to m_rrxn. m_rrxn is a vector of maps. m_rrxn has a length @@ -889,9 +887,7 @@ void InterfaceKinetics::installReagents(const ReactionData& r) nsFlt = r.pstoich[n]; ns = (size_t) nsFlt; if ((doublereal) ns != nsFlt) { - if (ns < 1) { - ns = 1; - } + ns = std::max(ns, 1); } /* * Add to m_prxn. m_prxn is a vector of maps. m_prxn has a length @@ -1038,9 +1034,7 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) if (exists) { if (!m_phaseExists[iphase]) { m_phaseExistsCheck--; - if (m_phaseExistsCheck < 0) { - m_phaseExistsCheck = 0; - } + m_phaseExistsCheck = std::max(m_phaseExistsCheck, 0); m_phaseExists[iphase] = true; } m_phaseIsStable[iphase] = true; diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index 6af6a0c38..98461a2db 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -122,9 +122,7 @@ vector_int ReactionPathDiagram::reactions() Path* p; for (i = 0; i < npaths; i++) { p = path(i); - if (p->flow() > flmax) { - flmax = p->flow(); - } + flmax = std::max(p->flow(), flmax); } m_rxns.clear(); for (i = 0; i < npaths; i++) { @@ -170,9 +168,7 @@ void ReactionPathDiagram::findMajorPaths(doublereal athreshold, size_t lda, k1 = m_speciesNumber[n]; k2 = m_speciesNumber[m]; fl = fabs(netFlow(k1,k2)); - if (fl > netmax) { - netmax = fl; - } + netmax = std::max(fl, netmax); } } for (n = 0; n < nn; n++) { @@ -253,18 +249,13 @@ void ReactionPathDiagram::exportToDot(ostream& s) if (flx < 0.0) { flx = -flx; } - if (flx > flmax) { - flmax = flx; - } + flmax = std::max(flx, flmax); } } } else { flmax = scale; } - - if (flmax < 1.e-10) { - flmax = 1.e-10; - } + flmax = std::max(flmax, 1e-10); // loop over all unique pairs of nodes @@ -346,9 +337,7 @@ void ReactionPathDiagram::exportToDot(ostream& s) else { for (size_t i = 0; i < nPaths(); i++) { p = path(i); - if (p->flow() > flmax) { - flmax = p->flow(); - } + flmax = std::max(p->flow(), flmax); } for (size_t i = 0; i < nPaths(); i++) { @@ -434,9 +423,7 @@ void ReactionPathDiagram::linkNodes(size_t k1, size_t k2, size_t rxn, } ff->addReaction(rxn, value, legend); m_rxns[rxn] = 1; - if (ff->flow() > m_flxmax) { - m_flxmax = ff->flow(); - } + m_flxmax = std::max(ff->flow(), m_flxmax); } std::vector ReactionPathDiagram::species() diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index 712e1972e..23233b0ee 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -154,10 +154,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, doublereal t_real = 0.0, update_norm = 1.0E6; bool do_time = false, not_converged = true; - - if (m_ioflag > 1) { - m_ioflag = 1; - } + m_ioflag = std::min(m_ioflag, 1); /* * Set the initial value of the do_time parameter @@ -665,10 +662,7 @@ static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, i *label = int(i); } } - - if (damp < 1.0e-2) { - damp = 1.0e-2; - } + damp = std::max(damp, 1e-2); /* * Only allow the damping parameter to increase by a factor of three each * iteration. Heuristic to avoid oscillations in the value of damp diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index 599acce64..3601d5b6e 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -369,9 +369,7 @@ doublereal BandMatrix::oneNorm() const for (int i = j - ku; i <= j + kl; i++) { sum += fabs(colP[kl + ku + i - j]); } - if (sum > value) { - value = sum; - } + value = std::max(sum, value); } return value; } @@ -380,15 +378,11 @@ size_t BandMatrix::checkRows(doublereal& valueSmall) const { valueSmall = 1.0E300; size_t iSmall = npos; - double vv; for (int i = 0; i < (int) m_n; i++) { double valueS = 0.0; for (int j = i - (int) m_kl; j <= i + (int) m_ku; j++) { if (j >= 0 && j < (int) m_n) { - vv = fabs(value(i,j)); - if (vv > valueS) { - valueS = vv; - } + valueS = std::max(fabs(value(i,j)), valueS); } } if (valueS < valueSmall) { @@ -406,15 +400,11 @@ size_t BandMatrix::checkColumns(doublereal& valueSmall) const { valueSmall = 1.0E300; size_t jSmall = npos; - double vv; for (int j = 0; j < (int) m_n; j++) { double valueS = 0.0; for (int i = j - (int) m_ku; i <= j + (int) m_kl; i++) { if (i >= 0 && i < (int) m_n) { - vv = fabs(value(i,j)); - if (vv > valueS) { - valueS = vv; - } + valueS = std::max(fabs(value(i,j)), valueS); } } if (valueS < valueSmall) { diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index 136cf4d28..a250ea4ab 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -813,9 +813,7 @@ void NonlinearSolver::calcSolnToResNormVector() if (resNormOld > 0.0) { m_ScaleSolnNormToResNorm = resNormOld; } - if (m_ScaleSolnNormToResNorm < 1.0E-8) { - m_ScaleSolnNormToResNorm = 1.0E-8; - } + m_ScaleSolnNormToResNorm = std::max(m_ScaleSolnNormToResNorm, 1e-8); // Recalculate the residual weights now that we know the value of m_ScaleSolnNormToResNorm computeResidWts(); diff --git a/src/numerics/SquareMatrix.cpp b/src/numerics/SquareMatrix.cpp index 4074b160d..d030e4164 100644 --- a/src/numerics/SquareMatrix.cpp +++ b/src/numerics/SquareMatrix.cpp @@ -351,9 +351,7 @@ size_t SquareMatrix::checkRows(doublereal& valueSmall) const for (size_t i = 0; i < m_nrows; i++) { double valueS = 0.0; for (size_t j = 0; j < m_nrows; j++) { - if (fabs(value(i,j)) > valueS) { - valueS = fabs(value(i,j)); - } + valueS = std::max(fabs(value(i,j)), valueS); } if (valueS < valueSmall) { iSmall = i; @@ -370,9 +368,7 @@ size_t SquareMatrix::checkColumns(doublereal& valueSmall) const for (size_t j = 0; j < m_nrows; j++) { double valueS = 0.0; for (size_t i = 0; i < m_nrows; i++) { - if (fabs(value(i,j)) > valueS) { - valueS = fabs(value(i,j)); - } + valueS = std::max(fabs(value(i,j)), valueS); } if (valueS < valueSmall) { jSmall = j; diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 3cc7743cb..0cd943a09 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -134,9 +134,7 @@ doublereal norm_square(const doublereal* x, for (j = 0; j < np; j++) { f = step[nv*j + n]/ewt; sum += f*f; - if (f*f > f2max) { - f2max = f*f; - } + f2max = std::max(f*f, f2max); } } return sum; diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 1e7daa61f..a95f7df34 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -338,9 +338,7 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x, if (m == 100) { dt *= 1.5; } - if (dt > m_tmax) { - dt = m_tmax; - } + dt = std::min(dt, m_tmax); } // No solution could be found with this time step. diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index c8fd24af0..16ee7f1a4 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -290,9 +290,7 @@ void Sim1D::solve(int loglevel, bool refine_grid) } else { nsteps = m_steps[istep]; } - if (dt > m_tmax) { - dt = m_tmax; - } + dt = std::min(dt, m_tmax); } } if (loglevel > 0) { diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index 7bf78a993..7a9a8e4de 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -789,9 +789,7 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg, r[k+1] = m_work[k + ioffset] * m_sphase->size(k) * rs0; r[k+1] -= rdt*(x[k+1] - prevSoln(k+1,0)); diag[k+1] = 1; - if (x[k+1] > maxx) { - maxx = x[k+1]; - } + maxx = std::max(x[k+1], maxx); } r[1] = 1.0 - sum; diag[1] = 0; diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index c16add358..b4d92be29 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -1292,10 +1292,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const m_IionicMolality += m_molalities[k] * z_k * z_k; } m_IionicMolality /= 2.0; - - if (m_IionicMolality > m_maxIionicStrength) { - m_IionicMolality = m_maxIionicStrength; - } + m_IionicMolality = std::min(m_IionicMolality, m_maxIionicStrength); /* * Calculate the stoichiometric ionic charge @@ -1313,10 +1310,7 @@ void DebyeHuckel::s_update_lnMolalityActCoeff() const } } m_IionicMolalityStoich /= 2.0; - - if (m_IionicMolalityStoich > m_maxIionicStrength) { - m_IionicMolalityStoich = m_maxIionicStrength; - } + m_IionicMolalityStoich = std::min(m_IionicMolalityStoich, m_maxIionicStrength); /* * Possibly update the stored value of the diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 4482c58d1..7883f1e2c 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -1383,10 +1383,7 @@ void HMWSoln::calcMolalitiesCropped() const for (size_t k = 0; k < m_kk; k++) { m_molalitiesCropped[k] = m_molalities[k]; - Itmp = m_molalities[k] * charge(k) * charge(k); - if (Itmp > Imax) { - Imax = Itmp; - } + Imax = std::max(m_molalities[k] * charge(k) * charge(k), Imax); } int cropMethod = 1; diff --git a/src/thermo/MolalityVPSSTP.cpp b/src/thermo/MolalityVPSSTP.cpp index 8c43a0b6b..581a72185 100644 --- a/src/thermo/MolalityVPSSTP.cpp +++ b/src/thermo/MolalityVPSSTP.cpp @@ -131,10 +131,7 @@ doublereal MolalityVPSSTP::moleFSolventMin() const void MolalityVPSSTP::calcMolalities() const { getMoleFractions(DATA_PTR(m_molalities)); - double xmolSolvent = m_molalities[m_indexSolvent]; - if (xmolSolvent < m_xmolSolventMIN) { - xmolSolvent = m_xmolSolventMIN; - } + double xmolSolvent = std::max(m_molalities[m_indexSolvent], m_xmolSolventMIN); double denomInv = 1.0/ (m_Mnaught * xmolSolvent); for (size_t k = 0; k < m_kk; k++) { m_molalities[k] *= denomInv; @@ -294,10 +291,7 @@ void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const { getMolalityActivityCoefficients(ac); AssertThrow(m_indexSolvent==0, "MolalityVPSSTP::getActivityCoefficients"); - double xmolSolvent = moleFraction(m_indexSolvent); - if (xmolSolvent < m_xmolSolventMIN) { - xmolSolvent = m_xmolSolventMIN; - } + double xmolSolvent = std::max(moleFraction(m_indexSolvent), m_xmolSolventMIN); for (size_t k = 1; k < m_kk; k++) { ac[k] /= xmolSolvent; } diff --git a/src/thermo/NasaThermo.cpp b/src/thermo/NasaThermo.cpp index fffdacd2a..a90e25d22 100644 --- a/src/thermo/NasaThermo.cpp +++ b/src/thermo/NasaThermo.cpp @@ -90,12 +90,8 @@ void NasaThermo::install(const std::string& name, size_t index, int type, m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid, ref_pressure, &clow[0])); - if (tlow > m_tlow_max) { - m_tlow_max = tlow; - } - if (thigh < m_thigh_min) { - m_thigh_min = thigh; - } + m_tlow_max = std::max(tlow, m_tlow_max); + m_thigh_min = std::min(thigh, m_thigh_min); if (m_tlow.size() < index + 1) { m_tlow.resize(index + 1, tlow); m_thigh.resize(index + 1, thigh); diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp index c3c1018b3..e652fbd54 100644 --- a/src/thermo/PDSS_HKFT.cpp +++ b/src/thermo/PDSS_HKFT.cpp @@ -953,9 +953,7 @@ doublereal PDSS_HKFT::f(const doublereal temp, const doublereal pres, const int if (TC < 155.0) { return 0.0; } - if (TC > 355.0) { - TC = 355.0; - } + TC = std::min(TC, 355.0); if (presBar > 1000.) { return 0.0; } diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp index b571d9bbc..02c800cbf 100644 --- a/src/thermo/RedlichKisterVPSSTP.cpp +++ b/src/thermo/RedlichKisterVPSSTP.cpp @@ -731,11 +731,7 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) * Get the string containing all of the values */ ctml::getFloatArray(xmlChild, hParams, true, "toSI", "excessEnthalpy"); - size_t nParamsFound = hParams.size(); - if (nParamsFound > Npoly) { - Npoly = nParamsFound; - } - + Npoly = std::max(hParams.size(), Npoly); } if (nodeName == "excessentropy") { @@ -743,10 +739,7 @@ void RedlichKisterVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies) * Get the string containing all of the values */ ctml::getFloatArray(xmlChild, sParams, true, "toSI", "excessEntropy"); - size_t nParamsFound = sParams.size(); - if (nParamsFound > Npoly) { - Npoly = nParamsFound; - } + Npoly = std::max(sParams.size(), Npoly); } } hParams.resize(Npoly, 0.0); diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index abfb3206f..7f175bce3 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -1025,11 +1025,7 @@ doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGu double btmp; calculateAB(TKelvin, atmp, btmp); - doublereal pres = presGuess; - double pp = psatEst(TKelvin); - if (pres < pp) { - pres = pp; - } + doublereal pres = std::max(psatEst(TKelvin), presGuess); double Vroot[3]; bool foundLiq = false; diff --git a/src/thermo/ShomateThermo.h b/src/thermo/ShomateThermo.h index 693c9c2ee..814694b3c 100644 --- a/src/thermo/ShomateThermo.h +++ b/src/thermo/ShomateThermo.h @@ -170,13 +170,8 @@ public: refPressure, chigh)); m_low[igrp-1].push_back(ShomatePoly(index, tlow, tmid, refPressure, clow)); - if (tlow > m_tlow_max) { - m_tlow_max = tlow; - } - if (thigh < m_thigh_min) { - m_thigh_min = thigh; - } - + m_tlow_max = std::max(m_tlow_max, tlow); + m_thigh_min = std::min(m_thigh_min, thigh); if (m_tlow.size() < index + 1) { m_tlow.resize(index + 1, tlow); m_thigh.resize(index + 1, thigh); diff --git a/src/thermo/SingleSpeciesTP.cpp b/src/thermo/SingleSpeciesTP.cpp index 8da09153f..cbdcef659 100644 --- a/src/thermo/SingleSpeciesTP.cpp +++ b/src/thermo/SingleSpeciesTP.cpp @@ -243,12 +243,7 @@ void SingleSpeciesTP::setState_HP(doublereal h, doublereal p, doublereal dt; setPressure(p); for (int n = 0; n < 50; n++) { - dt = (h - enthalpy_mass())/cp_mass(); - if (dt > 100.0) { - dt = 100.0; - } else if (dt < -100.0) { - dt = -100.0; - } + dt = clip((h - enthalpy_mass())/cp_mass(), -100.0, 100.0); setState_TP(temperature() + dt, p); if (fabs(dt) < tol) { return; @@ -267,12 +262,7 @@ void SingleSpeciesTP::setState_UV(doublereal u, doublereal v, setDensity(1.0/v); } for (int n = 0; n < 50; n++) { - dt = (u - intEnergy_mass())/cv_mass(); - if (dt > 100.0) { - dt = 100.0; - } else if (dt < -100.0) { - dt = -100.0; - } + dt = clip((u - intEnergy_mass())/cv_mass(), -100.0, 100.0); setTemperature(temperature() + dt); if (fabs(dt) < tol) { return; @@ -289,12 +279,7 @@ void SingleSpeciesTP::setState_SP(doublereal s, doublereal p, doublereal dt; setPressure(p); for (int n = 0; n < 50; n++) { - dt = (s - entropy_mass())*temperature()/cp_mass(); - if (dt > 100.0) { - dt = 100.0; - } else if (dt < -100.0) { - dt = -100.0; - } + dt = clip((s - entropy_mass())*temperature()/cp_mass(), -100.0, 100.0); setState_TP(temperature() + dt, p); if (fabs(dt) < tol) { return; @@ -313,12 +298,7 @@ void SingleSpeciesTP::setState_SV(doublereal s, doublereal v, setDensity(1.0/v); } for (int n = 0; n < 50; n++) { - dt = (s - entropy_mass())*temperature()/cv_mass(); - if (dt > 100.0) { - dt = 100.0; - } else if (dt < -100.0) { - dt = -100.0; - } + dt = clip((s - entropy_mass())*temperature()/cv_mass(), -100.0, 100.0); setTemperature(temperature() + dt); if (fabs(dt) < tol) { return; diff --git a/src/thermo/VPSSMgr.cpp b/src/thermo/VPSSMgr.cpp index 36dd03a18..334c8a51f 100644 --- a/src/thermo/VPSSMgr.cpp +++ b/src/thermo/VPSSMgr.cpp @@ -393,14 +393,8 @@ void VPSSMgr::initThermoXML(XML_Node& phaseNode, const std::string& id) m_p0 = kPDSS->refPressure(); for (size_t i = 0; i < m_kk; i++) { const PDSS* kPDSS = m_vptp_ptr->providePDSS(i); - doublereal mint = kPDSS->minTemp(); - if (mint > m_minTemp) { - m_minTemp = mint; - } - mint = kPDSS->maxTemp(); - if (mint < m_maxTemp) { - m_maxTemp = mint; - } + m_minTemp = std::max(m_minTemp, kPDSS->minTemp()); + m_maxTemp = std::min(m_maxTemp, kPDSS->maxTemp()); } } diff --git a/src/thermo/VPSSMgr_General.cpp b/src/thermo/VPSSMgr_General.cpp index 85a1eef63..5ca631e34 100644 --- a/src/thermo/VPSSMgr_General.cpp +++ b/src/thermo/VPSSMgr_General.cpp @@ -216,19 +216,9 @@ VPSSMgr_General::createInstallPDSS(size_t k, const XML_Node& speciesNode, m_PDSS_ptrs.resize(k+1, 0); } m_PDSS_ptrs[k] = kPDSS; - if ((k+1) >= m_kk) { - m_kk = k+1; - } - - doublereal minTemp = kPDSS->minTemp(); - if (minTemp > m_minTemp) { - m_minTemp = minTemp; - } - - doublereal maxTemp = kPDSS->maxTemp(); - if (maxTemp < m_maxTemp) { - m_maxTemp = maxTemp; - } + m_kk = std::max(m_kk, k+1); + m_minTemp = std::max(m_minTemp, kPDSS->minTemp()); + m_maxTemp = std::min(m_maxTemp, kPDSS->maxTemp()); doublereal p0 = kPDSS->refPressure(); if (k == 0) { diff --git a/src/thermo/WaterPropsIAPWSphi.cpp b/src/thermo/WaterPropsIAPWSphi.cpp index 5e996d99f..6333697b9 100644 --- a/src/thermo/WaterPropsIAPWSphi.cpp +++ b/src/thermo/WaterPropsIAPWSphi.cpp @@ -13,6 +13,7 @@ #include #include +#include namespace Cantera { @@ -1066,9 +1067,7 @@ doublereal WaterPropsIAPWSphi::dfind(doublereal p_red, doublereal tau, double if (n < 10) { dpdx = dpddelta * 1.1; } - if (dpdx < 0.001) { - dpdx = 0.001; - } + dpdx = std::max(dpdx, 0.001); /* * Formulate the update to reduced density using diff --git a/src/transport/MMCollisionInt.cpp b/src/transport/MMCollisionInt.cpp index 404b36b6d..274e73502 100644 --- a/src/transport/MMCollisionInt.cpp +++ b/src/transport/MMCollisionInt.cpp @@ -291,27 +291,21 @@ void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, i c.size(), DATA_PTR(c)); } m_o22poly.push_back(c); - if (rmserr > e22) { - e22 = rmserr; - } + e22 = std::max(e22, rmserr); rmserr = fitDelta(1, i, DeltaDegree, DATA_PTR(c)); m_apoly.push_back(c); if (DEBUG_MODE_ENABLED && log_level > 3) m_xml->XML_writeVector(*logfile, indent, "astar", c.size(), DATA_PTR(c)); - if (rmserr > ea) { - ea = rmserr; - } + ea = std::max(ea, rmserr); rmserr = fitDelta(2, i, DeltaDegree, DATA_PTR(c)); m_bpoly.push_back(c); if (DEBUG_MODE_ENABLED && log_level > 3) m_xml->XML_writeVector(*logfile, indent, "bstar", c.size(), DATA_PTR(c)); - if (rmserr > eb) { - eb = rmserr; - } + eb = std::max(eb, rmserr); rmserr = fitDelta(3, i, DeltaDegree, DATA_PTR(c)); m_cpoly.push_back(c); @@ -319,9 +313,7 @@ void MMCollisionInt::init(XML_Writer* xml, doublereal tsmin, doublereal tsmax, i m_xml->XML_writeVector(*logfile, indent, "cstar", c.size(), DATA_PTR(c)); } - if (rmserr > ec) { - ec = rmserr; - } + ec = std::max(ec, rmserr); if (DEBUG_MODE_ENABLED && log_level > 3) { m_xml->XML_close(*logfile, "dstar_fit"); @@ -372,10 +364,7 @@ doublereal MMCollisionInt::omega22(double ts, double deltastar) break; } int i1, i2; - i1 = i - 1; - if (i1 < 0) { - i1 = 0; - } + i1 = std::max(i - 1, 0); i2 = i1+3; if (i2 > 36) { i2 = 36; @@ -400,10 +389,7 @@ doublereal MMCollisionInt::astar(double ts, double deltastar) break; } int i1, i2; - i1 = i - 1; - if (i1 < 0) { - i1 = 0; - } + i1 = std::max(i - 1, 0); i2 = i1+3; if (i2 > 36) { i2 = 36; @@ -428,10 +414,7 @@ doublereal MMCollisionInt::bstar(double ts, double deltastar) break; } int i1, i2; - i1 = i - 1; - if (i1 < 0) { - i1 = 0; - } + i1 = std::max(i - 1, 0); i2 = i1+3; if (i2 > 36) { i2 = 36; @@ -456,10 +439,7 @@ doublereal MMCollisionInt::cstar(double ts, double deltastar) break; } int i1, i2; - i1 = i - 1; - if (i1 < 0) { - i1 = 0; - } + i1 = std::max(i - 1,0); i2 = i1+3; if (i2 > 36) { i2 = 36; diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index 380d237be..0c4dd4d6c 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -447,7 +447,7 @@ void TransportFactory::setupMM(std::ostream& flog, const std::vector tstar_max) { - tstar_max = ts2; - } + tstar_min = std::min(tstar_min, Boltzmann * tr.tmin/epsilon(i,j)); + tstar_max = std::max(tstar_max, Boltzmann * tr.tmax/epsilon(i,j)); // the effective dipole moment for (i,j) collisions tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j)); @@ -1300,12 +1294,8 @@ void TransportFactory::fitProperties(GasTransportParams& tr, } err = fit - val; relerr = err/val; - if (fabs(err) > mxerr) { - mxerr = fabs(err); - } - if (fabs(relerr) > mxrelerr) { - mxrelerr = fabs(relerr); - } + mxerr = std::max(mxerr, fabs(err)); + mxrelerr = std::max(mxrelerr, fabs(relerr)); } // evaluate max fit errors for conductivity @@ -1320,12 +1310,8 @@ void TransportFactory::fitProperties(GasTransportParams& tr, } err = fit - val; relerr = err/val; - if (fabs(err) > mxerr_cond) { - mxerr_cond = fabs(err); - } - if (fabs(relerr) > mxrelerr_cond) { - mxrelerr_cond = fabs(relerr); - } + mxerr_cond = std::max(mxerr_cond, fabs(err)); + mxrelerr_cond = std::max(mxrelerr_cond, fabs(relerr)); } tr.visccoeffs.push_back(c); tr.condcoeffs.push_back(c2); @@ -1428,12 +1414,8 @@ void TransportFactory::fitProperties(GasTransportParams& tr, } err = fit - val; relerr = err/val; - if (fabs(err) > mxerr) { - mxerr = fabs(err); - } - if (fabs(relerr) > mxrelerr) { - mxrelerr = fabs(relerr); - } + mxerr = std::max(mxerr, fabs(err)); + mxrelerr = std::max(mxrelerr, fabs(relerr)); } tr.diffcoeffs.push_back(c); if (DEBUG_MODE_ENABLED && tr.log_level >= 2 && m_verbose) { diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 373dba12d..c02454047 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -68,9 +68,7 @@ void ConstPressureReactor::initialize(doublereal t0) for (size_t m = 0; m < m_nwalls; m++) { if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); - if (nt > maxnt) { - maxnt = nt; - } + maxnt = std::max(maxnt, nt); if (m_wall[m]->kinetics(m_lr[m])) { if (&m_kin->thermo(0) != &m_wall[m]->kinetics(m_lr[m])->thermo(0)) { diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index a0f40f3a1..7bcf3d2b3 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -80,9 +80,7 @@ void IdealGasConstPressureReactor::initialize(doublereal t0) for (size_t m = 0; m < m_nwalls; m++) { if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); - if (nt > maxnt) { - maxnt = nt; - } + maxnt = std::max(maxnt, nt); if (m_wall[m]->kinetics(m_lr[m])) { if (&m_kin->thermo(0) != &m_wall[m]->kinetics(m_lr[m])->thermo(0)) { diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 40d1b5b87..91f59010e 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -83,9 +83,7 @@ void IdealGasReactor::initialize(doublereal t0) m_wall[m]->initialize(); if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); - if (nt > maxnt) { - maxnt = nt; - } + maxnt = std::max(maxnt, nt); if (m_wall[m]->kinetics(m_lr[m])) { if (&m_kin->thermo(0) != &m_wall[m]->kinetics(m_lr[m])->thermo(0)) { diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 33f24afb8..5e0739327 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -81,9 +81,7 @@ void Reactor::initialize(doublereal t0) m_wall[m]->initialize(); if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); - if (nt > maxnt) { - maxnt = nt; - } + maxnt = std::max(maxnt, nt); if (m_wall[m]->kinetics(m_lr[m])) { if (&m_kin->thermo(0) != &m_wall[m]->kinetics(m_lr[m])->thermo(0)) { diff --git a/test_problems/cathermo/HMW_dupl_test/HMW_dupl_test.cpp b/test_problems/cathermo/HMW_dupl_test/HMW_dupl_test.cpp index 3d6bceb90..230ac8a8f 100644 --- a/test_problems/cathermo/HMW_dupl_test/HMW_dupl_test.cpp +++ b/test_problems/cathermo/HMW_dupl_test/HMW_dupl_test.cpp @@ -152,14 +152,7 @@ int main(int argc, char** argv) /* * Make sure we are at the saturation pressure or above. */ - - double psat = HMW->satPressure(T); - - pres = OneAtm; - if (psat > pres) { - pres = psat; - } - + pres = std::max(HMW->satPressure(T), OneAtm); HMW->setState_TPM(T, pres, moll); diff --git a/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp b/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp index f56a2917e..a5cfef5e8 100644 --- a/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp +++ b/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp @@ -133,14 +133,7 @@ int main(int argc, char** argv) /* * Make sure we are at the saturation pressure or above. */ - - double psat = HMW->satPressure(T); - - pres = OneAtm; - if (psat > pres) { - pres = psat; - } - + pres = std::max(HMW->satPressure(T), OneAtm); HMW->setState_TPM(T, pres, moll); diff --git a/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp b/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp index 58584e479..4153e1042 100644 --- a/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp +++ b/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp @@ -125,11 +125,7 @@ int main(int argc, char** argv) double T = TTable.T[i]; double RT = GasConstant * T; - double psat = HMW->satPressure(T); - pres = OneAtm; - if (psat > pres) { - pres = psat; - } + pres = std::max(HMW->satPressure(T), OneAtm); HMW->setState_TPM(T, pres, moll); solid->setState_TP(T, pres); diff --git a/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp b/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp index 18e6977b3..90682c8c8 100644 --- a/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp +++ b/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp @@ -133,14 +133,7 @@ int main(int argc, char** argv) /* * Make sure we are at the saturation pressure or above. */ - - double psat = HMW->satPressure(T); - - pres = OneAtm; - if (psat > pres) { - pres = psat; - } - + pres = std::max(HMW->satPressure(T), OneAtm); HMW->setState_TPM(T, pres, moll); diff --git a/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp b/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp index bd93bc852..8d674274b 100644 --- a/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp +++ b/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp @@ -129,14 +129,7 @@ int main(int argc, char** argv) /* * Make sure we are at the saturation pressure or above. */ - - double psat = HMW->satPressure(T); - - pres = OneAtm; - if (psat > pres) { - pres = psat; - } - + pres = std::max(HMW->satPressure(T), OneAtm); HMW->setState_TPM(T, pres, moll); diff --git a/test_problems/surfSolverTest/surfaceSolver.cpp b/test_problems/surfSolverTest/surfaceSolver.cpp index a44fdcc5d..df1ddc005 100644 --- a/test_problems/surfSolverTest/surfaceSolver.cpp +++ b/test_problems/surfSolverTest/surfaceSolver.cpp @@ -369,11 +369,7 @@ int main(int argc, char** argv) */ double pres = gasTP->pressure(); gasTP->getMoleFractions(x); - double tmp = 0.3 * x[0]; - double tmp2 = 0.3 * x[1]; - if (tmp2 < tmp) { - tmp = tmp2; - } + double tmp = 0.3 * std::min(x[0], x[1]); x[0] += tmp; x[1] -= tmp; gasTP->setState_PX(pres, x); diff --git a/test_problems/surfSolverTest/surfaceSolver2.cpp b/test_problems/surfSolverTest/surfaceSolver2.cpp index 2f92138dc..e44e034d6 100644 --- a/test_problems/surfSolverTest/surfaceSolver2.cpp +++ b/test_problems/surfSolverTest/surfaceSolver2.cpp @@ -410,11 +410,7 @@ int main(int argc, char** argv) */ double pres = gasTP->pressure(); gasTP->getMoleFractions(x); - double tmp = 0.3 * x[0]; - double tmp2 = 0.3 * x[1]; - if (tmp2 < tmp) { - tmp = tmp2; - } + double tmp = 0.3 * std::min(x[0], x[1]); x[0] += tmp; x[1] -= tmp; gasTP->setState_PX(pres, x);