From 801334c84ea37778bc9c540c82483d0ac03ebd8c Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 31 Jul 2015 20:48:45 -0400 Subject: [PATCH] Merge nested if statements --- src/base/xml.cpp | 94 +++---- src/equil/ChemEquil.cpp | 82 +++--- src/equil/MultiPhaseEquil.cpp | 20 +- src/equil/vcs_Gibbs.cpp | 8 +- src/equil/vcs_MultiPhaseEquil.cpp | 24 +- src/equil/vcs_VolPhase.cpp | 59 ++-- src/equil/vcs_elem.cpp | 148 +++++----- src/equil/vcs_elem_rearrange.cpp | 8 +- src/equil/vcs_inest.cpp | 20 +- src/equil/vcs_phaseStability.cpp | 121 +++----- src/equil/vcs_prep.cpp | 6 +- src/equil/vcs_report.cpp | 7 +- src/equil/vcs_rxnadj.cpp | 18 +- src/equil/vcs_setMolesLinProg.cpp | 16 +- src/equil/vcs_solve.cpp | 12 +- src/equil/vcs_solve_TP.cpp | 265 ++++++++---------- src/kinetics/InterfaceKinetics.cpp | 56 ++-- src/kinetics/ReactionPath.cpp | 24 +- src/kinetics/importKinetics.cpp | 10 +- src/kinetics/solveSP.cpp | 12 +- src/numerics/IDA_Solver.cpp | 10 +- src/numerics/ResidJacEval.cpp | 8 +- src/numerics/RootFind.cpp | 94 +++---- src/oneD/MultiNewton.cpp | 14 +- src/thermo/DebyeHuckel.cpp | 60 ++-- src/thermo/HMWSoln.cpp | 128 ++++----- src/thermo/HMWSoln_input.cpp | 98 +++---- src/thermo/IdealMolalSoln.cpp | 8 +- src/thermo/IdealSolidSolnPhase.cpp | 8 +- src/thermo/IonsFromNeutralVPSSTP.cpp | 32 +-- src/thermo/LatticePhase.cpp | 6 +- src/thermo/MetalSHEelectrons.cpp | 8 +- src/thermo/MineralEQ3.cpp | 8 +- src/thermo/MixtureFugacityTP.cpp | 78 +++--- src/thermo/MolalityVPSSTP.cpp | 16 +- src/thermo/Mu0Poly.cpp | 8 +- src/thermo/PDSS_Water.cpp | 8 +- src/thermo/RedlichKwongMFTP.cpp | 18 +- src/thermo/ThermoFactory.cpp | 20 +- src/thermo/ThermoPhase.cpp | 6 +- src/thermo/VPSSMgrFactory.cpp | 27 +- src/thermo/VPSSMgr_IdealGas.cpp | 10 +- src/thermo/WaterProps.cpp | 14 +- src/transport/LTPspecies.cpp | 6 +- .../cathermo/HMW_graph_GvI/HMW_graph_GvI.cpp | 10 +- 45 files changed, 694 insertions(+), 1019 deletions(-) diff --git a/src/base/xml.cpp b/src/base/xml.cpp index f6930162c..2c0b97755 100644 --- a/src/base/xml.cpp +++ b/src/base/xml.cpp @@ -173,11 +173,9 @@ int XML_Reader::findQuotedString(const std::string& s, std::string& rstring) con ilocStart = iloc2; qtype = q2; } - if (iloc1 != string::npos) { - if (iloc1 < ilocStart) { - ilocStart = iloc1; - qtype = q1; - } + if (iloc1 != string::npos && iloc1 < ilocStart) { + ilocStart = iloc1; + qtype = q1; } if (qtype == ' ') { return 0; @@ -349,11 +347,9 @@ XML_Node& XML_Node::operator=(const XML_Node& right) { if (&right != this) { for (size_t i = 0; i < m_children.size(); i++) { - if (m_children[i]) { - if (m_children[i]->parent() == this) { - delete m_children[i]; - m_children[i] = 0; - } + if (m_children[i] && m_children[i]->parent() == this) { + delete m_children[i]; + m_children[i] = 0; } } m_children.resize(0); @@ -368,11 +364,9 @@ XML_Node::~XML_Node() writelog("XML_Node::~XML_Node: deleted a locked XML_Node: "+name()); } for (size_t i = 0; i < m_children.size(); i++) { - if (m_children[i]) { - if (m_children[i]->parent() == this) { - delete m_children[i]; - m_children[i] = 0; - } + if (m_children[i] && m_children[i]->parent() == this) { + delete m_children[i]; + m_children[i] = 0; } } } @@ -380,11 +374,9 @@ XML_Node::~XML_Node() void XML_Node::clear() { for (size_t i = 0; i < m_children.size(); i++) { - if (m_children[i]) { - if (m_children[i]->parent() == this) { - delete m_children[i]; - m_children[i] = 0; - } + if (m_children[i] && m_children[i]->parent() == this) { + delete m_children[i]; + m_children[i] = 0; } } m_value.clear(); @@ -593,10 +585,8 @@ bool XML_Node::isComment() const void XML_Node::_require(const std::string& a, const std::string& v) const { - if (hasAttrib(a)) { - if (attrib(a) == v) { - return; - } + if (hasAttrib(a) && attrib(a) == v) { + return; } string msg="XML_Node "+name()+" is required to have an attribute named " + a + " with the value \"" + v +"\", but instead the value is \"" + attrib(a); @@ -609,10 +599,8 @@ XML_Node* XML_Node::findNameID(const std::string& nameTarget, XML_Node* scResult = 0; XML_Node* sc; std::string idattrib = id(); - if (name() == nameTarget) { - if (idTarget == "" || idTarget == idattrib) { - return const_cast(this); - } + if (name() == nameTarget && (idTarget == "" || idTarget == idattrib)) { + return const_cast(this); } for (size_t n = 0; n < m_children.size(); n++) { sc = m_children[n]; @@ -645,12 +633,8 @@ XML_Node* XML_Node::findNameIDIndex(const std::string& nameTarget, std::string ii = attrib("index"); std::string index_s = int2str(index_i); int iMax = -1000000; - if (name() == nameTarget) { - if (idTarget == "" || idTarget == idattrib) { - if (index_s == ii) { - return const_cast(this); - } - } + if (name() == nameTarget && (idTarget == "" || idTarget == idattrib) && index_s == ii) { + return const_cast(this); } for (size_t n = 0; n < m_children.size(); n++) { sc = m_children[n]; @@ -658,10 +642,8 @@ XML_Node* XML_Node::findNameIDIndex(const std::string& nameTarget, ii = sc->attrib("index"); int indexR = atoi(ii.c_str()); idattrib = sc->id(); - if (idTarget == idattrib || idTarget == "") { - if (index_s == ii) { - return sc; - } + if ((idTarget == idattrib || idTarget == "") && index_s == ii) { + return sc; } if (indexR > iMax) { scResult = sc; @@ -674,10 +656,8 @@ XML_Node* XML_Node::findNameIDIndex(const std::string& nameTarget, XML_Node* XML_Node::findID(const std::string& id_, const int depth) const { - if (hasAttrib("id")) { - if (attrib("id") == id_) { - return const_cast(this); - } + if (hasAttrib("id") && attrib("id") == id_) { + return const_cast(this); } if (depth > 0) { XML_Node* r = 0; @@ -694,10 +674,8 @@ XML_Node* XML_Node::findID(const std::string& id_, const int depth) const XML_Node* XML_Node::findByAttr(const std::string& attr, const std::string& val, int depth) const { - if (hasAttrib(attr)) { - if (attrib(attr) == val) { - return const_cast(this); - } + if (hasAttrib(attr) && attrib(attr) == val) { + return const_cast(this); } if (depth > 0) { XML_Node* r = 0; @@ -829,25 +807,17 @@ void XML_Node::copyUnion(XML_Node* const node_dest) const for (size_t idc = 0; idc < ndc; idc++) { XML_Node* dcc = vsc[idc]; if (dcc->name() == sc->name()) { - if (sc->hasAttrib("id")) { - if (sc->attrib("id") != dcc->attrib("id")) { - break; - } + if (sc->hasAttrib("id") && sc->attrib("id") != dcc->attrib("id")) { + break; } - if (sc->hasAttrib("name")) { - if (sc->attrib("name") != dcc->attrib("name")) { - break; - } + if (sc->hasAttrib("name") && sc->attrib("name") != dcc->attrib("name")) { + break; } - if (sc->hasAttrib("model")) { - if (sc->attrib("model") != dcc->attrib("model")) { - break; - } + if (sc->hasAttrib("model") && sc->attrib("model") != dcc->attrib("model")) { + break; } - if (sc->hasAttrib("title")) { - if (sc->attrib("title") != dcc->attrib("title")) { - break; - } + if (sc->hasAttrib("title") && sc->attrib("title") != dcc->attrib("title")) { + break; } dc = vsc[idc]; } diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index 748729a80..30ea03698 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -755,10 +755,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, } } // Delta Damping - if (m == mm) { - if (fabs(res_trial[mm]) > 0.2) { - fctr = std::min(fctr, 0.2/fabs(res_trial[mm])); - } + if (m == mm && fabs(res_trial[mm]) > 0.2) { + fctr = std::min(fctr, 0.2/fabs(res_trial[mm])); } } if (fctr != 1.0 && DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { @@ -1057,10 +1055,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } for (m = 0; m < m_mm; m++) { - if (m != m_eloc) { - if (elMoles[m] <= options.absElemTol) { - x[m] = -200.; - } + if (m != m_eloc && elMoles[m] <= options.absElemTol) { + x[m] = -200.; } } @@ -1124,10 +1120,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, } } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { - if (!normalStep) { - writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); - } + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && !normalStep) { + writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); } if (!normalStep) { beta = 1.0; @@ -1135,16 +1129,14 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, for (im = 0; im < m_mm; im++) { m = m_orderVectorElements[im]; resid[m] = 0.0; - if (im < m_nComponents) { - if (elMoles[m] > 0.001 * elMolesTotal) { - if (eMolesCalc[m] > 1000. * elMoles[m]) { - resid[m] = -0.5; - resid[m_mm] -= 0.5; - } - if (1000 * eMolesCalc[m] < elMoles[m]) { - resid[m] = 0.5; - resid[m_mm] += 0.5; - } + if (im < m_nComponents && elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + resid[m] = -0.5; + resid[m_mm] -= 0.5; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + resid[m] = 0.5; + resid[m_mm] += 0.5; } } } @@ -1193,21 +1185,19 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, size_t kMSp2 = npos; int nSpeciesWithElem = 0; for (k = 0; k < m_kk; k++) { - if (n_i_calc[k] > nCutoff) { - if (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++) { - if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { - lumpSum[m] = 0; - break; - } + 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++) { + if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { + lumpSum[m] = 0; + break; } - } else { - kMSp = k; } + } else { + kMSp = k; } } } @@ -1443,10 +1433,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, beta = std::min(beta, -1.0 / resid[m]); } } - if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { - if (beta != 1.0) { - writelogf("(it %d) Beta = %g\n", iter, beta); - } + if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0 && beta != 1.0) { + writelogf("(it %d) Beta = %g\n", iter, beta); } } /* @@ -1498,17 +1486,13 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) double maxNegVal = -1.0; if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) { for (k = 0; k < m_kk; k++) { - if (nAtoms(k,m_eloc) > 0.0) { - if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { - maxPosVal = m_molefractions[k]; - maxPosEloc = k; - } + if (nAtoms(k,m_eloc) > 0.0 && m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { + maxPosVal = m_molefractions[k]; + maxPosEloc = k; } - if (nAtoms(k,m_eloc) < 0.0) { - if (m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) { - maxNegVal = m_molefractions[k]; - maxNegEloc = k; - } + if (nAtoms(k,m_eloc) < 0.0 && m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) { + maxNegVal = m_molefractions[k]; + maxNegEloc = k; } } } diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index d2634a95c..968c51973 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -41,13 +41,11 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_ // it from the calculation. Electrons are a special case, // since a species can have a negative number of 'atoms' // of electrons (positive ions). - if (m_mix->elementMoles(m) <= 0.0) { - if (m != m_eloc) { - m_incl_element[m] = 0; - for (k = 0; k < m_nsp_mix; k++) { - if (m_mix->nAtoms(k,m) != 0.0) { - m_incl_species[k] = 0; - } + if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) { + m_incl_element[m] = 0; + for (k = 0; k < m_nsp_mix; k++) { + if (m_mix->nAtoms(k,m) != 0.0) { + m_incl_species[k] = 0; } } } @@ -348,11 +346,9 @@ void MultiPhaseEquil::getComponents(const std::vector& order) doublereal maxmoles = -999.0; size_t kmax = 0; for (k = m+1; k < nColumns; k++) { - if (m_A(m,k) != 0.0) { - if (fabs(m_moles[m_order[k]]) > maxmoles) { - kmax = k; - maxmoles = fabs(m_moles[m_order[k]]); - } + if (m_A(m,k) != 0.0 && fabs(m_moles[m_order[k]]) > maxmoles) { + kmax = k; + maxmoles = fabs(m_moles[m_order[k]]); } } diff --git a/src/equil/vcs_Gibbs.cpp b/src/equil/vcs_Gibbs.cpp index 6602a0d0b..7c48ebe58 100644 --- a/src/equil/vcs_Gibbs.cpp +++ b/src/equil/vcs_Gibbs.cpp @@ -44,11 +44,9 @@ double VCS_SOLVE::vcs_GibbsPhase(size_t iphase, const double* const w, double g = 0.0; double phaseMols = 0.0; for (size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { - if (m_phaseID[kspec] == iphase) { - if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - g += w[kspec] * fe[kspec]; - phaseMols += w[kspec]; - } + if (m_phaseID[kspec] == iphase && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + g += w[kspec] * fe[kspec]; + phaseMols += w[kspec]; } } diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 3b9cca048..75aee1f77 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -321,13 +321,11 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, Tlow = Tnow; Slow = Snow; } else { - if (Slow > Starget) { - if (Snow < Slow) { - Thigh = Tlow; - Shigh = Slow; - Tlow = Tnow; - Slow = Snow; - } + if (Slow > Starget && Snow < Slow) { + Thigh = Tlow; + Shigh = Slow; + Tlow = Tnow; + Slow = Snow; } } } else { @@ -573,10 +571,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, } plogf("------------------------------------------" "-------------------\n"); - if (printLvl > 2) { - if (m_vsolve.m_timing_print_lvl > 0) { - plogf("Total time = %12.6e seconds\n", te); - } + if (printLvl > 2 && m_vsolve.m_timing_print_lvl > 0) { + plogf("Total time = %12.6e seconds\n", te); } } return iSuccess; @@ -1362,10 +1358,8 @@ int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int } plogf("------------------------------------------" "-------------------\n"); - if (printLvl > 2) { - if (m_vsolve.m_timing_print_lvl > 0) { - plogf("Total time = %12.6e seconds\n", te); - } + if (printLvl > 2 && m_vsolve.m_timing_print_lvl > 0) { + plogf("Total time = %12.6e seconds\n", te); } } return iStable; diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 605b324ed..82c4112c7 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -493,10 +493,8 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc, * then we have a valid state. If the phase went away, it would * be a valid starting point for F_k's. So, save the state. */ - if (stateCalc == VCS_STATECALC_OLD) { - if (v_totalMoles > 0.0) { - creationMoleNumbers_ = Xmol_; - } + if (stateCalc == VCS_STATECALC_OLD && v_totalMoles > 0.0) { + creationMoleNumbers_ = Xmol_; } /* @@ -528,12 +526,9 @@ void vcs_VolPhase::setMolesFromVCSCheck(const int vcsStateStatus, void vcs_VolPhase::updateFromVCS_MoleNumbers(const int vcsStateStatus) { - if (!m_UpToDate || (vcsStateStatus != m_vcsStateStatus)) { - if (vcsStateStatus == VCS_STATECALC_OLD || vcsStateStatus == VCS_STATECALC_NEW) { - if (m_owningSolverObject) { - setMolesFromVCS(vcsStateStatus); - } - } + if ((!m_UpToDate || vcsStateStatus != m_vcsStateStatus) && m_owningSolverObject && + (vcsStateStatus == VCS_STATECALC_OLD || vcsStateStatus == VCS_STATECALC_NEW)) { + setMolesFromVCS(vcsStateStatus); } } @@ -591,10 +586,8 @@ double vcs_VolPhase::electricPotential() const void vcs_VolPhase::setState_TP(const double temp, const double pres) { - if (Temp_ == temp) { - if (Pres_ == pres) { - return; - } + if (Temp_ == temp && Pres_ == pres) { + return; } TP_ptr->setElectricPotential(m_phi); TP_ptr->setState_TP(temp, pres); @@ -906,10 +899,8 @@ void vcs_VolPhase::setPhiVarIndex(size_t phiVarIndex) { m_phiVarIndex = phiVarIndex; m_speciesUnknownType[m_phiVarIndex] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE; - if (m_singleSpecies) { - if (m_phiVarIndex == 0) { - m_existence = VCS_PHASE_EXIST_ALWAYS; - } + if (m_singleSpecies && m_phiVarIndex == 0) { + m_existence = VCS_PHASE_EXIST_ALWAYS; } } @@ -935,20 +926,14 @@ void vcs_VolPhase::setExistence(const int existence) } } } else if (DEBUG_MODE_ENABLED && m_totalMolesInert == 0.0) { - if (v_totalMoles == 0.0) { - if (!m_singleSpecies || m_phiVarIndex != 0) { - throw CanteraError("vcs_VolPhase::setExistence", - "setting true existence for phase with no moles"); - } + if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) { + throw CanteraError("vcs_VolPhase::setExistence", + "setting true existence for phase with no moles"); } } - if (DEBUG_MODE_ENABLED && m_singleSpecies) { - if (m_phiVarIndex == 0) { - if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) { - throw CanteraError("vcs_VolPhase::setExistence", - "Trying to set existence of an electron phase to false"); - } - } + if (DEBUG_MODE_ENABLED && m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) { + throw CanteraError("vcs_VolPhase::setExistence", + "Trying to set existence of an electron phase to false"); } m_existence = existence; } @@ -1040,10 +1025,8 @@ static bool hasChargedSpecies(const ThermoPhase* const tPhase) static bool chargeNeutralityElement(const ThermoPhase* const tPhase) { int hasCharge = hasChargedSpecies(tPhase); - if (tPhase->chargeNeutralityNecessary()) { - if (hasCharge) { - return true; - } + if (tPhase->chargeNeutralityNecessary() && hasCharge) { + return true; } return false; } @@ -1155,11 +1138,9 @@ size_t vcs_VolPhase::transferElementsFM(const ThermoPhase* const tPhase) * The logic isn't set in stone, and is just for a particular type * of problem that I'm solving first. */ - if (ns == 1) { - if (tPhase->charge(0) != 0.0) { - m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE; - setPhiVarIndex(0); - } + if (ns == 1 && tPhase->charge(0) != 0.0) { + m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE; + setPhiVarIndex(0); } return ne; diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index d0530acb4..8f5e3d0e7 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -32,54 +32,52 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) * Require 12 digits of accuracy on non-zero constraints. */ for (size_t i = 0; i < top; ++i) { - if (m_elementActive[i]) { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) { + if (m_elementActive[i] && fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > fabs(m_elemAbundancesGoal[i]) * 1.0e-12) { + /* + * This logic is for charge neutrality condition + */ + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY && + m_elemAbundancesGoal[i] != 0.0) { + throw CanteraError("VCS_SOLVE::vcs_elabcheck", + "Problem with charge neutrality condition"); + } + if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { + double scale = VCS_DELETE_MINORSPECIES_CUTOFF; /* - * This logic is for charge neutrality condition + * Find out if the constraint is a multisign constraint. + * If it is, then we have to worry about roundoff error + * in the addition of terms. We are limited to 13 + * digits of finite arithmetic accuracy. */ - if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY && - m_elemAbundancesGoal[i] != 0.0) { - throw CanteraError("VCS_SOLVE::vcs_elabcheck", - "Problem with charge neutrality condition"); - } - if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { - double scale = VCS_DELETE_MINORSPECIES_CUTOFF; - /* - * Find out if the constraint is a multisign constraint. - * If it is, then we have to worry about roundoff error - * in the addition of terms. We are limited to 13 - * digits of finite arithmetic accuracy. - */ - bool multisign = false; - for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { - double eval = m_formulaMatrix(kspec,i); - if (eval < 0.0) { - multisign = true; - } - if (eval != 0.0) { - scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec])); - } + bool multisign = false; + for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { + double eval = m_formulaMatrix(kspec,i); + if (eval < 0.0) { + multisign = true; } - if (multisign) { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { - return false; - } - } else { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { - return false; - } + if (eval != 0.0) { + scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec])); + } + } + if (multisign) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { + return false; } } else { - /* - * For normal element balances, we require absolute compliance - * even for ridiculously small numbers. - */ - if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { - return false; - } else { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { return false; } } + } else { + /* + * For normal element balances, we require absolute compliance + * even for ridiculously small numbers. + */ + if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { + return false; + } else { + return false; + } } } } @@ -91,10 +89,8 @@ void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) for (size_t j = 0; j < m_numElemConstraints; ++j) { elemAbundPhase[j] = 0.0; for (size_t i = 0; i < m_numSpeciesTot; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (m_phaseID[i] == iphase) { - elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; - } + if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_phaseID[i] == iphase) { + elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; } } } @@ -385,21 +381,17 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { - if (m_elemAbundances[i] > 0.0) { - if (m_formulaMatrix(kspec,i) < 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i) ; - m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); - vcs_elab(); - break; - } + if (m_elemAbundances[i] > 0.0 && m_formulaMatrix(kspec,i) < 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i) ; + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); + vcs_elab(); + break; } - if (m_elemAbundances[i] < 0.0) { - if (m_formulaMatrix(kspec,i) > 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i); - m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); - vcs_elab(); - break; - } + if (m_elemAbundances[i] < 0.0 && m_formulaMatrix(kspec,i) > 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i); + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); + vcs_elab(); + break; } } } @@ -420,38 +412,30 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) bool useZeroed = true; for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (dev < 0.0) { - if (m_formulaMatrix(kspec,i) < 0.0) { - if (m_molNumSpecies_old[kspec] > 0.0) { - useZeroed = false; - } + if (m_formulaMatrix(kspec,i) < 0.0 && m_molNumSpecies_old[kspec] > 0.0) { + useZeroed = false; } } else { - if (m_formulaMatrix(kspec,i) > 0.0) { - if (m_molNumSpecies_old[kspec] > 0.0) { - useZeroed = false; - } + if (m_formulaMatrix(kspec,i) > 0.0 && m_molNumSpecies_old[kspec] > 0.0) { + useZeroed = false; } } } for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) { - if (dev < 0.0) { - if (m_formulaMatrix(kspec,i) < 0.0) { - double delta = dev / m_formulaMatrix(kspec,i) ; - m_molNumSpecies_old[kspec] += delta; - m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); - vcs_elab(); - break; - } + if (dev < 0.0 && m_formulaMatrix(kspec,i) < 0.0) { + double delta = dev / m_formulaMatrix(kspec,i) ; + m_molNumSpecies_old[kspec] += delta; + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); + vcs_elab(); + break; } - if (dev > 0.0) { - if (m_formulaMatrix(kspec,i) > 0.0) { - double delta = dev / m_formulaMatrix(kspec,i) ; - m_molNumSpecies_old[kspec] += delta; - m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); - vcs_elab(); - break; - } + if (dev > 0.0 && m_formulaMatrix(kspec,i) > 0.0) { + double delta = dev / m_formulaMatrix(kspec,i) ; + m_molNumSpecies_old[kspec] += delta; + m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0); + vcs_elab(); + break; } } } diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 6e6cef0b8..ade50444b 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -68,11 +68,9 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, */ k = m_numElemConstraints; for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) { - if (m_elementActive[ielem]) { - if (aw[ielem] != test) { - k = ielem; - break; - } + if (m_elementActive[ielem] && aw[ielem] != test) { + k = ielem; + break; } } if (k == m_numElemConstraints) { diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 7a972bc05..583566e24 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -212,10 +212,9 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, /* *********************************************************** */ double par = 0.5; for (size_t kspec = 0; kspec < m_numComponents; ++kspec) { - if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) { - par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]; - } + if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && + par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) { + par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]; } } par = 1. / par; @@ -239,10 +238,9 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, } } for (size_t kspec = m_numComponents; kspec < nspecies; ++kspec) { - if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (m_deltaMolNumSpecies[kspec] != 0.0) { - m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par; - } + if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && + m_deltaMolNumSpecies[kspec] != 0.0) { + m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par; } } /* @@ -265,10 +263,8 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, if (s == 0.0) { break; } - if (s < 0.0) { - if (ikl == 0) { - break; - } + if (s < 0.0 && ikl == 0) { + break; } /* ***************************************** */ /* *** TRY HALF STEP SIZE ****************** */ diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 3ddfaa342..788182427 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -70,10 +70,8 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const foundJrxn = true; // We can do the reaction if all other reactant components have positive mole fractions for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) { - if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0) { - if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { - foundJrxn = false; - } + if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { + foundJrxn = false; } } if (foundJrxn) { @@ -89,10 +87,8 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const } // We can do the backwards reaction if all of the product components species are positive for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) { - if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0) { - if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { - foundJrxn = false; - } + if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { + foundJrxn = false; } } if (foundJrxn) { @@ -125,23 +121,18 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() * The logic below calculates zeroedComponentLinkedPhasePops */ for (size_t j = 0; j < m_numComponents; j++) { - if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - if (m_molNumSpecies_old[j] <= 0.0) { - std::vector &jList = zeroedComponentLinkedPhasePops[j]; - size_t iph = m_phaseID[j]; - jList.push_back(iph); - for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) { - size_t kspec = irxn + m_numComponents; - iph = m_phaseID[kspec]; - vcs_VolPhase* Vphase = m_VolPhaseList[iph]; - int existence = Vphase->exists(); - if (existence < 0) { - if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) { - if (std::find(jList.begin(), jList.end(), iph) != jList.end()) { - jList.push_back(iph); - } - } - } + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS && m_molNumSpecies_old[j] <= 0.0) { + std::vector &jList = zeroedComponentLinkedPhasePops[j]; + size_t iph = m_phaseID[j]; + jList.push_back(iph); + for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) { + size_t kspec = irxn + m_numComponents; + iph = m_phaseID[kspec]; + vcs_VolPhase* Vphase = m_VolPhaseList[iph]; + int existence = Vphase->exists(); + if (existence < 0 && m_stoichCoeffRxnMatrix(j,irxn) > 0.0 && + std::find(jList.begin(), jList.end(), iph) != jList.end()) { + jList.push_back(iph); } } } @@ -168,26 +159,20 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() size_t kspec = Vphase->spGlobalIndexVCS(k); size_t irxn = kspec - m_numComponents; for (size_t j = 0; j < m_numComponents; j++) { - if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - if (m_molNumSpecies_old[j] <= 0.0) { - if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) { - bool foundPos = false; - for (size_t kk = 0; kk < nsp; kk++) { - size_t kkspec = Vphase->spGlobalIndexVCS(kk); - if (kkspec >= m_numComponents) { - size_t iirxn = kkspec - m_numComponents; - if (m_stoichCoeffRxnMatrix(j,iirxn) > 0.0) { - foundPos = true; - } - } - } - if (!foundPos) { - if (std::find(iphList.begin(), iphList.end(), j) != iphList.end()) { - iphList.push_back(j); - } + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS && m_molNumSpecies_old[j] <= 0.0 && m_stoichCoeffRxnMatrix(j,irxn) < 0.0) { + bool foundPos = false; + for (size_t kk = 0; kk < nsp; kk++) { + size_t kkspec = Vphase->spGlobalIndexVCS(kk); + if (kkspec >= m_numComponents) { + size_t iirxn = kkspec - m_numComponents; + if (m_stoichCoeffRxnMatrix(j,iirxn) > 0.0) { + foundPos = true; } } } + if (!foundPos && std::find(iphList.begin(), iphList.end(), j) != iphList.end()) { + iphList.push_back(j); + } } } } @@ -350,18 +335,14 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) if (Vphase->m_singleSpecies) { double s = 0.0; for (size_t j = 0; j < m_numComponents; ++j) { - if (!m_SSPhase[j]) { - if (m_molNumSpecies_old[j] > 0.0) { - s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; - } + if (!m_SSPhase[j] && m_molNumSpecies_old[j] > 0.0) { + s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; } } for (size_t j = 0; j < m_numPhases; j++) { Vphase = m_VolPhaseList[j]; - if (! Vphase->m_singleSpecies) { - if (m_tPhaseMoles_old[j] > 0.0) { - s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; - } + if (! Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) { + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } if (s != 0.0) { @@ -379,15 +360,13 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) */ for (size_t j = 0; j < m_numComponents; ++j) { double stoicC = m_stoichCoeffRxnMatrix(j,irxn); - if (stoicC != 0.0) { - if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec]; - if (negChangeComp > m_molNumSpecies_old[j]) { - if (m_molNumSpecies_old[j] > 0.0) { - m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC; - } else { - m_deltaMolNumSpecies[kspec] = 0.0; - } + if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec]; + if (negChangeComp > m_molNumSpecies_old[j]) { + if (m_molNumSpecies_old[j] > 0.0) { + m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC; + } else { + m_deltaMolNumSpecies[kspec] = 0.0; } } } @@ -421,10 +400,8 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) irxn = kspec - m_numComponents; for (size_t j = 0; j < m_numComponents; ++j) { double stoicC = m_stoichCoeffRxnMatrix(j,irxn); - if (stoicC != 0.0) { - if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - molNumSpecies_tmp[j] += stoicC * delmol; - } + if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + molNumSpecies_tmp[j] += stoicC * delmol; } } } @@ -456,10 +433,8 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) // Here we create a damp > 1 to account for this possibility. // We adjust upwards to make sure that a component in an existing multispecies // phase is modified by a factor of 1/1000. - if (ratioComp > 1.0E-30) { - if (ratioComp < 0.001) { - damp = 0.001 / ratioComp; - } + if (ratioComp > 1.0E-30 && ratioComp < 0.001) { + damp = 0.001 / ratioComp; } if (damp <= 1.0E-6) { return 3; @@ -721,15 +696,11 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]), 1.0E-8/fabs(delFrac[k])); } - if (delFrac[k] < 0.0) { - if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) { - damp = fracDelta_old[k] / (2.0 * -delFrac[k]); - } + if (delFrac[k] < 0.0 && 2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) { + damp = fracDelta_old[k] / (2.0 * -delFrac[k]); } - if (delFrac[k] > 0.0) { - if (2.0 * damp * delFrac[k] > fracDelta_old[k]) { - damp = fracDelta_old[k] / (2.0 * delFrac[k]); - } + if (delFrac[k] > 0.0 && 2.0 * damp * delFrac[k] > fracDelta_old[k]) { + damp = fracDelta_old[k] / (2.0 * delFrac[k]); } } damp = std::max(damp, 0.000001); diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index f993d47a4..8bc6f5195 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -32,10 +32,8 @@ void VCS_SOLVE::vcs_SSPhase() if (TPhInertMoles[iph] > 0.0) { Vphase->setExistence(2); } - if (numPhSpecies[iph] <= 1) { - if (TPhInertMoles[iph] == 0.0) { - Vphase->m_singleSpecies = true; - } + if (numPhSpecies[iph] <= 1 && TPhInertMoles[iph] == 0.0) { + Vphase->m_singleSpecies = true; } } diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index f53a175cd..04e8d91ef 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -219,10 +219,9 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("%-12.12s |",VPhase->PhaseName.c_str()); plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale); totalMoles += m_tPhaseMoles_old[iphase]; - if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles()) { - if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) { - throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem"); - } + if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() && + !vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) { + throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem"); } vcs_elabPhase(iphase, &gaPhase[0]); for (size_t j = 0; j < m_numElemConstraints; j++) { diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 1bbd55cc9..ef2e14bf3 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -166,18 +166,14 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) s = 1.0 / m_molNumSpecies_old[kspec]; } for (size_t j = 0; j < m_numComponents; ++j) { - if (!m_SSPhase[j]) { - if (m_molNumSpecies_old[j] > 0.0) { - s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; - } + if (!m_SSPhase[j] && m_molNumSpecies_old[j] > 0.0) { + s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j]; } } for (size_t j = 0; j < m_numPhases; j++) { vcs_VolPhase* Vphase = m_VolPhaseList[j]; - if (!Vphase->m_singleSpecies) { - if (m_tPhaseMoles_old[j] > 0.0) { - s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; - } + if (!Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) { + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } if (s != 0.0) { @@ -455,10 +451,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg() } } for (size_t j = 0; j < m_numPhases; j++) { - if (!m_VolPhaseList[j]->m_singleSpecies) { - if (m_tPhaseMoles_old[j] > 0.0) { - s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; - } + if (!m_VolPhaseList[j]->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) { + s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j]; } } if (s != 0.0) { diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 7399d049e..5c69b08cc 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -128,13 +128,11 @@ int VCS_SOLVE::vcs_setMolesLinProg() delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu); // if a component has nearly zero moles, redo // with a new set of components - if (!redo) { - if (delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { - plogf(" --- Component too small: %s\n", m_speciesName[jcomp].c_str()); - } - redo = true; + if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) { + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" --- Component too small: %s\n", m_speciesName[jcomp].c_str()); } + redo = true; } dxi_min = std::min(dxi_min, delta_xi); } @@ -153,10 +151,8 @@ int VCS_SOLVE::vcs_setMolesLinProg() } m_molNumSpecies_old[jcomp] += sc_irxn[jcomp] * dsLocal; m_molNumSpecies_old[jcomp] = max(0.0, m_molNumSpecies_old[jcomp]); - if (full) { - if (m_molNumSpecies_old[jcomp] < 1.0E-60) { - redo = true; - } + if (full && m_molNumSpecies_old[jcomp] < 1.0E-60) { + redo = true; } } } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 7903b7ca5..1b859488a 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -484,10 +484,8 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) if (pub->gai.size() != 0) { for (size_t i = 0; i < nelements; i++) { m_elemAbundancesGoal[i] = pub->gai[i]; - if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO) { - if (m_elemAbundancesGoal[i] < 1.0E-10) { - m_elemAbundancesGoal[i] = 0.0; - } + if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) { + m_elemAbundancesGoal[i] = 0.0; } } } else { @@ -501,10 +499,8 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } } - if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) { - if (m_elemAbundancesGoal[j] < 1.0E-10 * sum) { - m_elemAbundancesGoal[j] = 0.0; - } + if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10 * sum) { + m_elemAbundancesGoal[j] = 0.0; } } } else { diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index af6b25502..2fd3eb77d 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -908,12 +908,10 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, if (DEBUG_MODE_ENABLED) { doPhaseDeleteKspec = kspec; - if (m_debug_print_lvl >= 2) { - if (m_speciesStatus[kspec] >= 0) { - plogf(" --- SS species changed to zeroedss: "); - plogf("%-12s", m_speciesName[kspec].c_str()); - plogendl(); - } + if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) { + plogf(" --- SS species changed to zeroedss: "); + plogf("%-12s", m_speciesName[kspec].c_str()); + plogendl(); } } m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; @@ -1387,26 +1385,22 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, size_t kspec = m_indexRxnToSpecies[irxn]; int speciesType = vcs_species_type(kspec); if (speciesType < VCS_SPECIES_MINOR) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { - if (m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { - plogf(" --- major/minor species is now zeroed out: %s\n", - m_speciesName[kspec].c_str()); - } + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { + plogf(" --- major/minor species is now zeroed out: %s\n", + m_speciesName[kspec].c_str()); } ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MINOR) { - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { - if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { - if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) { - plogf(" --- Noncomponent turned from major to minor: "); - } else if (kspec < m_numComponents) { - plogf(" --- Component turned into a minor species: "); - } else { - plogf(" --- Zeroed Species turned into a " - "minor species: "); - } - plogf("%s\n", m_speciesName[kspec].c_str()); + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2 && m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { + if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) { + plogf(" --- Noncomponent turned from major to minor: "); + } else if (kspec < m_numComponents) { + plogf(" --- Component turned into a minor species: "); + } else { + plogf(" --- Zeroed Species turned into a " + "minor species: "); } + plogf("%s\n", m_speciesName[kspec].c_str()); } ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MAJOR) { @@ -1776,11 +1770,9 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) * If the component has a zero concentration and is a reactant * in the formation reaction, then dx == 0.0, and we just return. */ - if (m_molNumSpecies_old[j] <= 0.0) { - if (sc_irxn[j] < 0.0) { - *delta_ptr = 0.0; - return 0; - } + if (m_molNumSpecies_old[j] <= 0.0 && sc_irxn[j] < 0.0) { + *delta_ptr = 0.0; + return 0; } } /* @@ -1878,23 +1870,18 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec) * Check to see whether we have just annihilated a multispecies phase. * If it is extinct, call the delete_multiphase() function. */ - if (! m_SSPhase[klast]) { - if (Vphase->exists() != VCS_PHASE_EXIST_ALWAYS) { - bool stillExists = false; - for (size_t k = 0; k < m_numSpeciesRdc; k++) { - if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (m_phaseID[k] == iph) { - if (m_molNumSpecies_old[k] > 0.0) { - stillExists = true; - break; - } - } - } - } - if (!stillExists) { - vcs_delete_multiphase(iph); + if (! m_SSPhase[klast] && Vphase->exists() != VCS_PHASE_EXIST_ALWAYS) { + bool stillExists = false; + for (size_t k = 0; k < m_numSpeciesRdc; k++) { + if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && + m_phaseID[k] == iph && m_molNumSpecies_old[k] > 0.0) { + stillExists = true; + break; } } + if (!stillExists) { + vcs_delete_multiphase(iph); + } } /* * When the total number of noncomponent species is zero, we @@ -1940,10 +1927,8 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) if (Vphase->exists() == VCS_PHASE_EXIST_NO) { Vphase->setExistence(VCS_PHASE_EXIST_YES); for (size_t k = 0; k < m_numSpeciesTot; k++) { - if (m_phaseID[k] == iph) { - if (m_speciesStatus[k] != VCS_SPECIES_DELETED) { - m_speciesStatus[k] = VCS_SPECIES_MINOR; - } + if (m_phaseID[k] == iph && m_speciesStatus[k] != VCS_SPECIES_DELETED) { + m_speciesStatus[k] = VCS_SPECIES_MINOR; } } } @@ -2310,18 +2295,16 @@ size_t VCS_SOLVE::vcs_add_all_deleted() for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { size_t kspec = m_indexRxnToSpecies[irxn]; size_t iph = m_phaseID[kspec]; - if (m_tPhaseMoles_old[iph] > 0.0) { - if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) { - if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > - VCS_DELETE_MINORSPECIES_CUTOFF) || - (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { - retn++; - if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { - plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", - m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], - m_deltaGRxn_old[irxn]); - plogendl(); - } + if (m_tPhaseMoles_old[iph] > 0.0 && fabs(m_deltaGRxn_old[irxn]) > m_tolmin) { + if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > + VCS_DELETE_MINORSPECIES_CUTOFF) || + (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { + retn++; + if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { + plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", + m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], + m_deltaGRxn_old[irxn]); + plogendl(); } } } @@ -2587,35 +2570,31 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], int minNonZeroes = 100000; int nonZeroesKspec = 0; for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) { - if (aw[kspec] >= 0.0) { - if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - maxConcPossKspec = 1.0E10; - nonZeroesKspec = 0; - for (size_t j = 0; j < m_numElemConstraints; ++j) { - if (m_elementActive[j]) { - if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { - double nu = m_formulaMatrix(kspec,j); - if (nu != 0.0) { - nonZeroesKspec++; - maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec); - } - } + if (aw[kspec] >= 0.0 && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + maxConcPossKspec = 1.0E10; + nonZeroesKspec = 0; + for (size_t j = 0; j < m_numElemConstraints; ++j) { + if (m_elementActive[j] && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + double nu = m_formulaMatrix(kspec,j); + if (nu != 0.0) { + nonZeroesKspec++; + maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec); } } - if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) { - if (nonZeroesKspec <= minNonZeroes) { - if (kfound == npos || nonZeroesKspec < minNonZeroes) { + } + if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) { + if (nonZeroesKspec <= minNonZeroes) { + if (kfound == npos || nonZeroesKspec < minNonZeroes) { + kfound = kspec; + } else { + // ok we are sitting pretty equal here decide on the raw ss Gibbs energy + if (m_SSfeSpecies[kspec] <= m_SSfeSpecies[kfound]) { kfound = kspec; - } else { - // ok we are sitting pretty equal here decide on the raw ss Gibbs energy - if (m_SSfeSpecies[kspec] <= m_SSfeSpecies[kfound]) { - kfound = kspec; - } } } - minNonZeroes = std::min(minNonZeroes, nonZeroesKspec); - maxConcPoss = std::max(maxConcPoss, maxConcPossKspec); } + minNonZeroes = std::min(minNonZeroes, nonZeroesKspec); + maxConcPoss = std::max(maxConcPoss, maxConcPossKspec); } } } @@ -2812,17 +2791,13 @@ L_END_LOOP: juse = npos; jlose = npos; for (size_t j = 0; j < m_numElemConstraints; j++) { - if (!m_elementActive[j]) { - if (!strcmp(m_elementName[j].c_str(), "E")) { - juse = j; - } + if (!m_elementActive[j] && !strcmp(m_elementName[j].c_str(), "E")) { + juse = j; } } for (size_t j = 0; j < m_numElemConstraints; j++) { - if (m_elementActive[j]) { - if (!strncmp((m_elementName[j]).c_str(), "cn_", 3)) { - jlose = j; - } + if (m_elementActive[j] && !strncmp((m_elementName[j]).c_str(), "cn_", 3)) { + jlose = j; } } for (k = 0; k < m_numSpeciesTot; k++) { @@ -3029,10 +3004,8 @@ size_t VCS_SOLVE::vcs_basisOptMax(const double* const molNum, const size_t j, bool doSwap = false; if (m_SSPhase[j]) { doSwap = (molNum[i] * m_spSize[i]) > big; - if (!m_SSPhase[i]) { - if (doSwap) { - doSwap = molNum[i] > (molNum[largest] * 1.001); - } + if (!m_SSPhase[i] && doSwap) { + doSwap = molNum[i] > (molNum[largest] * 1.001); } } else { if (m_SSPhase[i]) { @@ -3066,10 +3039,8 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const // ---------- Treat zeroed out species first ---------------- if (m_molNumSpecies_old[kspec] <= 0.0) { - if (m_tPhaseMoles_old[iph] <= 0.0) { - if (!m_SSPhase[kspec]) { - return VCS_SPECIES_ZEROEDMS; - } + if (m_tPhaseMoles_old[iph] <= 0.0 && !m_SSPhase[kspec]) { + return VCS_SPECIES_ZEROEDMS; } /* @@ -3145,21 +3116,19 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const } } - if (irxn >= 0) { - if (m_deltaGRxn_old[irxn] >= 0.0) { - /* - * We are here when the species is or should remain zeroed out - */ - if (m_SSPhase[kspec]) { - return VCS_SPECIES_ZEROEDSS; + if (irxn >= 0 && m_deltaGRxn_old[irxn] >= 0.0) { + /* + * We are here when the species is or should remain zeroed out + */ + if (m_SSPhase[kspec]) { + return VCS_SPECIES_ZEROEDSS; + } else { + if (phaseExist >= VCS_PHASE_EXIST_YES) { + return VCS_SPECIES_ACTIVEBUTZERO; + } else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) { + return VCS_SPECIES_ZEROEDPHASE; } else { - if (phaseExist >= VCS_PHASE_EXIST_YES) { - return VCS_SPECIES_ACTIVEBUTZERO; - } else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) { - return VCS_SPECIES_ZEROEDPHASE; - } else { - return VCS_SPECIES_ZEROEDMS; - } + return VCS_SPECIES_ZEROEDMS; } } } @@ -3226,12 +3195,8 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const } else { double szAdj = m_scSize[irxn] * std::sqrt((double)m_numRxnTot); for (size_t k = 0; k < m_numComponents; ++k) { - if (!m_SSPhase[k]) { - if (m_stoichCoeffRxnMatrix(k,irxn) != 0.0) { - if (m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) { - return VCS_SPECIES_MAJOR; - } - } + if (!m_SSPhase[k] && m_stoichCoeffRxnMatrix(k,irxn) != 0.0 && m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) { + return VCS_SPECIES_MAJOR; } } } @@ -3689,10 +3654,8 @@ void VCS_SOLVE::check_tmoles() const double m_tPhaseMoles_old_a = TPhInertMoles[i]; for (size_t k = 0; k < m_numSpeciesTot; k++) { - if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM) { - if (m_phaseID[k] == i) { - m_tPhaseMoles_old_a += m_molNumSpecies_old[k]; - } + if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM && m_phaseID[k] == i) { + m_tPhaseMoles_old_a += m_molNumSpecies_old[k]; } } sum += m_tPhaseMoles_old_a; @@ -3784,10 +3747,8 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() } } } - if (kspec >= m_numComponents) { - if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { - ++m_numRxnMinorZeroed; - } + if (kspec >= m_numComponents && m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { + ++m_numRxnMinorZeroed; } } if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) { @@ -4163,15 +4124,13 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, bool zeroedPhase = true; for (size_t irxn = 0; irxn < irxnl; ++irxn) { size_t kspec = m_indexRxnToSpecies[irxn]; - if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (m_phaseID[kspec] == iphase) { - if (m_molNumSpecies_old[kspec] > 0.0) { - zeroedPhase = false; - } - deltaGRxn[irxn] = feSpecies[kspec]; - for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { - deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp]; - } + if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_phaseID[kspec] == iphase) { + if (m_molNumSpecies_old[kspec] > 0.0) { + zeroedPhase = false; + } + deltaGRxn[irxn] = feSpecies[kspec]; + for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) { + deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp]; } } } @@ -4214,24 +4173,22 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, * All of the dg[]'s will be equal. If dg[] is negative, then * the phase will come back into existence. */ - if (alterZeroedPhases) { - if (zeroedPhase) { - double phaseDG = 1.0; - for (size_t irxn = 0; irxn < irxnl; ++irxn) { - size_t kspec = m_indexRxnToSpecies[irxn]; - if (m_phaseID[kspec] == iphase) { - deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0); - phaseDG -= exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec]; - } + if (alterZeroedPhases && zeroedPhase) { + double phaseDG = 1.0; + for (size_t irxn = 0; irxn < irxnl; ++irxn) { + size_t kspec = m_indexRxnToSpecies[irxn]; + if (m_phaseID[kspec] == iphase) { + deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0); + phaseDG -= exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec]; } - /* - * Overwrite the individual dg's with the phase DG. - */ - for (size_t irxn = 0; irxn < irxnl; ++irxn) { - size_t kspec = m_indexRxnToSpecies[irxn]; - if (m_phaseID[kspec] == iphase) { - deltaGRxn[irxn] = 1.0 - phaseDG; - } + } + /* + * Overwrite the individual dg's with the phase DG. + */ + for (size_t irxn = 0; irxn < irxnl; ++irxn) { + size_t kspec = m_indexRxnToSpecies[irxn]; + if (m_phaseID[kspec] == iphase) { + deltaGRxn[irxn] = 1.0 - phaseDG; } } } @@ -4373,10 +4330,8 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) dx = std::min(dx, - 0.3333* m_molNumSpecies_old[j] / sc_irxn[j]); } } - if (m_molNumSpecies_old[j] <= 0.0) { - if (sc_irxn[j] < 0.0) { - dx = 0.0; - } + if (m_molNumSpecies_old[j] <= 0.0 && sc_irxn[j] < 0.0) { + dx = 0.0; } } } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index de20d2366..9dbbb3fce 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -510,52 +510,40 @@ void InterfaceKinetics::updateROP() for (size_t j = 0; j != nReactions(); ++j) { if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) { for (size_t p = 0; p < nPhases(); p++) { - if (m_rxnPhaseIsProduct[j][p]) { - if (! m_phaseExists[p]) { - m_ropnet[j] = 0.0; - m_ropr[j] = m_ropf[j]; - if (m_ropf[j] > 0.0) { - for (size_t rp = 0; rp < nPhases(); rp++) { - if (m_rxnPhaseIsReactant[j][rp]) { - if (! m_phaseExists[rp]) { - m_ropnet[j] = 0.0; - m_ropr[j] = m_ropf[j] = 0.0; - } - } + if (m_rxnPhaseIsProduct[j][p] && !m_phaseExists[p]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j]; + if (m_ropf[j] > 0.0) { + for (size_t rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsReactant[j][rp] && !m_phaseExists[rp]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j] = 0.0; } } } } - if (m_rxnPhaseIsReactant[j][p]) { - if (! m_phaseIsStable[p]) { - m_ropnet[j] = 0.0; - m_ropr[j] = m_ropf[j]; - } + if (m_rxnPhaseIsReactant[j][p] && !m_phaseIsStable[p]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j]; } } } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) { for (size_t p = 0; p < nPhases(); p++) { - if (m_rxnPhaseIsReactant[j][p]) { - if (! m_phaseExists[p]) { - m_ropnet[j] = 0.0; - m_ropf[j] = m_ropr[j]; - if (m_ropf[j] > 0.0) { - for (size_t rp = 0; rp < nPhases(); rp++) { - if (m_rxnPhaseIsProduct[j][rp]) { - if (! m_phaseExists[rp]) { - m_ropnet[j] = 0.0; - m_ropf[j] = m_ropr[j] = 0.0; - } - } + if (m_rxnPhaseIsReactant[j][p] && !m_phaseExists[p]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j]; + if (m_ropf[j] > 0.0) { + for (size_t rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsProduct[j][rp] && !m_phaseExists[rp]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j] = 0.0; } } } } - if (m_rxnPhaseIsProduct[j][p]) { - if (! m_phaseIsStable[p]) { - m_ropnet[j] = 0.0; - m_ropf[j] = m_ropr[j]; - } + if (m_rxnPhaseIsProduct[j][p] && !m_phaseIsStable[p]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j]; } } } diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index 01b6e5e5e..86f1a0516 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -259,10 +259,8 @@ void ReactionPathDiagram::exportToDot(ostream& s) for (i2 = i1+1; i2 < nNodes(); i2++) { k2 = m_speciesNumber[i2]; flx = netFlow(k1, k2); - if (m_local != npos) { - if (k1 != m_local && k2 != m_local) { - flx = 0.0; - } + if (m_local != npos && k1 != m_local && k2 != m_local) { + flx = 0.0; } if (flx != 0.0) { // set beginning and end of the path based on the @@ -817,16 +815,14 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, (m_atoms(kkr,m) < m_elatoms(m, i))) { map >& g = m_transfer[i]; if (g.empty()) { - if (!warn[i]) { - if (!quiet) { - output << endl; - output << "*************** REACTION IGNORED ***************" << endl; - output << "Warning: no rule to determine partitioning of " << element - << endl << " in reaction " << s.reactionString(i) << "." << endl - << "*************** REACTION IGNORED **************" << endl; - output << endl; - warn[i] = 1; - } + if (!warn[i] && !quiet) { + output << endl; + output << "*************** REACTION IGNORED ***************" << endl; + output << "Warning: no rule to determine partitioning of " << element + << endl << " in reaction " << s.reactionString(i) << "." << endl + << "*************** REACTION IGNORED **************" << endl; + output << endl; + warn[i] = 1; } f = 0.0; } else { diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index e3a79c813..d25bba29e 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -152,12 +152,10 @@ bool importKinetics(const XML_Node& phase, std::vector th, string owning_phase = phase["id"]; bool check_for_duplicates = false; - if (phase.parent()) { - if (phase.parent()->hasChild("validate")) { - const XML_Node& d = phase.parent()->child("validate"); - if (d["reactions"] == "yes") { - check_for_duplicates = true; - } + if (phase.parent() && phase.parent()->hasChild("validate")) { + const XML_Node& d = phase.parent()->child("validate"); + if (d["reactions"] == "yes") { + check_for_duplicates = true; } } diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index 2f8f44b91..c3490fc3c 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -357,13 +357,11 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, * End Newton's method. If not converged, print error message and * recalculate sdot's at equal site fractions. */ - if (not_converged) { - if (m_ioflag) { - printf("#$#$#$# Error in solveSP $#$#$#$ \n"); - printf("Newton iter on surface species did not converge, " - "update_norm = %e \n", update_norm); - printf("Continuing anyway\n"); - } + if (not_converged && m_ioflag) { + printf("#$#$#$# Error in solveSP $#$#$#$ \n"); + printf("Newton iter on surface species did not converge, " + "update_norm = %e \n", update_norm); + printf("Continuing anyway\n"); } /* diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index 3199fdaa9..c81100338 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -269,12 +269,10 @@ doublereal IDA_Solver::getCurrentStepFromIDA() void IDA_Solver::setJacobianType(int formJac) { m_formJac = formJac; - if (m_ida_mem) { - if (m_formJac == 1) { - int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian); - if (flag != IDA_SUCCESS) { - throw IDA_Err("IDADlsSetDenseJacFn failed."); - } + if (m_ida_mem && m_formJac == 1) { + int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian); + if (flag != IDA_SUCCESS) { + throw IDA_Err("IDADlsSetDenseJacFn failed."); } } } diff --git a/src/numerics/ResidJacEval.cpp b/src/numerics/ResidJacEval.cpp index 4b64d364e..82ffb617e 100644 --- a/src/numerics/ResidJacEval.cpp +++ b/src/numerics/ResidJacEval.cpp @@ -116,11 +116,9 @@ void ResidJacEval::calcSolnScales(const doublereal t, const doublereal* const ysolnOld, doublereal* const ysolnScales) { - if (ysolnScales) { - if (ysolnScales[0] == 0.0) { - for (int i = 0; i < neq_; i++) { - ysolnScales[i] = 1.0; - } + if (ysolnScales && ysolnScales[0] == 0.0) { + for (int i = 0; i < neq_; i++) { + ysolnScales[i] = 1.0; } } } diff --git a/src/numerics/RootFind.cpp b/src/numerics/RootFind.cpp index 7dae79792..9975cd41b 100644 --- a/src/numerics/RootFind.cpp +++ b/src/numerics/RootFind.cpp @@ -522,10 +522,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun * the user has said that it is ok to take */ doublereal xDelMax = 1.5 * fabs(x2 - x1); - if (specifiedDeltaXnorm_) { - if (0.5 * DeltaXnorm_ > xDelMax) { - xDelMax = 0.5 *DeltaXnorm_ ; - } + if (specifiedDeltaXnorm_ && 0.5 * DeltaXnorm_ > xDelMax) { + xDelMax = 0.5 *DeltaXnorm_ ; } if (fabs(xDelMax) < fabs(xnew - x2)) { xnew = x2 + sign(xnew-x2) * xDelMax; @@ -590,10 +588,8 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun } } } - if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_) { - if (xorig != xnew) { - fprintf(fp, " | xstraddle = %-11.5E", xnew); - } + if (DEBUG_MODE_ENABLED && printLvl >= 3 && writeLogAllowed_ && xorig != xnew) { + fprintf(fp, " | xstraddle = %-11.5E", xnew); } } @@ -601,17 +597,15 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun * Enforce a minimum stepsize if we haven't found a straddle. */ deltaXnew = xnew - x2; - if (fabs(deltaXnew) < 1.2 * delXMeaningful(xnew)) { - if (!foundStraddle) { - sgn = 1.0; - if (x2 > xnew) { - sgn = -1.0; - } - deltaXnew = 1.2 * delXMeaningful(xnew) * sgn; - rfT.reasoning += "Enforcing minimum stepsize from " + fp2str(xnew - x2) + - " to " + fp2str(deltaXnew); - xnew = x2 + deltaXnew; + if (fabs(deltaXnew) < 1.2 * delXMeaningful(xnew) && !foundStraddle) { + sgn = 1.0; + if (x2 > xnew) { + sgn = -1.0; } + deltaXnew = 1.2 * delXMeaningful(xnew) * sgn; + rfT.reasoning += "Enforcing minimum stepsize from " + fp2str(xnew - x2) + + " to " + fp2str(deltaXnew); + xnew = x2 + deltaXnew; } /* @@ -861,42 +855,38 @@ int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& fun /* * Check for excess convergence in the x coordinate */ - if (!converged) { - if (foundStraddle) { - doublereal denom = fabs(x1 - x2); - if (denom < 1.0E-200) { - retn = ROOTFIND_FAILEDCONVERGENCE; - converged = true; - rfT.reasoning += "ConvergenceFZero but X1X2Identical"; - } - if (theSame(x2, x1, 1.0E-2)) { - converged = true; - rfT.reasoning += " ConvergenceF and XSame"; - retn = ROOTFIND_SUCCESS; - } + if (!converged && foundStraddle) { + doublereal denom = fabs(x1 - x2); + if (denom < 1.0E-200) { + retn = ROOTFIND_FAILEDCONVERGENCE; + converged = true; + rfT.reasoning += "ConvergenceFZero but X1X2Identical"; + } + if (theSame(x2, x1, 1.0E-2)) { + converged = true; + rfT.reasoning += " ConvergenceF and XSame"; + retn = ROOTFIND_SUCCESS; } } } else { /* * We are here when F is not converged, but we may want to end anyway */ - if (!converged) { - if (foundStraddle) { - doublereal denom = fabs(x1 - x2); - if (denom < 1.0E-200) { - retn = ROOTFIND_FAILEDCONVERGENCE; - converged = true; - rfT.reasoning += "FNotConverged but X1X2Identical"; - } - /* - * The premise here is that if x1 and x2 get close to one another, - * then the accuracy of the calculation gets destroyed. - */ - if (theSame(x2, x1, 1.0E-5)) { - converged = true; - retn = ROOTFIND_SUCCESS_XCONVERGENCEONLY; - rfT.reasoning += "FNotConverged but XSame"; - } + if (!converged && foundStraddle) { + doublereal denom = fabs(x1 - x2); + if (denom < 1.0E-200) { + retn = ROOTFIND_FAILEDCONVERGENCE; + converged = true; + rfT.reasoning += "FNotConverged but X1X2Identical"; + } + /* + * The premise here is that if x1 and x2 get close to one another, + * then the accuracy of the calculation gets destroyed. + */ + if (theSame(x2, x1, 1.0E-5)) { + converged = true; + retn = ROOTFIND_SUCCESS_XCONVERGENCEONLY; + rfT.reasoning += "FNotConverged but XSame"; } } } @@ -961,11 +951,9 @@ done: std::swap(f1, f2); std::swap(x1, x2); *xbest = x2; - if (fabs(fnew) < fabs(f1)) { - if (f1 * fnew > 0.0) { - std::swap(f1, fnew); - std::swap(x1, xnew); - } + if (fabs(fnew) < fabs(f1) && f1 * fnew > 0.0) { + std::swap(f1, fnew); + std::swap(x1, xnew); } rfT.its = its; diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 4f1f948ef..486e8e6ea 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -53,14 +53,12 @@ doublereal bound_step(const doublereal* x, const doublereal* step, for (j = 0; j < np; j++) { val = x[index(m,j)]; - if (loglevel > 0) { - if (val > above + 1.0e-12 || val < below - 1.0e-12) { - sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n", - int2str(r.domainIndex()).c_str(), - r.componentName(m).c_str(), int2str(j).c_str(), - val, below, above); - writelog(string("\nERROR: solution out of bounds.\n")+buf); - } + if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) { + sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n", + int2str(r.domainIndex()).c_str(), + r.componentName(m).c_str(), int2str(j).c_str(), + val, below, above); + writelog(string("\nERROR: solution out of bounds.\n")+buf); } newval = val + step[index(m,j)]; diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp index b35b3a4ce..f72ad45a1 100644 --- a/src/thermo/DebyeHuckel.cpp +++ b/src/thermo/DebyeHuckel.cpp @@ -932,18 +932,16 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_) /* * Now look at the activity coefficient database */ - if (acNodePtr) { - if (acNodePtr->hasChild("stoichIsMods")) { - XML_Node& sIsNode = acNodePtr->child("stoichIsMods"); - map msIs; - getMap(sIsNode, msIs); - for (map::const_iterator _b = msIs.begin(); - _b != msIs.end(); - ++_b) { - size_t kk = speciesIndex(_b->first); - double val = fpValue(_b->second); - m_speciesCharge_Stoich[kk] = val; - } + if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) { + XML_Node& sIsNode = acNodePtr->child("stoichIsMods"); + map msIs; + getMap(sIsNode, msIs); + for (map::const_iterator _b = msIs.begin(); + _b != msIs.end(); + ++_b) { + size_t kk = speciesIndex(_b->first); + double val = fpValue(_b->second); + m_speciesCharge_Stoich[kk] = val; } } } @@ -978,33 +976,29 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_) for (size_t k = 0; k < m_kk; k++) { std::string kname = speciesName(k); const XML_Node* spPtr = xspecies[k]; - if (spPtr) { - if (spPtr->hasChild("electrolyteSpeciesType")) { - std::string est = getChildValue(*spPtr, "electrolyteSpeciesType"); - if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) { - throw CanteraError("DebyeHuckel:initThermoXML", - "Bad electrolyte type: " + est); - } + if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) { + std::string est = getChildValue(*spPtr, "electrolyteSpeciesType"); + if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) { + throw CanteraError("DebyeHuckel:initThermoXML", + "Bad electrolyte type: " + est); } } } /* * Then look at the phase thermo specification */ - if (acNodePtr) { - if (acNodePtr->hasChild("electrolyteSpeciesType")) { - XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType"); - map msEST; - getMap(ESTNode, msEST); - for (map::const_iterator _b = msEST.begin(); - _b != msEST.end(); - ++_b) { - size_t kk = speciesIndex(_b->first); - std::string est = _b->second; - if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) { - throw CanteraError("DebyeHuckel:initThermoXML", - "Bad electrolyte type: " + est); - } + if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) { + XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType"); + map msEST; + getMap(ESTNode, msEST); + for (map::const_iterator _b = msEST.begin(); + _b != msEST.end(); + ++_b) { + size_t kk = speciesIndex(_b->first); + std::string est = _b->second; + if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) { + throw CanteraError("DebyeHuckel:initThermoXML", + "Bad electrolyte type: " + est); } } } diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 421ac87ba..fc93828d4 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -2004,12 +2004,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (psi_ijk[n] != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), - molality[j]*molality[k]*psi_ijk[n]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && psi_ijk[n] != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); } } } @@ -2020,12 +2018,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all cations if (j != i) { sum2 += molality[j]*(2.0*Phi[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j] * Phi[counterIJ])!= 0.0) { - std::string snj = speciesName(j) + ":"; - printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(), - molality[j]*(2.0*Phi[counterIJ])); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * Phi[counterIJ])!= 0.0) { + std::string snj = speciesName(j) + ":"; + printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(), + molality[j]*(2.0*Phi[counterIJ])); } } for (size_t k = 1; k < m_kk; k++) { @@ -2033,12 +2029,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // two inner sums over anions n = k + j * m_kk + i * m_kk * m_kk; sum2 += molality[j]*molality[k]*psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (psi_ijk[n] != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), - molality[j]*molality[k]*psi_ijk[n]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && psi_ijk[n] != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); } /* * Find the counterIJ for the j,k interaction @@ -2047,12 +2041,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const size_t counterIJ2 = m_CounterIJ[n]; sum4 += (fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), - fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), + fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); } } } @@ -2063,12 +2055,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ if (charge(j) == 0) { sum5 += molality[j]*2.0*m_Lambda_nj(j,i); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { - std::string snj = speciesName(j) + ":"; - printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), - molality[j]*2.0*m_Lambda_nj(j,i)); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { + std::string snj = speciesName(j) + ":"; + printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), + molality[j]*2.0*m_Lambda_nj(j,i)); } /* * Zeta interaction term @@ -2150,12 +2140,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) > 0) { n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (psi_ijk[n] != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), - molality[j]*molality[k]*psi_ijk[n]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && psi_ijk[n] != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); } } } @@ -2169,12 +2157,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // sum over all anions if (j != i) { sum2 += molality[j]*(2.0*Phi[counterIJ]); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j] * Phi[counterIJ])!= 0.0) { - std::string snj = speciesName(j) + ":"; - printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(), - molality[j]*(2.0*Phi[counterIJ])); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j] * Phi[counterIJ])!= 0.0) { + std::string snj = speciesName(j) + ":"; + printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(), + molality[j]*(2.0*Phi[counterIJ])); } } for (size_t k = 1; k < m_kk; k++) { @@ -2182,12 +2168,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const // two inner sums over cations n = k + j * m_kk + i * m_kk * m_kk; sum2 += molality[j]*molality[k]*psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (psi_ijk[n] != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), - molality[j]*molality[k]*psi_ijk[n]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && psi_ijk[n] != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); } /* * Find the counterIJ for the symmetric binary interaction @@ -2196,12 +2180,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const size_t counterIJ2 = m_CounterIJ[n]; sum4 += fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), - fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(), + fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]); } } } @@ -2212,12 +2194,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const */ if (charge(j) == 0.0) { sum5 += molality[j]*2.0*m_Lambda_nj(j,i); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { - std::string snj = speciesName(j) + ":"; - printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), - molality[j]*2.0*m_Lambda_nj(j,i)); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) { + std::string snj = speciesName(j) + ":"; + printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(), + molality[j]*2.0*m_Lambda_nj(j,i)); } /* * Zeta interaction term @@ -2263,12 +2243,10 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const double sum3 = 0.0; for (size_t j = 1; j < m_kk; j++) { sum1 += molality[j]*2.0*m_Lambda_nj(i,j); - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (m_Lambda_nj(i,j) != 0.0) { - std::string snj = speciesName(j) + ":"; - printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(), - molality[j]*2.0*m_Lambda_nj(i,j)); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && m_Lambda_nj(i,j) != 0.0) { + std::string snj = speciesName(j) + ":"; + printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(), + molality[j]*2.0*m_Lambda_nj(i,j)); } /* * Zeta term -> we piggyback on the psi term @@ -2278,23 +2256,19 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const if (charge(k) < 0.0) { size_t n = k + j * m_kk + i * m_kk * m_kk; sum3 += molality[j]*molality[k]*psi_ijk[n]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (psi_ijk[n] != 0.0) { - std::string snj = speciesName(j) + "," + speciesName(k) + ":"; - printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), - molality[j]*molality[k]*psi_ijk[n]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && psi_ijk[n] != 0.0) { + std::string snj = speciesName(j) + "," + speciesName(k) + ":"; + printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(), + molality[j]*molality[k]*psi_ijk[n]); } } } } } double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]; - if (DEBUG_MODE_ENABLED && m_debugCalc) { - if (m_Mu_nnn[i] != 0.0) { - printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n", - 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]); - } + if (DEBUG_MODE_ENABLED && m_debugCalc && m_Mu_nnn[i] != 0.0) { + printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n", + 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]); } m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3; gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]); diff --git a/src/thermo/HMWSoln_input.cpp b/src/thermo/HMWSoln_input.cpp index ab261a380..2da4115c3 100644 --- a/src/thermo/HMWSoln_input.cpp +++ b/src/thermo/HMWSoln_input.cpp @@ -465,11 +465,9 @@ void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt) stemp = xmlChild.value(); double old = m_Theta_ij[counter]; m_Theta_ij[counter] = fpValueCheck(stemp); - if (old != 0.0) { - if (old != m_Theta_ij[counter]) { - throw CanteraError("HMWSoln::readXMLPsiCommonCation", - "conflicting values"); - } + if (old != 0.0 && old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLPsiCommonCation", + "conflicting values"); } } if (nodeName == "psi") { @@ -601,11 +599,9 @@ void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt) stemp = xmlChild.value(); double old = m_Theta_ij[counter]; m_Theta_ij[counter] = fpValueCheck(stemp); - if (old != 0.0) { - if (old != m_Theta_ij[counter]) { - throw CanteraError("HMWSoln::readXMLPsiCommonAnion", - "conflicting values"); - } + if (old != 0.0 && old != m_Theta_ij[counter]) { + throw CanteraError("HMWSoln::readXMLPsiCommonAnion", + "conflicting values"); } } if (nodeName == "psi") { @@ -1414,19 +1410,17 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) /* * Now look at the activity coefficient database */ - if (acNodePtr) { - if (acNodePtr->hasChild("stoichIsMods")) { - XML_Node& sIsNode = acNodePtr->child("stoichIsMods"); - map msIs; - getMap(sIsNode, msIs); - for (map::const_iterator _b = msIs.begin(); - _b != msIs.end(); - ++_b) { - size_t kk = speciesIndex(_b->first); - if (kk != npos) { - double val = fpValue(_b->second); - m_speciesCharge_Stoich[kk] = val; - } + if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) { + XML_Node& sIsNode = acNodePtr->child("stoichIsMods"); + map msIs; + getMap(sIsNode, msIs); + for (map::const_iterator _b = msIs.begin(); + _b != msIs.end(); + ++_b) { + size_t kk = speciesIndex(_b->first); + if (kk != npos) { + double val = fpValue(_b->second); + m_speciesCharge_Stoich[kk] = val; } } } @@ -1496,34 +1490,30 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) std::vector xspecies = speciesData(); for (size_t k = 0; k < m_kk; k++) { const XML_Node* spPtr = xspecies[k]; - if (spPtr) { - if (spPtr->hasChild("electrolyteSpeciesType")) { - string est = getChildValue(*spPtr, "electrolyteSpeciesType"); - if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) { - throw CanteraError("HMWSoln::initThermoXML", - "Bad electrolyte type: " + est); - } + if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) { + string est = getChildValue(*spPtr, "electrolyteSpeciesType"); + if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) { + throw CanteraError("HMWSoln::initThermoXML", + "Bad electrolyte type: " + est); } } } /* * Then look at the phase thermo specification */ - if (acNodePtr) { - if (acNodePtr->hasChild("electrolyteSpeciesType")) { - XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType"); - map msEST; - getMap(ESTNode, msEST); - for (map::const_iterator _b = msEST.begin(); - _b != msEST.end(); - ++_b) { - size_t kk = speciesIndex(_b->first); - if (kk != npos) { - string est = _b->second; - if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) { - throw CanteraError("HMWSoln::initThermoXML", - "Bad electrolyte type: " + est); - } + if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) { + XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType"); + map msEST; + getMap(ESTNode, msEST); + for (map::const_iterator _b = msEST.begin(); + _b != msEST.end(); + ++_b) { + size_t kk = speciesIndex(_b->first); + if (kk != npos) { + string est = _b->second; + if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) { + throw CanteraError("HMWSoln::initThermoXML", + "Bad electrolyte type: " + est); } } } @@ -1585,16 +1575,14 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) } } } - if (notDone) { - if (kMaxC != npos) { - if (mf[kMaxC] > (1.1 * sum / charge(kMaxC))) { - mf[kMaxC] -= sum / charge(kMaxC); - mf[0] += sum / charge(kMaxC); - } else { - mf[kMaxC] *= 0.5; - mf[0] += mf[kMaxC]; - notDone = true; - } + if (notDone && kMaxC != npos) { + if (mf[kMaxC] > (1.1 * sum / charge(kMaxC))) { + mf[kMaxC] -= sum / charge(kMaxC); + mf[0] += sum / charge(kMaxC); + } else { + mf[kMaxC] *= 0.5; + mf[0] += mf[kMaxC]; + notDone = true; } } } diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp index 0d831fbc3..943d8cc86 100644 --- a/src/thermo/IdealMolalSoln.cpp +++ b/src/thermo/IdealMolalSoln.cpp @@ -460,11 +460,9 @@ void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) */ initThermo(); - if (id_.size() > 0) { - if (phaseNode.id() != id_) { - throw CanteraError("IdealMolalSoln::initThermo", - "phasenode and Id are incompatible"); - } + if (id_.size() > 0 && phaseNode.id() != id_) { + throw CanteraError("IdealMolalSoln::initThermo", + "phasenode and Id are incompatible"); } /* diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp index b974ca172..7093c32b4 100644 --- a/src/thermo/IdealSolidSolnPhase.cpp +++ b/src/thermo/IdealSolidSolnPhase.cpp @@ -478,11 +478,9 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_) { - if (id_.size() > 0) { - if (phaseNode.id() != id_) { - throw CanteraError("IdealSolidSolnPhase::initThermoXML", - "phasenode and Id are incompatible"); - } + if (id_.size() > 0 && phaseNode.id() != id_) { + throw CanteraError("IdealSolidSolnPhase::initThermoXML", + "phasenode and Id are incompatible"); } /* diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index 6511ed788..44451d911 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -184,11 +184,9 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id_) { - if (id_.size() > 0) { - if (phaseNode.id() != id_) { - throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", - "phasenode and Id are incompatible"); - } + if (id_.size() > 0 && phaseNode.id() != id_) { + throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML", + "phasenode and Id are incompatible"); } /* @@ -769,16 +767,14 @@ static double factorOverlap(const std::vector& elnamesVN , { double fMax = 1.0E100; for (size_t mi = 0; mi < nElementsI; mi++) { - if (elnamesVI[mi] != "E") { - if (elemVectorI[mi] > 1.0E-13) { - double eiNum = elemVectorI[mi]; - for (size_t mn = 0; mn < nElementsN; mn++) { - if (elnamesVI[mi] == elnamesVN[mn]) { - if (elemVectorN[mn] <= 1.0E-13) { - return 0.0; - } - fMax = std::min(fMax, elemVectorN[mn]/eiNum); + if (elnamesVI[mi] != "E" && elemVectorI[mi] > 1.0E-13) { + double eiNum = elemVectorI[mi]; + for (size_t mn = 0; mn < nElementsN; mn++) { + if (elnamesVI[mi] == elnamesVN[mn]) { + if (elemVectorN[mn] <= 1.0E-13) { + return 0.0; } + fMax = std::min(fMax, elemVectorN[mn]/eiNum); } } } @@ -787,11 +783,9 @@ static double factorOverlap(const std::vector& elnamesVN , } void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) { - if (id_.size() > 0) { - if (phaseNode.id() != id_) { - throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", - "phasenode and Id are incompatible"); - } + if (id_.size() > 0 && phaseNode.id() != id_) { + throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML", + "phasenode and Id are incompatible"); } /* diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index c8da5fb45..bfacf251d 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -327,10 +327,8 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_) throw CanteraError(" LatticePhase::initThermoXML", "database problems"); } XML_Node* ss = s->findByName("standardState"); - if (ss) { - if (ss->findByName("molarVolume")) { - m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI"); - } + if (ss && ss->findByName("molarVolume")) { + m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI"); } } diff --git a/src/thermo/MetalSHEelectrons.cpp b/src/thermo/MetalSHEelectrons.cpp index 535830c7d..e7d11fc99 100644 --- a/src/thermo/MetalSHEelectrons.cpp +++ b/src/thermo/MetalSHEelectrons.cpp @@ -58,11 +58,9 @@ MetalSHEelectrons::MetalSHEelectrons(const std::string& infile, std::string id_) MetalSHEelectrons::MetalSHEelectrons(XML_Node& xmlphase, const std::string& id_) : xdef_(0) { - if (id_ != "") { - if (id_ != xmlphase["id"]) { - throw CanteraError("MetalSHEelectrons::MetalSHEelectrons", - "id's don't match"); - } + if (id_ != "" && id_ != xmlphase["id"]) { + throw CanteraError("MetalSHEelectrons::MetalSHEelectrons", + "id's don't match"); } if (xmlphase.child("thermo")["model"] != "MetalSHEelectrons") { throw CanteraError("MetalSHEelectrons::MetalSHEelectrons", diff --git a/src/thermo/MineralEQ3.cpp b/src/thermo/MineralEQ3.cpp index 80e416628..1a7c871c5 100644 --- a/src/thermo/MineralEQ3.cpp +++ b/src/thermo/MineralEQ3.cpp @@ -49,11 +49,9 @@ MineralEQ3::MineralEQ3(const std::string& infile, std::string id_) MineralEQ3::MineralEQ3(XML_Node& xmlphase, const std::string& id_) { - if (id_ != "") { - if (id_ != xmlphase["id"]) { - throw CanteraError("MineralEQ3::MineralEQ3", - "id's don't match"); - } + if (id_ != "" && id_ != xmlphase["id"]) { + throw CanteraError("MineralEQ3::MineralEQ3", + "id's don't match"); } std::string model = xmlphase.child("thermo")["model"]; if (model != "StoichSubstance" && model != "MineralEQ3") { diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index 4e6e79959..1fa2d66e1 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -580,26 +580,20 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa, * of 0.1 times the current volume */ double delMV = - (presBase - presPa) / dpdV; - if (!gasSide || delMV < 0.0) { - if (fabs(delMV) > 0.2 * molarVolBase) { - delMV = delMV / fabs(delMV) * 0.2 * molarVolBase; - } + if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) { + delMV = delMV / fabs(delMV) * 0.2 * molarVolBase; } /* * Only go 1/10 the way towards the spinodal at any one time. */ if (TKelvin < tcrit) { if (gasSide) { - if (delMV < 0.0) { - if (-delMV > 0.5 * (molarVolBase - molarVolSpinodal)) { - delMV = - 0.5 * (molarVolBase - molarVolSpinodal); - } + if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) { + delMV = - 0.5 * (molarVolBase - molarVolSpinodal); } } else { - if (delMV > 0.0) { - if (delMV > 0.5 * (molarVolSpinodal - molarVolBase)) { - delMV = 0.5 * (molarVolSpinodal - molarVolBase); - } + if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) { + delMV = 0.5 * (molarVolSpinodal - molarVolBase); } } } @@ -828,37 +822,35 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola } } - if (foundGas && foundLiquid) { - if (presGas != presLiquid) { - pres = 0.5 * (presLiquid + presGas); - bool goodLiq; - bool goodGas; - for (int i = 0; i < 50; i++) { - densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood); - if (densLiquid <= 0.0) { - goodLiq = false; - } else { - goodLiq = true; - RhoLiquidGood = densLiquid; - presLiquid = pres; - } - densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood); - if (densGas <= 0.0) { - goodGas = false; - } else { - goodGas = true; - RhoGasGood = densGas; - presGas = pres; - } - if (goodGas && goodLiq) { - break; - } - if (!goodLiq && !goodGas) { - pres = 0.5 * (pres + presLiquid); - } - if (goodLiq || goodGas) { - pres = 0.5 * (presLiquid + presGas); - } + if (foundGas && foundLiquid && presGas != presLiquid) { + pres = 0.5 * (presLiquid + presGas); + bool goodLiq; + bool goodGas; + for (int i = 0; i < 50; i++) { + densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood); + if (densLiquid <= 0.0) { + goodLiq = false; + } else { + goodLiq = true; + RhoLiquidGood = densLiquid; + presLiquid = pres; + } + densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood); + if (densGas <= 0.0) { + goodGas = false; + } else { + goodGas = true; + RhoGasGood = densGas; + presGas = pres; + } + if (goodGas && goodLiq) { + break; + } + if (!goodLiq && !goodGas) { + pres = 0.5 * (pres + presLiquid); + } + if (goodLiq || goodGas) { + pres = 0.5 * (presLiquid + presGas); } } } diff --git a/src/thermo/MolalityVPSSTP.cpp b/src/thermo/MolalityVPSSTP.cpp index f34ceaa7d..f041d18c6 100644 --- a/src/thermo/MolalityVPSSTP.cpp +++ b/src/thermo/MolalityVPSSTP.cpp @@ -203,17 +203,13 @@ void MolalityVPSSTP::setMolalitiesByName(const compositionMap& mMap) for (size_t k = 0; k < m_kk; k++) { double ch = charge(k); if (mf[k] > 0.0) { - if (ch > 0.0) { - if (ch * mf[k] > cPos) { - largePos = k; - cPos = ch * mf[k]; - } + if (ch > 0.0 && ch * mf[k] > cPos) { + largePos = k; + cPos = ch * mf[k]; } - if (ch < 0.0) { - if (fabs(ch) * mf[k] > cNeg) { - largeNeg = k; - cNeg = fabs(ch) * mf[k]; - } + if (ch < 0.0 && fabs(ch) * mf[k] > cNeg) { + largeNeg = k; + cNeg = fabs(ch) * mf[k]; } } sum += mf[k] * ch; diff --git a/src/thermo/Mu0Poly.cpp b/src/thermo/Mu0Poly.cpp index f2af5c3e8..30a1a66f5 100644 --- a/src/thermo/Mu0Poly.cpp +++ b/src/thermo/Mu0Poly.cpp @@ -180,11 +180,9 @@ void Mu0Poly::processCoeffs(const doublereal* coeffs) iT298 = i; ifound = true; } - if (i < nPoints - 1) { - if (coeffs[iindex+2] <= T1) { - throw CanteraError("Mu0Poly", - "Temperatures are not monotonic increasing"); - } + if (i < nPoints - 1 && coeffs[iindex+2] <= T1) { + throw CanteraError("Mu0Poly", + "Temperatures are not monotonic increasing"); } iindex += 2; } diff --git a/src/thermo/PDSS_Water.cpp b/src/thermo/PDSS_Water.cpp index 9cf3aea73..bafd0fc8b 100644 --- a/src/thermo/PDSS_Water.cpp +++ b/src/thermo/PDSS_Water.cpp @@ -316,11 +316,9 @@ void PDSS_Water::setPressure(doublereal p) // We are only putting the phase check here because of speed considerations. m_iState = m_sub.phaseState(true); - if (! m_allowGasPhase) { - if (m_iState != WATER_SUPERCRIT && m_iState != WATER_LIQUID && m_iState != WATER_UNSTABLELIQUID) { - throw CanteraError("PDSS_Water::setPressure", - "Water State isn't liquid or crit"); - } + if (!m_allowGasPhase && m_iState != WATER_SUPERCRIT && m_iState != WATER_LIQUID && m_iState != WATER_UNSTABLELIQUID) { + throw CanteraError("PDSS_Water::setPressure", + "Water State isn't liquid or crit"); } } diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 9a8bc2f60..763b5af23 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -1318,13 +1318,11 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; if (fabs(tmp) > 1.0E-4) { for (int j = 0; j < 3; j++) { - if (j != i) { - if (fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) { - writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + - fp2str(pres) + "): WARNING roots have merged: " + - fp2str(Vroot[i]) + ", " + fp2str(Vroot[j])); - writelogendl(); - } + if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) { + writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + + fp2str(pres) + "): WARNING roots have merged: " + + fp2str(Vroot[i]) + ", " + fp2str(Vroot[j])); + writelogendl(); } } } @@ -1400,10 +1398,8 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do } } } else { - if (nSolnValues == 2) { - if (delta > 0.0) { - nSolnValues = -2; - } + if (nSolnValues == 2 && delta > 0.0) { + nSolnValues = -2; } } return nSolnValues; diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 9bdd036de..da6279f04 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -377,13 +377,11 @@ void importPhase(XML_Node& phase, ThermoPhase* th) * present. ***************************************************************/ vector sparrays = phase.getChildren("speciesArray"); - if (ssConvention != cSS_CONVENTION_SLAVE) { - if (sparrays.empty()) { - throw CanteraError("importPhase", - "phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n" - + " There must be at least one speciesArray nodes " - "with one or more species"); - } + if (ssConvention != cSS_CONVENTION_SLAVE && sparrays.empty()) { + throw CanteraError("importPhase", + "phase, " + th->id() + ", has zero \"speciesArray\" XML nodes.\n" + + " There must be at least one speciesArray nodes " + "with one or more species"); } vector dbases; vector_int sprule(sparrays.size(),0); @@ -450,11 +448,9 @@ void importPhase(XML_Node& phase, ThermoPhase* th) } size_t nsp = spDataNodeList.size(); - if (ssConvention == cSS_CONVENTION_SLAVE) { - if (nsp > 0) { - throw CanteraError("importPhase()", "For Slave standard states, number of species must be zero: " - + int2str(nsp)); - } + if (ssConvention == cSS_CONVENTION_SLAVE && nsp > 0) { + throw CanteraError("importPhase()", "For Slave standard states, number of species must be zero: " + + int2str(nsp)); } for (size_t k = 0; k < nsp; k++) { XML_Node* s = spDataNodeList[k]; diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index 27662e5bd..ee596c8ce 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -625,10 +625,8 @@ void ThermoPhase::setSpeciesThermo(SpeciesThermo* spthermo) "Use of SpeciesThermo classes other than " "GeneralSpeciesThermo is deprecated."); } - if (m_spthermo) { - if (m_spthermo != spthermo) { - delete m_spthermo; - } + if (m_spthermo && m_spthermo != spthermo) { + delete m_spthermo; } m_spthermo = spthermo; } diff --git a/src/thermo/VPSSMgrFactory.cpp b/src/thermo/VPSSMgrFactory.cpp index 2fb641d53..31f8afb36 100644 --- a/src/thermo/VPSSMgrFactory.cpp +++ b/src/thermo/VPSSMgrFactory.cpp @@ -209,17 +209,15 @@ VPSSMgr* VPSSMgrFactory::newVPSSMgr(VPStandardStateTP* vp_ptr, // First look for any explicit instructions within the XML Database // for the standard state manager and the variable pressure // standard state manager - if (phaseNode_ptr) { - if (phaseNode_ptr->hasChild("thermo")) { - const XML_Node& thermoNode = phaseNode_ptr->child("thermo"); - if (thermoNode.hasChild("standardStateManager")) { - const XML_Node& ssNode = thermoNode.child("standardStateManager"); - ssManager = ssNode["model"]; - } - if (thermoNode.hasChild("variablePressureStandardStateManager")) { - const XML_Node& vpssNode = thermoNode.child("variablePressureStandardStateManager"); - vpssManager = vpssNode["model"]; - } + if (phaseNode_ptr && phaseNode_ptr->hasChild("thermo")) { + const XML_Node& thermoNode = phaseNode_ptr->child("thermo"); + if (thermoNode.hasChild("standardStateManager")) { + const XML_Node& ssNode = thermoNode.child("standardStateManager"); + ssManager = ssNode["model"]; + } + if (thermoNode.hasChild("variablePressureStandardStateManager")) { + const XML_Node& vpssNode = thermoNode.child("variablePressureStandardStateManager"); + vpssManager = vpssNode["model"]; } } @@ -270,10 +268,9 @@ VPSSMgr* VPSSMgrFactory::newVPSSMgr(VPStandardStateTP* vp_ptr, } } } - if (inasaCV || ishomateCV || isimpleCV) { - if (!inasaIG && !ishomateIG && !isimpleIG && !itpx && !ihptx && !iother) { - return new VPSSMgr_ConstVol(vp_ptr, spth); - } + if ((inasaCV || ishomateCV || isimpleCV) && + !inasaIG && !ishomateIG && !isimpleIG && !itpx && !ihptx && !iother) { + return new VPSSMgr_ConstVol(vp_ptr, spth); } return new VPSSMgr_General(vp_ptr, spth); } diff --git a/src/thermo/VPSSMgr_IdealGas.cpp b/src/thermo/VPSSMgr_IdealGas.cpp index 44e9defa4..05bec774a 100644 --- a/src/thermo/VPSSMgr_IdealGas.cpp +++ b/src/thermo/VPSSMgr_IdealGas.cpp @@ -88,12 +88,10 @@ PDSS* VPSSMgr_IdealGas::createInstallPDSS(size_t k, const XML_Node& speciesNode, const XML_Node* const phaseNode_ptr) { const XML_Node* ss = speciesNode.findByName("standardState"); - if (ss) { - if (ss->attrib("model") != "ideal_gas") { - throw CanteraError("VPSSMgr_IdealGas::createInstallPDSS", - "standardState model for species isn't " - "ideal_gas: " + speciesNode["name"]); - } + if (ss && ss->attrib("model") != "ideal_gas") { + throw CanteraError("VPSSMgr_IdealGas::createInstallPDSS", + "standardState model for species isn't " + "ideal_gas: " + speciesNode["name"]); } if (m_Vss.size() < k+1) { m_Vss.resize(k+1, 0.0); diff --git a/src/thermo/WaterProps.cpp b/src/thermo/WaterProps.cpp index 56b0b8103..3c01ccfa5 100644 --- a/src/thermo/WaterProps.cpp +++ b/src/thermo/WaterProps.cpp @@ -353,14 +353,12 @@ doublereal WaterProps::viscosityWater() const // Apply the near-critical point corrections if necessary doublereal mu2bar = 1.0; - if ((tbar >= 0.9970) && tbar <= 1.0082) { - if ((rhobar >= 0.755) && (rhobar <= 1.290)) { - doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho(); - drhodp *= presStar / rhoStar; - doublereal xsi = rhobar * drhodp; - if (xsi >= 21.93) { - mu2bar = 0.922 * std::pow(xsi, 0.0263); - } + if (tbar >= 0.9970 && tbar <= 1.0082 && rhobar >= 0.755 && rhobar <= 1.290) { + doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho(); + drhodp *= presStar / rhoStar; + doublereal xsi = rhobar * drhodp; + if (xsi >= 21.93) { + mu2bar = 0.922 * std::pow(xsi, 0.0263); } } diff --git a/src/transport/LTPspecies.cpp b/src/transport/LTPspecies.cpp index a1b3f2b69..a71d078ce 100644 --- a/src/transport/LTPspecies.cpp +++ b/src/transport/LTPspecies.cpp @@ -57,10 +57,8 @@ LTPspecies::LTPspecies(const XML_Node* const propNode, const std::string name, m_thermo(thermo), m_mixWeight(1.0) { - if (propNode) { - if (propNode->hasChild("mixtureWeighting")) { - m_mixWeight = getFloat(*propNode, "mixtureWeighting"); - } + if (propNode && propNode->hasChild("mixtureWeighting")) { + m_mixWeight = getFloat(*propNode, "mixtureWeighting"); } } diff --git a/test_problems/cathermo/HMW_graph_GvI/HMW_graph_GvI.cpp b/test_problems/cathermo/HMW_graph_GvI/HMW_graph_GvI.cpp index bf6a03ed5..cfc2d2020 100644 --- a/test_problems/cathermo/HMW_graph_GvI/HMW_graph_GvI.cpp +++ b/test_problems/cathermo/HMW_graph_GvI/HMW_graph_GvI.cpp @@ -78,12 +78,10 @@ int main(int argc, char** argv) + ISQRTbot*(1.0 - (double)i/(its - 1.0)); Is = ISQRT * ISQRT; - if (!doneSp) { - if (Is > 6.146) { - Is = 6.146; - doneSp = true; - i--; - } + if (!doneSp && Is > 6.146) { + Is = 6.146; + doneSp = true; + i--; } moll[i1] = Is; moll[i2] = Is;