{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3ec7988e",
   "metadata": {},
   "source": [
    "# Equilibrium calculations\n",
    "\n",
    "Different equilibrium calculations could be performed with the `yaeos` Excess\n",
    "Gibbs models.\n",
    "\n",
    "Let's instantiate an example model.\n",
    "\n",
    "## Example model\n",
    "As always, all the described methods are available for all Excess Gibbs models.\n",
    "But we are going to use an ethanol-hexane-water mixture modeled with UNIFAC\n",
    "as example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f23e7c6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaeos\n",
    "\n",
    "\n",
    "functional_groups = [\n",
    "    {1: 1, 2: 1, 14: 1},    # ethanol groups\n",
    "    {1: 2, 2: 4},           # n-hexane groups\n",
    "    {16: 1},                # water group\n",
    "]\n",
    "\n",
    "model = yaeos.UNIFACVLE(functional_groups)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c749a17",
   "metadata": {},
   "source": [
    "## Flash $\\gamma \\text{-} \\gamma$\n",
    "\n",
    "$$\n",
    "    y_{i} \\gamma_{i}^{y} = x_{i} \\gamma_{i}^{x}\n",
    "$$\n",
    "\n",
    "\n",
    "Liquid-liquid flash calculation with the $\\gamma \\text{-} \\gamma$ approach\n",
    "could be performed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7adf2ba3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'x': array([0.03080572, 0.96735068, 0.0018436 ]),\n",
       " 'y': array([0.31604605, 0.01086864, 0.67308532]),\n",
       " 'T': 303.15,\n",
       " 'beta': 0.5931639638170902}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "T = 30 + 273.15             # Temperature [K]\n",
    "n = np.array([5, 10, 10])   # Moles of each substance [mol]\n",
    "\n",
    "# IMPORTANT: flash calculations are performed specifying mole fractions\n",
    "z = n / n.sum()             # Mole fraction of eac substance [-]\n",
    "\n",
    "\n",
    "flash = model.flash_t(z, T)\n",
    "\n",
    "flash"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b8aca91",
   "metadata": {},
   "source": [
    "The result of the flash calculation is returned as a Python dictionary with the\n",
    "keys:\n",
    "\n",
    "- `x`: mole fractions of phase $x$\n",
    "- `y`: mole fractions of phase $y$\n",
    "- `T`: temperature in Kelvin that was specified\n",
    "- `beta`: mole fraction of phase $y$\n",
    "\n",
    "We can calculate the activities on each phase to check that are equal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a9e0b45a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.46710663 0.97593751 0.87010966]\n",
      "[0.46710662 0.9759375  0.87010961]\n"
     ]
    }
   ],
   "source": [
    "x = flash[\"x\"]\n",
    "y = flash[\"y\"]\n",
    "beta = flash[\"beta\"]\n",
    "\n",
    "ln_gammas_x = model.ln_gamma(x, T)\n",
    "ln_gammas_y = model.ln_gamma(y, T)\n",
    "\n",
    "gammas_x = np.exp(ln_gammas_x)\n",
    "gammas_y = np.exp(ln_gammas_y)\n",
    "\n",
    "activities_x = gammas_x * x\n",
    "activities_y = gammas_y * y\n",
    "\n",
    "print(activities_x)\n",
    "print(activities_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4769a61",
   "metadata": {},
   "source": [
    "As shown the activities in both phases are equals, the flash calculation was\n",
    "satisfactory.\n",
    "\n",
    "We can use the `beta` value to obtain the moles of each phase:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7614b8e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mixutre: ethanol-hexane-water\n",
      "Global mole vector:  [ 5 10 10]\n",
      "Phase y mole vector:  [4.68667812 0.16117207 9.9812489 ]\n",
      "Phase x mole vector:  [0.31332188 9.83882793 0.0187511 ]\n"
     ]
    }
   ],
   "source": [
    "nty = beta * n.sum()         # Total moles of phase y\n",
    "ntx = (1 - beta) * n.sum()   # total moles of phase x\n",
    "\n",
    "ny = y * nty  # moles vector of phase y\n",
    "nx = x * ntx  # moles vector of phase x\n",
    "\n",
    "print(\"Mixutre: ethanol-hexane-water\")\n",
    "print(\"Global mole vector: \", n)\n",
    "print(\"Phase y mole vector: \", ny)\n",
    "print(\"Phase x mole vector: \", nx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ff13629",
   "metadata": {},
   "source": [
    "There is an additional parameter to the `flash_t`, `k0`, which allows to\n",
    "provide an initialization to the flash calculation with:\n",
    "\n",
    "$$\n",
    "    k_{0} = \\frac{y}{x}\n",
    "$$\n",
    "\n",
    "The complete documentation of the `flash_t` can be accessed with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1c8577c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method flash_t in module yaeos.core:\n",
      "\n",
      "flash_t(z, temperature: float, k0=None) -> dict method of yaeos.models.excess_gibbs.unifac.UNIFACVLE instance\n",
      "    Two-phase split with specification of temperature and pressure.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    z : array_like\n",
      "        Global mole fractions\n",
      "    temperature : float\n",
      "        Temperature [K]\n",
      "    k0 : array_like, optional\n",
      "        Initial guess for the split, by default None (will use stability\n",
      "        analysis)\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    dict\n",
      "        Flash result dictionary with keys:\n",
      "            - x: \"phase 1\" mole fractions\n",
      "            - y: \"phase 2\" mole fractions\n",
      "            - T: temperature [K]\n",
      "            - beta: \"phase 2\" fraction\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(model.flash_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01edc10d",
   "metadata": {},
   "source": [
    "## Stability analysis\n",
    "\n",
    "It's possible to perform stability analysis by `tm` minimization proposed by\n",
    "Michelsen.\n",
    "\n",
    "Let's check with an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4b959248",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'w': array([2.07328795e-02, 1.47170347e-04, 9.79119950e-01]), 'tm': -0.5172683567506908}\n"
     ]
    }
   ],
   "source": [
    "T = 30 + 273.15             # Temperature [K]\n",
    "n = np.array([5, 10, 10])   # Moles vector [mol]\n",
    "\n",
    "z = n / n.sum()  # global mole fractions\n",
    "\n",
    "\n",
    "tm_min, tm_all_min = model.stability_analysis(z, T)\n",
    "\n",
    "print(tm_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f9d5148c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'tm': array([-0.51726836, -0.38462711, -0.51726836]), 'w': array([[2.07328795e-02, 1.47170347e-04, 9.79119950e-01],\n",
      "       [6.54275506e-03, 9.92256323e-01, 1.20092189e-03],\n",
      "       [2.07328782e-02, 1.47170324e-04, 9.79119951e-01]])}\n"
     ]
    }
   ],
   "source": [
    "print(tm_all_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebeff75d",
   "metadata": {},
   "source": [
    "Three minimization were performed, starting from near pure compositions of\n",
    "each compound.\n",
    "\n",
    "The first output `tm_min`, provide the value and composition of the lowest\n",
    "of the all found minimums.\n",
    "\n",
    "If the `tm` value of `tm_min` is negative, the phase of composition `z` is \n",
    "considered unstable. If you check the used composition, is the same used in the\n",
    "flash $\\gamma \\text{-} \\gamma$ tutorial. There, was shown that a flash\n",
    "calculation converges under the same conditions. The stability analysis in this\n",
    "case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d263a3ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.5172683567506908"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tm_min[\"tm\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68b8944f",
   "metadata": {},
   "source": [
    "has a negative value, confirming that the system was truly unstable.\n",
    "\n",
    "In the second output `tm_all_min` are listed all the found minimums. You can\n",
    "check that the `tm_min` is also listed there and repeated. Different\n",
    "initialization could converge to the same minimum.\n",
    "\n",
    "You can check in the minimums list that there is another minimum with value\n",
    "-0.38462711. We can use the composition of the two found minimums as\n",
    "initialization of a flash calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "65c4325b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'x': array([0.03080571, 0.96735068, 0.0018436 ]),\n",
       " 'y': array([0.31604605, 0.01086864, 0.67308532]),\n",
       " 'T': 303.15,\n",
       " 'beta': 0.5931639641540031}"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = tm_all_min[\"w\"][0]\n",
    "x = tm_all_min[\"w\"][1]\n",
    "\n",
    "k0 = y / x\n",
    "\n",
    "flash = model.flash_t(z, T, k0)\n",
    "\n",
    "flash"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "yaeos",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}