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.

[iurisegtovich]

thanks for the responses,

now, it is all clearer for me!

I have some follow up comments:

1. in the helmholtz expression, i think there may have been a typo

[was]

a = n * (np.log(v / (v - b)) - ak_mix / (b * SQRT2 * 2.0 * state.temperature)
            * np.log((v * (SQRT2 - 1.0) + b) / (v * (SQRT2 + 1.0) - b)))

[should be]

a = n * (np.log(v / (v - b)) - ak_mix / (b * SQRT2 * 2.0 * state.temperature)
            * np.log((v + (1.+ SQRT2) * b) / (v + (1-SQRT2 ) * b)))

then i can see that lnphi converges to zero (and phi converges to one) as volume increases or number of moles decreases

# parameters for propane
tc = SIArray1(369.96 * KELVIN)
pc = SIArray1(4250000.0 * PASCAL)
omega = np.array([0.153])
molar_weight = SIArray1(44.0962 * GRAM / MOL)

# create an instance of our python class and hand it over to rust
eos = EquationOfState.python(PyPengRobinson(tc, pc, omega, molar_weight))


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

gives

data type  :  <class 'builtins.Dual64'>
temperature:  300 + [0]ε
volume     :  24537000000000000000000000000 + [1]ε
moles      :  [602214076000000050000000 + [0]ε]
density    :  0.00002454310127562457 + [-0.0000000000000000000000000000000010002486561366334]ε
A/kT       :  -9847685637422656000000 + [0.00000040023059163927403]ε

data type  :  <class 'builtins.Dual64'>
temperature:  300 + [0]ε
volume     :  24537000000000000000000000000 + [0]ε
moles      :  [602214076000000050000000 + [1]ε]
density    :  0.00002454310127562457 + [0.00000000000000000000000000004075477849777887]ε
A/kT       :  -9847685637422656000000 + [-0.0326597209336492]ε

array([0.98391276])

[g-bauer]

Hi Iuri,

thanks for your interest in FeOs and your detailed look into the notebooks.

Your suggestions are very valuable for us to further improve these examples. As Philipp already mentioned, we should add more/clearer information about the units that are internally used in the State objects. We should also add a paragraph about the ideal gas part in the notebook to avoid confusion about the return value of the (rust) helmholtz_energy() method that by default includes the ideal gas contribution and thus is different from the print statement within the python method.

Concerning the typo, I'll have a look at it as soon as time permits - thanks for the comment.

[iurisegtovich]

i noticed that in the example core_user_defined_eos the helmholtz energy is typed differently
* np.log((SQRT2 - 1.0 + b * rho) / (SQRT2 + 1.0 - b * rho)))
(using rho instead of v)
but i think it has the same problem
the correct should be something similar to
* np.log( (1+(1+SQRT2)* b * rho) / (1 + (1-SQRT2) * b * rho)))

this form with \rho instead of v is similar to that of Smith,vanNess and Abbott
it is the term that they call "I"
(despite they only presenting expression for isobaric-residual-entropy, isobaric-residual-enthalpy, and isobaric-residual-gibbs, not isometric-residual-helmholtz )

[iurisegtovich]

hi,
any update on my question about the possible bug in the expression for Helmholtz energy of Peng Robinson?

[prehner]

Thank you for the reminder, I forgot that there is still this open point in here. As a quick check, I compared the both expressions and it turns out they reproduce the same pressures but not the same Helmholtz energies/chemical potentials/entropies. Maybe they can be rewritten so that this only becomes a question of different reference states or there is something significantly wrong here. I will have a closer look at that asap.

For the future, we should use this discussion forum more in its intended way (think of Stack Overflow) then we won't lose track of open discussions

[iurisegtovich]

Hi, thanks again for the responses,
it am finding it very interesting to use these notebooks.

i have another follow up comment:

i tried extending the "core_dual_numbers" to a binary mixture and testing the "ln_phi_pure" method there
( expecting to get the lnphi for each component if it was in a pure liquid )
then i noticed the "subset" method gets called, and then i got an error,
i think the state object on the subset was rejected by the init for having scalar quantities, instead of arrays of len 1

here is the test

N=1*MOL
P=10*BAR
T=180*KELVIN
x=np.array([.4,.6])

state1 = State(eos, temperature=T, pressure=P, total_moles=N, density_initialization='liquid',molefracs=x)

ans1=state1.ln_phi_pure()

subeos0 = EquationOfState.python( pyeos.subset([0]) )
subeos1 = EquationOfState.python( pyeos.subset([1]) )

state01=State(subeos0, temperature=T, pressure=P, total_moles=N, density_initialization='liquid',)
state11=State(subeos1, temperature=T, pressure=P, total_moles=N, density_initialization='liquid',)

ans2=state01.ln_phi_pure() #is of a pure liquid reference
ans3=state11.ln_phi_pure()

print(ans1,ans2,ans3)

here is the original subset code with prints and error

     def subset(self, i: [int]):
         """Return new equation of state containing a subset of all components."""
         if self.n > 1:
             tc = self.tc[i] 
             pc = self.pc[i]
             mw = self.mw[i]
             omega = self.omega[i]
             print('tc->',tc) #[369.96] -> list
             print('tc->',tc*KELVIN) #369.96 K -> scalar SI
             return PyPengRobinson(tc*KELVIN, pc*PASCAL, omega, mw*GRAM/MOL)
         else:
             return self        
    #PanicException: called `Result::unwrap()` on an `Err` value: PyErr { type: <class 'TypeError'>, value: TypeError("object of type 'si_units.SINumber' has no len()"), traceback: Some(<traceback object at 0x7fe208622e00>) }

here is a potential fix

    def subset(self, i: [int]):
        """Return new equation of state containing a subset of all components."""
        if self.n > 1:
            tc = SIArray1([tci*KELVIN for tci in self.tc[i]]) 
            pc = SIArray1([pci*PASCAL for pci in self.pc[i]]) 
            mw = SIArray1([mwi*GRAM/MOL for mwi in self.mw[i]]) 
            omega = self.omega[i]
            print('tc->',tc) #tc-> [369.96] K -> SI array of length 1
            return PyPengRobinson(tc, pc, omega, mw)
        else:
            return self

[g-bauer]

The focus of the dual number notebook is to explain and showcase (using the print statements) how the data types are changed from within Rust depending on the derivatives that are needed and that there is no special Python syntax needed to work with the generalized (hyper-) dual numbers.

That being said, I agree that the two notebooks (the one explaining dual numbers and the one explaining how to implement the Python class) should use a consistent implementation of the PR EoS. That's certainly something we have to improve.

[iurisegtovich]

regarding my last comment:

the statement
np.array([369.9 , 190.4 ])*KELVIN
gives an feos.si.SIArray1
[369.9, 190.4] K

but the statement
np.array([369.9 ])*KELVIN
gives an feos.si.SINumber
369.9K
which has no len method

i think this caused the subset method to fail in the ln_phi_pure test that i tried to make

[g-bauer]

That's correct. A numpy array of length one will be casted to a scalar when it's considered for multiplication so we have to use SIArray(). This is documented in the core_working_with_units example notebook in the "Arrays of quantities" section.

[iurisegtovich]

i think i understood the behavior,
but the thing is that the "subset" function provided in the example core_dual_numbers won't work as is,
because self.tc arrives as numpy array,

i suggested a fix in the above comment "here is a potential fix"

now i see that in the next example core_user_defined_eos
there is a different subset method that also works

    def subset(self, i: [int]):
        """Return new equation of state containing a subset of all components."""
        if self.n > 1:
            tc = self.tc[i] 
            pc = self.pc[i]
            mw = self.mw[i]
            omega = self.omega[i]
            return PyPengRobinson(
                SIArray1(tc * KELVIN), 
                SIArray1(pc * PASCAL), 
                omega, 
                SIArray1(mw * GRAM / MOL)
            )
        else:
            return self

[g-bauer]

Yes, that's what I was reffering to with my comment above. We will adjust the notebooks so that both implementations are identical to avoid unnecessary confusion.

[iurisegtovich]

uhm ok, got it!

[iurisegtovich]

one more thing, possible bug in the dual operator overloading:

while making changes to the mixing rule of the eos i got to a point where i have a variable A and a variable B
where the types are
A: <class 'numpy.ndarray'>
B: <class 'builtins.Dual2_64'>

(B is one of the state variables and A is a constant of the eos in array format, with len==1)

then B*A gives <class 'builtins.Dual2_64'>
but A*B gives <class 'numpy.ndarray'>

and the latter results in an error when calculating derivatives involved in lnphi

i believe the type of A*B should also be <class 'builtins.Dual2_64'>

seems like the __mul__ of duals, involved in B*A, is ok but
maybe the code of some __rmul__ dunder between <class 'builtins.Dual2_64'> and <class 'numpy.ndarray'> is missing or overriden

ultimately, the error is
TypeError('not implemented!')
when doing A/B
where A is <class 'builtins.Dual2_64'>
1.4212379489852138 + [-0]ε1 + [0]ε1²
and B is <class 'numpy.ndarray'>
[0.38640102878407845 + [-0.002934806836090203]ε1 + [0.00005394160087480313]ε1²]

looks like a Dual accidentally encapsulated into an array

[g-bauer]

We don't export the Python bindings for the dual numbers in FeOs. The above example uses the num_dual package (which defines the dual numbers) which exposes some dual number types and their constructors as well as some helper functions to calculate derivatives. You can find the repository and documentation here.

I am note sure how access to these constructors is needed within a helmholtz_energy method when implementing an equation of state. Wrapping these data types into numpy arrays is already possible.

Since we implemented dual numbers, the rust-numpy crate made it possible to define numpy arrays of custom objects from Rust which would allow us to implement operations with numpy arrays so that is definitely something we will work on in the future.

[iurisegtovich]

i believe access to the dual constructor is really not needed when implementing the EOS,
but with acess to the dual constructor i could change the __mul__ and __rmul__ functions to overcome said limitation of the binding library,
so i could have a sensible array operations instead of "not implemented" error

in case the helholtz energy involved some vectorizable operation like

f (state.temperature) * self.b[:]
(array parametrer) * (dual variable)

( currently i think f(t) * b[:] would work, but b[:] * f(t) would not )

[g-bauer]

Maybe you can try to come up with a proof of concept using the dual numbers and functions exported in num_dual? I am not sure if it's possible to overload methods defined on datatypes defined in extension modules.

[iurisegtovich]

sure!
https://colab.research.google.com/drive/1hbzNPa9VYRTL6IIH9xk9DXL_mJhLnYEW

we can't set attributes of some built-in/extension type as 'int' and 'numpy.ndarray',
but it worked with the Dual64

# int.__mul__=print #TypeError: can't set attributes of built-in/extension type 'int'
# np.ndarray.__mul__ = print #TypeError: can't set attributes of built-in/extension type 'numpy.ndarray'
Dual64.__mul__ = print #works
Dual2_64.from_re(123.) * Dual2_64.from_re(456.)
>>> 456 + [0]ε1 + [0]ε1²

in the colab notebook there is a complete example with Dual64 operating with ndarray of 1 or more dimensions

ans1 = array1*dual1
ans2=dual1*array1
ans3=array2*dual1
ans4=dual1*array2
ans1,ans2,ans3,ans4
#>>> (array([1 + [0]ε], dtype=object),
#>>> array([1 + [0]ε], dtype=object),
#>>> array([1 + [0]ε, 2 + [0]ε], dtype=object),
#>>> array([1 + [0]ε, 2 + [0]ε], dtype=object))

dual1*array22
#>>>array([[0.04320965668136245 + [0]ε, 0.08769116589324533 + [0]ε],
#>>>       [0.6426217879947215 + [0]ε, 0.10796352160245104 + [0]ε]],
#>>>      dtype=object)

array22*dual1  
#>>>array([[0.04320965668136245 + [0]ε, 0.08769116589324533 + [0]ε],
#>>>       [0.6426217879947215 + [0]ε, 0.10796352160245104 + [0]ε]],
#>>>      dtype=object)

[g-bauer]

FYI - we implemented your suggestions to improve working numpy arrays. These operations should now work with the latest versions of num_dual.

[iurisegtovich]

possible issue regarding the example 'core_user_defined_eos'
runs perfectly fine with the pip version
but fails with the latest version

s_cp = State.critical_point(eos)

print("Critical point")

print("temperature: ", s_cp.temperature)

print("density    : ", s_cp.mass_density())

print("pressure   : ", s_cp.pressure())

thread '<unnamed>' panicked at 'called `Result::unwrap()` on an `Err` value: PyErr { type: <class 'AttributeError'>, value: AttributeError("'builtins.StateHDDV2' object has no attribute 'moles'"), traceback: Some(<traceback object at 0x7f2084b0cd00>) }', feos-core/src/python/user_defined.rs:210:1
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

---------------------------------------------------------------------------
PanicException                            Traceback (most recent call last)
/tmp/ipykernel_78197/3722799551.py in <module>
----> 1 s_cp = State.critical_point(eos)
      2 print("Critical point")
      3 print("temperature: ", s_cp.temperature)
      4 print("density    : ", s_cp.mass_density())
      5 print("pressure   : ", s_cp.pressure())

PanicException: called `Result::unwrap()` on an `Err` value: PyErr { type: <class 'AttributeError'>, value: AttributeError("'builtins.StateHDDV2' object has no attribute 'moles'"), traceback: Some(<traceback object at 0x7f2084b0cd00>) }

latest in my case is

commit 54b3ecff5ab960e30bea39efa2a20246488c2611 (HEAD -> main, origin/main, origin/HEAD)
Author: Philipp Rehner <69816385+prehner@users.noreply.github.com>
Date:   Mon Jan 2 13:44:02 2023 +0100

    Add DFT benchmark (#103)

[g-bauer]

That's an issue/bug we are currently working on and that will be fixed in the next version that we release. 😄 You can keep up with the implementation in this issue.

[iurisegtovich]

is this example (from core_working_with_units) working as expected?

# Mathematical operations work on arrays just like they work with scalar values.

# compute the norm of a velocity vector

velocity = np.array([3.5, 2.3, -1.3]) * METER / SECOND

speed = np.linalg.norm(velocity)

speed

>>> [3.5, 2.3, 1.3] m/s

in comparison with

# compute the norm of a velocity vector

velocity = np.array([3.5, 2.3, -1.3]) # METER / SECOND

speed = np.linalg.norm(velocity) * METER / SECOND

speed

>>> 4.3852m/s

[prehner]

It is not. Thanks for bringing it to our attention!

[iurisegtovich]

hi, this also does not seem right to me:

# parameters for propane
tc = SIArray1(369.96 * KELVIN)
pc = SIArray1(4250000.0 * PASCAL)
omega = np.array([0.153])
molar_weight = SIArray1(44.0962 * GRAM / MOL)

# create an instance of our python class and hand it over to rust
pyeos=PyPengRobinson(tc, pc, omega, molar_weight)
eos = EquationOfState.python(pyeos)

# - - - 

N=1*MOL
V=24.537*LITER
T=300*KELVIN
state = State(eos, temperature=T, volume=V, total_moles=N)
state.c_v(Contributions.IdealGas) , state.c_p(Contributions.IdealGas) 
>>>(-0  J/mol/K, 8.31446261815324  J/mol/K)

Cp for ideal pressure gas is Cv+R
but Cv for ideal gas should not be zero.
see for example the simple correlations of smith van ness and abbott in the appendices

[g-bauer]

Currently, there is no way for users to specify the ideal gas contribution for equations of state that are implemented in Python but it's something we will add in the future. The tracking issue is here and we recently started discussing how the whole ideal gas part can be made more transparent to users here.

It's generally possible to enable the use of a contribution that is already implemented in Rust (Joback) or define a custom one (in Python) but it's not something we implemented yet.

As you have recognized, currently for Python implementations, we set the heat capacity to zero.

[iurisegtovich]

i had seen in the linked issue #38 that it was not possible to specify the ideal gas part from the python interface.
but as i had seen the ideal gas part of the helmholtz energy being calculated before,
i thought there was gonna be a default ideal gas part of the cv as well

now i see that the relation for ideal gas helmholtz as function of T is such that gives constant (nonzero) ideal gas entropy, zero ideal gas internal energy and zero ideal gas cv,
so it is missing calorimetric representation,
but it is all consistent .

#propane

N=1*MOL

V=24.537*LITER

T=299*KELVIN

state = State(eos, temperature=T, volume=V, total_moles=N)

a1=state.helmholtz_energy(Contributions.IdealGas)

print( state.c_v(Contributions.IdealGas) , state.helmholtz_energy(Contributions.IdealGas) )

​
T=301*KELVIN

state = State(eos, temperature=T, volume=V, total_moles=N)

a2=state.helmholtz_energy(Contributions.IdealGas)

print( state.c_v(Contributions.IdealGas) , state.helmholtz_energy(Contributions.IdealGas) )

​
T=300*KELVIN

state = State(eos, temperature=T, volume=V, total_moles=N)

​

print( state.entropy(Contributions.IdealGas), '<-- returned constant for any T' )

​

print( (a2-a1)/2 )

​

print( state.internal_energy(Contributions.IdealGas) )

>>> -0  J/mol/K -28.875370770871704 kJ
>>> -0  J/mol/K -29.068517063653456 kJ
>>> 96.57314639087527  J/K <-- returned constant for any T
>>> -96.57314639087599  J
>>> 0  J

is the Joback contribution (in Rust) also added as a ideal gas helmholtz energy modification, in order for the ideal gas entropy, internal energy and Cv remain consistent?

[g-bauer]

For the ideal gas contribution, one has to add a model for the de Broglie wavelength. We provide an implementation using the model of Joback and Reid which is used e.g. with PC-SAFT.

With the de Broglie wavelength provided, the Helmholtz energy of the ideal gas contribution can readily be calculated - and by taking the respective derivatives also entropy, internal energy, etc in a consistent manner.

You can see the Rust trait (interface) for the ideal gas here. And here is the actual implementation for Joback.

[iurisegtovich]

nice, thanks for the responses;
also thanks for the source code links!

[iurisegtovich]

would you like me to open a new discussion specific to this check? i can provide a more thorough example regarding fugacity coefficient trend with both forms

…

On Mon, Mar 20, 2023, 09:10 Philipp Rehner ***@***.***> wrote: Thank you for the reminder, I forgot that there is still this open point in here. As a quick check, I compared the both expressions and it turns out they reproduce the same pressures but not the same Helmholtz energies/chemical potentials/entropies. Maybe they can be rewritten so that this only becomes a question of different reference states or there is something significantly wrong here. I will have a closer look at that asap. For the future, we should use this discussion forum more in its intended way (think of Stack Overflow) then we won't lose track of open discussions — Reply to this email directly, view it on GitHub <#106 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACTMSVLV3IP7IG625RLHDLTW5BCMDANCNFSM6AAAAAATZFXHEI> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[prehner]

Since we are now talking about a potential bug in the code can you open an issue on that specifically?

[iurisegtovich]

sure, i will do it later today

…

On Mon, Mar 20, 2023, 09:26 Philipp Rehner ***@***.***> wrote: Since we are now talking about a potential bug in the code can you open an issue on that specifically? — Reply to this email directly, view it on GitHub <#106 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACTMSVPGT6TSBLA5R6JTCMLW5BEHJANCNFSM6AAAAAATZFXHEI> . You are receiving this because you authored the thread.Message ID: ***@***.***>