core_dual_numbers example - confusion with units and references

[iurisegtovich]

hi,
I was reading and trying the first example 'core_dual_numbers'
and got confused with some aspects:
(some of which apply to the example "core_user_defined_eos.ipynb as well")

see attachment notebook examples.zip

in the rev1 attachment i highlight some points of confusion
in the rev2 attachment i try to refactor the Helmholtz function of the EOS to S.I. and Reid/Poling/Prausnitz notation, and try to get the same results as before for the dimensionless results (a_kt)

main issues:

1. results for a_kt of the pyfunc and a_kt of the "state.helmholtz_energy()/(KB300KELVIN)" method does not seem to match - what does that mean?

see values in the attachment

2. the reference of the chemical potential, and its relation with fugacity coefficient is not clear for me,
   also lnphi and lnphi_pure do not seem to match, what does that mean?
   (in my opinion the chemical potential / lnphi is a very motivating application for autodiffing)

see values in the attachment

3. i observed that the pyobject STATE, in the notebook body, contains quantities with units,
   but the type of the object state that arrives at the Helmholtz function of the class is of the type "<class 'builtins.StateF'>"
   and contains floats with implicit units - that also caused me some confusion when trying to take advantage of the units system inside the Helmholtz function, perhaps rename the argument of the "def helmholtz_energy(self, state):" to
   "def helmholtz_energy(self, state_floats):" / warn in the docstring of the class in the example that the received object contains only floats and array of floats (mol fracs)

suggestions:

1. the example represents a bulk system (no confinement), but the values used for the example are of few angstrom**3 and few molecules size

please see suggestion of macroscopic values

2. the example uses values inside the helmholtz energy not in S.I., and with nonstandard notation for the a and a/(bt) or a/(brt) term

please see suggestion of my attempt at an eos reprogrammed in S.I. and with notation from reid/poling/prausnitz
where a has R**2 instead of R and helmholtz uses a term like a/(brt) instead of a/(bt)

thanks for your attention

[prehner]

Hi, thank you for the feedback on the notebooks.

1. The difference between the value printed within the Python function and the result of state.helmholtz_energy is that the first one is the residual Helmholtz energy (that can be computed from an eos like PR) whereas feos adds an ideal gas contribution. You can call

state.helmholtz_energy(Contributions.ResidualNvt) / (KB * 300*KELVIN)

to verify that the output is indeed the same

data type  :  <class 'float'>
temperature:  300.0
volume     :  40744.0
moles      :  [1.0]
density    :  2.4543491066169253e-05
A/kT       :  [5.06008546]

5.060085461999782

2. For the chemical potential, again the contributions are important. If you want to manually evaluate the fugacity coefficient

lnphi = state.chemical_potential(Contributions.ResidualNvt)/RGAS/T - np.log(Z)

will give you the correct result.

ln_phi and ln_phi_pure give the same result but only for stable liquid phases, e.g.,

N=300*MOL
V=24.537*LITER
T=300*KELVIN
state = State(eos, temperature=T, volume=V, total_moles=N)
state.ln_phi_pure(), state.ln_phi()

(array([3.0549047]), array([3.0549047]))

The point of the ln_phi_pure() function is basically to help with the calculation of activity coefficients. We should reflect that more clearly in the name actually.

3. The data type of the state that is passed changes depending on the derivative that is being calculated. For a simple evaluation of the Helmholtz energy it is a state with floats, for first derivatives (e.g., pressure) it is a state with dual numbers. To include units together with the data types for automatic differentiation would be extremely cumbersome. The units are chosen to feel natural when implementing molecular (SAFT-type) equations of state. For the implementation of cubics it is indeed somewhat unergonomic. Just to stress the point: for the implementation of models (basically a one time task) the model needs to be expressed in reduced (molecular) units. While using the model everything is expressed very explicitly and comfortably in SI units.

regarding your suggestions:

1. We have either awkardly small inputs (volumes in ANGSTROM, total_moles in # of molecules rather than moles) or awkwardly large outputs in the example print statements in the Helmholtz energy. In this notebook, the focus is on how the data type in the Helmholtz energy function depends on the derivative, therefore, the focus was to have a nice output.
2. Yes, I guess I addressed that already. To make the example more readable, it would indeed be helpful to give the variables more meaningful names.