diff --git a/include/cantera/thermo/IonsFromNeutralVPSSTP.h b/include/cantera/thermo/IonsFromNeutralVPSSTP.h index 76e6bcaed..bcc431118 100644 --- a/include/cantera/thermo/IonsFromNeutralVPSSTP.h +++ b/include/cantera/thermo/IonsFromNeutralVPSSTP.h @@ -824,6 +824,7 @@ public: ThermoPhase* neutralMoleculePhase_; private: + GibbsExcessVPSSTP *geThermo; //! If true then we own the underlying neutral Molecule Phase /*! diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index 7147beef0..9cacb95a7 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -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) diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index b9ac0c4e4..d424699c3 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -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(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(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(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(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(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(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(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(neutralMoleculePhase_); if (!geThermo) { throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed"); } diff --git a/src/thermo/MargulesVPSSTP.cpp b/src/thermo/MargulesVPSSTP.cpp index ae2e79d96..d88579939 100644 --- a/src/thermo/MargulesVPSSTP.cpp +++ b/src/thermo/MargulesVPSSTP.cpp @@ -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; } } //==================================================================================================================== diff --git a/src/transport/LiquidTransport.cpp b/src/transport/LiquidTransport.cpp index 437685a1c..718c0fe3b 100644 --- a/src/transport/LiquidTransport.cpp +++ b/src/transport/LiquidTransport.cpp @@ -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 m_thermo; //MargulesVPSSTP * tempMarg = dynamic_cast (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) {