From 0d241fd64aa454d04969fc43641cce74539af0c3 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 13 Nov 2016 22:58:15 -0500 Subject: [PATCH] Use toSI instead of explicit constant for unit conversions --- src/thermo/MineralEQ3.cpp | 8 ++--- src/thermo/PDSS_HKFT.cpp | 46 ++++++++++++++--------------- src/thermo/SpeciesThermoFactory.cpp | 14 ++++----- test/kinetics/pdep.cpp | 2 +- 4 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/thermo/MineralEQ3.cpp b/src/thermo/MineralEQ3.cpp index 44b171035..7afda21f8 100644 --- a/src/thermo/MineralEQ3.cpp +++ b/src/thermo/MineralEQ3.cpp @@ -253,15 +253,15 @@ void MineralEQ3::convertDGFormation() } } // Ok, now do the calculation. Convert to joules kmol-1 - doublereal dg = m_deltaG_formation_pr_tr * 4.184 * 1.0E3; + doublereal dg = m_deltaG_formation_pr_tr * toSI("cal/gmol"); //! Store the result into an internal variable. m_Mu0_pr_tr = dg + totalSum; - double Hcalc = m_Mu0_pr_tr + 298.15 * m_Entrop_pr_tr * 4184.0; - double DHjmol = m_deltaH_formation_pr_tr * 4184.0; + double Hcalc = m_Mu0_pr_tr + 298.15 * m_Entrop_pr_tr * toSI("cal/gmol"); + double DHjmol = m_deltaH_formation_pr_tr * toSI("kal/gmol"); // If the discrepancy is greater than 100 cal gmol-1, print an error - if (fabs(Hcalc -DHjmol) > 10.* 1.0E6 * 4.184) { + if (fabs(Hcalc -DHjmol) > 100 * toSI("cal/gmol")) { throw CanteraError("installMinEQ3asShomateThermoFromXML()", "DHjmol is not consistent with G and S: {} vs {}", Hcalc, DHjmol); diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp index b220a35a1..d20a6ec2b 100644 --- a/src/thermo/PDSS_HKFT.cpp +++ b/src/thermo/PDSS_HKFT.cpp @@ -203,7 +203,7 @@ doublereal PDSS_HKFT::enthalpy_mole() const doublereal PDSS_HKFT::enthalpy_mole2() const { - double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * 1.0E3 * 4.184; + double enthTRPR = m_Mu0_tr_pr + 298.15 * m_Entrop_tr_pr * toSI("cal/gmol"); return deltaH() + enthTRPR; } @@ -214,7 +214,7 @@ doublereal PDSS_HKFT::intEnergy_mole() const doublereal PDSS_HKFT::entropy_mole() const { - return m_Entrop_tr_pr * 1.0E3 * 4.184 + deltaS(); + return m_Entrop_tr_pr * toSI("cal/gmol") + deltaS(); } doublereal PDSS_HKFT::gibbs_mole() const @@ -276,7 +276,7 @@ doublereal PDSS_HKFT::cp_mole() const doublereal Cp_calgmol = c1term + c2term + a3term + a4term + yterm + xterm + otterm + rterm; // Convert to Joules / kmol - doublereal Cp = Cp_calgmol * 1.0E3 * 4.184; + doublereal Cp = Cp_calgmol * toSI("cal/gmol"); return Cp; } @@ -321,7 +321,7 @@ doublereal PDSS_HKFT::molarVolume() const doublereal molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm; // Convert to m**3 / kmol from (cal/gmol/Pa) - return molVol_calgmolPascal * 4.184 * 1.0E3; + return molVol_calgmolPascal * toSI("cal/gmol"); } doublereal PDSS_HKFT::density() const @@ -403,23 +403,23 @@ void PDSS_HKFT::initThermo() // Ok, we have mu. Let's check it against the input value // of DH_F to see that we have some internal consistency - doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184); - doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184; + doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol")); + doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol"); // If the discrepancy is greater than 100 cal gmol-1, print // an error and exit. - if (fabs(Hcalc -DHjmol) > 100.* 1.0E3 * 4.184) { + if (fabs(Hcalc -DHjmol) > 100.* toSI("cal/gmol")) { std::string sname = m_tp->speciesName(m_spindex); if (s_InputInconsistencyErrorExit) { throw CanteraError("PDSS_HKFT::initThermo()", "For {}, DHjmol is" " not consistent with G and S: {} vs {} cal gmol-1", - sname, Hcalc/4.184E3, m_deltaH_formation_tr_pr); + sname, Hcalc/toSI("cal/gmol"), m_deltaH_formation_tr_pr); } else { writelog("PDSS_HKFT::initThermo() WARNING: DHjmol for {} is not" " consistent with G and S: calculated {} vs input {} cal gmol-1", - sname, Hcalc/4.184E3, m_deltaH_formation_tr_pr); - writelog(" : continuing with consistent DHjmol = {}", Hcalc/4.184E3); - m_deltaH_formation_tr_pr = Hcalc / (1.0E3 * 4.184); + sname, Hcalc/toSI("cal/gmol"), m_deltaH_formation_tr_pr); + writelog(" : continuing with consistent DHjmol = {}", Hcalc/toSI("cal/gmol")); + m_deltaH_formation_tr_pr = Hcalc / toSI("cal/gmol"); } } doublereal nu = 166027; @@ -571,25 +571,25 @@ void PDSS_HKFT::constructPDSSXML(VPStandardStateTP* tp, size_t spindex, if (hasDHO == 0) { m_charge_j = m_tp->charge(m_spindex); convertDGFormation(); - doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184); - m_deltaH_formation_tr_pr = Hcalc / (1.0E3 * 4.184); + doublereal Hcalc = m_Mu0_tr_pr + 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol")); + m_deltaH_formation_tr_pr = Hcalc / toSI("cal/gmol"); } if (hasDGO == 0) { - doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184; - m_Mu0_tr_pr = DHjmol - 298.15 * (m_Entrop_tr_pr * 1.0E3 * 4.184); - m_deltaG_formation_tr_pr = m_Mu0_tr_pr / (1.0E3 * 4.184); + doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol"); + m_Mu0_tr_pr = DHjmol - 298.15 * (m_Entrop_tr_pr * toSI("cal/gmol")); + m_deltaG_formation_tr_pr = m_Mu0_tr_pr / toSI("cal/gmol"); double tmp = m_Mu0_tr_pr; m_charge_j = m_tp->charge(m_spindex); convertDGFormation(); double totalSum = m_Mu0_tr_pr - tmp; m_Mu0_tr_pr = tmp; - m_deltaG_formation_tr_pr = (m_Mu0_tr_pr - totalSum)/ (1.0E3 * 4.184); + m_deltaG_formation_tr_pr = (m_Mu0_tr_pr - totalSum)/ toSI("cal/gmol"); } if (hasSO == 0) { m_charge_j = m_tp->charge(m_spindex); convertDGFormation(); - doublereal DHjmol = m_deltaH_formation_tr_pr * 1.0E3 * 4.184; - m_Entrop_tr_pr = (DHjmol - m_Mu0_tr_pr) / (298.15 * 1.0E3 * 4.184); + doublereal DHjmol = m_deltaH_formation_tr_pr * toSI("cal/gmol"); + m_Entrop_tr_pr = (DHjmol - m_Mu0_tr_pr) / (298.15 * toSI("cal/gmol")); } } @@ -668,7 +668,7 @@ doublereal PDSS_HKFT::deltaH() const + yterm + yrterm + wterm + wrterm + otterm + otrterm; // Convert to Joules / kmol - return deltaH_calgmol * 1.0E3 * 4.184; + return deltaH_calgmol * toSI("cal/gmol"); } doublereal PDSS_HKFT::deltaG() const @@ -702,7 +702,7 @@ doublereal PDSS_HKFT::deltaG() const doublereal deltaG_calgmol = sterm + c1term + a1term + a2term + c2term + a3term + a4term + wterm + wrterm + yterm; // Convert to Joules / kmol - return deltaG_calgmol * 1.0E3 * 4.184; + return deltaG_calgmol * toSI("cal/gmol"); } doublereal PDSS_HKFT::deltaS() const @@ -743,7 +743,7 @@ doublereal PDSS_HKFT::deltaS() const doublereal deltaS_calgmol = c1term + c2term + a3term + a4term + wterm + wrterm + otterm + otrterm; // Convert to Joules / kmol - return deltaS_calgmol * 1.0E3 * 4.184; + return deltaS_calgmol * toSI("cal/gmol"); } doublereal PDSS_HKFT::ag(const doublereal temp, const int ifunc) const @@ -898,7 +898,7 @@ void PDSS_HKFT::convertDGFormation() totalSum -= m_charge_j * LookupGe("H"); } // Ok, now do the calculation. Convert to joules kmol-1 - doublereal dg = m_deltaG_formation_tr_pr * 4.184 * 1.0E3; + doublereal dg = m_deltaG_formation_tr_pr * toSI("cal/gmol"); //! Store the result into an internal variable. m_Mu0_tr_pr = dg + totalSum; } diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 3674a8c98..207512a3c 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -166,11 +166,11 @@ SpeciesThermoInterpType* newShomateForMineralEQ3(const XML_Node& MinEQ3node) doublereal a = getFloat(MinEQ3node, "a", "toSI") / toSI("cal/gmol/K"); doublereal b = getFloat(MinEQ3node, "b", "toSI") / toSI("cal/gmol/K2"); doublereal c = getFloat(MinEQ3node, "c", "toSI") / toSI("cal-K/gmol"); - doublereal dg = deltaG_formation_pr_tr * 4.184 * 1.0E3; - doublereal DHjmol = deltaH_formation_pr_tr * 1.0E3 * 4.184; - doublereal fac = DHjmol - dg - 298.15 * Entrop_pr_tr * 1.0E3 * 4.184; + doublereal dg = deltaG_formation_pr_tr * toSI("cal/gmol"); + doublereal DHjmol = deltaH_formation_pr_tr * toSI("cal/gmol"); + doublereal fac = DHjmol - dg - 298.15 * Entrop_pr_tr * toSI("cal/gmol"); doublereal Mu0_tr_pr = fac + dg; - doublereal e = Entrop_pr_tr * 1.0E3 * 4.184; + doublereal e = Entrop_pr_tr * toSI("cal/gmol"); doublereal Hcalc = Mu0_tr_pr + 298.15 * e; // Now calculate the shomate polynomials @@ -181,11 +181,11 @@ SpeciesThermoInterpType* newShomateForMineralEQ3(const XML_Node& MinEQ3node) // Cp = As + Bs * t + Cs * t*t + Ds * t*t*t + Es / (t*t) // where // t = temperature(Kelvin) / 1000 - double As = a * 4.184; - double Bs = b * 4.184 * 1000.; + double As = a * toSI("cal"); + double Bs = b * toSI("cal") * 1000.; double Cs = 0.0; double Ds = 0.0; - double Es = c * 4.184 / (1.0E6); + double Es = c * toSI("cal") / (1.0E6); double t = 298.15 / 1000.; double H298smFs = As * t + Bs * t * t / 2.0 - Es / t; diff --git a/test/kinetics/pdep.cpp b/test/kinetics/pdep.cpp index 683831b32..b14940b9c 100644 --- a/test/kinetics/pdep.cpp +++ b/test/kinetics/pdep.cpp @@ -41,7 +41,7 @@ protected: void set_TP(double T, double P) { T_ = T; - RT_ = GasConstant / 4184.0 * T; + RT_ = GasConst_cal_mol_K * T; P_ = P; thermo_->setState_TP(T_, P_); }