Skip to content

core_dual_numbers example - confusion with units and references #105

New issue
New issue
Closed
Closed
core_dual_numbers example - confusion with units and references#105

Description

@iurisegtovich
iurisegtovich
opened on Jan 12, 2023

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

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

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

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

Activity

  1. prehner commented on Jan 12, 2023

    @prehner
    prehner
    on Jan 12, 2023
    Contributor

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

    1. 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.
  2. feos-org
    locked and limited conversation to collaborators on Jan 12, 2023
  3. prehner
    converted this issue into a discussion #106 on Jan 12, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions