diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp index e2ed1351a..f7c9a79f6 100644 --- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -637,7 +637,7 @@ namespace Cantera { dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k]; } } - + //==================================================================================================================== // Get the array of log concentration-like derivatives of the // log activity coefficients /* @@ -665,19 +665,28 @@ namespace Cantera { dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k]; } } - - // This is temporary. We will get rid of this + //==================================================================================================================== + void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) { + s_update_lnActCoeff(); + s_update_dlnActCoeff_dlnN(); + double *data = & dlnActCoeffdlnN_(0,0); + for (int k = 0; k < m_kk; k++) { + for (int m = 0; m < m_kk; m++) { + dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m]; + } + } + } + //==================================================================================================================== void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp) { double p = pressure(); IonsFromNeutralVPSSTP::setState_TP(temp, p); } - - // This is temporary. We will get rid of this + //==================================================================================================================== void IonsFromNeutralVPSSTP::setPressure(doublereal p) { double t = temperature(); IonsFromNeutralVPSSTP::setState_TP(t, p); } - + //==================================================================================================================== // Set the temperature (K) and pressure (Pa) /* * Setting the pressure may involve the solution of a nonlinear equation. @@ -741,7 +750,7 @@ namespace Cantera { } } - + //==================================================================================================================== // Calculate neutral molecule mole fractions /* * This routine calculates the neutral molecule mole @@ -869,8 +878,8 @@ namespace Cantera { } } - -// Calculate neutral molecule mole fractions + //==================================================================================================================== + // Calculate neutral molecule mole fractions /* * This routine calculates the neutral molecule mole * fraction given the vector of ion mole fractions, @@ -1211,6 +1220,7 @@ namespace Cantera { dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_); + dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0); } //==================================================================================================================== //! Return the factor overlap @@ -1734,11 +1744,9 @@ namespace Cantera { /* * Get the activity coefficients of the neutral molecules */ - GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); - + GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); if (!geThermo) { - - return; + throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed"); } int nsp_ge = geThermo->nSpecies(); geThermo->getdlnActCoeffdlnN(nsp_ge, &(dlnActCoeffdlnN_NeutralMolecule_(0,0))); @@ -1750,13 +1758,7 @@ namespace Cantera { // Do the cation list for (k = 0; k < (int) cationList_.size(); k++) { - for (m = 0; m < (int) cationList_.size(); m++) { - //! Get the id for the next cation - //icat = cationList_[k]; - //jNeut = fm_invert_ionForNeutral[icat]; - //fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk]; - //lnActCoeff_Scaled_[icat] = log(gammaNeutralMolecule_[jNeut])/fmij; - + for (m = 0; m < (int) cationList_.size(); m++) { kcat = cationList_[k]; kNeut = fm_invert_ionForNeutral[kcat]; @@ -1767,14 +1769,13 @@ namespace Cantera { mNeut = fm_invert_ionForNeutral[mcat]; mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk]; - dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij; + dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij; - - for (m = 0; m < numPassThroughSpecies_; m++) { - mcat = passThroughList_[m]; - mNeut = fm_invert_ionForNeutral[mcat]; - dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij; - } + } + for (m = 0; m < numPassThroughSpecies_; m++) { + mcat = passThroughList_[m]; + mNeut = fm_invert_ionForNeutral[mcat]; + dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij; } } @@ -1819,8 +1820,6 @@ namespace Cantera { throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type"); break; } - - } //==================================================================================================================== } diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h index 787b305df..d0c62c20a 100644 --- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.h +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.h @@ -458,6 +458,27 @@ namespace Cantera { */ virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const; + //! Get the array of derivatives of the ln activity coefficients with respect to the ln species mole numbers + /*! + * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a + * log of a species mole number (with all other species mole numbers held constant) + * + * units = 1 / kmol + * + * dlnActCoeffdlnN[ ld * k + m] will contain the derivative of log act_coeff for the mth + * species with respect to the number of moles of the kth species. + * + * \f[ + * \frac{d \ln(\gamma_m) }{d \ln( n_k ) }\Bigg|_{n_i} + * \f] + * + * @param ld Number of rows in the matrix + * @param dlnActCoeffdlnN Output vector of derivatives of the + * log Activity Coefficients. length = m_kk * m_kk + */ + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) ; + + //! Get the Salt Dissociation Coefficients //! Returns the vector of dissociation coefficients and vector of charges /*! diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index b80b2aa4c..93adb082c 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -829,7 +829,7 @@ namespace Cantera { virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const; - //! Get the array of derivatives of the log activity coefficients with respect to the ln species mole numbers + //! Get the array of derivatives of the ln activity coefficients with respect to the ln species mole numbers /*! * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a * log of a species mole number (with all other species mole numbers held constant) diff --git a/Cantera/src/thermo/MixtureFugacityTP.cpp b/Cantera/src/thermo/MixtureFugacityTP.cpp index 21e032fa6..901daa9c9 100644 --- a/Cantera/src/thermo/MixtureFugacityTP.cpp +++ b/Cantera/src/thermo/MixtureFugacityTP.cpp @@ -1008,7 +1008,7 @@ namespace Cantera { * WATER_SUPERCRIT above the critical temperature */ int MixtureFugacityTP::phaseState(bool checkState) const { - int state; + int state = iState_; if (checkState) { double t = temperature(); double tcrit = critTemperature(); diff --git a/Cantera/src/thermo/MixtureFugacityTP.h b/Cantera/src/thermo/MixtureFugacityTP.h index f5d603334..5aa9f7c3a 100644 --- a/Cantera/src/thermo/MixtureFugacityTP.h +++ b/Cantera/src/thermo/MixtureFugacityTP.h @@ -753,7 +753,7 @@ protected: * @return returns the estimated saturation pressure at the given temperature */ virtual doublereal psatEst(doublereal TKelvin) const; - + public: //! Estimate for the molar volume of the liquid /*! * Note: this is only used as a starting guess for later routines that actually calculate an @@ -770,7 +770,7 @@ protected: */ virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const; - protected: + public: //! Calculates the density given the temperature and the pressure and a guess at the density. /*! * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this. @@ -797,7 +797,7 @@ protected: */ virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested, doublereal rhoguess); - + protected: //! Utility routine in the calculation of the saturation pressure /*! * Private routine diff --git a/Cantera/src/thermo/RedlichKwongMFTP.cpp b/Cantera/src/thermo/RedlichKwongMFTP.cpp index ec1627ba7..4b225585a 100644 --- a/Cantera/src/thermo/RedlichKwongMFTP.cpp +++ b/Cantera/src/thermo/RedlichKwongMFTP.cpp @@ -342,7 +342,7 @@ namespace Cantera { double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current)); - if (fabs(pp -m_Pcurrent) > 1.0E-5 * m_Pcurrent) { + if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) { throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe"); } #endif @@ -1306,7 +1306,10 @@ namespace Cantera { */ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess) { - //setTemperature(TKelvin); + /* + * It's necessary to set the temperature so that m_a_current is set correctly. + */ + setTemperature(TKelvin); double tcrit = critTemperature(); doublereal mmw = meanMolecularWeight(); double densBase = 0.0; @@ -1331,7 +1334,7 @@ namespace Cantera { } } - + doublereal volguess = mmw / rhoguess; NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_); diff --git a/Cantera/src/thermo/RedlichKwongMFTP.h b/Cantera/src/thermo/RedlichKwongMFTP.h index 8728fad52..979957d96 100644 --- a/Cantera/src/thermo/RedlichKwongMFTP.h +++ b/Cantera/src/thermo/RedlichKwongMFTP.h @@ -185,7 +185,7 @@ namespace Cantera { */ virtual doublereal isothermalCompressibility() const; - protected: + protected: /** * Calculate the density of the mixture using the partial * molar volumes and mole fractions as input @@ -211,6 +211,7 @@ namespace Cantera { */ virtual void calcDensity(); + protected: //! Set the temperature (K) /*! * Overwritten setTemperature(double) from State.h. This @@ -608,7 +609,7 @@ namespace Cantera { * @return Returns the change in enthalpy in units of J kmol-1. */ virtual doublereal hresid() const; - + public: //! Estimate for the molar volume of the liquid /*! * Note: this is only used as a starting guess for later routines that actually calculate an @@ -624,7 +625,7 @@ namespace Cantera { */ virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const; - protected: + public: //! Calculates the density given the temperature and the pressure and a guess at the density. /*! * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this.