Closed iurisegtovich closed 1 year ago
Hi, thank you for the feedback on the notebooks.
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
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.
regarding your suggestions:
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")
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?
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)
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
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
thanks for your attention