examples/numerical example/Fig4 exact vs approximate medium entropy production.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%reset -f\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy as sp\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plot_settings as ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PySTFP import PySTFP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load analytical results, which have been calculated using the notebook\n",
    "# analytical solution.ipynb\n",
    "\n",
    "analytical_solution = {}\n",
    "\n",
    "with h5py.File('analytical_solution.h5','r') as f:\n",
    "    #\n",
    "    for key in f.keys():\n",
    "        analytical_solution[key] = f[key][()]\n",
    "\n",
    "#analytical_dictionary.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xP = analytical_solution['xP']\n",
    "P = analytical_solution['P']\n",
    "\n",
    "a_xP = analytical_solution['a_xP']\n",
    "D_xP = analytical_solution['D_xP']\n",
    "\n",
    "a_x0 = analytical_solution['a_x0']\n",
    "D_x0 = analytical_solution['D_x0']\n",
    "\n",
    "dt = analytical_solution['dt']\n",
    "x0 = analytical_solution['x0']\n",
    "\n",
    "\n",
    "j = a_xP[0] * P[0] \\\n",
    "            - D_xP[1] * P[0] \\\n",
    "            - D_xP[0] * P[1]\n",
    "\n",
    "Sm_dot_integrand = j/D_xP[0] * ( a_xP[0] - D_xP[1] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sm_dot = np.zeros_like(dt,dtype=float)\n",
    "\n",
    "for i,current_integrand in enumerate(Sm_dot_integrand):\n",
    "    Sm_dot[i] = np.trapz(current_integrand,xP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sm_dot_perturbative = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate medium entropy production rate from Gaussian propagator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for the Gaussian propagator\n",
    "\n",
    "\\begin{align}\n",
    "P_{\\mathrm{GP}}(\\Delta x)\n",
    "&=\n",
    "\\frac{1}{\\sqrt{4 \\pi D(x_0) \\Delta t}}\n",
    "\\exp\\left[\n",
    "- \\frac{ \\Delta t}{4D(x_0)} \n",
    "\\left( \\frac{ \\Delta x}{\\Delta t} - a(x_0)\\right)^2\n",
    "    \\right]\n",
    "\\end{align}\n",
    "\n",
    "it holds that \n",
    "\n",
    "\\begin{align}\n",
    "\\partial_x P_{\\mathrm{GP}}(\\Delta x)\n",
    "&=\n",
    "-2\n",
    " \\frac{ 1}{4D(x_0)} \n",
    "  \\left( \\frac{ \\Delta x}{\\Delta t} - a(x_0)\\right) P_{\\mathrm{GP}}(\\Delta x)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sm_dot_Gaussian = np.zeros_like(dt,dtype=float)\n",
    "\n",
    "for i,current_dt in enumerate(dt):\n",
    "        #\n",
    "        P_Gaussian = 1/np.sqrt(4*np.pi*D_x0[0]*current_dt)\n",
    "        P_Gaussian *= np.exp( - current_dt/(4*D_x0[0]) \\\n",
    "                * ( ( xP - x0)/current_dt - a_x0[0])**2\n",
    "                )\n",
    "        #\n",
    "        dP_Gaussian = - 2/(4*D_x0[0]) \\\n",
    "                * ( ( xP - x0)/current_dt - a_x0[0]) \\\n",
    "                * P_Gaussian\n",
    "        #\n",
    "        j_Gaussian = a_xP[0] * P_Gaussian \\\n",
    "                - D_xP[1] * P_Gaussian \\\n",
    "                - D_xP[0] * dP_Gaussian\n",
    "        #\n",
    "        integrand_Gaussian = j_Gaussian/D_xP[0] * ( a_xP[0] - D_xP[1] )\n",
    "        #\n",
    "        Sm_dot_Gaussian[i] = np.trapz(integrand_Gaussian,xP)\n",
    "\n",
    "\n",
    "Sm_dot_perturbative.append( Sm_dot_Gaussian )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate perturbative medium entropy production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = PySTFP()\n",
    "\n",
    "orders_for_perturbative_results = [1,4]\n",
    "\n",
    "for order in orders_for_perturbative_results:\n",
    "\t#\n",
    "\tcurrent_Sm_dot = p.get_dS_m_dt_from_derivatives(\n",
    "\t\t\t\tD_derivatives=D_x0, # numpy array, \n",
    "\t\t\t\t# D_derivatives[i] in units of self.D0 / self.L**i\n",
    "\t\t\t\ta_derivatives=a_x0, # numpy array,\n",
    "\t\t\t\t# a_derivatives[i] in units of self.L**(1-i)/self.T\n",
    "\t\t\t\torder=order, # order in epsilon for dimensionless (!) moment\n",
    "\t\t\t\t# note that if dimensionless_position == True, then \"order\"\n",
    "\t\t\t\t# is the order in physical time\n",
    "\t\t\t\tdimensionless_time=False,\n",
    "\t\t\t\t\t)\n",
    "\t#\n",
    "\tSm_dot_perturbative.append( current_Sm_dot(dt) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot medium entropy production rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axes = plt.subplots(3,1,figsize=(6,8))\n",
    "fig.subplots_adjust(hspace=0.00001,wspace=0.6)\n",
    "\n",
    "# Subplots:\n",
    "# axes[0]: Medium entropy production rate\n",
    "# axes[1]: Relative error for perturbative results\n",
    "# axes[2]: Running exponent of relative error\n",
    "\n",
    "\n",
    "# plot vertical line for the breakdown time in all subplots\n",
    "for i in range(3):\n",
    "    #\n",
    "    ps.add_vertical_breakdown_line(ax = axes[i],\n",
    "                                t_breakdown=analytical_solution['t_breakdown'],\n",
    "                                )\n",
    "\n",
    "\n",
    "###########################################\n",
    "# axes[0]: Medium entropy production rate #\n",
    "###########################################\n",
    "ax = axes[0]\n",
    "ax.set_title('Medium entropy production rate',y=1.03)\n",
    "\n",
    "# analytical\n",
    "ax.plot(\n",
    "    dt,\n",
    "    Sm_dot,\n",
    "    color=ps.analytical['color'],\n",
    "    label=ps.analytical['label'],\n",
    "    lw=ps.analytical['lw'],\n",
    "    )\n",
    "\n",
    "# perturbative\n",
    "for i, current_Sm_dot in enumerate(Sm_dot_perturbative):\n",
    "    #\n",
    "    ax.plot(\n",
    "        dt,\n",
    "        current_Sm_dot,\n",
    "        color=ps.perturbative['color'][i],\n",
    "        dashes=ps.perturbative['dashes'][i],\n",
    "        lw=ps.perturbative['lw'][i],\n",
    "        label=ps.perturbative['label'][i],\n",
    "    )\n",
    "\n",
    "\n",
    "####################################################\n",
    "# axes[1]: Relative error for perturbative results #\n",
    "####################################################\n",
    "ax = axes[1]\n",
    "\n",
    "# plot horizontal line\n",
    "ax.axhline(ps.horizontal['value'],\n",
    "            color=ps.horizontal['color'],\n",
    "            label=ps.horizontal['label'].format(ps.horizontal['value']),\n",
    "            lw=ps.horizontal['lw'])\n",
    "\n",
    "for i, current_Sm_dot in enumerate(Sm_dot_perturbative):\n",
    "    #\n",
    "    E_rel = np.fabs((current_Sm_dot-Sm_dot)/Sm_dot)\n",
    "    #\n",
    "    ax.plot(\n",
    "        dt,\n",
    "        E_rel,\n",
    "        color=ps.perturbative['color'][i],\n",
    "        dashes=ps.perturbative['dashes'][i],\n",
    "        lw=ps.perturbative['lw'][i],\n",
    "        label=ps.perturbative['label'][i],\n",
    "    )  \n",
    "    #\n",
    "    index = np.argmin(np.fabs( E_rel - ps.horizontal['value']) )\n",
    "    print('{0}: dt = {1:3.3f}, \\tE_rel = {2:3.3f}'.format(\n",
    "                        i,dt[index],ps.horizontal['value'])\n",
    "                        )\n",
    "\n",
    "\n",
    "###############################################\n",
    "# axes[2]: Running exponent of relative error #\n",
    "###############################################\n",
    "ax = axes[2]\n",
    "\n",
    "#\n",
    "for i, current_Sm_dot in enumerate(Sm_dot_perturbative):\n",
    "    #\n",
    "    diff = np.fabs((current_Sm_dot-Sm_dot)/Sm_dot)\n",
    "    #\n",
    "    local_exponent = np.log(diff[2:]) - np.log(diff[:-2])\n",
    "    local_exponent /= np.log(dt[2:]) - np.log(dt[:-2])\n",
    "    #\n",
    "    ax.plot(\n",
    "        dt[1:-1],\n",
    "        local_exponent,\n",
    "        color=ps.perturbative['color'][i],\n",
    "        dashes=ps.perturbative['dashes'][i],\n",
    "        lw=ps.perturbative['lw'][i],\n",
    "        label=ps.perturbative['label'][i],\n",
    "    )  \n",
    "\n",
    "####################################################\n",
    "# Set scaling, limits, labels etc for all subplots #\n",
    "####################################################\n",
    "xlims = [1e-3,1e0]\n",
    "ylims = {0:[-0.03,0.22],\n",
    "        1:[1e-9,1e1],\n",
    "        2:[0.8,5.5]}\n",
    "\n",
    "yticks = {0:[0.,0.1,0.2],\n",
    "        1:[1e-8,1e-6,1e-4,1e-2,1e0],\n",
    "        2:[1,2,3,4,5]}\n",
    "\n",
    "ylabels = {0:r'$\\dot{S}^{\\,\\mathrm{m}}\\cdot T$',\n",
    "            1:r'$E_{\\mathrm{rel}}$',\n",
    "            2:r'$\\kappa$'}\n",
    "\n",
    "for i in range(3):\n",
    "    ax = axes[i]\n",
    "    #\n",
    "    ax.set_xscale('log')\n",
    "    ax.set_xlim(*xlims)\n",
    "    ax.set_ylim(ylims[i])\n",
    "    if i == 1:\n",
    "        ax.set_yscale('log')\n",
    "    #\n",
    "    if i < 2:\n",
    "        ax.set_xticks([])\n",
    "    ax.set_yticks(yticks[i])\n",
    "    #\n",
    "    if i == 2:\n",
    "        ax.set_xlabel(r'$\\Delta t/T$')\n",
    "    ax.set_ylabel(ylabels[i])\n",
    "axes[1].set_zorder(-1)\n",
    "axes[0].legend(loc='lower left',\n",
    "        bbox_to_anchor=(0.,-0.2),\n",
    "        framealpha=0.99\n",
    "        )\n",
    "\n",
    "#################################################\n",
    "# Add enumeration (a), (b), (c) to the subplots #\n",
    "#################################################\n",
    "enumeration = [ '(a)','(b)','(c)']\n",
    "for i,e in enumerate(enumeration):\n",
    "    ax = axes[i]\n",
    "    #\n",
    "    xlims = ax.get_xlim()\n",
    "    dx = xlims[1] - xlims[0]\n",
    "    ylims = ax.get_ylim()\n",
    "    dy = ylims[1] - ylims[0]\n",
    "    #\n",
    "    ax.text(-0.2,\n",
    "            .95,\n",
    "            e,\n",
    "        horizontalalignment='center',\n",
    "        verticalalignment='center',\n",
    "        transform = ax.transAxes,\n",
    "                fontsize=25)\n",
    "\n",
    "plt.show()\n",
    "fig.savefig('Fig4_exact_vs_approximate_medium_entropy_production_rate.pdf',\n",
    "               bbox_inches='tight')\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "snowflake",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}