From a56402df4eb7444676fc6163bc02e047bfef8bc8 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 4 Nov 2022 18:03:48 -0400
Subject: [PATCH] Add I&ECR examples to docs
---
doc/source/examples/IECR examples.ipynb | 198 ++++++++++++++++++++++++
doc/source/examples/index.rst | 8 +
doc/source/index.rst | 1 +
3 files changed, 207 insertions(+)
create mode 100644 doc/source/examples/IECR examples.ipynb
create mode 100644 doc/source/examples/index.rst
diff --git a/doc/source/examples/IECR examples.ipynb b/doc/source/examples/IECR examples.ipynb
new file mode 100644
index 0000000..65a6de2
--- /dev/null
+++ b/doc/source/examples/IECR examples.ipynb
@@ -0,0 +1,198 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "1e108f86",
+ "metadata": {},
+ "source": [
+ "# The teqp paper in I&ECR \n",
+ "\n",
+ "A few minor changes have been made:\n",
+ " \n",
+ "* The ``get_splus`` method requires the molar concentrations to be a numpy array (to avoid copies) (as of version 0.14.0\n",
+ "* The top-level methods ``teqp.xxx`` have been deprecated, and the methods attached to the instance are preferred"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9539423f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import timeit, numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "plt.style.use('classic')\n",
+ "import teqp\n",
+ "\n",
+ "def build_models():\n",
+ " Tc_K, pc_Pa, acentric = 647.096, 22064000.0, 0.3442920843\n",
+ "\n",
+ " water = {\n",
+ " \"a0i / Pa m^6/mol^2\": 0.12277 , \"bi / m^3/mol\": 0.000014515, \"c1\": 0.67359, \n",
+ " \"Tc / K\": 647.096, \"epsABi / J/mol\": 16655.0, \"betaABi\": 0.0692, \"class\": \"4C\"\n",
+ " }\n",
+ " j = {\"cubic\": \"SRK\", \"pures\": [water], \"R_gas / J/mol/K\": 8.3144598}\n",
+ "\n",
+ " datapath = teqp.get_datapath()\n",
+ " def get_PCSAFT():\n",
+ " c = teqp.SAFTCoeffs()\n",
+ " # Values from https://doi.org/10.1016/j.fluid.2017.11.015, \n",
+ " # but association contribution is ignored\n",
+ " c.name = 'Water'\n",
+ " c.m = 2.5472\n",
+ " c.sigma_Angstrom = 2.1054\n",
+ " c.epsilon_over_k = 138.63\n",
+ " return teqp.PCSAFTEOS(coeffs=[c])\n",
+ " \n",
+ " return [\n",
+ " ('vdW', teqp.vdWEOS([Tc_K], [pc_Pa])),\n",
+ " ('PR', teqp.canonical_PR([Tc_K], [pc_Pa], [acentric])),\n",
+ " ('SRK', teqp.canonical_SRK([Tc_K], [pc_Pa], [acentric])),\n",
+ " ('PCSAFT', get_PCSAFT()),\n",
+ " ('CPA', teqp.CPAfactory(j)),\n",
+ " ('IAPWS', teqp.build_multifluid_model([\"Water\"], datapath))\n",
+ " ]\n",
+ "\n",
+ "fig, ax = plt.subplots(1,1,figsize=(7,6))\n",
+ "T = 700 # K\n",
+ "rhovec = np.geomspace(0.1, 30e3, 10000) # mol/m^3; critical density is 17873.8... mol/m^3\n",
+ "tic = timeit.default_timer()\n",
+ "for abbrv, model in build_models():\n",
+ " splus = np.array([model.get_splus(T, np.array([rho])) for rho in rhovec])\n",
+ " plt.plot(rhovec, splus, label=abbrv, lw = 1.5 if abbrv=='IAPWS' else 1)\n",
+ "elap = timeit.default_timer()-tic\n",
+ "plt.axvline(17873.8, dashes=[2,2])\n",
+ "plt.legend(loc='best')\n",
+ "plt.gca().set(xlabel=r'$\\rho$ / mol/m$^3$', ylabel=r'$s^+\\equiv (s^{\\rm ig}(T,\\rho)-s(T,\\rho))/R$')\n",
+ "plt.tight_layout(pad=0.2)\n",
+ "plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
+ "plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7)\n",
+ "plt.savefig('splus_water_700K.pdf')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fb28d259",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import json, timeit\n",
+ "import pandas, numpy as np, matplotlib.pyplot as plt\n",
+ "plt.style.use('classic')\n",
+ "import teqp\n",
+ "\n",
+ "Tc_K = [190.564, 154.581]\n",
+ "pc_Pa = [4599200, 5042800]\n",
+ "acentric = [0.011, 0.022]\n",
+ "model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
+ "fig, ax = plt.subplots(1,1,figsize=(7, 6), subplot_kw=dict(projection='3d'))\n",
+ "tic = timeit.default_timer()\n",
+ "for ifluid in [0,1]:\n",
+ " model0 = teqp.canonical_PR([Tc_K[ifluid]], [pc_Pa[ifluid]], [acentric[ifluid]])\n",
+ " for T in np.linspace(190, 120, 50):\n",
+ " if T > Tc_K[ifluid]: continue\n",
+ " [rhoL, rhoV] = model0.superanc_rhoLV(T)\n",
+ " rhovecL = np.array([0.0, 0.0]); rhovecL[ifluid] = rhoL\n",
+ " rhovecV = np.array([0.0, 0.0]); rhovecV[ifluid] = rhoV\n",
+ " opt = teqp.TVLEOptions(); opt.calc_criticality = True\n",
+ " df = pandas.DataFrame(model.trace_VLE_isotherm_binary(T, rhovecL, rhovecV, opt))\n",
+ " df['too_critical'] = df.apply(\n",
+ " lambda row: (abs(row['crit. conditions L'][0]) < 5e-8), axis=1)\n",
+ " first_too_critical = np.argmax(df['too_critical'])\n",
+ " df = df.iloc[0:(first_too_critical if first_too_critical else len(df))]\n",
+ " line, = ax.plot(xs=df['T / K'], ys=df['xL_0 / mole frac.'], zs=df['pL / Pa']/1e6, \n",
+ " lw=0.2, color='k')\n",
+ " ax.plot(xs=df['T / K'], ys=df['xV_0 / mole frac.'], zs=df['pL / Pa']/1e6, \n",
+ " dashes=[2,2], color=line.get_color(), lw=0.2)\n",
+ "elap = timeit.default_timer()-tic\n",
+ "ax.view_init(elev=10., azim=130)\n",
+ "ax.set(xlabel='$T$ / K', ylabel='$x_1$ / mole frac.', zlabel='$p$ / MPa')\n",
+ "fig.text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
+ "fig.text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7)\n",
+ "plt.tight_layout(pad=0.2)\n",
+ "plt.savefig('PR_VLE_trace.pdf')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cbe191dd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import timeit\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt \n",
+ "plt.style.use('classic')\n",
+ "import pandas\n",
+ "import teqp\n",
+ "\n",
+ "def get_critical_curve(ipure):\n",
+ " \"\"\" Return curve as pandas DataFrame \"\"\"\n",
+ " names = ['Nitrogen', 'Ethane']\n",
+ " model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
+ " T0 = model.get_Tcvec()[ipure]\n",
+ " rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)\n",
+ " rho0[1-ipure] = 0\n",
+ " o = teqp.TCABOptions() \n",
+ " o.init_dt = 1.0 # step in the parameter\n",
+ " o.rel_err = 1e-8\n",
+ " o.abs_err = 1e-5\n",
+ " o.integration_order = 5\n",
+ " o.calc_stability = True\n",
+ " o.polish = True\n",
+ " curveJSON = model.trace_critical_arclength_binary(T0, rho0, '', o)\n",
+ " df = pandas.DataFrame(curveJSON)\n",
+ " rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
+ " df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot\n",
+ " return df\n",
+ "\n",
+ "if __name__ == '__main__':\n",
+ " fig, ax = plt.subplots(1,1,figsize=(7, 6))\n",
+ " tic = timeit.default_timer()\n",
+ " for ipure in [1,0]:\n",
+ " df = get_critical_curve(ipure)\n",
+ " first_unstable = np.argmax(~df['locally stable'])\n",
+ " df = df.iloc[0:(first_unstable if first_unstable else len(df))]\n",
+ " line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
+ " plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', \n",
+ " color=line.get_color())\n",
+ "\n",
+ " elap = timeit.default_timer()-tic\n",
+ " plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa',\n",
+ " xlim=(100, 350), ylim=(1, 1e3))\n",
+ " plt.yscale('log')\n",
+ " plt.tight_layout(pad=0.2)\n",
+ " plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
+ " plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7)\n",
+ " plt.savefig('N2_ethane_critical.pdf')\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.10.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/examples/index.rst b/doc/source/examples/index.rst
new file mode 100644
index 0000000..0be3570
--- /dev/null
+++ b/doc/source/examples/index.rst
@@ -0,0 +1,8 @@
+Examples
+==========
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Contents:
+
+ IECR_examples
\ No newline at end of file
diff --git a/doc/source/index.rst b/doc/source/index.rst
index c537fab..24647d1 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -10,6 +10,7 @@ Welcome to teqp's documentation!
derivs/index
algorithms/index
api/modules
+ examples/index
Indices and tables
==================
--
GitLab