On Dual Numbers in FeOs¶
In FeOs, we use generalized dual numbers to compute partial derivatives of the Helmholtz energy. In this notebook, we take a closer look at how that works.
We will make use of the EquationOfState.python() feature: we use a simple version of the Peng-Robinson equation of state implemented as Python class which is registered in FeOs as equation of state. If you want to learn how to implement an equation of state as Python class in conjunction with FeOs, take a look at the example implementation in the examples section. In short - a class that implements a helmholtz_energy function that takes a StateD (or any state object, where
the internal data types can be any generalized dual number) and returns the Helmholtz energy (plus some minor functions), can be used with FeOs. We use the equation of state implemented in Python simply because we can add print statements and inspect the data types at runtime.
Dual Numbers¶
We won’t go into much detail regarding the theory of generalized dual numbers. Suffice it to say that generalized dual numbers are mathematical constructs consisting of one real value and one or more non-real values (just like complex numbers) that enable evaluation of arithmetic operations such that the operation and the derivative(s) can simultaneously be computed. We call these generalized dual numbers because they can have a different number of non-real values. For example a dual number has a real and a single non-real value and can be used to compute first (partial) derivatives, while a hyper-dual number has a real and three non-real values and can be used to compute first and second (partial) derivatives. We can - in principle - create numbers that allow computation of an arbitrarily high derivative at the cost of execution speed.
Similar to numerical differentiation and complex step differentiation, a derivative is computed by defining a “step size”. For complex step and dual number differentiation this step is introduced in the non-real part. Derivatives of dual numbers, however, are exact (to machine precision) and independent of the step size used. Hence, we use unity as step size.
For example, when we feed a temperature as Dual64 data type
print(temperature)
300 + [1]ε
into a function, its return value will be a Dual64 as well. The non-real part of the function result is the derivative with respect to temperature.
Do we need to create dual numbers by hand?¶
No. If you use already implemented equations of state, you will not even recognize that dual numbers are used. However, if you use Python or Rust to implement you own equation of state, it is useful to know how FeOs creates and uses dual numbers. It’s easier to look at an example than to talk about it. Below, you find the implementation of the Peng-Robinson equation of state where we added some print statement to the helmholtz_energy function to keep track of the values and data
types of the thermodynamic state that enters the Helmholtz energy computation.
[1]:
from feos.eos import *
from feos.si import *
import numpy as np
SQRT2 = np.sqrt(2)
class PyPengRobinson:
def __init__(self, critical_temperature, critical_pressure, acentric_factor, molar_weight, delta_ij=None):
"""Peng-Robinson Equation of State
Parameters
----------
critical_temperature : SIArray1
critical temperature of each component.
critical_pressure : SIArray1
critical pressure of each component.
acentric_factor : np.array[float]
acentric factor of each component (dimensionless).
molar_weight: SIArray1
molar weight of each component.
delta_ij : np.array[[float]], optional
binary parameters. Shape=[n, n], n = number of components.
defaults to zero for all binary interactions.
Raises
------
ValueError: if the input values have incompatible sizes.
"""
self.n = len(critical_temperature)
if len(set((len(critical_temperature), len(critical_pressure), len(acentric_factor)))) != 1:
raise ValueError("Input parameters must all have the same lenght.")
if self.n == 1:
self.tc = (critical_temperature / KELVIN)[0]
self.pc = (critical_pressure / PASCAL)[0]
self.omega = acentric_factor[0]
self.mw = (molar_weight / GRAM * MOL)[0]
else:
self.tc = critical_temperature / KELVIN
self.pc = critical_pressure / PASCAL
self.omega = acentric_factor
self.mw = molar_weight / GRAM * MOL
self.a_r = 0.45724 * critical_temperature**2 * RGAS / critical_pressure / ANGSTROM**3 / NAV / KELVIN
self.b = 0.07780 * critical_temperature * RGAS / critical_pressure / ANGSTROM**3 / NAV
self.kappa = 0.37464 + (1.54226 - 0.26992 * acentric_factor) * acentric_factor
self.delta_ij = np.zeros((self.n, self.n)) if delta_ij is None else delta_ij
def helmholtz_energy(self, state):
"""Return Helmholtz energy.
Parameters
----------
state : StateHD
The thermodynamic state.
Returns
-------
helmholtz_energy: float | any dual number
The return type depends on the input types.
"""
n = np.sum(state.moles)
x = state.molefracs
tr = 1.0 / self.tc * state.temperature
ak = ((1.0 - np.sqrt(tr)) * self.kappa + 1.0)**2 * self.a_r
ak_mix = 0.0
if self.n > 1:
for i in range(self.n):
for j in range(self.n):
ak_mix += np.sqrt(ak[i] * ak[j]) * (x[i] * x[j] * (1.0 - self.delta_ij[i, j]))
else:
ak_mix = ak
b = np.sum(x * self.b)
v = 1.0 / state.density
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)))
# some print statements to inspect data types
print()
print("data type : ", type(state.temperature))
print("temperature: ", state.temperature)
print("volume : ", state.volume)
print("moles : ", state.moles)
print("density : ", state.density)
print("A/kT : ", a)
return a
def components(self) -> int:
"""Number of components."""
return self.n
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(tc*KELVIN, pc*PASCAL, omega, mw*GRAM/MOL)
else:
return self
def molar_weight(self) -> SIArray1:
if isinstance(self.mw, float):
return np.array([self.mw]) * GRAM / MOL
else:
return self.mw * GRAM / MOL
def max_density(self, moles:[float]) -> float:
b = np.sum(moles * self.b) / np.sum(moles);
return 0.9 / b
[2]:
# 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))
Properties and dual numbers¶
After initializing the equation of state for a single substance (in the above cell), let us now build a thermodynamic state. For the sake of this example, we use the natural variables of the Helmholtz energy (\(\mathbf{N}\), \(V\), \(T\)) as input so that no volume or temperature has to be iterated.
[3]:
state = State(eos, temperature=300*KELVIN, volume=40744*ANGSTROM**3, total_moles=1/NAV)
When computing a thermodynamic property for a State, FeOs internally generates a new object which contains the same information as State but with dual numbers as data types. Which dual numbers (dual, hyper-dual, etc.) are used depends on the property to compute or rather which partial derivatives are needed. FeOs then modifies the non-real parts of those properties (temperature, volume or amount of substance) so that the correct derivatives are computed and feeds this
“generalized dual state” object into the function that computes the Helmholtz energy. The result will have the same data type as the state’s variables and the needed derivative(s) are extracted and returned.
FeOs does all of that independent from the implementation of the equation of state. All you have to do - if you implement a Helmholtz energy function - is to write your code knowing that the state’s properties are generalized dual numbers. In Python, your can write code just like “regular” Python because dual numbers implement all arithmetic opartions you’d typically use.
Let’s take a look at an example. First, we compute the molar_helmholtz_energy. Since no derivatives are needed to compute the Helmholtz energy, the data types of the State’s properties are regular floating point numbers (float).
If you are in a running noteboook, execute the cell below and check the output. If you run the cell a second time, you will notice that the information about data types and values aren’t printed to the screen anymore. That’s because we cache already computed derivatives for a given state. A second call to molar_helmholtz_energy simply pulls the already computed value from cache and returns it without entering the helmholtz_energy function.
[4]:
state.helmholtz_energy() / (KB * 300*KELVIN)
data type : <class 'float'>
temperature: 300.0
volume : 40744.0
moles : [1.0]
density : 2.4543491066169253e-05
A/kT : [5.06008546]
[4]:
-6.554978406564657
Next, we compute the pressure, for which we compute the first partial derivative of the Helmholtz energy with respect to the volume.
\(p = -\left(\frac{\partial A}{\partial V}\right)_{\mathbf{n}, T}\)
If you execute the cell below, you’ll notice that the data types are now dual numbers (Dual64, a dual number with 64bit floating point numbers) with a single non-real part. This non-real or dual part of the volume is unity while all other non-real parts are zero, which means that the first partial derivative with respect to volume will be computed.
FeOs modifies the dual parts of the temperature, volume or amount of substance (for a component) when computing derivatives. Note that the density is calculated as \(\rho = \frac{\sum n_i}{V}\) (where \(n_i\) is the amount of substance of component \(i\)) and hence also has a dual part which is different from zero. It is also different from unity - the value is calculated according to the artihmetic operation of division between two dual numbers.
[5]:
state.pressure()
data type : <class 'builtins.PyDual64'>
temperature: 300 + [0]ε
volume : 40744 + [1]ε
moles : [1 + [0]ε]
density : 0.000024543491066169253 + [-0.00000000060238295371513]ε
A/kT : 5.060085461999783 + [0.00000040024326944247174]ε
[5]:
Similar to the pressure, the entropy is computed by taking the first partial derivative with respect to temperature. Here, the dual part of the density is zero because both dual parts of the volume and the amount of substance are zero, too.
[6]:
state.molar_entropy()
data type : <class 'builtins.PyDual64'>
temperature: 300 + [1]ε
volume : 40744 + [0]ε
moles : [1 + [0]ε]
density : 0.000024543491066169253 + [0]ε
A/kT : 5.060085461999783 + [-0.025513116823522065]ε
[6]:
Let’s take a look at a more involved property. For the Joule-Thomson coefficient, we need multiple second partial derivatives:
\(\mu_{JT}=\left(\frac{\partial T}{\partial p}\right)_{H,N_i} = -\frac{1}{C_p} \left(V + T \left(\frac{\partial p}{\partial T}\right)_{\mathbf{n}, V} \left(\frac{\partial p}{\partial V}\right)_{\mathbf{n}, T}^{-1} \right)\)
We need three evaluations of the Helmholtz energy (that are actually occuring in the computation of \(C_p\)):
\(\left(\frac{\partial p}{\partial T}\right)_{\mathbf{n}, V} \rightarrow -\left(\frac{\partial A}{\partial T \partial V}\right)_\mathbf{n}\): second partial derivative w.r.t volume and temperature,
\(\left(\frac{\partial p}{\partial V}\right)_{\mathbf{n}, T} \rightarrow -\left(\frac{\partial^2 A}{\partial V^2}\right)_{\mathbf{n}, T}\): second derivative w.r.t volume,
\(\left(\frac{\partial^2 A}{\partial T^2}\right)_{\mathbf{n}, V}\): 2nd partial derivative w.r.t temperature.
Since we compute second partial derivatives, the data type is now HyperDual64 (one real, 3 non-real values).
[7]:
state.joule_thomson()
data type : <class 'builtins.PyHyperDual64'>
temperature: 300 + [0]ε1 + [1]ε2 + [0]ε1ε2
volume : 40744 + [1]ε1 + [0]ε2 + [0]ε1ε2
moles : [1 + [0]ε1 + [0]ε2 + [0]ε1ε2]
density : 0.000024543491066169253 + [-0.00000000060238295371513]ε1 + [0]ε2 + [0]ε1ε2
A/kT : 5.060085461999783 + [0.00000040024326944247174]ε1 + [-0.025513116823522065]ε2 + [-0.0000000023037315630585307]ε1ε2
data type : <class 'builtins.PyDual2_64'>
temperature: 300 + [0]ε1 + [0]ε1²
volume : 40744 + [1]ε1 + [0]ε1²
moles : [1 + [0]ε1 + [0]ε1²]
density : 0.000024543491066169253 + [-0.00000000060238295371513]ε1 + [0.00000000000002956916128583988]ε1²
A/kT : 5.060085461999783 + [0.00000040024326944247174]ε1 + [-0.000000000019592452372041545]ε1²
data type : <class 'builtins.PyDual2_64'>
temperature: 300 + [1]ε1 + [0]ε1²
volume : 40744 + [0]ε1 + [0]ε1²
moles : [1 + [0]ε1 + [0]ε1²]
density : 0.000024543491066169253 + [0]ε1 + [0]ε1²
A/kT : 5.060085461999783 + [-0.025513116823522065]ε1 + [0.0001919137873859282]ε1²
[7]:
Summary¶
FeOsuses generalized dual numbers which enable the computation of higher order partial derivatives of the Helmholtz energy without the need of implementing these derivatives analytically.When a property is computed,
FeOscreates a new thermodynamic state with the correct data types according to the partial derivatives needed.If a property was computed before, it is pulled from the state’s cache instead of re-evaluating the Helmholtz energy.