{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# On Dual Numbers in `FeOs`\n", "\n", "In `FeOs`, we use [generalized dual numbers](https://www.frontiersin.org/articles/10.3389/fceng.2021.758090/full) to compute partial derivatives of the Helmholtz energy.\n", "In this notebook, we take a closer look at how that works.\n", "\n", "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. \n", "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.\n", "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)\n", "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.\n", "\n", "## Dual Numbers\n", "\n", "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.\n", "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.\n", "We can - in principle - create numbers that allow computation of an arbitrarily high derivative at the cost of execution speed.\n", "\n", "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.\n", "\n", "For example, when we feed a temperature as `Dual64` data type\n", "\n", "```python\n", "print(temperature)\n", "300 + [1]ε\n", "```\n", "\n", "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.\n", "\n", "## Do we need to create dual numbers by hand?\n", "\n", "No. If you use already implemented equations of state, you will not even recognize that dual numbers are used.\n", "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.\n", "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." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from feos.eos import *\n", "from feos.si import *\n", "import numpy as np\n", " \n", "SQRT2 = np.sqrt(2)\n", "\n", "class PyPengRobinson: \n", " def __init__(self, critical_temperature, critical_pressure, acentric_factor, molar_weight, delta_ij=None):\n", " \"\"\"Peng-Robinson Equation of State\n", " \n", " Parameters\n", " ----------\n", " critical_temperature : SIArray1\n", " critical temperature of each component.\n", " critical_pressure : SIArray1\n", " critical pressure of each component.\n", " acentric_factor : np.array[float] \n", " acentric factor of each component (dimensionless).\n", " molar_weight: SIArray1\n", " molar weight of each component.\n", " delta_ij : np.array[[float]], optional\n", " binary parameters. Shape=[n, n], n = number of components.\n", " defaults to zero for all binary interactions.\n", " \n", " Raises\n", " ------\n", " ValueError: if the input values have incompatible sizes.\n", " \"\"\"\n", " self.n = len(critical_temperature)\n", " if len(set((len(critical_temperature), len(critical_pressure), len(acentric_factor)))) != 1:\n", " raise ValueError(\"Input parameters must all have the same lenght.\")\n", " \n", " if self.n == 1:\n", " self.tc = (critical_temperature / KELVIN)[0]\n", " self.pc = (critical_pressure / PASCAL)[0]\n", " self.omega = acentric_factor[0]\n", " self.mw = (molar_weight / GRAM * MOL)[0]\n", " else:\n", " self.tc = critical_temperature / KELVIN\n", " self.pc = critical_pressure / PASCAL\n", " self.omega = acentric_factor\n", " self.mw = molar_weight / GRAM * MOL\n", " \n", " self.a_r = 0.45724 * critical_temperature**2 * RGAS / critical_pressure / ANGSTROM**3 / NAV / KELVIN\n", " self.b = 0.07780 * critical_temperature * RGAS / critical_pressure / ANGSTROM**3 / NAV\n", " self.kappa = 0.37464 + (1.54226 - 0.26992 * acentric_factor) * acentric_factor\n", " self.delta_ij = np.zeros((self.n, self.n)) if delta_ij is None else delta_ij\n", " \n", " def helmholtz_energy(self, state):\n", " \"\"\"Return Helmholtz energy.\n", " \n", " Parameters\n", " ----------\n", " state : StateHD\n", " The thermodynamic state.\n", " \n", " Returns\n", " -------\n", " helmholtz_energy: float | any dual number\n", " The return type depends on the input types.\n", " \"\"\"\n", " n = np.sum(state.moles)\n", " x = state.molefracs\n", " tr = 1.0 / self.tc * state.temperature\n", " ak = ((1.0 - np.sqrt(tr)) * self.kappa + 1.0)**2 * self.a_r\n", " ak_mix = 0.0\n", " if self.n > 1:\n", " for i in range(self.n):\n", " for j in range(self.n):\n", " ak_mix += np.sqrt(ak[i] * ak[j]) * (x[i] * x[j] * (1.0 - self.delta_ij[i, j]))\n", " else:\n", " ak_mix = ak\n", " b = np.sum(x * self.b)\n", " v = 1.0 / state.density\n", " a = n * (np.log(v / (v - b)) - ak_mix / (b * SQRT2 * 2.0 * state.temperature)\n", " * np.log((v * (SQRT2 - 1.0) + b) / (v * (SQRT2 + 1.0) - b)))\n", " \n", " # some print statements to inspect data types\n", " print()\n", " print(\"data type : \", type(state.temperature))\n", " print(\"temperature: \", state.temperature)\n", " print(\"volume : \", state.volume)\n", " print(\"moles : \", state.moles)\n", " print(\"density : \", state.density)\n", " print(\"A/kT : \", a)\n", " return a\n", " \n", " def components(self) -> int: \n", " \"\"\"Number of components.\"\"\"\n", " return self.n\n", " \n", " def subset(self, i: [int]):\n", " \"\"\"Return new equation of state containing a subset of all components.\"\"\"\n", " if self.n > 1:\n", " tc = self.tc[i] \n", " pc = self.pc[i]\n", " mw = self.mw[i]\n", " omega = self.omega[i]\n", " return PyPengRobinson(tc*KELVIN, pc*PASCAL, omega, mw*GRAM/MOL)\n", " else:\n", " return self\n", " \n", " def molar_weight(self) -> SIArray1:\n", " if isinstance(self.mw, float):\n", " return np.array([self.mw]) * GRAM / MOL\n", " else:\n", " return self.mw * GRAM / MOL\n", " \n", " def max_density(self, moles:[float]) -> float:\n", " b = np.sum(moles * self.b) / np.sum(moles);\n", " return 0.9 / b " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# parameters for propane\n", "tc = SIArray1(369.96 * KELVIN)\n", "pc = SIArray1(4250000.0 * PASCAL)\n", "omega = np.array([0.153])\n", "molar_weight = SIArray1(44.0962 * GRAM / MOL)\n", "\n", "# create an instance of our python class and hand it over to rust\n", "eos = EquationOfState.python(PyPengRobinson(tc, pc, omega, molar_weight))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Properties and dual numbers\n", "\n", "After initializing the equation of state for a single substance (in the above cell), let us now build a thermodynamic state.\n", "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "state = State(eos, temperature=300*KELVIN, volume=40744*ANGSTROM**3, total_moles=1/NAV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "Which dual numbers (dual, hyper-dual, etc.) are used depends on the property to compute or rather which partial derivatives are needed.\n", "`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.\n", "The result will have the same data type as the state's variables and the needed derivative(s) are extracted and returned.\n", "\n", "`FeOs` does all of that independent from the implementation of the equation of state.\n", "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.\n", "In Python, your can write code just like \"regular\" Python because dual numbers implement all arithmetic opartions you'd typically use." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at an example. First, we compute the `molar_helmholtz_energy`.\n", "Since *no derivatives* are needed to compute the Helmholtz energy, the data types of the `State`'s properties are regular floating point numbers (`float`).\n", "\n", "If you are in a running noteboook, execute the cell below and check the output.\n", "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.\n", "That's because we *cache* already computed derivatives for a given state.\n", "A second call to `molar_helmholtz_energy` simply pulls the already computed value from cache and returns it without entering the `helmholtz_energy` function." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "data type : \n", "temperature: 300.0\n", "volume : 40744.0\n", "moles : [1.0]\n", "density : 2.4543491066169253e-05\n", "A/kT : [5.06008546]\n" ] }, { "data": { "text/plain": [ "-6.554978406564657" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state.helmholtz_energy() / (KB * 300*KELVIN)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we compute the `pressure`, for which we compute the first partial derivative of the Helmholtz energy with respect to the volume.\n", "\n", "$p = -\\left(\\frac{\\partial A}{\\partial V}\\right)_{\\mathbf{n}, T}$\n", "\n", "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.\n", "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.\n", "\n", "`FeOs` modifies the dual parts of the temperature, volume or amount of substance (for a component) when computing derivatives.\n", "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.\n", "It is also different from unity - the value is calculated according to the artihmetic operation of division between two dual numbers." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "data type : \n", "temperature: 300 + [0]ε\n", "volume : 40744 + [1]ε\n", "moles : [1 + [0]ε]\n", "density : 0.000024543491066169253 + [-0.00000000060238295371513]ε\n", "A/kT : 5.060085461999783 + [0.00000040024326944247174]ε\n" ] }, { "data": { "text/latex": [ "$100\\,\\mathrm{kPa}$" ], "text/plain": [ "100.00005278190909 kPa" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state.pressure()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similar to the `pressure`, the `entropy` is computed by taking the first partial derivative with respect to temperature.\n", "Here, the dual part of the density is zero because both dual parts of the volume and the amount of substance are zero, too." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "data type : \n", "temperature: 300 + [1]ε\n", "volume : 40744 + [0]ε\n", "moles : [1 + [0]ε]\n", "density : 0.000024543491066169253 + [0]ε\n", "A/kT : 5.060085461999783 + [-0.025513116823522065]ε\n" ] }, { "data": { "text/latex": [ "$118.14\\,\\mathrm{\\frac{ J}{molK}}$" ], "text/plain": [ "118.13947975470876 J/mol/K" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state.molar_entropy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at a more involved property.\n", "For the Joule-Thomson coefficient, we need multiple second partial derivatives:\n", "\n", "$\\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)$\n", "\n", "We need three evaluations of the Helmholtz energy (that are actually occuring in the computation of $C_p$):\n", "\n", "1. $\\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,\n", "2. $\\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,\n", "3. $\\left(\\frac{\\partial^2 A}{\\partial T^2}\\right)_{\\mathbf{n}, V}$: 2nd partial derivative w.r.t temperature.\n", "\n", "Since we compute second partial derivatives, the data type is now `HyperDual64` (one real, 3 non-real values)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "data type : \n", "temperature: 300 + [0]ε1 + [1]ε2 + [0]ε1ε2\n", "volume : 40744 + [1]ε1 + [0]ε2 + [0]ε1ε2\n", "moles : [1 + [0]ε1 + [0]ε2 + [0]ε1ε2]\n", "density : 0.000024543491066169253 + [-0.00000000060238295371513]ε1 + [0]ε2 + [0]ε1ε2\n", "A/kT : 5.060085461999783 + [0.00000040024326944247174]ε1 + [-0.025513116823522065]ε2 + [-0.0000000023037315630585307]ε1ε2\n", "\n", "data type : \n", "temperature: 300 + [0]ε1 + [0]ε1²\n", "volume : 40744 + [1]ε1 + [0]ε1²\n", "moles : [1 + [0]ε1 + [0]ε1²]\n", "density : 0.000024543491066169253 + [-0.00000000060238295371513]ε1 + [0.00000000000002956916128583988]ε1²\n", "A/kT : 5.060085461999783 + [0.00000040024326944247174]ε1 + [-0.000000000019592452372041545]ε1²\n", "\n", "data type : \n", "temperature: 300 + [1]ε1 + [0]ε1²\n", "volume : 40744 + [0]ε1 + [0]ε1²\n", "moles : [1 + [0]ε1 + [0]ε1²]\n", "density : 0.000024543491066169253 + [0]ε1 + [0]ε1²\n", "A/kT : 5.060085461999783 + [-0.025513116823522065]ε1 + [0.0001919137873859282]ε1²\n" ] }, { "data": { "text/latex": [ "$-1.4939\\times10^{-4}\\,\\mathrm{\\frac{ms^{2}K}{kg}}$" ], "text/plain": [ "-1.4938695195180555e-4 m kg^-1 s^2 K" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state.joule_thomson()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "- `FeOs` uses generalized dual numbers which enable the computation of higher order partial derivatives of the Helmholtz energy without the need of implementing these derivatives analytically.\n", "- When a property is computed, `FeOs` creates a new thermodynamic state with the correct data types according to the partial derivatives needed.\n", "- If a property was computed before, it is pulled from the state's cache instead of re-evaluating the Helmholtz energy." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 4 }