diff --git a/Cantera/src/thermo/DebyeHuckel.cpp b/Cantera/src/thermo/DebyeHuckel.cpp index 94822538c..fbc348f44 100644 --- a/Cantera/src/thermo/DebyeHuckel.cpp +++ b/Cantera/src/thermo/DebyeHuckel.cpp @@ -32,7 +32,6 @@ namespace Cantera { MolalityVPSSTP(), m_formDH(DHFORM_DILUTE_LIMIT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(30.0), m_useHelgesonFixedForm(false), @@ -62,7 +61,6 @@ namespace Cantera { MolalityVPSSTP(), m_formDH(DHFORM_DILUTE_LIMIT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(30.0), m_useHelgesonFixedForm(false), @@ -85,7 +83,6 @@ namespace Cantera { MolalityVPSSTP(), m_formDH(DHFORM_DILUTE_LIMIT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(3.0), m_useHelgesonFixedForm(false), @@ -132,7 +129,6 @@ namespace Cantera { MolalityVPSSTP::operator=(b); m_formDH = b.m_formDH; m_formGC = b.m_formGC; - m_Pcurrent = b.m_Pcurrent; m_Aionic = b.m_Aionic; m_npActCoeff = b.m_npActCoeff; m_IionicMolality = b.m_IionicMolality; @@ -172,7 +168,7 @@ namespace Cantera { } - /** + /* * ~DebyeHuckel(): (virtual) * * Destructor for DebyeHuckel. Release objects that @@ -187,7 +183,7 @@ namespace Cantera { } } - /** + /* * duplMyselfAsThermoPhase(): * * This routine operates at the ThermoPhase level to @@ -199,7 +195,7 @@ namespace Cantera { return (ThermoPhase *) mtp; } - /** + /* * Equation of state type flag. The base class returns * zero. Subclasses should define this to return a unique * non-zero value. Constants defined for this purpose are @@ -227,7 +223,7 @@ namespace Cantera { // // -------- Molar Thermodynamic Properties of the Solution --------------- // - /** + /* * Molar enthalpy of the solution. Units: J/kmol. */ doublereal DebyeHuckel::enthalpy_mole() const { @@ -235,7 +231,7 @@ namespace Cantera { return mean_X(DATA_PTR(m_tmpV)); } - /** + /* * Molar internal energy of the solution. Units: J/kmol. * * This is calculated from the soln enthalpy and then @@ -249,7 +245,7 @@ namespace Cantera { return uu; } - /** + /* * Molar soln entropy at constant pressure. Units: J/kmol/K. * * This is calculated from the partial molar entropies. @@ -259,13 +255,13 @@ namespace Cantera { return mean_X(DATA_PTR(m_tmpV)); } - /// Molar Gibbs function. Units: J/kmol. + // Molar Gibbs function. Units: J/kmol. doublereal DebyeHuckel::gibbs_mole() const { getChemPotentials(DATA_PTR(m_tmpV)); return mean_X(DATA_PTR(m_tmpV)); } - /** + /* * Molar heat capacity at constant pressure. Units: J/kmol/K. * * Returns the solution heat capacition at constant pressure. @@ -289,7 +285,7 @@ namespace Cantera { // ------- Mechanical Equation of State Properties ------------------------ // - /** + /* * Pressure. Units: Pa. * For this incompressible system, we return the internally storred * independent value of the pressure. @@ -298,22 +294,29 @@ namespace Cantera { return m_Pcurrent; } - /** + /* * Set the pressure at constant temperature. Units: Pa. * This method sets a constant within the object. * The mass density is not a function of pressure. */ void DebyeHuckel::setPressure(doublereal p) { + #ifdef DEBUG_MODE //printf("setPressure: %g\n", p); #endif - double temp = temperature(); - if (m_waterSS) { - /* - * Call the water SS and set it's internal state - */ - m_waterSS->setTempPressure(temp, p); + /* + * Store the current pressure + */ + m_Pcurrent = p; + /* + * update the standard state thermo + * -> This involves calling the water function and setting the pressure + */ + _updateStandardStateThermo(); + + if (m_waterSS) { + /* * Store the internal density of the water SS. * Note, we would have to do this for all other @@ -321,10 +324,6 @@ namespace Cantera { */ m_densWaterSS = m_waterSS->density(); } - /* - * Store the current pressure - */ - m_Pcurrent = p; /* * Calculate all of the other standard volumes * -> note these are constant for now @@ -361,7 +360,7 @@ namespace Cantera { } - /** + /* * The isothermal compressibility. Units: 1/Pa. * The isothermal compressibility is defined as * \f[ @@ -377,7 +376,7 @@ namespace Cantera { return 0.0; } - /** + /* * The thermal expansion coefficient. Units: 1/K. * The thermal expansion coefficient is defined as * @@ -418,7 +417,7 @@ namespace Cantera { } } - /** + /* * Overwritten setMolarDensity() function is necessary because the * density is not an indendent variable. * @@ -441,9 +440,7 @@ namespace Cantera { * the value propagates to underlying objects. */ void DebyeHuckel::setTemperature(doublereal temp) { - if (m_waterSS) { - m_waterSS->setTemperature(temp); - } + _updateStandardStateThermo(); State::setTemperature(temp); } @@ -452,7 +449,7 @@ namespace Cantera { // ------- Activities and Activity Concentrations // - /** + /* * This method returns an array of generalized concentrations * \f$ C_k\f$ that are defined such that * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ @@ -495,7 +492,7 @@ namespace Cantera { return 1.0 / mvSolvent; } - /** + /* * Returns the natural logarithm of the standard * concentration of the kth species */ @@ -504,7 +501,7 @@ namespace Cantera { return log(c_solvent); } - /** + /* * Returns the units of the standard and general concentrations * Note they have the same units, as their divisor is * defined to be equal to the activity of the kth species @@ -538,7 +535,7 @@ namespace Cantera { } - /** + /* * Get the array of non-dimensional activities at * the current solution temperature, pressure, and * solution concentration. @@ -546,6 +543,7 @@ namespace Cantera { * */ void DebyeHuckel::getActivities(doublereal* ac) const { + _updateStandardStateThermo(); /* * Update the molality array, m_molalities() * This requires an update due to mole fractions @@ -561,7 +559,7 @@ namespace Cantera { exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent; } - /** + /* * getMolalityActivityCoefficients() (virtual, const) * * Get the array of non-dimensional Molality based @@ -574,7 +572,7 @@ namespace Cantera { */ void DebyeHuckel:: getMolalityActivityCoefficients(doublereal* acMolality) const { - + _updateStandardStateThermo(); A_Debye_TP(-1.0, -1.0); s_update_lnMolalityActCoeff(); copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality); @@ -586,7 +584,7 @@ namespace Cantera { // // ------ Partial Molar Properties of the Solution ----------------- // - /** + /* * Get the species chemical potentials. Units: J/kmol. * * This function returns a vector of chemical potentials of the @@ -633,7 +631,7 @@ namespace Cantera { } - /** + /* * Returns an array of partial molar enthalpies for the species * in the mixture. * Units (J/kmol) @@ -756,7 +754,7 @@ namespace Cantera { } } - /** + /* * getPartialMolarVolumes() (virtual, const) * * returns an array of partial molar volumes of the species @@ -838,7 +836,7 @@ namespace Cantera { * in the Solution ------------------ */ - /** + /* * getStandardChemPotentials() (virtual, const) * * @@ -855,6 +853,7 @@ namespace Cantera { * units = J / kmol */ void DebyeHuckel::getStandardChemPotentials(doublereal* mu) const { + _updateStandardStateThermo(); getGibbs_ref(mu); doublereal pref; doublereal delta_p; @@ -868,7 +867,7 @@ namespace Cantera { } } - /** + /* * Get the nondimensional gibbs function for the species * standard states at the current T and P of the solution. * @@ -884,6 +883,7 @@ namespace Cantera { * standard state gibbs function for species k. */ void DebyeHuckel::getGibbs_RT(doublereal* grt) const { + _updateStandardStateThermo(); getPureGibbs(grt); doublereal invRT = 1.0 / _RT(); for (int k = 0; k < m_kk; k++) { @@ -891,7 +891,7 @@ namespace Cantera { } } - /** + /* * * getPureGibbs() * @@ -924,6 +924,7 @@ namespace Cantera { */ void DebyeHuckel:: getEnthalpy_RT(doublereal* hrt) const { + _updateStandardStateThermo(); getEnthalpy_RT_ref(hrt); doublereal pref; doublereal delta_p; @@ -956,6 +957,7 @@ namespace Cantera { */ void DebyeHuckel:: getEntropy_R(doublereal* sr) const { + _updateStandardStateThermo(); getEntropy_R_ref(sr); if (m_waterSS) { sr[0] = m_waterSS->entropy_mole(); @@ -982,6 +984,7 @@ namespace Cantera { * constant pressure heat capacity for species k. */ void DebyeHuckel::getCp_R(doublereal* cpr) const { + _updateStandardStateThermo(); getCp_R_ref(cpr); if (m_waterSS) { cpr[0] = m_waterSS->cp_mole(); @@ -996,6 +999,7 @@ namespace Cantera { * units = m^3 / kmol */ void DebyeHuckel::getStandardVolumes(doublereal *vol) const { + _updateStandardStateThermo(); copy(m_speciesSize.begin(), m_speciesSize.end(), vol); if (m_waterSS) { @@ -1003,7 +1007,6 @@ namespace Cantera { vol[0] = molecularWeight(0)/dd; } } - /* * ------ Thermodynamic Values for the Species Reference States --- @@ -1815,7 +1818,7 @@ namespace Cantera { return dAdT; } - /** + /* * d2A_DebyedT2_TP() (virtual) * * Returns the 2nd derivative of the A_Debye parameter with @@ -1849,7 +1852,7 @@ namespace Cantera { return d2AdT2; } - /** + /* * dA_DebyedP_TP() (virtual) * * Returns the derivative of the A_Debye parameter with @@ -2692,6 +2695,32 @@ namespace Cantera { } } + /* + * Updates the standard state thermodynamic functions at the current T and P of the solution. + * + * @internal + * + * This function gets called for every call to functions in this + * class. It checks to see whether the temperature or pressure has changed and + * thus the ss thermodynamics functions for all of the species + * must be recalculated. + * + * + * Note, this will throw an error. It must be reimplemented in derived classes. + */ + void DebyeHuckel::_updateStandardStateThermo(doublereal pnow) const { + _updateRefStateThermo(); + doublereal tnow = temperature(); + if (pnow == -1.0) { + pnow = m_Pcurrent; + } + if (m_tlast != tnow || m_plast != pnow) { + if (m_waterSS) { + m_waterSS->setTempPressure(tnow, pnow); + } + m_tlast = tnow; + m_plast = pnow; + } + } + } - - diff --git a/Cantera/src/thermo/DebyeHuckel.h b/Cantera/src/thermo/DebyeHuckel.h index 5205e7176..7873c368f 100644 --- a/Cantera/src/thermo/DebyeHuckel.h +++ b/Cantera/src/thermo/DebyeHuckel.h @@ -21,65 +21,6 @@ namespace Cantera { - - /*! - - * - * Major Parameters: - * - * m_formDH = Form of the Debye-Huckel expression - * - * DHFORM_DILUTE_LIMIT = 0 - * - * This form assumes a dilute limit to DH, and is mainly - * for informational purposes: - * - * ln(gamma_k) = -z_k**2 * alpha * sqrt(I) - * - * where I = 1/2 sum_k( molality_k * z_k**2) - * - * DHFORM_BDOT_AK = 1 - * - * This form assumes Bethke's format for the DH coefficient - * - * ln(gamma_k) = -z_k**2 * alpha * sqrt(I) / (1 + B * a_k * sqrt(I)) - * + bdot_k * I - * - * (note, this particular form where a_k can differ in - * multielectrolyte - * solutions has problems wrt a gibbs-duhem analysis. However - * we include it here because there is a lot of data fit to it) - * - * DHFORM_BDOT_AUNIFORM = 2 - * - * This form assumes Bethke's format for the DH coefficient - * - * ln(gamma_k) = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) - * + bdot_k * I - * - * The value of a is determined at the beginning of the - * calculation, and not changed. - * - * DHFORM_BETAIJ = 3 - * - * This form assumes a linear expansion in a virial coefficient form - * It is used extensively in Newmann's book, and is the beginning of - * more complex treatments for stronger electrolytes, like Pitzer - * and HMW treatments. - * - * ln(gamma_k) = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) - * + 2* sum_j (beta_jk m_j) - * - * DHFORM_PITZER_BETAIJ = 4 - * - * This form assumes an activity coefficient formulation consistent - * with a truncated form of Pitzer's formulation. - * - * ln(gamma_k) = -z_k**2 * alpha * sqrt(I) / (1 + B * a * sqrt(I)) - * -2 * z_k**2 * alpha * ln(1 + B * a * sqrt(I)) / (B * a) - * + 2 * sum_j (beta_jk m_j) - * - */ /*! * @name Formats for the Activity Coefficients * @@ -318,25 +259,57 @@ namespace Cantera { *
* %Application within %Kinetics Managers *
- * The standard concentration is equal to 1.0. This means that the - * kinetics operator works on an (activities basis). Since this - * is a stoichiometric substance, this means that the concentration - * of this phase drops out of kinetics expressions. + * For the time being, we have set the standard concentration for all species in + * this phase equal to the default concentration of the solvent at 298 K and 1 atm. + * This means that the + * kinetics operator essentially works on an activities basis, with units specified + * as if it were on a concentration basis. + * + * For example, a bulk-phase binary reaction between liquid species j and k, producing + * a new liquid species l would have the + * following equation for its rate of progress variable, \f$ R^1 \f$, which has + * units of kmol m-3 s-1. * - * An example of a reaction using this is a sticking coefficient - * reaction of a substance in an ideal gas phase on a surface with a bulk phase - * species in this phase. In this case, the rate of progress for this - * reaction, \f$ R_s \f$, may be expressed via the following equation: * \f[ - * R_s = k_s C_{gas} + * R^1 = k^1 C_j^a C_k^a = k^1 (C_o a_j) (C_o a_k) * \f] - * where the units for \f$ R_s \f$ are kmol m-2 s-1. \f$ C_{gas} \f$ has units - * of kmol m-3. Therefore, the kinetic rate constant, \f$ k_s \f$, has - * units of m s-1. Nowhere does the concentration of the bulk phase - * appear in the rate constant expression, since it's a stoichiometric - * phase and the activity is always equal to 1.0. + * where + * \f[ + * C_j^a = C_o a_j \quad and \quad C_k^a = C_o a_k + * \f] + * + * \f$ C_j^a \f$ is the activity concentration of species j, and + * \f$ C_k^a \f$ is the activity concentration of species k. \f$ C_o \f$ + * is the concentration of water at 298 K and 1 atm. \f$ a_j \f$ is + * the activity of species j at the current temperature and pressure + * and concentration of the liquid phase. \f$k^1 \f$ has units of m3 kmol-1 s-1. * + * The reverse rate constant can then be obtained from the law of microscopic reversibility + * and the equilibrium expression for the system. + * + * \f[ + * \frac{a_j a_k}{ a_l} = K^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} ) + * \f] + * + * \f$ K^{o,1} \f$ is the dimensionless form of the equilibrium constant. + * + * \f[ + * R^{-1} = k^{-1} C_l^a = k^{-1} (C_o a_l) + * \f] + * + * where + * + * \f[ + * k^{-1} = k^1 K^{o,1} C_o + * \f] + * + * \f$k^{-1} \f$ has units of s-1. + * + * Note, this treatment may be modified in the future, as events dictate. + * + *
* Instantiation of the Class + *
* * The constructor for this phase is NOT located in the default ThermoFactory * for %Cantera. However, a new %StoichSubstanceSSTP may be created by @@ -397,8 +370,8 @@ namespace Cantera { @endverbatim * - * The model attribute, "StoichSubstanceSSTP", on the thermo element identifies the phase as being - * a StoichSubstanceSSTP object. + * The model attribute, "StoichSubstanceSSTP", on the thermo element identifies the phase as + * being a StoichSubstanceSSTP object. * */ class DebyeHuckel : public MolalityVPSSTP { @@ -662,8 +635,8 @@ namespace Cantera { * @{ */ - /** - * This method returns an array of generalized concentrations + //! This method returns an array of generalized concentrations + /*! * \f$ C_k\f$ that are defined such that * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ * is a standard concentration @@ -698,7 +671,8 @@ namespace Cantera { * @param k Optional parameter indicating the species. The default * is to assume this refers to species 0. * @return - * Returns the standard Concentration in units of m3 kmol-1. + * Returns the standard Concentration in units of + * m3 kmol-1. */ virtual doublereal standardConcentration(int k=0) const; @@ -1357,6 +1331,24 @@ namespace Cantera { double _osmoticCoeffHelgesonFixedForm() const; double _lnactivityWaterHelgesonFixedForm() const; + protected: + + //! Updates the standard state thermodynamic functions at the current T and P of the solution. + /*! + * @internal + * + * This function gets called for every call to functions in this + * class. It checks to see whether the temperature or pressure has changed and + * thus the ss thermodynamics functions for all of the species + * must be recalculated. + * + * + * Note, this will throw an error. It must be reimplemented in derived classes. + */ + virtual void _updateStandardStateThermo(doublereal pres = -1.0) const; + + + //@} @@ -1408,11 +1400,6 @@ namespace Cantera { */ int m_formGC; - /** - * Current pressure in Pascal - */ - double m_Pcurrent; - //! Vector containing the electrolyte species type /*! * The possible types are: diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index 1af0b911e..bccc13695 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -31,7 +31,6 @@ namespace Cantera { m_formPitzer(PITZERFORM_BASE), m_formPitzerTemp(PITZER_TEMP_CONSTANT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(100.0), m_TempPitzerRef(298.15), @@ -61,7 +60,6 @@ namespace Cantera { m_formPitzer(PITZERFORM_BASE), m_formPitzerTemp(PITZER_TEMP_CONSTANT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(100.0), m_TempPitzerRef(298.15), @@ -85,7 +83,6 @@ namespace Cantera { m_formPitzer(PITZERFORM_BASE), m_formPitzerTemp(PITZER_TEMP_CONSTANT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(100.0), m_TempPitzerRef(298.15), @@ -133,7 +130,6 @@ namespace Cantera { m_formPitzer = b.m_formPitzer; m_formPitzerTemp = b.m_formPitzerTemp; m_formGC = b.m_formGC; - m_Pcurrent = b.m_Pcurrent; m_Aionic = b.m_Aionic; m_IionicMolality = b.m_IionicMolality; m_maxIionicStrength = b.m_maxIionicStrength; @@ -265,7 +261,6 @@ namespace Cantera { m_formPitzer(PITZERFORM_BASE), m_formPitzerTemp(PITZER_TEMP_CONSTANT), m_formGC(2), - m_Pcurrent(OneAtm), m_IionicMolality(0.0), m_maxIionicStrength(30.0), m_TempPitzerRef(298.15), diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index 3693e25f7..c5bc0e237 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -996,13 +996,6 @@ namespace Cantera { * bimolecular rxns which have units of m-3 kmol-1 s-1.) */ int m_formGC; - - /** - * Current pressure in Pascal. This is now the independent variable - * as it must be for multicomponent solutions. - */ - double m_Pcurrent; - vector_int m_electrolyteSpeciesType; diff --git a/Cantera/src/thermo/IdealMolalSoln.cpp b/Cantera/src/thermo/IdealMolalSoln.cpp index 5bde94fbd..29dec3aee 100644 --- a/Cantera/src/thermo/IdealMolalSoln.cpp +++ b/Cantera/src/thermo/IdealMolalSoln.cpp @@ -23,7 +23,6 @@ namespace Cantera { */ IdealMolalSoln::IdealMolalSoln() : MolalityVPSSTP(), - m_Pcurrent(OneAtm), m_formGC(2) { } @@ -55,7 +54,6 @@ namespace Cantera { if (&b != this) { MolalityVPSSTP::operator=(b); m_speciesMolarVolume = b.m_speciesMolarVolume; - m_Pcurrent = b.m_Pcurrent; m_formGC = b.m_formGC; m_expg0_RT = b.m_expg0_RT; m_pe = b.m_pe; diff --git a/Cantera/src/thermo/IdealMolalSoln.h b/Cantera/src/thermo/IdealMolalSoln.h index 687dc75ca..ddc93d665 100644 --- a/Cantera/src/thermo/IdealMolalSoln.h +++ b/Cantera/src/thermo/IdealMolalSoln.h @@ -1003,13 +1003,6 @@ namespace Cantera { */ array_fp m_speciesMolarVolume; - /*! - * Current pressure in Pascal. - * - * This is an independent variable in the problem. - */ - double m_Pcurrent; - /** * The standard concentrations can have three different forms * depending on the value of the member attribute m_formGC, which diff --git a/Cantera/src/thermo/SingleSpeciesTP.h b/Cantera/src/thermo/SingleSpeciesTP.h index 0617cb4db..7df40963c 100644 --- a/Cantera/src/thermo/SingleSpeciesTP.h +++ b/Cantera/src/thermo/SingleSpeciesTP.h @@ -670,9 +670,17 @@ namespace Cantera { doublereal m_tmin; //! Upper value of the temperature for which reference thermo is valid doublereal m_tmax; - //! Current value of the pressure (Pascals) + + //! The current pressure of the solution (Pa) + /*! + * It gets initialized to 1 atm. + */ doublereal m_press; - //! Value of the reference pressure (Pascals) + + /*! + * Reference pressure (Pa) must be the same for all species + * - defaults to 1 atm. + */ doublereal m_p0; //! Last temperature used to evaluate the thermodynamic polynomial. diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index 27de64bc9..b8e3893be 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -30,8 +30,13 @@ namespace Cantera { */ VPStandardStateTP::VPStandardStateTP() : ThermoPhase(), + m_Pcurrent(OneAtm), m_tlast(-1.0), - m_plast(-1.0) + m_tlast_ref(-1.0), + m_plast(-1.0), + m_p0(OneAtm), + m_useTmpRefStateStorage(true), + m_useTmpStandardStateStorage(false) { } @@ -46,8 +51,13 @@ namespace Cantera { */ VPStandardStateTP::VPStandardStateTP(const VPStandardStateTP &b) : ThermoPhase(), + m_Pcurrent(OneAtm), m_tlast(-1.0), - m_plast(-1.0) + m_tlast_ref(-1.0), + m_plast(-1.0), + m_p0(OneAtm), + m_useTmpRefStateStorage(true), + m_useTmpStandardStateStorage(false) { *this = b; } @@ -69,13 +79,17 @@ namespace Cantera { /* * However, we have to handle data that we own. */ + m_Pcurrent = b.m_Pcurrent; m_tlast = b.m_tlast; + m_tlast_ref = b.m_tlast_ref; m_plast = b.m_plast; + m_p0 = b.m_p0; + m_useTmpRefStateStorage = b.m_useTmpRefStateStorage; m_h0_RT = b.m_h0_RT; m_cp0_R = b.m_cp0_R; m_g0_RT = b.m_g0_RT; m_s0_R = b.m_s0_R; - m_V0 = b.m_V0; + m_useTmpStandardStateStorage = b.m_useTmpStandardStateStorage; m_hss_RT = b.m_hss_RT; m_cpss_R = b.m_cpss_R; m_gss_RT = b.m_gss_RT; @@ -148,18 +162,33 @@ namespace Cantera { } void VPStandardStateTP::getEnthalpy_RT(doublereal* hrt) const { - _updateStandardStateThermo(); - copy(m_hss_RT.begin(), m_hss_RT.end(), hrt); + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_hss_RT.begin(), m_hss_RT.end(), hrt); + } else { + err("getEnthalpy_RT ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); + } } void VPStandardStateTP::getEntropy_R(doublereal* srt) const { - _updateStandardStateThermo(); - copy(m_sss_R.begin(), m_sss_R.end(), srt); + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_sss_R.begin(), m_sss_R.end(), srt); + } else { + err("getEntropy_R ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); + } } void VPStandardStateTP::getGibbs_RT(doublereal* grt) const { - _updateStandardStateThermo(); - copy(m_gss_RT.begin(), m_gss_RT.end(), grt); + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_gss_RT.begin(), m_gss_RT.end(), grt); + } else { + err("getGibbs_RT ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); + } } void VPStandardStateTP::getPureGibbs(doublereal* g) const { @@ -171,23 +200,38 @@ namespace Cantera { } void VPStandardStateTP::getIntEnergy_RT(doublereal* urt) const { - _updateStandardStateThermo(); - copy(m_hss_RT.begin(), m_hss_RT.end(), urt); - doublereal RT = _RT(); - doublereal tmp = pressure() / RT; - for (int k = 0; k < m_kk; k++) { - urt[k] -= tmp * m_Vss[k]; + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_hss_RT.begin(), m_hss_RT.end(), urt); + doublereal RT = _RT(); + doublereal tmp = pressure() / RT; + for (int k = 0; k < m_kk; k++) { + urt[k] -= tmp * m_Vss[k]; + } + } else { + err("getIntEnergy_RT ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); } } void VPStandardStateTP::getCp_R(doublereal* cpr) const { - _updateStandardStateThermo(); - copy(m_cpss_R.begin(), m_cpss_R.end(), cpr); + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_cpss_R.begin(), m_cpss_R.end(), cpr); + } else { + err("getCp_R ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); + } } void VPStandardStateTP::getStandardVolumes(doublereal *vol) const { - _updateStandardStateThermo(); - copy(m_Vss.begin(), m_Vss.end(), vol); + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(); + copy(m_Vss.begin(), m_Vss.end(), vol); + } else { + err("getStandardVolumes ERROR: Must be overwritten in child classes"); + _updateStandardStateThermo(); + } } /* @@ -200,16 +244,23 @@ namespace Cantera { * the reference pressure for the species. */ void VPStandardStateTP::getEnthalpy_RT_ref(doublereal *hrt) const { - /* - * Call the function that makes sure the local copy of the - * species reference thermo functions are up to date for the - * current temperature. - */ - _updateRefStateThermo(); - /* - * Copy the enthalpy function into return vector. - */ - copy(m_h0_RT.begin(), m_h0_RT.end(), hrt); + if (m_useTmpRefStateStorage) { + /* + * Call the function that makes sure the local copy of the + * species reference thermo functions are up to date for the + * current temperature. + */ + _updateRefStateThermo(); + /* + * Copy the enthalpy function into return vector. + */ + copy(m_h0_RT.begin(), m_h0_RT.end(), hrt); + } else if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(m_p0); + copy(m_hss_RT.begin(), m_hss_RT.end(), hrt); + } else { + err("getEnthalpy_RT_ref() ERROR: not handled"); + } } /* @@ -218,16 +269,23 @@ namespace Cantera { * of the solution and the reference pressure for the species. */ void VPStandardStateTP::getGibbs_RT_ref(doublereal *grt) const { - /* - * Call the function that makes sure the local copy of - * the species reference thermo functions are up to date - * for the current temperature. - */ - _updateRefStateThermo(); - /* - * Copy the gibbs function into return vector. - */ - copy(m_g0_RT.begin(), m_g0_RT.end(), grt); + if (m_useTmpRefStateStorage) { + /* + * Call the function that makes sure the local copy of + * the species reference thermo functions are up to date + * for the current temperature. + */ + _updateRefStateThermo(); + /* + * Copy the gibbs function into return vector. + */ + copy(m_g0_RT.begin(), m_g0_RT.end(), grt); + } else if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(m_p0); + copy(m_gss_RT.begin(), m_gss_RT.end(), grt); + } else { + err("getGibbs_RT_ref() ERROR: not handled"); + } } /* @@ -253,16 +311,23 @@ namespace Cantera { * of the solution and the reference pressure for the species. */ void VPStandardStateTP::getEntropy_R_ref(doublereal *er) const { - /* - * Call the function that makes sure the local copy of - * the species reference thermo functions are up to date - * for the current temperature. - */ - _updateRefStateThermo(); - /* - * Copy the gibbs function into return vector. - */ - copy(m_s0_R.begin(), m_s0_R.end(), er); + if (m_useTmpRefStateStorage) { + /* + * Call the function that makes sure the local copy of + * the species reference thermo functions are up to date + * for the current temperature. + */ + _updateRefStateThermo(); + /* + * Copy the gibbs function into return vector. + */ + copy(m_s0_R.begin(), m_s0_R.end(), er); + } else if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(m_p0); + copy(m_sss_R.begin(), m_sss_R.end(), er); + } else { + err("getEntropy_R_ref() ERROR: not handled"); + } } /* @@ -272,17 +337,39 @@ namespace Cantera { * and reference pressure for the species. */ void VPStandardStateTP::getCp_R_ref(doublereal *cpr) const { - /* - * Call the function that makes sure the local copy of - * the species reference thermo functions are up to date - * for the current temperature. - */ - _updateRefStateThermo(); - /* - * Copy the gibbs function into return vector. - */ - copy(m_cp0_R.begin(), m_cp0_R.end(), cpr); + if (m_useTmpRefStateStorage) { + /* + * Call the function that makes sure the local copy of + * the species reference thermo functions are up to date + * for the current temperature. + */ + _updateRefStateThermo(); + /* + * Copy the gibbs function into return vector. + */ + copy(m_cp0_R.begin(), m_cp0_R.end(), cpr); + } else if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(m_p0); + copy(m_cpss_R.begin(), m_cpss_R.end(), cpr); + } else { + err("getCp_R_ref() ERROR: not handled"); + } } + + /* + * Get the molar volumes of the species reference states at the current + * T and P_ref of the solution. + * + * units = m^3 / kmol + */ + void VPStandardStateTP::getStandardVolumes_ref(doublereal *vol) const { + if (m_useTmpStandardStateStorage) { + _updateStandardStateThermo(m_p0); + copy(m_Vss.begin(), m_Vss.end(), vol); + } else { + err("getStandardVolumes_ref() ERROR: not handled"); + } + } /* * Perform initializations after all species have been @@ -300,16 +387,19 @@ namespace Cantera { void VPStandardStateTP::initLengths() { m_kk = nSpecies(); int leng = m_kk; - m_h0_RT.resize(leng); - m_g0_RT.resize(leng); - m_cp0_R.resize(leng); - m_s0_R.resize(leng); - m_V0.resize(leng); - m_hss_RT.resize(leng); - m_gss_RT.resize(leng); - m_cpss_R.resize(leng); - m_sss_R.resize(leng); - m_Vss.resize(leng); + if (m_useTmpRefStateStorage){ + m_h0_RT.resize(leng); + m_g0_RT.resize(leng); + m_cp0_R.resize(leng); + m_s0_R.resize(leng); + } + if (m_useTmpStandardStateStorage) { + m_hss_RT.resize(leng); + m_gss_RT.resize(leng); + m_cpss_R.resize(leng); + m_sss_R.resize(leng); + m_Vss.resize(leng); + } } /* @@ -337,21 +427,34 @@ namespace Cantera { /* * void _updateRefStateThermo() (protected, virtual, const) * - * This function gets called for every call to functions in this - * class. It checks to see whether the temperature has changed and + * This function checks to see whether the temperature has changed and * thus the reference thermodynamics functions for all of the species * must be recalculated. + * It must be called for every reference state function evaluation, + * if m_useTmpRefStateStorage is set to true. * If the temperature has changed, the species thermo manager is called - * to recalculate G, Cp, H, and S at the current temperature. + * to recalculate the following internal arrays at the current temperature and at + * the reference pressure: + * + * - m_h0_RT + * - m_g0_RT + * - m_s0_R + * - m_cp0_R + * + * This function may be reimplemented in child objects. However, it doesn't + * necessarily have to be, if the species thermo manager can carry + * out the full calculation. */ void VPStandardStateTP::_updateRefStateThermo() const { - doublereal tnow = temperature(); - if (m_tlast != tnow) { - m_spthermo->update(tnow, DATA_PTR(m_cp0_R), DATA_PTR(m_h0_RT), - DATA_PTR(m_s0_R)); - m_tlast = tnow; - for (int k = 0; k < m_kk; k++) { - m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + if (m_spthermo) { + doublereal tnow = temperature(); + if (m_tlast_ref != tnow) { + m_spthermo->update(tnow, DATA_PTR(m_cp0_R), DATA_PTR(m_h0_RT), + DATA_PTR(m_s0_R)); + m_tlast_ref = tnow; + for (int k = 0; k < m_kk; k++) { + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } } } } @@ -359,16 +462,31 @@ namespace Cantera { /* * void _updateStandardStateThermo() (protected, virtual, const) * - * This function gets called for every call to functions in this - * class. It checks to see whether the temperature has changed and + * If m_useTmpStandardStateStorage is true, + * This function must be called for every call to functions in this + * class that need standard state properties. + * Child classes may require that it be called even if m_useTmpStandardStateStorage + * is not true. + * It checks to see whether the temperature has changed and * thus the ss thermodynamics functions for all of the species * must be recalculated. + * + * This */ - void VPStandardStateTP::_updateStandardStateThermo() const { + void VPStandardStateTP::_updateStandardStateThermo(doublereal pnow) const { + _updateRefStateThermo(); doublereal tnow = temperature(); - doublereal pnow = pressure(); + if (pnow == -1.0) { + pnow = pressure(); + } if (m_tlast != tnow || m_plast != pnow) { - err("getStandardVolumes"); + err("_updateStandardStateThermo ERROR: Must be overwritten in child classes"); + /* + * Redo objects that need reevaluation. + */ + for (int k = 0; k < m_kk; k++) { + m_g0_RT[k] = m_g0_RT[k]; + } m_tlast = tnow; m_plast = pnow; } @@ -376,5 +494,3 @@ namespace Cantera { } - - diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index 51b3dfb49..1e9965d5a 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -34,15 +34,53 @@ namespace Cantera { * a variable pressure standard state for species. * * Several concepts are introduced. The first concept is there are temporary - * variables for holding the species standard values of Cp, H, S, and V at the + * variables for holding the species standard state values + * of Cp, H, S, G, and V at the * last temperature and pressure called. These functions are not recalculated * if a new call is made using the previous temperature and pressure. * * There are also temporary - * variables for holding the species reference-state values of Cp, H, S, and V at the + * variables for holding the species reference-state values of Cp, H, S, and G at the * last temperature and reference pressure called. These functions are not recalculated * if a new call is made using the previous temperature. * + * To support the above functionality, pressure and temperature variables, + * m_plast and m_tlast, are kept which store the last pressure and temperature + * used in the evaluation of standard state properties. An optional utility is provided + * to store the results from the last temperature and pressure standard + * state calculation and use it on subsequent calculations, if the temperature + * and pressure are unchanged. + * + * If #m_useTmpRefStateStorage is set to true, then the following internal + * arrays, containing information about the reference arrays, + * are calculated and kept up to date at every call. + * + * - #m_h0_RT + * - #m_g0_RT + * - #m_s0_R + * - #m_cp0_R + * + * The virtual function #_updateRefStateThermo() is supplied to do this + * and may be reimplemented in child routines. A default implementation + * based on the speciesThermo class is supplied in this base class. + * #_updateStandardStateThermo() is called whenever a reference state property is needed. + * + * When #m_useTmpStandardStateStorage is true, then the following + * internal arrays, containing information on the standard state properties + * are calculated and kept up to date. + * + * - #m_hss_RT; + * - #m_cpss_R; + * - #m_gss_RT; + * - #m_sss_R; + * - #m_Vss + * + * The virtual function #_updateStandardStateThermo() is supplied to do this + * and must be reimplemented in child routines, when #m_useTmpStandardStateStorage is true. + * It may be optionally reimplemented in child routines if + * #m_useTmpStandardStateStorage is false. + * #_updateStandardStateThermo() is called whenever a standard state property is needed. + * * This class is usually used for nearly incompressible phases. For those phases, it * makes sense to change the equation of state independent variable from density to pressure. * @@ -99,19 +137,21 @@ namespace Cantera { /// @name Partial Molar Properties of the Solution (VPStandardStateTP) //@{ - /** - * Get the array of non-dimensional species chemical potentials - * These are partial molar Gibbs free energies. + + //! Get the array of non-dimensional species chemical potentials + //! These are partial molar Gibbs free energies. + /*! * \f$ \mu_k / \hat R T \f$. * Units: unitless * * We close the loop on this function, here, calling - * getChemPotentials() and then dividing by RT. + * getChemPotentials() and then dividing by RT. No need for child + * classes to handle. * * @param mu Output vector of non-dimensional species chemical potentials * Length: m_kk. */ - virtual void getChemPotentials_RT(doublereal* mu) const; + void getChemPotentials_RT(doublereal* mu) const; //@} @@ -166,15 +206,20 @@ namespace Cantera { */ virtual void getGibbs_RT(doublereal* grt) const; - /** - * Get the nondimensional Gibbs functions for the standard - * state of the species at the current T and P. + + //! Get the nondimensional Gibbs functions for the standard + //! state of the species at the current T and P. + /*! + * (Note resolved at this level) * * @param gpure Output vector of standard state * Gibbs free energies. length = m_kk. * units are J/kmol. + * + * @todo This could be eliminated. It doesn't fit into the current + * naming convention. */ - virtual void getPureGibbs(doublereal* gpure) const; + void getPureGibbs(doublereal* gpure) const; /** * Returns the vector of nondimensional @@ -224,24 +269,33 @@ namespace Cantera { /*! * @internal * - * This function gets called for every call to functions in this + * If m_useTmpStandardStateStorage is true, + * this function must be called for every call to functions in this * class. It checks to see whether the temperature or pressure has changed and * thus the ss thermodynamics functions for all of the species * must be recalculated. * - * This function is responsible for updating the following internal members: + * This function is responsible for updating the following internal members, + * when m_useTmpStandardStateStorage is true. * - * m_hss_RT; - * m_cpss_R; - * m_gss_RT; - * m_sss_R; - * m_Vss + * - m_hss_RT; + * - m_cpss_R; + * - m_gss_RT; + * - m_sss_R; + * - m_Vss + * + * If m_useTmpStandardStateStorage is not true, this function may be + * required to be called by child classes to update internal member data. * * Note, this will throw an error. It must be reimplemented in derived classes. + * + * @param pres Pressure at which to carry out the calculation. + * The default is to use the current pressure, storred in m_Pcurrent. */ - virtual void _updateStandardStateThermo() const; + virtual void _updateStandardStateThermo(doublereal pres = -1.0) const; public: + //@} /// @name Thermodynamic Values for the Species Reference States (VPStandardStateTP) /*! @@ -298,7 +352,7 @@ namespace Cantera { */ virtual void getEntropy_R_ref(doublereal *er) const; - /** + /*! * Returns the vector of nondimensional * constant pressure heat capacities of the reference state * at the current temperature of the solution @@ -310,15 +364,15 @@ namespace Cantera { */ virtual void getCp_R_ref(doublereal *cprt) const; - //! Recalculate the Reference state thermo functions + //! Get the molar volumes of the species reference states at the current + //! T and P_ref of the solution. /*! - * This function checks to see whether the temperature has changed and - * thus the reference thermodynamics functions for all of the species - * must be recalculated. - * If the temperature has changed, the species thermo manager is called - * to recalculate G, Cp, H, and S at the current temperature and at - * the reference pressure. + * units = m^3 / kmol + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. */ + virtual void getStandardVolumes_ref(doublereal *vol) const; protected: @@ -327,12 +381,24 @@ namespace Cantera { * This function checks to see whether the temperature has changed and * thus the reference thermodynamics functions for all of the species * must be recalculated. + * It must be called for every reference state function evaluation, + * if m_useTmpRefStateStorage is set to true. * If the temperature has changed, the species thermo manager is called - * to recalculate G, Cp, H, and S at the current temperature and at - * the reference pressure. + * to recalculate the following internal arrays at the current temperature and at + * the reference pressure: + * + * - m_h0_RT + * - m_g0_RT + * - m_s0_R + * - m_cp0_R + * + * This function may be reimplemented in child objects. However, it doesn't + * necessarily have to be, if the species thermo manager can carry + * out the full calculation. */ virtual void _updateRefStateThermo() const; + //@} @@ -416,13 +482,34 @@ namespace Cantera { //@} protected: + + //! The current pressure of the solution (Pa) + /*! + * It gets initialized to 1 atm. + */ + mutable doublereal m_Pcurrent; //! The last temperature at which the reference thermodynamic properties were calculated at. mutable doublereal m_tlast; + //! The last temperature at which the reference thermodynamic properties were calculated at. + mutable doublereal m_tlast_ref; + //! The last pressure at which the Standard State thermodynamic properties were calculated at. mutable doublereal m_plast; + /*! + * Reference pressure (Pa) must be the same for all species + * - defaults to 1 atm. + */ + doublereal m_p0; + + /*! + * boolean indicating whether temporary reference state storage is used + * -> default is true + */ + bool m_useTmpRefStateStorage; + /*! * Vector containing the species reference enthalpies at T = m_tlast * and P = p_ref. @@ -447,11 +534,11 @@ namespace Cantera { */ mutable vector_fp m_s0_R; - /** - * Vector containing the species reference volumes - * at T = m_tlast and P = p_ref - */ - mutable vector_fp m_V0; + /*! + * boolean indicating whether temporary standard state storage is used + * -> default is false + */ + bool m_useTmpStandardStateStorage; /** * Vector containing the species Standard State enthalpies at T = m_tlast @@ -482,6 +569,8 @@ namespace Cantera { * at T = m_tlast and P = m_plast */ mutable vector_fp m_Vss; + + private: diff --git a/test_problems/cathermo/DH_graph_1/runtest b/test_problems/cathermo/DH_graph_1/runtest index 60de953e0..79ef10d1b 100755 --- a/test_problems/cathermo/DH_graph_1/runtest +++ b/test_problems/cathermo/DH_graph_1/runtest @@ -51,7 +51,6 @@ fi /bin/rm -f test.out test.diff DH_NaCl_acommon.csv echo 'Testing the DH dilute act calculation - act vs I' -/bin/rm DH_NaCl_acommon.csv $prog DH_NaCl_acommon.xml > DH_NaCl_acommon.csv retnStat=$? echo 'Making a comparison with the good saved solution: '