Make Matlab species thermo properties consistent and useful

[speth]

The documentation for the Matlab function enthalpies_RT says it returns the "standard-state" species enthalpies. However, it ultimately calls the C++ method ThermoPhase::getEnthalpy_RT_ref, which returns what Cantera calls the "reference state" enthalpies. This represents the confluence of several problems:

• The name of the function does not indicate that it is returning values for the standard state rather than the current state
• No function exists in the Matlab toolbox to return the actual species enthalpies
• The terminology of "reference state" vs "standard state" is inconsistent between the Matlab toolbox and the C++ base, with the Matlab terminology (I think) more closely matching the more common thermodynamic convention.

The same issue applies for gibbs_RT, entropies_R, and Cp_R, except that for these, the documentation doesn't even say that the values returned are for the reference state.

[decaluwe]

Hi Ray,

I've been thinking about your third bullet, some more:

The complex part is that there are two relevant "standard states:"

1. A standard state for reference species properties. This is the one that most people think about, and what you identify as the "common thermodynamic convention": species properties are readily calculated for a given temperature and reference pressure, assuming that X_k = 1.0 for a given species. This could also be called the "reference state" in that is it literally the state for which you have known species properties which you can reference.

2. However, there is also a phase "standard state." This distinct from and cannot be equal to the species reference state, as it refers to the mixture properties at a given state (including composition), and there is obviously no state for which X_k = 1 for all species. This state is frequently required for derivations many thermodynamic equations of state. It is a state at which phase properties (e.g. mixture-averaged properties) can be readily calculated. For a non-ideal EoS, for example, departure functions can be calculated by integrating from the "standard state" to the actual state.

So the two are distinct and necessarily so. The situation is easily confused because of the simplicity of the ideal gas EoS, where no "standard phase state" is actually required. One can calculate all thermodynamic properties quite readily, directly from the reference species properties, and mixture rules do not require any complex mathematics. So it's very easy (and quite common) to use the two terms synonymously.

That said, do we need both (or either) "standard state" and "reference species property" calculations featured in the codebase? You can't really avoid having the reference data loaded and stored as such, but I'm somewhat skeptical that we need the standard state data explicitly calculated or stored.

1. The whole point of the standard state is that its properties are readily calculated. So why not just plug these expressions into the departure function calculations? Explicitly storing both reference and standard information is needlessly confusing.
2. I don't really see the "standard state" values being anything that one would with to access for subsequent calculations, beyond what is already done in the class itself.
3. Also, explicitly calculating the standard state info separate from the departure functions is perhaps dangerous, if someone uses the current class as a basis for a new derived class (i.e. the two cannot be modified independently).

As I'm typing this, an additional confusion just reared its ugly head: the standardConcentrations used in the kinetic calculations are also wrapped up in this confusion. More specifically, the species activity concentrations are gamma_k*X_k/vbar, where gamma_k is the activity coefficient and vbar the mixture-averaged molar volume.

The form of gamma_k comes from the EoS: when one writes the chemical potential for a species k, it can be represented as:
mu_k = mu_k_o + RT ln(a_k) = mu_k_o + RT ln(gamma_k*X_k/vbar/C_k_o)

where C_k_o are the standardConcentrations, and mu_k_o are the species chemical potentials at the standard concentrations. There doesn't really seem to be any clear "correct" choice as to whether this is a "standard" state or a "reference" state:

1. On the one hand, the activity coefficients come from the EoS, so it would make sense that the same "standard state" used to derive the EoS is that used to calculate the activity coefficients.
2. On the other hand, this is clearly a "species property" calculation, and the mu_k_o will obviously involve the "reference" species data (regardless of whether they perform additional calculations).
3. The difference is all semantics, in this case: I just checked the gamma_k calculations for the RedlichKwongMFTP and the (in development) PengRobinsonMFTP classes, and the choice of this standard state has no effect on the calculated activity coefficient. So the choice of C_k_o simply affects the value of mu_k_o.

But obviously, it would be nice to get some clarity, here. I just don't yet know what that should be. I'm leaning toward just having the "reference" state being the only one referenced in the code, which seems to be your preference, as well. I'm just still a little wary that this may have unintended consequences or cause other unanticipated confusions.

[speth]

A status update for the "experimental" Matlab toolbox: The ThermoPhase.enthalpies_RT property ultimately calls ThermoPhase::getEnthalpy_RT_ref, as before. However, there are no properties to replace the old gibbs_RT, entropies_R, or Cp_R methods. I think this needs to be addressed in order to complete Cantera/enhancements#177, whether or not we ever settle on a consistent definition of what "standard state" means.

[ischoegl]

For the sake of consistency, this issue would be resolved by providing the same interface for MATLAB as for Python, correct? In the case of MATLAB's ThermoPhase.enthalpies_RT, Python uses ThermoPhase.standard_enthalpies_RT (while calling different C++ functions on the backend, i.e. ThermoPhase::getEnthalpy_RT_ref and ThermoPhase::getEnthalpy_RT).

The distinctions of the corresponding docstrings aren't necessarily clear to me with

cantera/include/cantera/thermo/ThermoPhase.h

Lines 944 to 951 in 9d0e54a

	//! Returns the vector of nondimensional enthalpies of the reference state
	//! at the current temperature of the solution and the reference pressure
	//! for the species.
	/*!
	* @param hrt Output vector containing the nondimensional reference
	* state enthalpies. Length: m_kk.
	*/
	virtual void getEnthalpy_RT_ref(double* hrt) const {

vs.

cantera/include/cantera/thermo/ThermoPhase.h

Lines 866 to 872 in 9d0e54a

	//! Get the nondimensional Enthalpy functions for the species at their
	//! standard states at the current *T* and *P* of the solution.
	/*!
	* @param hrt Output vector of nondimensional standard state enthalpies.
	* Length: m_kk.
	*/
	virtual void getEnthalpy_RT(double* hrt) const {

with experimental MATLAB

cantera/interfaces/matlab_experimental/Base/ThermoPhase.m

Lines 122 to 133 in 9d0e54a

	% Non-dimensional enthalpies ::
	%
	% >> v = tp.enthalpies_RT
	%
	% :param tp:
	% Instance of class :mat:class:`ThermoPhase` (or another
	% object that derives from ThermoPhase).
	% :return:
	% Vector of standard-state species enthalpies constant
	% and T is the temperature. For gaseous species, these
	% values are ideal gas enthalpies.
	enthalpies_RT

and Python

cantera/interfaces/cython/cantera/thermo.pyx

Lines 1664 to 1668 in 9d0e54a

	property standard_enthalpies_RT:
	"""
	Array of nondimensional species standard-state enthalpies at the
	current temperature and pressure.
	"""

The legacy MATLAB toolbox stated

%     Vector of standard-state species enthalpies
%     divided by RT, where R is the universal gas
%     constant and T is the temperature. For gaseous species, these
%     values are ideal gas enthalpies.

which indicates some copy-pasta issue in the new toolbox docstring.

Overall, is this just about one using a reference pressure while the other isn't? I.e. MATLAB references data that are typically provided in textbook appendices (e.g. Appendix A, Glassman 3rd ed)?

What is the relationship to what is calculated here? Is this what MATLAB uses?

cantera/include/cantera/thermo/SpeciesThermoInterpType.h

Lines 191 to 205 in 9d0e54a

	//! Compute the reference-state property of one species
	/*!
	* Given temperature T in K, this method updates the values of the non-
	* dimensional heat capacity at constant pressure, enthalpy, and entropy, at
	* the reference pressure, of the species.
	*
	* @param temp Temperature (Kelvin)
	* @param cp_R Vector of Dimensionless heat capacities. (length m_kk).
	* @param h_RT Vector of Dimensionless enthalpies. (length m_kk).
	* @param s_R Vector of Dimensionless entropies. (length m_kk).
	*/
	virtual void updatePropertiesTemp(const double temp,
	double* cp_R,
	double* h_RT,
	double* s_R) const;

The other properties mentioned on top would be ThermoPhase.standard_gibbs_RT, ThermoPhase.standard_entropies_R and ThermoPhase.standard_cp_R in the Python API.

Given that we have a new experimental MATLAB toolbox, there's a window of opportunity to make things consistent.

PS: one alternative would be to switch the Python API to the MATLAB convention by dropping the standard_ prefix?

PPS: this is code from (now removed) PDSS_IdealGas.cpp, which may shed some light on this? I.e. the ss presumably refers to "standard-state"?

void PDSS_IdealGas::setTemperature(double temp)
{
    m_temp = temp;
    m_spthermo->updatePropertiesTemp(temp, &m_cp0_R, &m_h0_RT, &m_s0_R);
    m_g0_RT = m_h0_RT - m_s0_R;
    m_V0 = GasConstant * m_temp / m_p0;
    m_hss_RT = m_h0_RT;
    m_cpss_R = m_cp0_R;
    m_sss_R = m_s0_R - log(m_pres/m_p0);
    m_gss_RT = m_hss_RT - m_sss_R;
    m_Vss = GasConstant * m_temp / m_pres;
}

[ischoegl]

Alright. Some further digging cleared things up (especially in conjunction with @decaluwe's comments). With definitions

cantera/include/cantera/thermo/PDSS.h

Lines 457 to 467 in 9d0e54a

	protected:
	double m_h0_RT; //!< Reference state enthalpy divided by RT
	double m_cp0_R; //!< Reference state heat capacity divided by R
	double m_s0_R; //!< Reference state entropy divided by R
	double m_g0_RT; //!< Reference state Gibbs free energy divided by RT
	double m_V0; //!< Reference state molar volume (m^3/kmol)
	double m_hss_RT; //!< Standard state enthalpy divided by RT
	double m_cpss_R; //!< Standard state heat capacity divided by R
	double m_sss_R; //!< Standard state entropy divided by R
	double m_gss_RT; //!< Standard state Gibbs free energy divided by RT
	double m_Vss; //!< Standard State molar volume (m^3/kmol)

and, similarly

cantera/include/cantera/thermo/VPStandardStateTP.h

Lines 259 to 261 in 9d0e54a

	//! Vector containing the species reference enthalpies at T = m_tlast
	//! and P = p_ref.
	mutable vector<double> m_h0_RT;

and

cantera/include/cantera/thermo/VPStandardStateTP.h

Lines 278 to 280 in 9d0e54a

	//! Vector containing the species Standard State enthalpies at T = m_tlast
	//! and P = m_plast.
	mutable vector<double> m_hss_RT;

as well as

cantera/src/thermo/VPStandardStateTP.cpp

Lines 46 to 50 in 9d0e54a

	void VPStandardStateTP::getEnthalpy_RT(double* hrt) const
	{
	updateStandardStateThermo();
	std::copy(m_hss_RT.begin(), m_hss_RT.end(), hrt);
	}

and

cantera/src/thermo/VPStandardStateTP.cpp

Lines 99 to 103 in 9d0e54a

	void VPStandardStateTP::getEnthalpy_RT_ref(double* hrt) const
	{
	updateStandardStateThermo();
	std::copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
	}

everything follows quite nicely.

Re-reading @decaluwe's comments, I tend to agree that it may not make sense to provide an API for standard-state quantities, which I believe (?) I read from

The situation is easily confused because of the simplicity of the ideal gas EoS, where no "standard phase state" is actually required. One can calculate all thermodynamic properties quite readily, directly from the reference species properties, and mixture rules do not require any complex mathematics. So it's very easy (and quite common) to use the two terms synonymously.

[...] do we need both (or either) "standard state" and "reference species property" calculations featured in the codebase? You can't really avoid having the reference data loaded and stored as such, but I'm somewhat skeptical that we need the standard state data explicitly calculated or stored.

and

I don't really see the "standard state" values being anything that one would with to access for subsequent calculations, beyond what is already done in the class itself.

The crux is that standard_* quantities are used in one Python sample (coverage_dependent_surf.py) as well as several tests. Better documentation with consistent nomenclature, in combination with consistent Python and MATLAB API's may be the best we can do?

PS: here's further context (sorry if I am documenting this here as I am reading up to this)

cantera/include/cantera/thermo/ThermoPhase.h

Lines 178 to 200 in 9d0e54a

	* ## Calculating and accessing thermodynamic properties
	*
	* The calculation of thermodynamic functions within %ThermoPhase is broken down roughly
	* into two or more steps. First, the standard state properties of all of the species
	* are calculated at the current temperature and at either the current pressure or at a
	* reference pressure. If the calculation is carried out at a reference pressure instead
	* of at the current pressure the calculation is called a "reference state properties"
	* calculation, just to make the distinction (even though it may be considered to be a
	* fixed-pressure standard-state calculation). The next step is to adjust the reference
	* state calculation to the current pressure. The thermodynamic functions then are
	* considered to be at the standard state of each species. Lastly the mixing
	* contributions are added to arrive at the thermodynamic functions for the solution.
	*
	* The %ThermoPhase class provides interfaces to thermodynamic properties calculated for
	* the reference state of each species, the standard state values for each species, the
	* thermodynamic functions for solution values, both on a per mole of solution basis
	* (such as ThermoPhase::enthalpy_mole()), on a per kg of solution basis, and on a
	* partial molar basis for each species (such as
	* ThermoPhase::getPartialMolarEnthalpies). At each level, functions for the enthalpy,
	* entropy, Gibbs free energy, internal energy, and volume are provided. So, 5 levels
	* (reference state, standard state, partial molar, per mole of solution, and per mass
	* of solution) and 5 functions multiplied together makes 25 possible functions. That's
	* why %ThermoPhase is such a large class.

In conclusion, I'd be in favor to not expose the first two levels (reference state and standard state) at advanced APIs, i.e. deprecate standard_* quantities in Python, and introduce partial_molar_* quantities in MATLAB.

[speth]

I'm 👍 on leaving the "reference state" properties out of the high-level APIs (and perhaps, they should even be protected members of the C++ ThermoPhase class) as they are an implementation detail.

The "standard state" properties are not as important as the "partial molar" ones, but I think there are some cases where they can be useful to look at. Of course, this runs into the perpetual issue of Cantera's definition of "standard state" being different from the ones generally used in thermodynamics textbooks.

The other source of variability in here that doesn't really help is the mix of dimensional versus nondimensional quantities, for example providing nondimensional reference state but dimensional partial molar properties. Likewise, here, I think it would be best to regard the nondimensional properties as implementation details. It's probably not worth immediately removing such accessors, but at least we can avoid increasing the proliferation of these functions.