From 584abf0d104e40c8fd05738ca3679d96fb898ad5 Mon Sep 17 00:00:00 2001 From: Victor Brunini Date: Wed, 12 Mar 2014 20:37:04 +0000 Subject: [PATCH] Performance improvements for HMWSoln and MaskellSolidSolnPhase. HMWSoln: avoid recomputing lambdas if ionic strength hasn't changed. MaskellSolidSolnPhase: Avoid recomputing activity coefficients if r and T haven't changed. --- include/cantera/thermo/HMWSoln.h | 3 ++ .../cantera/thermo/MaskellSolidSolnPhase.h | 25 ++++++---- src/thermo/HMWSoln.cpp | 24 +++++++--- src/thermo/MaskellSolidSolnPhase.cpp | 46 ++++++++++--------- test/thermo/MaskellSolidSolnPhase_Test.cpp | 4 +- 5 files changed, 62 insertions(+), 40 deletions(-) diff --git a/include/cantera/thermo/HMWSoln.h b/include/cantera/thermo/HMWSoln.h index df300e686..8ca3281b8 100644 --- a/include/cantera/thermo/HMWSoln.h +++ b/include/cantera/thermo/HMWSoln.h @@ -3272,6 +3272,9 @@ private: */ void calc_lambdas(double is) const; + //! Track the last ionic strength lambdas were calculated at to avoid unnecessarily recalculating them + mutable double last_is; + /** * Calculate etheta and etheta_prime * diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h index c59a21ce9..f35219148 100644 --- a/include/cantera/thermo/MaskellSolidSolnPhase.h +++ b/include/cantera/thermo/MaskellSolidSolnPhase.h @@ -65,17 +65,17 @@ public: //! Return the standard concentration for the kth species /*! * The standard concentration \f$ C^0_k \f$ used to normalize the - * generalized concentration. In many cases, this quantity will be the - * same for all species in a phase. However, for this case, we will return - * a distinct concentration for each species. This is the inverse of the - * species molar volume. Units for the standard concentration are kmol - * m-3. + * generalized concentration. * - * @param k Species number: this is a require parameter, - * a change from the ThermoPhase base class, where it was - * an optional parameter. + * @param k Species number: this is an optional parameter, */ - virtual doublereal standardConcentration(size_t k) const; + virtual doublereal standardConcentration(size_t k=0) const { return 1.0; } + + //! Natural logarithm of the standard concentration of the kth species. + /*! + * @param k index of the species (defaults to zero) + */ + virtual doublereal logStandardConc(size_t k=0) const { return 0.0; } //! @name Molar Thermodynamic Properties of the Solution //! @{ @@ -294,7 +294,12 @@ protected: * Function to call through to m_spthermo->update and fill m_h0_RT, * m_cp0_R, m_g0_RT, m_s0_R. */ - void _updateThermo() const; + bool _updateThermo() const; + + //! Vector containing the last computed activity coefficients at T = m_tlast and r = last_r + mutable std::vector last_ac; + //! Last value of r used to update activity coeffs + mutable doublereal last_r; //! Vector containing the species reference enthalpies at T = m_tlast mutable vector_fp m_h0_RT; diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 10a826f4b..c3dff31e8 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -73,6 +73,7 @@ HMWSoln::HMWSoln() : CROP_ln_gamma_o_max(3.0), CROP_ln_gamma_k_min(-5.0), CROP_ln_gamma_k_max(15.0), + last_is(-1.0), m_debugCalc(0) { for (size_t i = 0; i < 17; i++) { @@ -126,6 +127,7 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) : CROP_ln_gamma_o_max(3.0), CROP_ln_gamma_k_min(-5.0), CROP_ln_gamma_k_max(15.0), + last_is(-1.0), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -180,6 +182,7 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) : CROP_ln_gamma_o_max(3.0), CROP_ln_gamma_k_min(-5.0), CROP_ln_gamma_k_max(15.0), + last_is(-1.0), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -234,6 +237,7 @@ HMWSoln::HMWSoln(const HMWSoln& b) : CROP_ln_gamma_o_max(3.0), CROP_ln_gamma_k_min(-5.0), CROP_ln_gamma_k_max(15.0), + last_is(-1.0), m_debugCalc(0) { /* @@ -447,6 +451,7 @@ HMWSoln::HMWSoln(int testProb) : CROP_ln_gamma_o_max(3.0), CROP_ln_gamma_k_min(-5.0), CROP_ln_gamma_k_max(15.0), + last_is(-1.0), m_debugCalc(0) { if (testProb != 1) { @@ -1845,7 +1850,10 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const } } - + const double twoT = 2.0 * T; + const double invT = 1.0 / T; + const double invT2 = invT * invT; + const double twoinvT3 = 2.0 * invT * invT2; for (i = 1; i < m_kk; i++) { for (j = 1; j < m_kk; j++) { for (size_t k = 1; k < m_kk; k++) { @@ -1868,14 +1876,14 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const + Psi_coeff[4]*tln; m_Psi_ijk_L[n] = Psi_coeff[1] - + Psi_coeff[2]*2.0*T - - Psi_coeff[3]/(T*T) - + Psi_coeff[4]/T; + + Psi_coeff[2]*twoT + - Psi_coeff[3]*invT2 + + Psi_coeff[4]*invT; m_Psi_ijk_LL[n] = Psi_coeff[2]*2.0 - + 2.0*Psi_coeff[3]/(T*T*T) - - Psi_coeff[4]/(T*T); + + Psi_coeff[3]*twoinvT3 + - Psi_coeff[4]*invT2; } } } @@ -5456,6 +5464,9 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const void HMWSoln::calc_lambdas(double is) const { + const double tol = 1.e-12; + if( std::abs(is - last_is) < tol ) return; + double aphi, dj, jfunc, jprime, t, x, zprod; int i, ij, j; /* @@ -5513,6 +5524,7 @@ void HMWSoln::calc_lambdas(double is) const #endif } } + last_is = is; } void HMWSoln::calc_thetas(int z1, int z2, diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp index 3aee9aba0..67b71db2c 100644 --- a/src/thermo/MaskellSolidSolnPhase.cpp +++ b/src/thermo/MaskellSolidSolnPhase.cpp @@ -23,6 +23,8 @@ namespace Cantera MaskellSolidSolnPhase::MaskellSolidSolnPhase() : m_Pref(OneAtm), m_Pcurrent(OneAtm), + last_ac(2), + last_r(-1.0), m_h0_RT(2), m_cp0_R(2), m_g0_RT(2), @@ -36,6 +38,8 @@ MaskellSolidSolnPhase::MaskellSolidSolnPhase() : MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) : m_Pref(OneAtm), m_Pcurrent(OneAtm), + last_ac(2), + last_r(-1.0), m_h0_RT(2), m_cp0_R(2), m_g0_RT(2), @@ -64,23 +68,13 @@ ThermoPhase* MaskellSolidSolnPhase::duplMyselfAsThermoPhase() const void MaskellSolidSolnPhase:: getActivityConcentrations(doublereal* c) const { - std::vector pmv(m_kk); - getPartialMolarVolumes(&pmv[0]); getActivityCoefficients(c); for(unsigned sp=0; sp < m_kk; ++sp) { - c[sp] *= moleFraction(sp) / pmv[sp]; + c[sp] *= moleFraction(sp); } } -doublereal MaskellSolidSolnPhase::standardConcentration(size_t k) const -{ - std::vector pmv(m_kk); - getPartialMolarVolumes(&pmv[0]); - doublereal result = 1.0 / pmv[k]; - return result; -} - /******************************************************************** * Molar Thermodynamic Properties of the Solution ********************************************************************/ @@ -163,17 +157,23 @@ void MaskellSolidSolnPhase::setMolarDensity(const doublereal n) void MaskellSolidSolnPhase:: getActivityCoefficients(doublereal* ac) const { - _updateThermo(); + bool temp_changed = _updateThermo(); const doublereal r = moleFraction(product_species_index); - const doublereal pval = p(r); - const doublereal fmval = fm(r); - const doublereal rfm = r * fmval; - const doublereal RT = GasConstant * temperature(); - const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval)) / - (std::pow(1 - r - rfm, 1 + pval) * (1 - r)); - const doublereal B = pval * h_mixing / RT; - ac[product_species_index] = A * std::exp(B); - ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B); + const doublereal tol = 1.e-12; + if( temp_changed || std::abs(r - last_r) > tol ) + { + const doublereal pval = p(r); + const doublereal fmval = fm(r); + const doublereal rfm = r * fmval; + const doublereal RT = GasConstant * temperature(); + const doublereal A = (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval)) / + (std::pow(1 - r - rfm, 1 + pval) * (1 - r)); + const doublereal B = pval * h_mixing / RT; + last_ac[product_species_index] = A * std::exp(B); + last_ac[reactant_species_index] = 1 / (A * r * (1-r) ) * std::exp(-B); + last_r = r; + } + std::copy(last_ac.begin(), last_ac.end(), ac); } void MaskellSolidSolnPhase:: @@ -324,7 +324,7 @@ void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string VPStandardStateTP::initThermoXML(phaseNode, id_); } -void MaskellSolidSolnPhase::_updateThermo() const +bool MaskellSolidSolnPhase::_updateThermo() const { assert(m_kk == 2); doublereal tnow = temperature(); @@ -339,7 +339,9 @@ void MaskellSolidSolnPhase::_updateThermo() const m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; } m_tlast = tnow; + return true; } + return false; } doublereal MaskellSolidSolnPhase::s() const diff --git a/test/thermo/MaskellSolidSolnPhase_Test.cpp b/test/thermo/MaskellSolidSolnPhase_Test.cpp index be2cd2684..59bd3dbcc 100644 --- a/test/thermo/MaskellSolidSolnPhase_Test.cpp +++ b/test/thermo/MaskellSolidSolnPhase_Test.cpp @@ -123,8 +123,8 @@ TEST_F(MaskellSolidSolnPhase_Test, standardConcentrations) initializeTestPhaseWithXML(valid_file); ASSERT_TRUE(dynamic_cast(test_phase) != NULL); - EXPECT_DOUBLE_EQ(1.0 / 0.005, test_phase->standardConcentration(0)); - EXPECT_DOUBLE_EQ(1.0 / 0.01, test_phase->standardConcentration(1)); + EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(0)); + EXPECT_DOUBLE_EQ(1.0, test_phase->standardConcentration(1)); } TEST_F(MaskellSolidSolnPhase_Test, activityConcentrations)