{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ed8d1cd8-0c5b-4f4e-8668-ef653075ec02",
   "metadata": {},
   "source": [
    "# Superancillaries \n",
    "\n",
    "## Motivation\n",
    "\n",
    "VLE calculations for pure fluids require to solve the system of equations \n",
    "\n",
    "$$\n",
    "p(T,\\rho') = p(T,\\rho'')\n",
    "$$\n",
    "$$\n",
    "g(T,\\rho') = g(T,\\rho'')\n",
    "$$\n",
    "which is a complicated non-linear rootfinding problem. For a specified $T$, one must get guess values for $\\rho'(T)$ and $\\rho''(T)$ which are commonly obtained from ancillary functions that give \"good enough\" estimates of the densities from which a nonlinear rootfinding algorithm launches to solve for the co-existing densities.\n",
    "\n",
    "These calculations, while not very slow (order of microseconds per call), often, if not usually, represent a computational pinchpoint. So it would be nice to be able to use a numerical function that can represent the results of these iterative calculations so well that the iterative calculation itself can be avoided. The superancillary functions described here satisfy that goal.\n",
    "\n",
    "## Theory and Approach\n",
    "\n",
    "The development of superancillary functions have been laid out in a series of publications:\n",
    "\n",
    "* [Exceptionally reliable density-solving algorithms for multiparameter mixture models from Chebyshev expansion rootfinding](https://scholar.google.com/citations?view_op=view_citation&hl=en&user=WNn0e_4AAAAJ&cstart=20&pagesize=80&sortby=pubdate&citation_for_view=WNn0e_4AAAAJ:hMod-77fHWUC)\n",
    "* [Efficient and Precise Representation of Pure Fluid Phase Equilibria with Chebyshev Expansions](https://scholar.google.com/citations?view_op=view_citation&hl=en&user=WNn0e_4AAAAJ&cstart=20&pagesize=80&sortby=pubdate&citation_for_view=WNn0e_4AAAAJ:OU6Ihb5iCvQC)\n",
    "* [Superancillary Equations for the Multiparameter Equations of State in REFPROP 10.0](https://scholar.google.com/citations?view_op=view_citation&hl=en&user=WNn0e_4AAAAJ&sortby=pubdate&citation_for_view=WNn0e_4AAAAJ:AXPGKjj_ei8C)\n",
    "\n",
    "The term \"superancillary\" was coined by Ulrich Deiters to differentiate it from the ancillary functions that are commonly provided alongside reference equations of state.\n",
    "\n",
    "At their core, a superancillary function is constructed from a series of Chebyshev expansions. When considering the entire set of Chebyshev expansions, they span the entire range of the independent variable. In their current form, they support only 1D function approximation. Chebyshev expansions, and orthogonal polynomials more generally, are a well-studied numerical tool for function approximation and can permit function approximation to the level of numerical precision.\n",
    "\n",
    "To build the superancillary equations, one does vapor-liquid equilibrium (VLE) calculations at carefully selected temperatures and *does some math* to construct the Chebyshev expansion. If the expansion is not accurate enough, the domain is subdivided into two halves and the process is then carried out in each section.\n",
    "\n",
    "To ensure highly accurate and reliable results in the very near vicinity of the critical point, as well as at very pressures (e.g., near the triple point of propane), it is necessary to do the phase equilibrium calculations in extended precision arithmetic. The ``boost::multiprecision`` library is used in C++ to do the extended precision calculations, in concert with the new ``teqp`` library for the equation of state part. The code used to do this exercise is in a fork at https://github.com/CoolProp/fastchebpure and [the releases](https://github.com/CoolProp/fastchebpure/releases) contain the obtained functions.\n",
    "\n",
    "### Caveats:\n",
    "* When superancillaries are enabled, the \"critical point\" is the numerical one that is obtained by enforcing $\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_T = \\left(\\frac{\\partial^2 p}{\\partial \\rho^2}\\right)_T = 0$, rather than the one reported by the EOS developers. This is because the superancillaries are developed to try to fix the critical region as well as possible. The differences can be sometimes not too small: [An Analysis of the Critical Region of Multiparameter Equations of State](https://scholar.google.com/citations?view_op=view_citation&hl=en&user=WNn0e_4AAAAJ&sortby=pubdate&citation_for_view=WNn0e_4AAAAJ:1qzjygNMrQYC)\n",
    "* When pressure is provided, the inverse superancillary for $T(p)$ is used, which introduces a very small error\n",
    "* When superancillaries are used, there is an overhead per `AbstractState` instantiation on the order of 100 $\\mu$s; if you cannot tolerate that overhead, either a) use the low-level interface, where that overhead can be amortized, or define the environment variable `COOLPROP_DISABLE_SUPERANCILLARIES_ENTIRELY` to turn off the parsing of superancillary data structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab07d156-01ad-43c2-902c-6f45ee76dfa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import timeit\n",
    "import CoolProp.CoolProp as CP\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import urllib\n",
    "from zipfile import ZipFile\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bad6e729-92f7-42d2-a53e-ff1c11263af1",
   "metadata": {},
   "outputs": [],
   "source": [
    "AS = CP.AbstractState('HEOS','Water')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d7a3783-b7fa-471b-902b-539fe3f8769a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# See caveat noted above. The use of superancillary functions necessarily changes the critical point for the fluid\n",
    "# Not usually by an amount that will be meaningful in practical applications\n",
    "for superanc in [False, True]:\n",
    "    CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, superanc)\n",
    "    print(superanc, AS.T_critical())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58913d10-19c7-4b47-858c-60c7e9a66ec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The JSON data for the expansions can be accessed from CoolProp\n",
    "jSuper = json.loads(CP.get_fluid_param_string(\"WATER\", \"JSON\"))[0][\"EOS\"][0][\"SUPERANCILLARY\"]\n",
    "SA = CP.SuperAncillary(json.dumps(jSuper))\n",
    "\n",
    "AS = CP.AbstractState('HEOS','Water')\n",
    "Tt = AS.Ttriple()\n",
    "Tc = AS.T_critical()\n",
    "\n",
    "Ts = np.linspace(Tt, 647.0959999999873, 5*10**5)\n",
    "ps = np.zeros_like(Ts)\n",
    "SA.eval_sat_many(Ts, 'P', 0, ps)\n",
    "plt.plot(1/Ts, ps)\n",
    "plt.yscale('log'); plt.gca().set(xlabel='$1/T$ / 1/K', ylabel='$p$ / Pa');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418b1527-9d80-44ca-9461-5cf1cb5ad35d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# At the lower level, calling with a buffer of points\n",
    "tic = timeit.default_timer()\n",
    "SA.eval_sat_many(Ts, 'P', 0, ps)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(Ts)*1e6, 'μs/call')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54f6ad37-0813-4f8a-b7ba-8752b550587e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# With superancillaries disabled\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n",
    "QT_INPUTS = CP.QT_INPUTS\n",
    "tic = timeit.default_timer()\n",
    "for T_ in Ts:\n",
    "    AS.update(QT_INPUTS, 0, T_)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(Ts)*1e6, 'μs/call')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d254da7-0080-4549-9222-0f0b276fd6d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# With superancillaries enabled, three superancillary functions are evaluated\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n",
    "QT_INPUTS = CP.QT_INPUTS\n",
    "tic = timeit.default_timer()\n",
    "for T_ in Ts:\n",
    "    AS.update(QT_INPUTS, 0, T_)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(Ts)*1e6, 'μs/call')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8731b27-90f9-413a-8cfc-363d6feea0f1",
   "metadata": {},
   "source": [
    "If pressure is specified, the speedup is even greater because one must also iterate for the pressure of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72969532-5b98-4c44-adc9-204d16127e9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pssmall = np.geomspace(ps[0], ps[-1]*(1-1e-10), 10**4) # The full equilibrium is slow\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n",
    "PQ_INPUTS = CP.PQ_INPUTS\n",
    "tic = timeit.default_timer()\n",
    "for p_ in pssmall:\n",
    "    AS.update(PQ_INPUTS, p_, 0)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(pssmall)*1e6, 'μs/call')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960463cb-9574-4ab7-b3be-23a932470d30",
   "metadata": {},
   "outputs": [],
   "source": [
    "# With superancillaries enabled, one evaluates T(p) from an inverse function and then uses that T\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n",
    "PQ_INPUTS = CP.PQ_INPUTS\n",
    "tic = timeit.default_timer()\n",
    "for p_ in pssmall:\n",
    "    AS.update(PQ_INPUTS, p_, 0)\n",
    "toc = timeit.default_timer()\n",
    "print((toc-tic)/len(pssmall)*1e6, 'μs/call')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56dee520-1a03-45dc-bd2e-c8343da27e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cecc435c-89ad-47f5-a368-3daba910a294",
   "metadata": {},
   "source": [
    "## Other validation checks to confirm accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7635fe7-94f7-4415-bcb3-403c5efb5d42",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_hfg(T):\n",
    "    \"\"\" Latent heat \"\"\"\n",
    "    AS.update(CP.QT_INPUTS, 0, T)\n",
    "    return AS.saturated_vapor_keyed_output(CP.iHmolar)-AS.saturated_liquid_keyed_output(CP.iHmolar)\n",
    "    \n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n",
    "HFG_SA = np.array([get_hfg(T_) for T_ in Ts])\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n",
    "HFG_NON = np.array([get_hfg(T_) for T_ in Ts])\n",
    "\n",
    "plt.plot(Ts, np.abs(HFG_NON/HFG_SA-1))\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n",
    "plt.yscale('log')\n",
    "plt.ylabel(r'$|\\Delta h_{fg, non}/\\Delta h_{fg, sa} - 1|$')\n",
    "plt.xlabel('$T$ / K');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7eb7d11-f782-46b4-a245-07d05aa038f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sf(T):\n",
    "    \"\"\" liquid entropy \"\"\"\n",
    "    AS.update(CP.QT_INPUTS, 0, T)\n",
    "    return AS.saturated_liquid_keyed_output(CP.iSmolar)\n",
    "    \n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n",
    "SF_SA = np.array([get_sf(T_) for T_ in Ts])\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n",
    "SF_NON = np.array([get_sf(T_) for T_ in Ts])\n",
    "\n",
    "plt.plot(Ts, np.abs(SF_NON - SF_SA))\n",
    "plt.yscale('log')\n",
    "plt.ylabel(r\"$|s'_{non} - s'_{sa}|$ / J/mol\")\n",
    "plt.xlabel('$T$ / K');\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(Ts, SF_NON)\n",
    "plt.title(r'Note value goes towards zero at $T\\to 273.16$ K')\n",
    "plt.gca().set(ylabel=r\"$s'_{non}$ / J/mol\", xlabel='$T$ / K')\n",
    "\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f05b9de1-27d7-42ce-a35c-08cbb5b99d25",
   "metadata": {},
   "source": [
    "## Checks against extended precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933cb2e5-d31e-4fac-a1d9-b4aa82e079a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "outputversion = '2025.04.27'\n",
    "if not Path(f'{outputversion}.zip').exists():\n",
    "    print('Downloading the chebyshev output file to ', Path('.').absolute())\n",
    "    urllib.request.urlretrieve(f'https://github.com/CoolProp/fastchebpure/archive/refs/tags/{outputversion}.zip', f'{outputversion}.zip')\n",
    "\n",
    "    with ZipFile(f'{outputversion}.zip') as z:\n",
    "        print(Path(f'{outputversion}.zip').parent.absolute())\n",
    "        z.extractall(Path(f'{outputversion}.zip').parent.absolute())\n",
    "\n",
    "    with (Path('.') / f'fastchebpure-{outputversion}' / '.gitignore').open('w') as fp:\n",
    "        fp.write(\"*\")\n",
    "\n",
    "outputcheck = (Path('.') / f'fastchebpure-{outputversion}' / 'outputcheck').absolute()\n",
    "assert outputcheck.exists()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e26150cd-0e6a-477a-90c9-4aa11e076d7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import CoolProp.CoolProp as CP\n",
    "\n",
    "fluid = 'Water'\n",
    "\n",
    "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n",
    "\n",
    "AS = CP.AbstractState('HEOS', f\"{fluid}\")\n",
    "jSuper = json.loads(CP.get_fluid_param_string(f\"{fluid}\", \"JSON\"))[0]['EOS'][0]['SUPERANCILLARY']\n",
    "superanc = CP.SuperAncillary(json.dumps(jSuper))\n",
    "RPname = AS.fluid_param_string(\"REFPROP_name\")\n",
    "\n",
    "# Load extended precision calcs from the release on github\n",
    "chk = json.load(open(f'{outputcheck}/{fluid}_check.json'))\n",
    "df = pd.DataFrame(chk['data'])\n",
    "# df.info() # uncomment to see what fields are available\n",
    "\n",
    "Tcrit_num = AS.get_fluid_parameter_double(0, \"SUPERANC::Tcrit_num\")\n",
    "T = df['T / K'].to_numpy()\n",
    "Theta = (Tcrit_num-T)/Tcrit_num\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6,10))\n",
    "\n",
    "\n",
    "plt.sca(axes[0])\n",
    "rhoL_anc = np.zeros_like(T)\n",
    "superanc.eval_sat_many(T, 'D', 0, rhoL_anc)\n",
    "err = np.abs(df[\"rho'(mp) / mol/m^3\"]/rhoL_anc-1)\n",
    "plt.plot(Theta, err, 'o', label=r'$\\Upsilon$:SA')\n",
    "\n",
    "errCP = np.abs(df[\"rho'(mp) / mol/m^3\"]/CP.PropsSI('Dmolar', 'T', T, 'Q', 0, f'HEOS::{fluid}')-1)\n",
    "plt.plot(Theta, errCP, '^', label=r'$\\Upsilon$:HEOS')\n",
    "\n",
    "try:\n",
    "    errRP = np.abs(df[\"rho'(mp) / mol/m^3\"]/CP.PropsSI('Dmolar', 'T', T, 'Q', 0, f'REFPROP::{RPname}')-1)\n",
    "    plt.plot(Theta, errRP, 'x', label=r'$\\Upsilon$:REFPROP')\n",
    "except BaseException as BE:\n",
    "    print(BE)\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.ylabel(r\"$|\\rho_{{\\rm \\Upsilon}}'/\\rho_{{\\rm ep}}'-1|$\")\n",
    "plt.yscale('log')\n",
    "\n",
    "\n",
    "\n",
    "plt.sca(axes[1])\n",
    "rhoV_anc = np.zeros_like(T)\n",
    "superanc.eval_sat_many(T, 'D', 1, rhoV_anc)\n",
    "err = np.abs(df[\"rho''(mp) / mol/m^3\"]/rhoV_anc-1)\n",
    "plt.plot(Theta, err, 'o', label=r'$\\Upsilon$:SA')\n",
    "\n",
    "errCP = np.abs(df[\"rho''(mp) / mol/m^3\"]/CP.PropsSI('Dmolar', 'T', T, 'Q', 1, f'HEOS::{fluid}')-1)\n",
    "plt.plot(Theta, errCP, '^', label=r'$\\Upsilon$:HEOS')\n",
    "\n",
    "try:\n",
    "    errRP = np.abs(df[\"rho''(mp) / mol/m^3\"]/CP.PropsSI('Dmolar', 'T', T, 'Q', 1, f'REFPROP::{RPname}')-1)\n",
    "    plt.plot(Theta, errRP, 'x', label=r'$\\Upsilon$:REFPROP')\n",
    "except BaseException as BE:\n",
    "    print(BE)\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.ylabel(r\"$|\\rho_{{\\rm \\Upsilon}}''/\\rho_{{\\rm ep}}''-1|$\")\n",
    "plt.yscale('log')\n",
    "\n",
    "\n",
    "\n",
    "plt.sca(axes[2])\n",
    "p_anc = np.zeros_like(T)\n",
    "superanc.eval_sat_many(T, 'P', 1, p_anc)\n",
    "err = np.abs(df[\"p(mp) / Pa\"]/p_anc-1)\n",
    "plt.plot(Theta, err, 'o', label=r'$\\Upsilon$:SA')\n",
    "\n",
    "errCP = np.abs(df[\"p(mp) / Pa\"]/CP.PropsSI('P', 'T', T, 'Q', 1, f'HEOS::{fluid}')-1)\n",
    "plt.plot(Theta, errCP, '^', label=r'$\\Upsilon$:HEOS')\n",
    "\n",
    "try:\n",
    "    errRP = np.abs(df[\"p(mp) / Pa\"]/CP.PropsSI('P', 'T', T, 'Q', 1, f'REFPROP::{RPname}')-1)\n",
    "    plt.plot(Theta, errRP, 'x', label=r'$\\Upsilon$:REFPROP')\n",
    "except BaseException as BE:\n",
    "    print(BE)\n",
    "plt.legend(loc='best')\n",
    "\n",
    "# print(CP.PropsSI('gas_constant', 'T', T[0], 'Q', 1, f'HEOS::{fluid}'))\n",
    "# print(CP.PropsSI('gas_constant', 'T', T[0], 'Q', 1, f'REFPROP::{fluid}'))\n",
    "\n",
    "plt.ylabel(r\"$|p_{{\\rm \\Upsilon}}/p_{{\\rm ep}}-1|$\")\n",
    "plt.yscale('log')\n",
    "\n",
    "plt.sca(axes[2])\n",
    "plt.xlabel(r'$\\Theta=(T_{{\\rm crit,num}}-T)/T_{{\\rm crit,num}}$')\n",
    "plt.xscale('log')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.axhline(1e-12, dashes=[2,2])\n",
    "    ax.axvline(1e-6, dashes=[2,2])\n",
    "\n",
    "plt.tight_layout(pad=0.2)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
