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

Performance improvements focused on Stefan-Maxwell liquid transport

Browse source 
model for molten salts.

Summary:
- Division is expensive, replace repeated division with multiplication
  by stored 1/x.
- New and delete are expensive, make working vectors in some functions
  static so that they are not repeatedly instantiated and deleted.
- Reorder a few loops to reduce computation and hopefully cache misses.
This commit is contained in:
Victor Brunini 2012-12-13 22:14:36 +00:00
parent 4f41718e5b
commit d4d1093957
5 changed files with 111 additions and 108 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
include/cantera/thermo/IonsFromNeutralVPSSTP.h
View file

@ -824,6 +824,7 @@ public:
ThermoPhase* neutralMoleculePhase_;
private:
GibbsExcessVPSSTP *geThermo;
//! If true then we own the underlying neutral Molecule Phase
/*!

6
src/thermo/GibbsExcessVPSSTP.cpp
View file

@ -184,10 +184,9 @@ void GibbsExcessVPSSTP::setPressure(doublereal p)
void GibbsExcessVPSSTP::calcDensity()
{
doublereal* vbar = NULL;
vbar = new doublereal[m_kk];
static vector_fp vbar(m_kk);
// double *vbar = &m_pp[0];
getPartialMolarVolumes(vbar);
getPartialMolarVolumes(&vbar[0]);
doublereal vtotal = 0.0;
for (size_t i = 0; i < m_kk; i++) {
@ -195,7 +194,6 @@ void GibbsExcessVPSSTP::calcDensity()
}
doublereal dd = meanMolecularWeight() / vtotal;
Phase::setDensity(dd);
delete [] vbar;
}
void GibbsExcessVPSSTP::setState_TP(doublereal t, doublereal p)

30
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -47,6 +47,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(0),
geThermo(0),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
@ -101,6 +102,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile,
IOwnNThermoPhase_ = false;
}
initThermoFile(inputFile, id);
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
}
//====================================================================================================================
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot,
@ -124,6 +126,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot,
IOwnNThermoPhase_ = false;
}
importPhase(*findXMLPhase(&phaseRoot, id), this);
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
}
//====================================================================================================================
@ -144,6 +147,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP& b) :
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(0),
geThermo(0),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
@ -185,6 +189,7 @@ operator=(const IonsFromNeutralVPSSTP& b)
} else {
neutralMoleculePhase_ = b.neutralMoleculePhase_;
}
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
GibbsExcessVPSSTP::operator=(b);
@ -955,7 +960,7 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
{
doublereal fmij;
vector_fp y;
static vector_fp y;
y.resize(numNeutralMoleculeSpecies_,0.0);
doublereal sumy, sumdy;
@ -963,6 +968,7 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
//! Zero the vector we are trying to find.
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
y[k] = 0.0;
dy[k] = 0.0;
}
@ -985,8 +991,9 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
if (jNeut != npos) {
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
AssertTrace(fmij != 0.0);
dy[jNeut] += dx[icat] / fmij;
y[jNeut] += moleFractions_[icat] / fmij;
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
y[jNeut] += moleFractions_[icat] * temp;
}
}
@ -994,8 +1001,9 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
size_t icat = passThroughList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
dy[jNeut] += dx[icat] / fmij;
y[jNeut] += moleFractions_[icat] / fmij;
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
y[jNeut] += moleFractions_[icat] * temp;
}
#ifdef DEBUG_MODE_NOT
//check dy sum to zero
@ -1029,8 +1037,9 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const
sumy += y[k];
sumdy += dy[k];
}
sumy = 1.0 / sumy;
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
dy[k] = dy[k]/sumy - y[k]*sumdy/sumy/sumy;
dy[k] = dy[k] * sumy - y[k]*sumdy*sumy*sumy;
}
break;
@ -1474,7 +1483,6 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
/*
* Get the activity coefficients of the neutral molecules
*/
GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
if (!geThermo) {
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
@ -1483,8 +1491,8 @@ void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const double
}
size_t numNeutMolSpec = geThermo->nSpecies();
vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
vector_fp dX_NeutralMolecule(numNeutMolSpec);
static vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
static vector_fp dX_NeutralMolecule(numNeutMolSpec);
getNeutralMoleculeMoleGrads(DATA_PTR(dXds),DATA_PTR(dX_NeutralMolecule));
@ -1545,7 +1553,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
/*
* Get the activity coefficients of the neutral molecules
*/
GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
if (!geThermo) {
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
return;
@ -1604,7 +1611,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
/*
* Get the activity coefficients of the neutral molecules
*/
GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
if (!geThermo) {
dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
return;
@ -1663,7 +1669,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
/*
* Get the activity coefficients of the neutral molecules
*/
GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
if (!geThermo) {
dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
return;
@ -1726,7 +1731,6 @@ void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
/*
* Get the activity coefficients of the neutral molecules
*/
GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
if (!geThermo) {
throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
}

159
src/thermo/MargulesVPSSTP.cpp
View file

@ -479,25 +479,23 @@ void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
for (size_t iK = 0; iK < m_kk; iK++) {
delAK = 0;
delBK = 0;
for (size_t i = 0; i < numBinaryInteractions_; i++) {
}
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
const doublereal temp1 = g0 + g1 * XB;
const doublereal all = -1.0*XA*XB*temp1 - XA*XB*XB*g1;
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
for (size_t iK = 0; iK < m_kk; iK++) {
vbar[iK] += all;
// vbar[iK] += XA*XB*temp1+((delAK-XA)*XB+XA*(delBK-XB))*temp1+XB*XA*(delBK-XB)*g1;
}
vbar[iA] += XB * temp1;
vbar[iB] += XA * temp1 + XA*XB*g1;
}
}
@ -634,25 +632,26 @@ void MargulesVPSSTP::s_update_lnActCoeff() const
size_t iA, iB, iK, delAK, delBK;
double XA, XB, g0 , g1;
double T = temperature();
double RT = GasConstant*T;
lnActCoeff_Scaled_.assign(m_kk, 0.0);
double invRT = 1.0 / (GasConstant*T);
lnActCoeff_Scaled_.resize(m_kk);
for (iK = 0; iK < m_kk; iK++) {
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
lnActCoeff_Scaled_[iK] = 0.0;
}
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) * invRT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) * invRT;
XA = moleFractions_[iA];
XB = moleFractions_[iB];
const doublereal XAXB = XA * XB;
const doublereal g0g1XB = (g0 + g1 * XB);
const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
for (iK = 0; iK < m_kk; iK++) {
lnActCoeff_Scaled_[iK] += all;
}
lnActCoeff_Scaled_[iA] += XB * g0g1XB;
lnActCoeff_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
}
}
//===================================================================================================================
@ -667,29 +666,35 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dT() const
{
size_t iA, iB, iK, delAK, delBK;
doublereal XA, XB, g0, g1;
doublereal T = temperature();
doublereal RTT = GasConstant*T*T;
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
doublereal invT = 1.0 / temperature();
doublereal invRTT = 1.0 / (GasConstant)*invT*invT;
dlnActCoeffdT_Scaled_.resize(m_kk);
d2lnActCoeffdT2_Scaled_.resize(m_kk);
for (iK = 0; iK < m_kk; iK++) {
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = -m_HE_b_ij[i] / RTT;
g1 = -m_HE_c_ij[i] / RTT;
double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
dlnActCoeffdT_Scaled_[iK] += temp;
d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T;
dlnActCoeffdT_Scaled_[iK] = 0.0;
d2lnActCoeffdT2_Scaled_[iK] = 0.0;
}
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
g0 = -m_HE_b_ij[i] * invRTT;
g1 = -m_HE_c_ij[i] * invRTT;
const doublereal XAXB = XA * XB;
const doublereal g0g1XB = (g0 + g1 * XB);
const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
const doublereal mult = 2.0 * invT;
const doublereal dT2all = mult * all;
for (iK = 0; iK < m_kk; iK++) {
// double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
dlnActCoeffdT_Scaled_[iK] += all;
d2lnActCoeffdT2_Scaled_[iK] -= dT2all;
}
dlnActCoeffdT_Scaled_[iA] += XB * g0g1XB;
dlnActCoeffdT_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
d2lnActCoeffdT2_Scaled_[iA] -= mult * XB * g0g1XB;
d2lnActCoeffdT2_Scaled_[iB] -= mult * XA * g0g1XB + XAXB * g1;
}
}
//====================================================================================================================
@ -736,33 +741,27 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal*
for (iK = 0; iK < m_kk; iK++) {
dlnActCoeffds[iK] = 0.0;
}
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
delAK = 0;
delBK = 0;
if (iA==iK) {
delAK = 1;
} else if (iB==iK) {
delBK = 1;
}
XA = moleFractions_[iA];
XB = moleFractions_[iB];
dXA = dXds[iA];
dXB = dXds[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB
+ dlnActCoeffdT_Scaled_[iK]*dTds;
for (size_t i = 0; i < numBinaryInteractions_; i++) {
iA = m_pSpecies_A_ij[i];
iB = m_pSpecies_B_ij[i];
XA = moleFractions_[iA];
XB = moleFractions_[iB];
dXA = dXds[iA];
dXB = dXds[iB];
g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
const doublereal g02g1XB = g0 + 2*g1*XB;
const doublereal g2XAdXB = 2*g1*XA*dXB;
const doublereal all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB;
for (iK = 0; iK < m_kk; iK++) {
// dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB
// + dlnActCoeffdT_Scaled_[iK]*dTds;
dlnActCoeffds[iK] += all + dlnActCoeffdT_Scaled_[iK]*dTds;
}
dlnActCoeffds[iA] += dXB * g02g1XB;
dlnActCoeffds[iB] += dXA * g02g1XB + g2XAdXB;
}
}
//====================================================================================================================

23
src/transport/LiquidTransport.cpp
View file

@ -1523,7 +1523,7 @@ void LiquidTransport::updateHydrodynamicRadius_T()
void LiquidTransport::update_Grad_lnAC()
{
doublereal grad_T;
vector_fp grad_lnAC(m_nsp), grad_X(m_nsp);
static vector_fp grad_lnAC(m_nsp), grad_X(m_nsp);
// IonsFromNeutralVPSSTP * tempIons = dynamic_cast<IonsFromNeutralVPSSTP *> m_thermo;
//MargulesVPSSTP * tempMarg = dynamic_cast<MargulesVPSSTP *> (tempIons->neutralMoleculePhase_);
@ -1634,8 +1634,8 @@ void LiquidTransport::stefan_maxwell_solve()
* considerations involving species concentrations going to zero.
*
*/
for (size_t i = 0; i < m_nsp; i++) {
for (size_t a = 0; a < m_nDim; a++) {
for (size_t a = 0; a < m_nDim; a++) {
for (size_t i = 0; i < m_nsp; i++) {
m_Grad_mu[a*m_nsp + i] =
m_chargeSpecies[i] * Faraday * m_Grad_V[a]
//+ (m_volume_spec[i] - M[i]/dens_) * m_Grad_P[a]
@ -1648,8 +1648,8 @@ void LiquidTransport::stefan_maxwell_solve()
double mwSolvent = m_thermo->molecularWeight(iSolvent);
double mnaught = mwSolvent/ 1000.;
double lnmnaught = log(mnaught);
for (size_t i = 1; i < m_nsp; i++) {
for (size_t a = 0; a < m_nDim; a++) {
for (size_t a = 0; a < m_nDim; a++) {
for (size_t i = 1; i < m_nsp; i++) {
m_Grad_mu[a*m_nsp + i] -=
m_molefracs[i] * GasConstant * m_Grad_T[a] * lnmnaught;
}
@ -1662,6 +1662,7 @@ void LiquidTransport::stefan_maxwell_solve()
*/
double condSum1;
const doublereal invRT = 1.0 / (GasConstant * T);
switch (m_nDim) {
case 1: /* 1-D approximation */
@ -1685,7 +1686,7 @@ void LiquidTransport::stefan_maxwell_solve()
"Unknown reference velocity provided.");
}
for (size_t i = 1; i < m_nsp; i++) {
m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
m_B(i,0) = m_Grad_mu[i] * invRT;
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {
@ -1748,8 +1749,8 @@ void LiquidTransport::stefan_maxwell_solve()
"Unknown reference velocity provided.");
}
for (size_t i = 1; i < m_nsp; i++) {
m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
m_B(i,1) = m_Grad_mu[m_nsp + i] / (GasConstant * T);
m_B(i,0) = m_Grad_mu[i] * invRT;
m_B(i,1) = m_Grad_mu[m_nsp + i] * invRT;
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {
@ -1792,9 +1793,9 @@ void LiquidTransport::stefan_maxwell_solve()
"Unknown reference velocity provided.");
}
for (size_t i = 1; i < m_nsp; i++) {
m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
m_B(i,1) = m_Grad_mu[m_nsp + i] / (GasConstant * T);
m_B(i,2) = m_Grad_mu[2*m_nsp + i] / (GasConstant * T);
m_B(i,0) = m_Grad_mu[i] * invRT;
m_B(i,1) = m_Grad_mu[m_nsp + i] * invRT;
m_B(i,2) = m_Grad_mu[2*m_nsp + i] * invRT;
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {