research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Performance improvements for HMWSoln and MaskellSolidSolnPhase.

Browse source 
HMWSoln: avoid recomputing lambdas if ionic strength hasn't changed.

MaskellSolidSolnPhase: Avoid recomputing activity coefficients if
r and T haven't changed.
This commit is contained in:
Victor Brunini 2014-03-12 20:37:04 +00:00
parent ec6b65513b
commit 584abf0d10
5 changed files with 62 additions and 40 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
include/cantera/thermo/HMWSoln.h
View file

@ -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
*

25
include/cantera/thermo/MaskellSolidSolnPhase.h
View file

@ -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<SUP>-3</SUP>.
* 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<doublereal> 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;

24
src/thermo/HMWSoln.cpp
View file

@ -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,

46
src/thermo/MaskellSolidSolnPhase.cpp
View file

@ -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<doublereal> 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<doublereal> 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

4
test/thermo/MaskellSolidSolnPhase_Test.cpp
View file

@ -123,8 +123,8 @@ TEST_F(MaskellSolidSolnPhase_Test, standardConcentrations)
initializeTestPhaseWithXML(valid_file);
ASSERT_TRUE(dynamic_cast<MaskellSolidSolnPhase *>(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)