core_dual_numbers example - confusion with units and references about feos

[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.