From 59a800218f5b52cebf6ca5ae4b548737adf100be Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 14 Apr 2011 18:24:13 +0000 Subject: [PATCH] Changed getdlnActCoeffdlnN to be nonconst. added getdlnActCoeffdlnN_numderiv() function --- Cantera/src/thermo/GibbsExcessVPSSTP.h | 2 +- Cantera/src/thermo/MargulesVPSSTP.cpp | 2 +- Cantera/src/thermo/MargulesVPSSTP.h | 2 +- Cantera/src/thermo/MolalityVPSSTP.h | 4 +- Cantera/src/thermo/PhaseCombo_Interaction.cpp | 3 +- Cantera/src/thermo/PhaseCombo_Interaction.h | 2 +- Cantera/src/thermo/ThermoPhase.cpp | 85 ++++++++++++++++++- Cantera/src/thermo/ThermoPhase.h | 4 +- 8 files changed, 91 insertions(+), 13 deletions(-) diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.h b/Cantera/src/thermo/GibbsExcessVPSSTP.h index 68a4f3ad5..521f7f344 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.h +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.h @@ -323,7 +323,7 @@ namespace Cantera { * @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) const { + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) { err(" getdlnActCoeffdlnN: nonzero and nonimplemented"); } diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index b2c57cde8..15fc3edb2 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -1051,7 +1051,7 @@ namespace Cantera { } } //==================================================================================================================== - void MargulesVPSSTP::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) const { + void MargulesVPSSTP::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) { s_update_dlnActCoeff_dlnN(); double *data = & dlnActCoeffdlnN_(0,0); for (int k = 0; k < m_kk; k++) { diff --git a/Cantera/src/thermo/MargulesVPSSTP.h b/Cantera/src/thermo/MargulesVPSSTP.h index 8913c53c2..b80b2aa4c 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.h +++ b/Cantera/src/thermo/MargulesVPSSTP.h @@ -847,7 +847,7 @@ namespace Cantera { * @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) const; + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) ; //@} diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index 4bb286c64..5999b0add 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -812,8 +812,8 @@ namespace Cantera { * @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) const { - err(" getdlnActCoeffdlnN: nonzero and nonimplemented"); + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) { + getdlnActCoeffdlnN_numderiv(ld, dlnActCoeffdlnN); } //! returns a summary of the state of the phase as a string diff --git a/Cantera/src/thermo/PhaseCombo_Interaction.cpp b/Cantera/src/thermo/PhaseCombo_Interaction.cpp index ae891ec49..574d3477e 100644 --- a/Cantera/src/thermo/PhaseCombo_Interaction.cpp +++ b/Cantera/src/thermo/PhaseCombo_Interaction.cpp @@ -1110,7 +1110,7 @@ namespace Cantera { /* * HKM - Checked for Transition */ - void PhaseCombo_Interaction::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) const { + void PhaseCombo_Interaction::getdlnActCoeffdlnN(const int ld, doublereal *dlnActCoeffdlnN) { s_update_dlnActCoeff_dlnN(); double *data = & dlnActCoeffdlnN_(0,0); for (int k = 0; k < m_kk; k++) { @@ -1141,7 +1141,6 @@ namespace Cantera { m_pSpecies_A_ij.resize(num, -1); m_pSpecies_B_ij.resize(num, -1); - throw CanteraError("", "unimplemented"); } //==================================================================================================================== diff --git a/Cantera/src/thermo/PhaseCombo_Interaction.h b/Cantera/src/thermo/PhaseCombo_Interaction.h index 59ffb72e0..a14049990 100644 --- a/Cantera/src/thermo/PhaseCombo_Interaction.h +++ b/Cantera/src/thermo/PhaseCombo_Interaction.h @@ -846,7 +846,7 @@ namespace Cantera { * @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) const; + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN); //@} diff --git a/Cantera/src/thermo/ThermoPhase.cpp b/Cantera/src/thermo/ThermoPhase.cpp index 0acf3a7ca..b014eaf89 100644 --- a/Cantera/src/thermo/ThermoPhase.cpp +++ b/Cantera/src/thermo/ThermoPhase.cpp @@ -21,6 +21,7 @@ #endif #include "ThermoPhase.h" +#include "mdp_allo.h" #include //@{ @@ -1098,12 +1099,88 @@ namespace Cantera { * @param dlnActCoeffdN Output vector of derivatives of the * log Activity Coefficients. length = m_kk * m_kk */ - void ThermoPhase::getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) const { - for (int m = 0; m < m_kk; m++) { - for (int k = 0; k < m_kk; k++) { - dlnActCoeffdlnN[ld * k + m] = 0.0; + void ThermoPhase::getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN) { + + + for (int m = 0; m < m_kk; m++) { + for (int k = 0; k < m_kk; k++) { + dlnActCoeffdlnN[ld * k + m] = 0.0; + } } + return; + } + //==================================================================================================================== + void ThermoPhase::getdlnActCoeffdlnN_numderiv(const int ld, doublereal * const dlnActCoeffdlnN) { + + int k, j; + double deltaMoles_j = 0.0; + double pres = pressure(); + + /* + * Evaluate the current base activity coefficients if necessary + */ + std::vector ActCoeff_Base(m_kk); + getActivityCoefficients(DATA_PTR(ActCoeff_Base)); + std::vector Xmol_Base(m_kk); + getMoleFractions(DATA_PTR(Xmol_Base)); + + // Make copies of ActCoeff and Xmol_ for use in taking differences + std::vector ActCoeff(m_kk); + std::vector Xmol(m_kk); + double v_totalMoles = 1.0; + double TMoles_base = v_totalMoles; + + /* + * Loop over the columns species to be deltad + */ + for (j = 0; j < m_kk; j++) { + /* + * Calculate a value for the delta moles of species j + * -> NOte Xmol_[] and Tmoles are always positive or zero + * quantities. + * -> experience has shown that you always need to make the deltas greater than needed to + * change the other mole fractions in order to capture some effects. + */ + double moles_j_base = v_totalMoles * Xmol_Base[j]; + deltaMoles_j = 1.0E-7 * moles_j_base + v_totalMoles * 1.0E-13 + 1.0E-150; + /* + * Now, update the total moles in the phase and all of the + * mole fractions based on this. + */ + v_totalMoles = TMoles_base + deltaMoles_j; + for (k = 0; k < m_kk; k++) { + Xmol[k] = Xmol_Base[k] * TMoles_base / v_totalMoles; + } + Xmol[j] = (moles_j_base + deltaMoles_j) / v_totalMoles; + + /* + * Go get new values for the activity coefficients. + * -> Note this calls setState_PX(); + */ + setState_PX(pres, DATA_PTR(Xmol)); + getActivityCoefficients(DATA_PTR(ActCoeff)); + + /* + * Calculate the column of the matrix + */ + double * const lnActCoeffCol = dlnActCoeffdlnN + ld * j; + for (k = 0; k < m_kk; k++) { + lnActCoeffCol[k] = (2*moles_j_base + deltaMoles_j) *(ActCoeff[k] - ActCoeff_Base[k]) / + ((ActCoeff[k] + ActCoeff_Base[k]) * deltaMoles_j); + } + /* + * Revert to the base case Xmol_, v_totalMoles + */ + v_totalMoles = TMoles_base; + mdp::mdp_copy_dbl_1(DATA_PTR(Xmol), DATA_PTR(Xmol_Base), m_kk); } + /* + * Go get base values for the activity coefficients. + * -> Note this calls setState_TPX() again; + * -> Just wanted to make sure that cantera is in sync + * with VolPhase after this call. + */ + setState_PX(pres, DATA_PTR(Xmol_Base)); } //==================================================================================================================== /* diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index a27d999b8..d1b15ebbf 100644 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -2137,7 +2137,9 @@ namespace Cantera { * @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) const; + virtual void getdlnActCoeffdlnN(const int ld, doublereal * const dlnActCoeffdlnN); + + virtual void getdlnActCoeffdlnN_numderiv(const int ld, doublereal * const dlnActCoeffdlnN); /** * @}