From b6cf6b313996c5c6f8308e3b7dc2945c8833d73e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 3 Aug 2010 21:07:28 +0000 Subject: [PATCH] Fleshed out the gradient of actcoeff wrt mole number. --- Cantera/src/base/Array.h | 11 +++ Cantera/src/thermo/GibbsExcessVPSSTP.cpp | 22 +++++- Cantera/src/thermo/GibbsExcessVPSSTP.h | 19 +++-- Cantera/src/thermo/MargulesVPSSTP.cpp | 96 +++++++++++++++++++++--- Cantera/src/thermo/MargulesVPSSTP.h | 17 ++++- 5 files changed, 144 insertions(+), 21 deletions(-) diff --git a/Cantera/src/base/Array.h b/Cantera/src/base/Array.h index bde5f28c6..a0c637db5 100644 --- a/Cantera/src/base/Array.h +++ b/Cantera/src/base/Array.h @@ -223,6 +223,17 @@ namespace Cantera { for (; b != end(); ++b, ++xb, ++yb) *b = a*(*xb) + *yb; } + //! Set all of the entries to zero + inline void zero() { + int nn = m_nrows * m_ncols; + if (nn > 0) { + /* + * Using memset is the fastest way to zero a contiguous + * section of memory. + */ + (void) memset((void *) &m_data[0], 0, nn * sizeof(doublereal)); + } + } //! Allows setting elements using the syntax A(i,j) = x. /*! diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index fe665941a..a2cc35383 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -33,7 +33,15 @@ namespace Cantera { * */ GibbsExcessVPSSTP::GibbsExcessVPSSTP() : - VPStandardStateTP() + VPStandardStateTP(), + moleFractions_(0), + lnActCoeff_Scaled_(0), + dlnActCoeffdT_Scaled_(0), + d2lnActCoeffdT2_Scaled_(0), + dlnActCoeffdlnN_Scaled_(0), + dlnActCoeffdlnX_Scaled_(0), + dlnActCoeffdN_Scaled_(0,0), + m_pp(0) { } @@ -44,7 +52,15 @@ namespace Cantera { * has a working copy constructor */ GibbsExcessVPSSTP::GibbsExcessVPSSTP(const GibbsExcessVPSSTP &b) : - VPStandardStateTP() + VPStandardStateTP(), + moleFractions_(0), + lnActCoeff_Scaled_(0), + dlnActCoeffdT_Scaled_(0), + d2lnActCoeffdT2_Scaled_(0), + dlnActCoeffdlnN_Scaled_(0), + dlnActCoeffdlnX_Scaled_(0), + dlnActCoeffdN_Scaled_(0,0), + m_pp(0) { GibbsExcessVPSSTP::operator=(b); } @@ -69,6 +85,7 @@ namespace Cantera { d2lnActCoeffdT2_Scaled_ = b.d2lnActCoeffdT2_Scaled_; dlnActCoeffdlnX_Scaled_ = b.dlnActCoeffdlnX_Scaled_; dlnActCoeffdlnN_Scaled_ = b.dlnActCoeffdlnN_Scaled_; + dlnActCoeffdN_Scaled_ = b.dlnActCoeffdN_Scaled_; m_pp = b.m_pp; return *this; @@ -328,6 +345,7 @@ namespace Cantera { d2lnActCoeffdT2_Scaled_.resize(m_kk); dlnActCoeffdlnX_Scaled_.resize(m_kk); dlnActCoeffdlnN_Scaled_.resize(m_kk); + dlnActCoeffdN_Scaled_.resize(m_kk, m_kk); m_pp.resize(m_kk); } diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index 832c75c45..da716c5fa 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -542,6 +542,8 @@ namespace Cantera { protected: + // HKM get rid of _Scaled_ prefix + //! Storage for the current values of the mole fractions of the species /*! * This vector is kept up-to-date when the setState functions are called. @@ -553,29 +555,36 @@ namespace Cantera { mutable std::vector moleFractions_; //! Storage for the current values of the activity coefficients of the - //! species, divided by RT + //! species mutable std::vector lnActCoeff_Scaled_; //! Storage for the current derivative values of the //! gradients with respect to temperature of the - //! log of theactivity coefficients of the species + //! log of the activity coefficients of the species mutable std::vector dlnActCoeffdT_Scaled_; //! Storage for the current derivative values of the //! gradients with respect to temperature of the - //! log of theactivity coefficients of the species + //! log of the activity coefficients of the species mutable std::vector d2lnActCoeffdT2_Scaled_; //! Storage for the current derivative values of the //! gradients with respect to logarithm of the mole fraction of the - //! log of theactivity coefficients of the species + //! log of the activity coefficients of the species @deprecated mutable std::vector dlnActCoeffdlnN_Scaled_; //! Storage for the current derivative values of the //! gradients with respect to logarithm of the mole fraction of the - //! log of theactivity coefficients of the species + //! log of theactivity coefficients of the species @deprecated mutable std::vector dlnActCoeffdlnX_Scaled_; + //! Storage for the current derivative values of the gradients with respect to logarithm of the species mole number of the + //! log of the activity coefficients of the species + /*! + * dlnActCoeffdN_Scaled_(k, m) is the derivative of ln(gamma_k) wrt ln mole number of species m + */ + mutable Array2D dlnActCoeffdN_Scaled_; + //! Temporary storage space that is fair game mutable std::vector m_pp; diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index 8664a1ba9..140f71b1f 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -375,7 +375,7 @@ namespace Cantera { //return -1.0; return 0.0; } - + //==================================================================================================================== // Get the array of non-dimensional molar-based activity coefficients at // the current solution temperature, pressure, and solution concentration. /* @@ -667,7 +667,7 @@ namespace Cantera { // been identified. void MargulesVPSSTP::initLengths() { m_kk = nSpecies(); - + dlnActCoeffdN_Scaled_.resize(m_kk, m_kk); } /* @@ -759,7 +759,7 @@ namespace Cantera { double T = temperature(); double RT = GasConstant*T; fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); - for ( iK = 0; iK < m_kk; iK++ ){ + for (iK = 0; iK < m_kk; iK++) { XK = moleFractions_[iK]; for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -786,13 +786,12 @@ namespace Cantera { */ void MargulesVPSSTP::s_update_dlnActCoeff_dT() const { int iA, iB, iK, delAK, delBK; - doublereal XA, XB, XK, g0, g1; + doublereal XA, XB, g0, g1; doublereal T = temperature(); doublereal RTT = GasConstant*T*T; fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); - for ( iK = 0; iK < m_kk; iK++ ){ - XK = moleFractions_[iK]; + for (iK = 0; iK < m_kk; iK++) { for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; iB = m_pSpecies_B_ij[i]; @@ -923,6 +922,73 @@ namespace Cantera { } } + //==================================================================================================================== + // Update the derivative of the log of the activity coefficients wrt dlnN + /* + * This function will be called to update the internally stored gradients of the + * logarithm of the activity coefficients. These are used in the determination + * of the diffusion coefficients. + * + * he = X_A X_B(B + C X_B) + */ + void MargulesVPSSTP::s_update_dlnActCoeff_dN() const { + int iA, iB; + doublereal delAK, delBK; + double XA, XB, g0 , g1; + double T = temperature(); + double RT = GasConstant*T; + + doublereal delAM, delBM; + + dlnActCoeffdN_Scaled_.zero(); + + /* + * Loop over the activity coefficient gamma_k + */ + for (int iK = 0; iK < m_kk; iK++) { + for (int iM = 0; iM < m_kk; iM++) { + + for (int i = 0; i < numBinaryInteractions_; i++) { + + iA = m_pSpecies_A_ij[i]; + iB = m_pSpecies_B_ij[i]; + + delAK = 0.0; + delBK = 0.0; + delAM = 0.0; + delBM = 0.0; + if (iA==iK) delAK = 1.0; + else if (iB==iK) delBK = 1.0; + if (iA==iM) delAM = 1.0; + else if (iB==iM) delBM = 1.0; + + 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; + + double gfac = g0 + g1 * XB; + double gggg = (delBK - XB) * g1; + + // all values of dlnActCoeffdN_Scaled_(iK, iM) hare an additional divisor of n_total + + dlnActCoeffdN_Scaled_(iK, iM) += gfac * delAK * ( - XB + delBM); + + dlnActCoeffdN_Scaled_(iK, iM) += gfac * delBK * ( - XA + delAM); + + dlnActCoeffdN_Scaled_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM); + + dlnActCoeffdN_Scaled_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM); + + dlnActCoeffdN_Scaled_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM); + + dlnActCoeffdN_Scaled_(iK, iM) += - g1 * XA * XB * (- XB + delBM); + } + } + } + } + //==================================================================================================================== void MargulesVPSSTP::s_update_dlnActCoeff_dlnX() const { int iA, iB; @@ -950,21 +1016,31 @@ namespace Cantera { } } - + //==================================================================================================================== void MargulesVPSSTP::getdlnActCoeffdlnN(doublereal *dlnActCoeffdlnN) const { s_update_dlnActCoeff_dlnN(); for (int k = 0; k < m_kk; k++) { dlnActCoeffdlnN[k] = dlnActCoeffdlnN_Scaled_[k]; } } + //==================================================================================================================== void MargulesVPSSTP::getdlnActCoeffdlnX(doublereal *dlnActCoeffdlnX) const { s_update_dlnActCoeff_dlnX(); for (int k = 0; k < m_kk; k++) { dlnActCoeffdlnX[k] = dlnActCoeffdlnX_Scaled_[k]; } } - - + //==================================================================================================================== + void MargulesVPSSTP::getdlnActCoeffdN(const int ld, doublereal *dlnActCoeffdN) const { + s_update_dlnActCoeff_dN(); + double *data = & dlnActCoeffdN_Scaled_(0,0); + for (int k = 0; k < m_kk; k++) { + for (int m = 0; m < m_kk; m++) { + dlnActCoeffdN[ld * k + m] = data[m_kk * k + m]; + } + } + } + //==================================================================================================================== void MargulesVPSSTP::resizeNumInteractions(const int num) { numBinaryInteractions_ = num; m_HE_b_ij.resize(num, 0.0); @@ -984,7 +1060,7 @@ namespace Cantera { m_pSpecies_B_ij.resize(num, -1); } - + //==================================================================================================================== /* * Process an XML node called "binaryNeutralSpeciesParameters" diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index 9ebfd9aa1..b35719d33 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -113,7 +113,7 @@ namespace Cantera { * \f[ * R T \ln( \gamma_k )= \sum_i \left( \left( \delta_{Ai,k} X_{Bi} + \delta_{Bi,k} X_{Ai} - X_{Ai} X_{Bi} \right) * \left( g^E_{o,i} + g^E_{1,i} X_{Bi} \right) + - * \left( \delta_{Ai,k} - X_{Bi} \right) X_{Ai} X_{Bi} g^E_{1,i} \right) + * \left( \delta_{Bi,k} - X_{Bi} \right) X_{Ai} X_{Bi} g^E_{1,i} \right) * \f] * where * \f$ g^E_{o,i} = h_{o,i} - T s_{o,i} \f$ and \f$ g^E_{1,i} = h_{1,i} - T s_{1,i} \f$ @@ -834,9 +834,8 @@ namespace Cantera { * @param dlnActCoeffdN Output vector of derivatives of the * log Activity Coefficients. length = m_kk * m_kk */ - virtual void getdlnActCoeffdN(const int ld, doublereal * const dlnActCoeffdN) const { - err("getdlnActCoeffdN"); - } + virtual void getdlnActCoeffdN(const int ld, doublereal * const dlnActCoeffdN) const; + //@} private: @@ -898,6 +897,14 @@ namespace Cantera { */ void s_update_dlnActCoeff_dlnN() const; + //! Update the derivative of the log of the activity coefficients wrt log(moles_m) + /*! + * This function will be called to update the internally storred + * derivative of the natural logarithm of the activity coefficients + * wrt logarithm of the mole number of species + */ + void s_update_dlnActCoeff_dN() const; + private: //! Error function @@ -961,6 +968,8 @@ namespace Cantera { //! Entropy term for the quaternary mole fraction interaction of the //! excess gibbs free energy expression mutable vector_fp m_VSE_d_ij; + + //! vector of species indices representing species A in the interaction /*!