IllinoisGRMHD/doc/Tutorial-IllinoisGRMHD__inlined_functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Tutorial-IllinoisGRMHD: inlined_functions.C\n",
    "\n",
    "## Authors: Leo Werneck & Zach Etienne\n",
    "\n",
    "<font color='red'>**This module is currently under development**</font>\n",
    "\n",
    "## This notebook presents a series of inline functions that are used by major functions within IllinoisGRMHD.\n",
    "\n",
    "### Required and recommended citations:\n",
    "\n",
    "* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
    "* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
    "* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module is organized as follows\n",
    "\n",
    "0. [Step 0](#src_dir): **Source directory creation**\n",
    "1. [Step 1](#introduction): **Introduction**\n",
    "1. [Step 2](#find_cp_cm): **The `find_cp_cm()` function**\n",
    "1. [Step 3](#compute_v02): **The `compute_v02()` function**\n",
    "1. [Step 4](#font_fix__rhob_loop): **The `font_fix__rhob_loop()` function**\n",
    "1. [Step 5](#lower_4vector_output_spatial_part): **The `lower_4vector_output_spatial_part()` function**\n",
    "1. [Step 6](#impose_speed_limit_output_u0): **The `impose_speed_limit_output_u0()` function**\n",
    "1. [Step 7](#enforce_pressure_floor_ceiling): **The `enforce_pressure_floor_ceiling()` function**\n",
    "1. [Step 8](#compute_smallba_b2_and_u_i_over_u0_psi4): **The `compute_smallba_b2_and_u_i_over_u0_psi4()` function**\n",
    "1. [Step 9](#code_validation): **Code validation**\n",
    "1. [Step 10](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='src_dir'></a>\n",
    "\n",
    "# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
    "$$\\label{src_dir}$$\n",
    "\n",
    "We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory, if it does not exist yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Creation of the IllinoisGRMHD source directory\n",
    "# Step 0a: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "import os,sys\n",
    "nrpy_dir_path = os.path.join(\"..\",\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "# Step 0b: Load up cmdline_helper and create the directory\n",
    "import cmdline_helper as cmd\n",
    "IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
    "cmd.mkdir(IGM_src_dir_path)\n",
    "\n",
    "# Step 0c: Create the output file path\n",
    "outfile_path__inlined_functions__C = os.path.join(IGM_src_dir_path,\"inlined_functions.C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='introduction'></a>\n",
    "\n",
    "# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
    "$$\\label{introduction}$$\n",
    "\n",
    "In this tutorial notebook, we explain functions of `IllinoisGRMHD` which are called for various purposes. This means that this notebook does not have a specific \"theme\". We will cover functions whose purposes vary from a simple optimization when squaring numbers to computing minimum and maximum characteristic speeds at cell interfaces.\n",
    "\n",
    "We have tried our best to keep this tutorial module as independent from the others as possible. When new concepts appear, we offer useful references. The mathematical requirements of each function are also covered in great detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='find_cp_cm'></a>\n",
    "\n",
    "# Step 2: The `find_cp_cm()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{find_cp_cm}$$\n",
    "\n",
    "We will now explain the inlined function `find_cp_cm`. Keep in mind that this function depends on the function `compute_v02`, [which is implemented below](#compute_v02). This function is called with the objective of computing the minimum ($-$) and maximum ($+$) characteristic speeds at each cell interface, $c_{\\pm}^{r,l}$.\n",
    "\n",
    "We approximate the general GRMHD dispersion relation (eq. 27 of [Gammie & McKinney (2003)](https://arxiv.org/pdf/astro-ph/0301509.pdf)) by the simpler expression\n",
    "\n",
    "$$\n",
    "\\omega_{\\rm cm}^{2} = \\left[v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right)\\right]k_{\\rm cm}^{2}\\ ,\n",
    "$$\n",
    "\n",
    "where $\\omega_{\\rm cm}=-k_{\\mu}u^{\\mu}$ is the frequency and $k_{\\rm cm}^{2} = K_{\\mu}K^{\\mu}$ the wavenumber of an MHD wave mode in the frame comoving with the fluid, where $K_{\\mu}$ is defined as the projection of the wave vector $k^{\\nu}$ onto the direction normal to $u^{\\nu}$: $K_{\\mu} = \\left(g_{\\mu\\nu}+u_{\\mu}u_{\\nu}\\right)k^{\\nu}$. $c_{\\rm s}$ is the sound speed, and $v_{\\rm A}$ is the Alfvén speed, given by\n",
    "\n",
    "$$\n",
    "v_{\\rm A} = \\sqrt{\\frac{b^{2}}{\\rho_{b}h + b^{2}}}\\ .\n",
    "$$\n",
    "\n",
    "With these definitions, we may then solve the approximate dispersion relation above along direction $i$, noting that in the comoving frame $k_{\\mu} = \\left(-\\omega,k_{j}\\delta^{j}_{\\ i}\\right)$ and the wave (phase) velocity is $c_{\\pm} = \\left.\\omega\\middle/\\left(k_{j}\\delta^{j}_{\\ i}\\right)\\right.$. The dispersion can then be written as a quadratic equation for $c_{\\pm}$:\n",
    "\n",
    "$$\n",
    "ac_{\\pm}^{2} + bc_{\\pm} + c = 0\\ ,\n",
    "$$\n",
    "\n",
    "with\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "a &= \\left(1-v_{0}^{2}\\right)\\left(u^{0}\\right)^{2} - v_{0}^{2}g^{00}\\ ,\\\\\n",
    "b &= 2v_{0}^{2}g^{i0} - 2u^{i}u^{0}\\left(1-v^{2}_{0}\\right)\\ ,\\\\\n",
    "c &= \\left(1-v_{0}^{2}\\right)\\left(u^{i}\\right)^{2} - v_{0}^{2}g^{ii}\\ ,\\\\\n",
    "v_{0}^{2} &= v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right)\\ ,\\\\\n",
    "c_{\\rm s} &= \\left.\\left[\\frac{dP_{\\rm cold}}{d\\rho_{b}} + \\Gamma_{\\rm th}\\left(\\Gamma_{\\rm th}-1\\right)\\epsilon_{\\rm th}\\right]\\middle/h\\right.\\ ,\\\\\n",
    "c_{+} &= \\max\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ ,\\\\\n",
    "c_{-} &= \\min\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ .\n",
    "\\end{align}\n",
    "}\n",
    "$$\n",
    "\n",
    "For the implementation of $v_{0}^{2}$, please see [Step 4 below](#compute_v02)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "static inline void find_cp_cm(CCTK_REAL &cplus,CCTK_REAL &cminus,CCTK_REAL v02,CCTK_REAL u0,\n",
    "                       CCTK_REAL vi,CCTK_REAL ONE_OVER_LAPSE_SQUARED,CCTK_REAL shifti,CCTK_REAL psim4,CCTK_REAL gupii) {\n",
    "  // This computes phase speeds in the direction given by flux_dirn.\n",
    "  //  Note that we replace the full dispersion relation with a simpler\n",
    "  //     one, which overestimates the max. speeds by a factor of ~2.\n",
    "  //     See full discussion around Eqs. 49 and 50 in\n",
    "  //     http://arxiv.org/pdf/astro-ph/0503420.pdf .\n",
    "  //  What follows is a complete derivation of the quadratic we solve.\n",
    "  // wcm = (-k_0 u0 - k_x ux)\n",
    "  // kcm^2 = K_{\\mu} K^{\\mu},\n",
    "  // K_{\\mu} K^{\\mu} = (g_{\\mu a} + u_{\\mu} u_a) k^a * g^{\\mu b} [ (g_{c b} + u_c u_b) k^c ]\n",
    "  // --> g^{\\mu b} (g_{c b} + u_{c} u_{b}) k^c = (\\delta^{\\mu}_c + u_c u^{\\mu} ) k^c\n",
    "  //                 = (g_{\\mu a} + u_{\\mu} u_a) k^a * (\\delta^{\\mu}_c + u_c u^{\\mu} ) k^c\n",
    "  //                 =[(g_{\\mu a} + u_{\\mu} u_a) \\delta^{\\mu}_c + (g_{\\mu a} + u_{\\mu} u_a) u_c u^{\\mu} ] k^c k^a\n",
    "  //                 =[(g_{c a} + u_c u_a) + (u_c u_a -  u_a u_c] k^c k^a\n",
    "  //                 =(g_{c a} + u_c u_a) k^c k^a\n",
    "  //                 = k_a k^a + u^c u^a k_c k_a\n",
    "  // k^a = g^{\\mu a} k_{\\mu} = g^{0 a} k_0 + g^{x a} k_x\n",
    "  // k_a k^a = k_0 g^{0 0} k_0 + k_x k_0 g^{0 x} + g^{x 0} k_0 k_x + g^{x x} k_x k_x\n",
    "  //         = g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2\n",
    "  // u^c u^a k_c k_a = (u^0 k_0 + u^x k_x) (u^0 k_0 + u^x k_x) = (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2\n",
    "  // (k_0 u0)^2  + 2 k_x ux k_0 u0 + (k_x ux)^2 = v02 [ (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 + g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2]\n",
    "  // (1-v02) (u^0 k_0 + u^x k_x)^2 = v02 (g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2)\n",
    "  // (1-v02) (u^0 k_0/k_x + u^x)^2 = v02 (g^{00} (k_0/k_x)^2 + 2 g^{x0} k_0/k_x + g^{xx})\n",
    "  // (1-v02) (u^0 X + u^x)^2 = v02 (g^{00} X^2 + 2 g^{x0} X + g^{xx})\n",
    "  // (1-v02) (u0^2 X^2 + 2 ux u0 X + ux^2) = v02 (g^{00} X^2 + 2 g^{x0} X + g^{xx})\n",
    "  // X^2 ( (1-v02) u0^2 - v02 g^{00}) + X (2 ux u0 (1-v02) - 2 v02 g^{x0}) + (1-v02) ux^2 - v02 g^{xx}\n",
    "  // a = (1-v02) u0^2 - v02 g^{00} = (1-v02) u0^2 + v02/lapse^2 <-- VERIFIED\n",
    "  // b = 2 ux u0 (1-v02) - 2 v02 shiftx/lapse^2 <-- VERIFIED, X->-X, because X = -w/k_1, and we are solving for -X.\n",
    "  // c = (1-v02) ux^2 - v02 (gupxx*psim4 - (shiftx/lapse)^2) <-- VERIFIED\n",
    "  // v02 = v_A^2 + c_s^2 (1 - v_A^2)\n",
    "  CCTK_REAL u0_SQUARED=SQR(u0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by setting\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "a &= \\left(1-v_{0}^{2}\\right)\\left(u^{0}\\right)^{2} - v_{0}^{2}g^{00}\\\\\n",
    "b &= 2v_{0}^{2}g^{i0} - 2u^{i}u^{0}\\left(1-v^{2}_{0}\\right)\\\\\n",
    "c &= \\left(1-v_{0}^{2}\\right)\\left(u^{i}\\right)^{2} - v_{0}^{2}g^{ii}\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "  //Find cplus, cminus:\n",
    "  CCTK_REAL a = u0_SQUARED * (1.0-v02) + v02*ONE_OVER_LAPSE_SQUARED;\n",
    "  CCTK_REAL b = 2.0* ( shifti*ONE_OVER_LAPSE_SQUARED * v02 - u0_SQUARED * vi * (1.0-v02) );\n",
    "  CCTK_REAL c = u0_SQUARED*SQR(vi) * (1.0-v02) - v02 * ( psim4*gupii -\n",
    "                                                         SQR(shifti)*ONE_OVER_LAPSE_SQUARED);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we find the minimum ($-$) and maximum ($+$) characteristic speeds\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "c_{+} &= \\max\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ ,\\\\\n",
    "c_{-} &= \\min\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ .\n",
    "\\end{align}\n",
    "}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $IGM_src_dir_path/inlined_functions.C\n",
    "\n",
    "  CCTK_REAL detm = b*b - 4.0*a*c;\n",
    "  //ORIGINAL LINE OF CODE:\n",
    "  //if(detm < 0.0) detm = 0.0;\n",
    "  //New line of code (without the if() statement) has the same effect:\n",
    "  detm = sqrt(0.5*(detm + fabs(detm))); /* Based on very nice suggestion from Roland Haas */\n",
    "\n",
    "  cplus = 0.5*(detm-b)/a;\n",
    "  cminus = -0.5*(detm+b)/a;\n",
    "  if (cplus < cminus) {\n",
    "    CCTK_REAL cp = cminus;\n",
    "    cminus = cplus;\n",
    "    cplus = cp;\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compute_v02'></a>\n",
    "\n",
    "# Step 3: The `compute_v02()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{compute_v02}$$\n",
    "\n",
    "This function is used to evaluate $v_{0}^{2}$, a quantity necessary for the computation of the minimum and maximum characteristic speeds at each cell interface, $c_{\\pm}^{r,l}$. For more information on this procedure, please see the [implementation of the `find_cp_cm` function in Step 3](#find_cp_cm).\n",
    "\n",
    "We start with the sound speed:\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "c_{\\rm s} = \\left.\\left[\\frac{dP_{\\rm cold}}{d\\rho_{b}} + \\Gamma_{\\rm th}\\left(\\Gamma_{\\rm th}-1\\right)\\epsilon_{\\rm th}\\right]\\middle/h\\right.\n",
    "}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "static inline void compute_v02(CCTK_REAL dPcold_drho,CCTK_REAL Gamma_th,CCTK_REAL eps_th,CCTK_REAL h,CCTK_REAL *smallb,CCTK_REAL *U,  CCTK_REAL &v02L) {\n",
    "\n",
    "  if(U[RHOB]<=0) { v02L=1.0; return; }\n",
    "\n",
    "  /* c_s = sound speed = (dP_c/drho + \\Gamma(\\Gamma-1) \\epsilon_th)/h */\n",
    "  CCTK_REAL c_s_squared  = (dPcold_drho + Gamma_th*(Gamma_th-1.0)*eps_th)/(h);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we compute the square of the Alfén speed, $v_{\\rm A}$, which is given by\n",
    "\n",
    "$$\n",
    "\\boxed{v_{\\rm A}^{2} = \\frac{b^{2}}{\\rho_{b}h + b^{2}}}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "  /* v_A = Alfven speed = sqrt( b^2/(rho0 h + b^2) ) */\n",
    "  CCTK_REAL v_A_squared = smallb[SMALLB2]/(smallb[SMALLB2] + U[RHOB]*(h));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, $v_{0}$ is related to the sound speed and the Alfén speed via\n",
    "\n",
    "$$\n",
    "\\boxed{v_{0}^{2} = v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right)}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "  v02L = v_A_squared + c_s_squared*(1.0-v_A_squared);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='font_fix__rhob_loop'></a>\n",
    "\n",
    "# Step 4: The `font_fix__rhob_loop()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{font_fix__rhob_loop}$$\n",
    "\n",
    "This function implements the main loop inside the [font_fix__hybrid_EOS()](Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.ipynb#font_fix_hybrid_eos) function (see [Font *et al.*](https://arxiv.org/pdf/gr-qc/9811015.pdf)).\n",
    "\n",
    "We now perform the following iterative process, which is described in detail in [Appendix A of Zachariah *et al.* (2012)](https://arxiv.org/pdf/1112.0568.pdf). We refer the reader to Eqs. (A60), (A61), and (A62).\n",
    "\n",
    "1. Store the previously computed values of $W_{n}$, $S_{{\\rm fluid},n}^{2}$, and $\\rho_{n}$\n",
    "2. Compute $h = 1 + \\epsilon_{\\rm cold} + P_{\\rm cold}/\\rho_{n}$\n",
    "3. Set\n",
    "$$\n",
    "\\boxed{\\rho_{n+1} = \\psi^{-6}\\rho_{\\star}\\left(1 + \\frac{S_{{\\rm fluid},n}^{2}}{\\rho_{\\star} h_{n}}\\right)^{-1/2}}\n",
    "$$\n",
    "\n",
    "4. For a given value of $n$, perform steps 1 (for $\\rho$), 2 and 3 until $\\left|\\rho_{n+1}-\\rho_{n}\\right| < \\rho_{n+1}\\epsilon$, where $\\epsilon$ is a user given tolerance\n",
    "5. After convergence is obtained, update:\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "h_{n+1} &= 1 + \\epsilon_{\\rm cold} + P_{\\rm cold}/\\rho_{n+1}\\\\\n",
    "W_{n+1} &= \\psi^{-6}\\sqrt{\\tilde{S}^{2}_{{\\rm fluid},n} + \\rho_{\\star}^{2} h_{n+1}^{2}}\\\\\n",
    "S_{{\\rm fluid},n+1}^{2} &= \\frac{W^{2}_{n+1}\\left(\\tilde{S}\\cdot\\tilde{S}\\right) + \\left(\\bar{B}\\cdot\\tilde{S}\\right)^{2}\\left(\\bar{B}^{2} + 2W_{n+1}\\right)}{\\left(W_{n+1} + \\bar{B}^{2}\\right)^{2}}\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$\n",
    "6. Repeat steps 1 through 5 until $\\left|W_{n+1}-W_{n}\\right| < W_{n+1}\\epsilon$ *and* $\\left|S^{2}_{{\\rm fluid},n+1}-S^{2}_{{\\rm fluid},n}\\right| < S^{2}_{{\\rm fluid},n+1}\\epsilon$ *or* we reach the maximum number of iterations\n",
    "7. If font fix fails, increase the tolerance and try again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "/* Function    : font_fix__rhob_loop()\n",
    " * Authors     : Leo Werneck\n",
    " * Description : Determines rhob using the font fix prescription\n",
    " * Dependencies: find_polytropic_K_and_Gamma_index()\n",
    " *             : compute_P_cold__eps_cold()\n",
    " * Reference   : Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf]\n",
    " *\n",
    " * Inputs      : maxits          - maximum number of iterations allowed\n",
    " *             : tol             - font fix tolerance\n",
    " *             : W               - See eq. (A26)\n",
    " *             : Sf2             - S_{fluid}^{2}, see eq. (A24)\n",
    " *             : Psim6           - This is equal to sqrt(\\gamma)\n",
    " *             : sdots           - \\tilde{S}_{\\mu}\\tilde{S}^{\\mu}\n",
    " *             : BbardotS2       - (\\bar{B}^{\\mu}S_{\\mu})^{2},\n",
    " *             : B2bar           - \\bar{B}^{2}, see eq. (A28)\n",
    " *             : CONSERVS        - Array of conservative variables\n",
    " *             : eos             - Struct of EOS parameters\n",
    " *             : rhob_in         - Initial value of rhob\n",
    " *             : rhob_out        - Output variable\n",
    " *\n",
    " * Outputs     : rhob_out        - Updated value of rhob\n",
    " *             : return value: 0 - Font fix worked\n",
    " *             : return value: 1 - Font fix failed\n",
    " */\n",
    "inline int font_fix__rhob_loop( int maxits, CCTK_REAL tol,\n",
    "                                CCTK_REAL W, CCTK_REAL Sf2, CCTK_REAL Psim6, CCTK_REAL sdots, CCTK_REAL BbardotS2, CCTK_REAL B2bar,\n",
    "                                CCTK_REAL *CONSERVS,\n",
    "                                eos_struct eos, CCTK_REAL rhob_in, CCTK_REAL &rhob_out ) {\n",
    "\n",
    "  /* Declare basic variables */\n",
    "  bool fontcheck=true;\n",
    "  int itcount = 0, j0, j1;\n",
    "  CCTK_REAL W0, Sf20, rhob0, rhob1, h, P_cold, eps_cold;\n",
    "\n",
    "  //////////////////////\n",
    "  // OUTER LOOP START //\n",
    "  //////////////////////\n",
    "  while(fontcheck && itcount < maxits) {\n",
    "\n",
    "    /* Set variables to their input values */\n",
    "    itcount++;\n",
    "    W0    = W;\n",
    "    Sf20  = Sf2;\n",
    "    rhob1 = rhob_in;\n",
    "\n",
    "    /* Based on rhob_in (i.e. rhob1), determine the\n",
    "     * polytropic index j1\n",
    "     */\n",
    "    j1 = find_polytropic_K_and_Gamma_index(eos,rhob1);\n",
    "\n",
    "    //////////////////////\n",
    "    // INNER LOOP START //\n",
    "    //////////////////////\n",
    "    do {\n",
    "\n",
    "      /* Set rhob0/j0 to be equal to the rhob/j used\n",
    "       * in the previous iteration, i.e. rhob1/j1.\n",
    "       */\n",
    "      rhob0 = rhob1;\n",
    "      j0    = j1;\n",
    "\n",
    "      /* Compute h using h_cold and our polytropic EOS\n",
    "       * .------------------------------------------.\n",
    "       * | h = h_cold = 1 + eps_cold + P_cold/rhob. |\n",
    "       * .------------------------------------------.\n",
    "       */\n",
    "      compute_P_cold__eps_cold(eos,rhob0, P_cold, eps_cold);\n",
    "      h = 1.0 + eps_cold + P_cold/rhob0;\n",
    "\n",
    "      /* Update rhob using eq. (A62) in Etienne et al. (2011)\n",
    "       *          https://arxiv.org/pdf/1112.0568.pdf\n",
    "       * .---------------------------------------------------------------------------.\n",
    "       * | rhob = rho_star * Psi^{-6} / sqrt( 1 + S_fluid^{2}/( (rho_star*h)^{2} ) ) |\n",
    "       * .---------------------------------------------------------------------------.\n",
    "       */\n",
    "      rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n",
    "\n",
    "      /* Update j1 */\n",
    "      j1 = find_polytropic_K_and_Gamma_index(eos,rhob1);\n",
    "\n",
    "    }  while( fabs(rhob1-rhob0) > rhob1*tol || j1 != j0);\n",
    "    //////////////////////\n",
    "    //  INNER LOOP END  //\n",
    "    //////////////////////\n",
    "\n",
    "    /* Output the last value of rhob */\n",
    "    rhob_out = rhob1;\n",
    "\n",
    "    /* Perform physical checks on the variables\n",
    "     * and output the last value of h obtained\n",
    "     */\n",
    "    compute_P_cold__eps_cold(eos,rhob_out, P_cold, eps_cold);\n",
    "    h = 1.0 + eps_cold + P_cold/rhob_out;\n",
    "\n",
    "    /* Set W based on eq. (A60) in Etienne et al. (2011)\n",
    "     *       https://arxiv.org/pdf/1112.0568.pdf\n",
    "     * .-------------------------------------------------------.\n",
    "     * | W = psi^{-6} * sqrt( S_fluid^{2} + (rho_star*h)^{2} ) |\n",
    "     * .-------------------------------------------------------.\n",
    "     */\n",
    "    W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n",
    "\n",
    "    /* Then update S_{fluid}^{2} using eq. (A61) in Etienne et al. (2011)\n",
    "     *           https://arxiv.org/pdf/1112.0568.pdf\n",
    "     * .---------------------------------------------------------------------------.\n",
    "     * | S_fluid^{2} = ( W^{2}*S^{2} + (B.S)^2*(B^{2} + 2W) )/( ( W + B^{2} )^{2} )|\n",
    "     * .---------------------------------------------------------------------------.\n",
    "     */\n",
    "    Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n",
    "\n",
    "    if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n",
    "\n",
    "  }\n",
    "  //////////////////////\n",
    "  //  OUTER LOOP END  //\n",
    "  //////////////////////\n",
    "\n",
    "  /* If the code converged before the max\n",
    "   * number of iterations were exceeded,\n",
    "   * return 0, otherwise return 1.\n",
    "   */\n",
    "  if(fontcheck || itcount >= maxits) {\n",
    "    return 1;\n",
    "  }\n",
    "  else {\n",
    "    return 0;\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lower_4vector_output_spatial_part'></a>\n",
    "\n",
    "# Step 5: The `lower_4vector_output_spatial_part()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{lower_4vector_output_spatial_part}$$\n",
    "\n",
    "This function is used to lower the indices of the spatial components of 4-vectors, $b^{\\mu}$. Consider\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "b_{i} &= g_{i\\mu}b^{\\mu}                                       \\\\\n",
    "      &= g_{i0}b^{0} + g_{ij}b^{j}                             \\\\\n",
    "      &= \\left(\\gamma_{ij}\\beta^{j}\\right)b^{0} + \\gamma_{ij}b^{j} \\\\\n",
    "      &= \\gamma_{ij}\\left(b^{j} + \\beta^{j}b^{0}\\right)\\ ,\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "or, using the conformal metric and each component separately\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "b_{x} &= \\psi^{4}\\left[\\bar{\\gamma}_{xx}\\left(b^{x} + \\beta^{x}b^{0}\\right)+\\bar{\\gamma}_{xy}\\left(b^{y} + \\beta^{y}b^{0}\\right)+\\bar{\\gamma}_{xz}\\left(b^{z} + \\beta^{z}b^{0}\\right)\\right]\\\\\n",
    "b_{y} &= \\psi^{4}\\left[\\bar{\\gamma}_{yx}\\left(b^{x} + \\beta^{x}b^{0}\\right)+\\bar{\\gamma}_{yy}\\left(b^{y} + \\beta^{y}b^{0}\\right)+\\bar{\\gamma}_{yz}\\left(b^{z} + \\beta^{z}b^{0}\\right)\\right]\\\\\n",
    "b_{z} &= \\psi^{4}\\left[\\bar{\\gamma}_{zx}\\left(b^{x} + \\beta^{x}b^{0}\\right)+\\bar{\\gamma}_{zy}\\left(b^{y} + \\beta^{y}b^{0}\\right)+\\bar{\\gamma}_{zz}\\left(b^{z} + \\beta^{z}b^{0}\\right)\\right]\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "// b_x = g_{\\mu x} b^{\\mu}\n",
    "//     = g_{t x} b^t + g_{i x} b^i\n",
    "//     = b^t gamma_{xj} beta^j + gamma_{ix} b^i\n",
    "//     = gamma_{xj} (b^j + beta^j b^t)\n",
    "static inline void lower_4vector_output_spatial_part(CCTK_REAL psi4,CCTK_REAL *METRIC,CCTK_REAL *smallb, CCTK_REAL *smallb_lower) {\n",
    "  smallb_lower[SMALLBX] = psi4*( METRIC[GXX]*(smallb[SMALLBX]+smallb[SMALLBT]*METRIC[SHIFTX]) + METRIC[GXY]*(smallb[SMALLBY]+smallb[SMALLBT]*METRIC[SHIFTY]) +\n",
    "                                 METRIC[GXZ]*(smallb[SMALLBZ]+smallb[SMALLBT]*METRIC[SHIFTZ]) );\n",
    "  smallb_lower[SMALLBY] = psi4*( METRIC[GXY]*(smallb[SMALLBX]+smallb[SMALLBT]*METRIC[SHIFTX]) + METRIC[GYY]*(smallb[SMALLBY]+smallb[SMALLBT]*METRIC[SHIFTY]) +\n",
    "                                 METRIC[GYZ]*(smallb[SMALLBZ]+smallb[SMALLBT]*METRIC[SHIFTZ]) );\n",
    "  smallb_lower[SMALLBZ] = psi4*( METRIC[GXZ]*(smallb[SMALLBX]+smallb[SMALLBT]*METRIC[SHIFTX]) + METRIC[GYZ]*(smallb[SMALLBY]+smallb[SMALLBT]*METRIC[SHIFTY]) +\n",
    "                                 METRIC[GZZ]*(smallb[SMALLBZ]+smallb[SMALLBT]*METRIC[SHIFTZ]) );\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='impose_speed_limit_output_u0'></a>\n",
    "\n",
    "# Step 6: The `impose_speed_limit_output_u0()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{impose_speed_limit_output_u0}$$\n",
    "\n",
    "We now call upon the `impose_speed_limit_output_u0()` function inside the `inlined_functions.C` code file of `IllinoisGRMHD`. The basic algorithm performed by this function is summarized here. We start by evaluating the quantity\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "{\\rm one\\_minus\\_one\\_over\\_alpha\\_u0\\_squared} \\equiv A \n",
    "&= \\gamma_{ij}\\left(\\frac{v^{i}+\\beta^{i}}{\\alpha}\\right)\\left(\\frac{v^{j}+\\beta^{j}}{\\alpha}\\right)\\\\\n",
    "&= \\frac{\\gamma_{ij}}{\\alpha^{2}}\\left[\\frac{\\gamma^{ik}u_{k}}{u^{0}} - \\beta^{i} + \\beta^{i}\\right]\\left[\\frac{\\gamma^{j\\ell}u_{\\ell}}{u^{0}} - \\beta^{j} + \\beta^{j}\\right]\\\\\n",
    "&=\\frac{\\gamma_{ij}u^{i}u^{j}}{\\left(\\alpha u^{0}\\right)^{2}}\\\\\n",
    "&=\\frac{\\left(\\alpha u^{0}\\right)^{2}-1}{\\left(\\alpha u^{0}\\right)^{2}}\\\\\n",
    "&=1 - \\frac{1}{\\left(\\alpha u^{0}\\right)^{2}}\\ \\\\\n",
    "\\implies \\boxed{A = 1 - \\frac{1}{\\left(\\alpha u^{0}\\right)^{2}}}\\ ,\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where when going from line 1 to 2 and from line 3 to 4 we have used Eqs. (53) and (56) from [Duez *et al.*](https://arxiv.org/pdf/astro-ph/0503420.pdf), respectively. Keep in mind that the equation we are going to implement below is\n",
    "\n",
    "$$\n",
    "\\boxed{{\\rm one\\_minus\\_one\\_over\\_alpha\\_u0\\_squared} = \\gamma_{ij}\\left(\\frac{v^{i}+\\beta^{i}}{\\alpha}\\right)\\left(\\frac{v^{j}+\\beta^{j}}{\\alpha}\\right)}\\ ,\n",
    "$$\n",
    "\n",
    "but it is important to know that this equation also equals $A$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "static inline void impose_speed_limit_output_u0(CCTK_REAL *METRIC,CCTK_REAL *U,CCTK_REAL psi4,CCTK_REAL ONE_OVER_LAPSE,output_stats &stats, CCTK_REAL &u0_out) {\n",
    "\n",
    "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "#endif\n",
    "\n",
    "\n",
    "  // Derivation of first equation:\n",
    "  // \\gamma_{ij} (v^i + \\beta^i)(v^j + \\beta^j)/(\\alpha)^2\n",
    "  //   = \\gamma_{ij} 1/(u^0)^2 ( \\gamma^{ik} u_k \\gamma^{jl} u_l /(\\alpha)^2 <- Using Eq. 53 of arXiv:astro-ph/0503420\n",
    "  //   = 1/(u^0 \\alpha)^2 u_j u_l \\gamma^{jl}  <- Since \\gamma_{ij} \\gamma^{ik} = \\delta^k_j\n",
    "  //   = 1/(u^0 \\alpha)^2 ( (u^0 \\alpha)^2 - 1 ) <- Using Eq. 56 of arXiv:astro-ph/0503420\n",
    "  //   = 1 - 1/(u^0 \\alpha)^2 <= 1\n",
    "  CCTK_REAL one_minus_one_over_alpha_u0_squared = psi4*(METRIC[GXX]* SQR(U[VX] + METRIC[SHIFTX]) +\n",
    "                                                        2.0*METRIC[GXY]*(U[VX] + METRIC[SHIFTX])*(U[VY] + METRIC[SHIFTY]) +\n",
    "                                                        2.0*METRIC[GXZ]*(U[VX] + METRIC[SHIFTX])*(U[VZ] + METRIC[SHIFTZ]) +\n",
    "                                                        METRIC[GYY]* SQR(U[VY] + METRIC[SHIFTY]) +\n",
    "                                                        2.0*METRIC[GYZ]*(U[VY] + METRIC[SHIFTY])*(U[VZ] + METRIC[SHIFTZ]) +\n",
    "                                                        METRIC[GZZ]* SQR(U[VZ] + METRIC[SHIFTZ]) )*SQR(ONE_OVER_LAPSE);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we construct the \"speed limit quantity\"\n",
    "\n",
    "$$\n",
    "{\\rm ONE\\_MINUS\\_ONE\\_OVER\\_GAMMA\\_SPEED\\_LIMIT\\_SQUARED} \\equiv B = 1-\\frac{1}{\\gamma^{2}_{\\rm speed\\ limit}}\\ .\n",
    "$$\n",
    "\n",
    "If $A > B$, then we construct the correction factor $C\\equiv A / B$, and adjust the velocities using\n",
    "\n",
    "$$\n",
    "\\boxed{v^{i} \\to \\left(v^{i}+\\beta^{i}\\right)C - \\beta^{i}}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "  /*** Limit velocity to GAMMA_SPEED_LIMIT ***/\n",
    "  const CCTK_REAL ONE_MINUS_ONE_OVER_GAMMA_SPEED_LIMIT_SQUARED = 1.0-1.0/SQR(GAMMA_SPEED_LIMIT);\n",
    "  if(one_minus_one_over_alpha_u0_squared > ONE_MINUS_ONE_OVER_GAMMA_SPEED_LIMIT_SQUARED) {\n",
    "    CCTK_REAL correction_fac = sqrt(ONE_MINUS_ONE_OVER_GAMMA_SPEED_LIMIT_SQUARED/one_minus_one_over_alpha_u0_squared);\n",
    "    U[VX] = (U[VX] + METRIC[SHIFTX])*correction_fac-METRIC[SHIFTX];\n",
    "    U[VY] = (U[VY] + METRIC[SHIFTY])*correction_fac-METRIC[SHIFTY];\n",
    "    U[VZ] = (U[VZ] + METRIC[SHIFTZ])*correction_fac-METRIC[SHIFTZ];\n",
    "    one_minus_one_over_alpha_u0_squared=ONE_MINUS_ONE_OVER_GAMMA_SPEED_LIMIT_SQUARED;\n",
    "    stats.failure_checker+=1000;\n",
    "  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, since $A$ is evaluated using the first line above, namely\n",
    "\n",
    "$$\n",
    "\\gamma_{ij}\\left(\\frac{v^{i}+\\beta^{i}}{\\alpha}\\right)\\left(\\frac{v^{j}+\\beta^{j}}{\\alpha}\\right) = A = 1 - \\frac{1}{\\left(\\alpha u^{0}\\right)^{2}}\\ ,\n",
    "$$\n",
    "\n",
    "we can then compute $u_{0}$ by simply doing\n",
    "\n",
    "$$\n",
    "\\boxed{u^{0} = \\frac{1}{\\alpha\\sqrt{1-A}}}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "  // A = 1.0-one_minus_one_over_alpha_u0_squared = 1-(1-1/(al u0)^2) = 1/(al u0)^2\n",
    "  // 1/sqrt(A) = al u0\n",
    "  //CCTK_REAL alpha_u0_minus_one = 1.0/sqrt(1.0-one_minus_one_over_alpha_u0_squared)-1.0;\n",
    "  //u0_out          = (alpha_u0_minus_one + 1.0)*ONE_OVER_LAPSE;\n",
    "  CCTK_REAL alpha_u0 = 1.0/sqrt(1.0-one_minus_one_over_alpha_u0_squared);\n",
    "  if(std::isnan(alpha_u0*ONE_OVER_LAPSE)) printf(\"BAD FOUND NAN U0 CALC: %.15e %.15e %.15e | %.15e %.15e\\n\",alpha_u0,ONE_OVER_LAPSE,one_minus_one_over_alpha_u0_squared,psi4, U[VX]);\n",
    "  u0_out = alpha_u0*ONE_OVER_LAPSE;\n",
    "}\n",
    "\n",
    "// The two lines of code below are written to reduce roundoff error and were in the above function. I don't think they reduce error.\n",
    "// one_over_alpha_u0  = sqrt(1.0-one_minus_one_over_alpha_u0_squared);\n",
    "/* Proof of following line: */\n",
    "/*   [        1-1/(alphau0)^2         ] / [ 1/(alphau0) (1 + 1/(alphau0)) ] */\n",
    "/* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [ 1/(alphau0) + 1/(alphau0)^2   ] */\n",
    "/* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [  (alphau0 + 1)/(alphau0)^2    ] */\n",
    "/* = [        (alphau0)^2 - 1)        ] / [        (alphau0 + 1)          ] */\n",
    "/*   [   (alphau0 + 1) (alphau0 - 1)  ] / [        (alphau0 + 1)          ] */\n",
    "/* =                              alphau0 - 1                               */\n",
    "//alpha_u0_minus_one = one_minus_one_over_alpha_u0_squared/one_over_alpha_u0/(1.0+one_over_alpha_u0);\n",
    "//u0_out = (alpha_u0_minus_one+1.0)*ONE_OVER_LAPSE;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='enforce_pressure_floor_ceiling'></a>\n",
    "\n",
    "# Step 7: The `enforce_pressure_floor_ceiling()` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{enforce_pressure_floor_ceiling}$$\n",
    "\n",
    "After the Newton-Raphson solver has successfully found a set of primitives, the primitives are checked for physicality, and if they are not in the physical range, they are minimally modified until they return to the physical range. First, if the velocity is found to be superluminal, the speed is reduced to `IllinoisGRMHD`’s default Lorentz factor limit, a procedure that we already explained above when we discussed the `impose_speed_limit_output_u0` function.\n",
    "\n",
    "Next, `IllinoisGRMHD` does not include any cooling mechanism, which means that for evolutions adopting a $\\Gamma$-law equation of state, the pressure should not physically drop below $P_{\\rm cold}$. So a pressure floor of $0.9P_{\\rm cold}$ is imposed. Increasing this floor to $P_{\\rm cold}$ exactly results in large central density drifts in TOV star evolutions.\n",
    "\n",
    "**NOTE**: Please keep in mind that the floor and ceiling values presented here were found ***empirically***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "static inline void enforce_pressure_floor_ceiling(output_stats &stats,CCTK_REAL kpoly,CCTK_REAL P_cold,CCTK_REAL Psi6,const CCTK_REAL Psi6threshold,CCTK_REAL rho_b,const CCTK_REAL rhobatm,  CCTK_REAL &P) {\n",
    "  CCTK_REAL P_min=0.9*P_cold;\n",
    "  if(P<P_min) {\n",
    "    stats.failure_checker+=10;\n",
    "    P=P_min;\n",
    "  }\n",
    "  //MAX(P,P_min);\n",
    "  //if(P < P_min) P=1.0*P_cold;\n",
    "\n",
    "  /* OLD: Discarded because lower limit is unphysical.\n",
    "     if(P <= 0.5*kpoly*P_cold) {\n",
    "     P=0.5*kpoly*P_cold;\n",
    "     }\n",
    "  */"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulations can crash in the other extreme, if $P/P_{\\rm cold}$ becomes too large. This typically only happens in very low-density regions or inside black holes.  So at densities $\\rho_{b}<100\\rho_{\\rm atm}$ or deep inside black hole horizons, a ceiling on $P$ of $100P_{\\rm cold}$ is enforced (see Appendix A of [Etienne *et al.* (2012)](https://arxiv.org/abs/1112.0568) for more details).\n",
    "\n",
    "We also introduce a parameter, $\\psi^{6}_{\\rm threshold}$, which determines whether the region under consideration is deep inside the BH horizon or not. For regions deep inside the BH horizon, defined by $\\sqrt{\\gamma} = \\psi^{6} > \\psi^{6}_{\\rm threshold}$, the primary goal is to keep the evolution stable and prevent inaccurate data from leaking out of the BH horizon. It was determined that in this situation, a better ceiling on $P$ is $10^{5}P_{\\rm cold}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "  //CCTK_REAL P_max = 10.0*P_cold;\n",
    "  CCTK_REAL P_max = 100.0*P_cold;\n",
    "  if(Psi6 > Psi6threshold) P_max = 1e5*P_cold; // <-- better than 10.\n",
    "\n",
    "  if((rho_b < 100.0*rhobatm || Psi6 > Psi6threshold) && P>P_max) {\n",
    "    P=P_max;\n",
    "    stats.failure_checker+=100;\n",
    "  }\n",
    "\n",
    "  /*\n",
    "    CCTK_REAL rho_horiz_cap = 1000.0*rhobatm;\n",
    "\n",
    "    //New density damping mechanism inside the horizon\n",
    "    if(Psi6 > Psi6threshold && rho_b>rho_horiz_cap) {\n",
    "    CCTK_REAL six_phi=log(Psi6);\n",
    "    CCTK_REAL six_phithreshold=log(Psi6threshold);\n",
    "    CCTK_REAL Psi6max_approx=350000;\n",
    "    rho_b = rho_horiz_cap+(rho_b-rho_horiz_cap)*exp(-200.0*SQR((six_phi-six_phithreshold)/log(Psi6max_approx)));\n",
    "    }\n",
    "  */\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compute_smallba_b2_and_u_i_over_u0_psi4'></a>\n",
    "\n",
    "# Step 8: The `compute_smallba_b2_and_u_i_over_u0_psi4` function \\[Back to [top](#toc)\\]\n",
    "$$\\label{compute_smallba_b2_and_u_i_over_u0_psi4}$$\n",
    "\n",
    "In this inlined function we will compute quantities related to the magnetic field measured in the comoving fluid frame, $b^{\\mu}$.\n",
    "\n",
    "We will need the following identities\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "v^{i}       &= \\frac{u^{i}}{u^{0}}\\ ,\\\\\n",
    "B^{0}_{(u)} &= \\frac{u_{i}B^{i}}{\\alpha}\\ ,\\\\\n",
    "B^{i}_{(u)} &= \\frac{1}{u^{0}}\\left(\\frac{B^{i}}{\\alpha} + u^{i}B^{0}_{(u)}\\right)\\ ,\\\\\n",
    "b^{\\mu}     &= \\frac{B^{\\mu}_{(u)}}{\\sqrt{4\\pi}}\\ .\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We start by setting the relation\n",
    "\n",
    "$$\n",
    "b^{0} = \\frac{u_{i}B^{i}}{\\alpha\\sqrt{4\\pi}} \\implies \\boxed{\\alpha\\sqrt{4\\pi}b^{0} = u_{i}B^{i}}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "static inline void compute_smallba_b2_and_u_i_over_u0_psi4(CCTK_REAL *METRIC,CCTK_REAL *METRIC_LAP_PSI4,CCTK_REAL *U,CCTK_REAL u0L,CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI,\n",
    "                                                    CCTK_REAL &u_x_over_u0_psi4,CCTK_REAL &u_y_over_u0_psi4,CCTK_REAL &u_z_over_u0_psi4,CCTK_REAL *smallb) {\n",
    "\n",
    "  // NOW COMPUTE b^{\\mu} and b^2 = b^{\\mu} b^{\\nu} g_{\\mu \\nu}\n",
    "  CCTK_REAL ONE_OVER_U0 = 1.0/u0L;\n",
    "  CCTK_REAL shiftx_plus_vx = (METRIC[SHIFTX]+U[VX]);\n",
    "  CCTK_REAL shifty_plus_vy = (METRIC[SHIFTY]+U[VY]);\n",
    "  CCTK_REAL shiftz_plus_vz = (METRIC[SHIFTZ]+U[VZ]);\n",
    "\n",
    "  // Eq. 56 in http://arxiv.org/pdf/astro-ph/0503420.pdf:\n",
    "  //  u_i = gamma_{ij} u^0 (v^j + beta^j), gamma_{ij} is the physical metric, and gamma_{ij} = Psi4 * METRIC[Gij], since METRIC[Gij] is the conformal metric.\n",
    "  u_x_over_u0_psi4 =  METRIC[GXX]*shiftx_plus_vx + METRIC[GXY]*shifty_plus_vy + METRIC[GXZ]*shiftz_plus_vz;\n",
    "  u_y_over_u0_psi4 =  METRIC[GXY]*shiftx_plus_vx + METRIC[GYY]*shifty_plus_vy + METRIC[GYZ]*shiftz_plus_vz;\n",
    "  u_z_over_u0_psi4 =  METRIC[GXZ]*shiftx_plus_vx + METRIC[GYZ]*shifty_plus_vy + METRIC[GZZ]*shiftz_plus_vz;\n",
    "\n",
    "  // Eqs. 23 and 31 in http://arxiv.org/pdf/astro-ph/0503420.pdf:\n",
    "  //   Compute alpha sqrt(4 pi) b^t = u_i B^i\n",
    "  CCTK_REAL alpha_sqrt_4pi_bt = ( u_x_over_u0_psi4*U[BX_CENTER] + u_y_over_u0_psi4*U[BY_CENTER] + u_z_over_u0_psi4*U[BZ_CENTER] ) * METRIC_LAP_PSI4[PSI4]*u0L;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we compute\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "b^{i} &= \\frac{B^{i}_{(u)}}{\\sqrt{4\\pi}}\\\\\n",
    "      &= \\frac{1}{u^{0}\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{\\alpha} + B^{0}_{(u)}u^{i}\\right)\\\\\n",
    "      &= \\frac{1}{u^{0}\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{\\alpha} + \\sqrt{4\\pi}b^{0}u^{i}\\right)\\\\\n",
    "      &= \\frac{1}{\\alpha\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{u^{0}} + \\alpha\\sqrt{4\\pi}b^{0}\\frac{u^{i}}{u^{0}}\\right)\\\\\n",
    "\\implies &\\boxed{b^{i} = \\frac{1}{\\alpha\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{u^{0}} + \\alpha\\sqrt{4\\pi}b^{0}v^{i}\\right)}\\ .\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "  // Eq. 24 in http://arxiv.org/pdf/astro-ph/0503420.pdf:\n",
    "  // b^i = B^i_u / sqrt(4 pi)\n",
    "  // b^i = ( B^i/alpha + B^0_u u^i ) / ( u^0 sqrt(4 pi) )\n",
    "  // b^i = ( B^i/alpha +  sqrt(4 pi) b^t u^i ) / ( u^0 sqrt(4 pi) )\n",
    "  // b^i = ( B^i +  alpha sqrt(4 pi) b^t u^i ) / ( alpha u^0 sqrt(4 pi) )\n",
    "  // b^i = ( B^i/u^0 +  alpha sqrt(4 pi) b^t u^i/u^0 ) / ( alpha sqrt(4 pi) )\n",
    "  // b^i = ( B^i/u^0 +  alpha sqrt(4 pi) b^t v^i ) / ( alpha sqrt(4 pi) )\n",
    "  smallb[SMALLBX] = (U[BX_CENTER]*ONE_OVER_U0 + U[VX]*alpha_sqrt_4pi_bt)*ONE_OVER_LAPSE_SQRT_4PI;\n",
    "  smallb[SMALLBY] = (U[BY_CENTER]*ONE_OVER_U0 + U[VY]*alpha_sqrt_4pi_bt)*ONE_OVER_LAPSE_SQRT_4PI;\n",
    "  smallb[SMALLBZ] = (U[BZ_CENTER]*ONE_OVER_U0 + U[VZ]*alpha_sqrt_4pi_bt)*ONE_OVER_LAPSE_SQRT_4PI;\n",
    "  // Eq. 23 in http://arxiv.org/pdf/astro-ph/0503420.pdf, with alpha sqrt (4 pi) b^2 = u_i B^i already computed above\n",
    "  smallb[SMALLBT] = alpha_sqrt_4pi_bt * ONE_OVER_LAPSE_SQRT_4PI;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we compute\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "b^{2} &= g_{\\mu\\nu}b^{\\mu}b^{\\nu}\\\\\n",
    "      &= g_{00}\\left(b^{0}\\right)^{2} + g_{ij}b^{i}b^{j} + 2g_{0i}b^{0}b^{i}\\\\\n",
    "      &= \\left(-\\alpha^{2} + \\gamma_{ij}\\beta^{i}\\beta^{j}\\right)\\left(b^{0}\\right)^{2} + \\gamma_{ij}b^{i}b^{j} + 2b^{0}\\gamma_{ij}\\beta^{j}b^{i}\\\\\n",
    "      &= -\\left(\\alpha b^{0}\\right)^{2} + \\gamma_{ij}\\left[b^{i}b^{j} + 2b^{0}b^{i}\\beta^{j} + \\left(b^{0}\\right)^{2}\\beta^{i}\\beta^{j}\\right]\\\\\n",
    "\\implies &\\boxed{b^{2} = -\\left(\\alpha b^{0}\\right)^{2} + \\gamma_{ij}\\left(b^{i} + b^{0}\\beta^{i}\\right)\\left(b^{j} + b^{0}\\beta^{j}\\right)}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/inlined_functions.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__inlined_functions__C\n",
    "\n",
    "\n",
    "  // b^2 = g_{\\mu \\nu} b^{\\mu} b^{\\nu}\n",
    "  //     = gtt bt^2 + gxx bx^2 + gyy by^2 + gzz bz^2 + 2 (gtx bt bx + gty bt by + gtz bt bz + gxy bx by + gxz bx bz + gyz by bz)\n",
    "  //     = (-al^2 + gamma_{ij} betai betaj) bt^2 + b^i b^j gamma_{ij} + 2 g_{t i} b^t b^i\n",
    "  //     = - (alpha b^t)^2 + (b^t)^2 gamma_{ij} beta^i beta^j + b^i b^j gamma_{ij} + 2 b^t g_{t i} b^i\n",
    "  //     = - (alpha b^t)^2 + (b^t)^2 gamma_{ij} beta^i beta^j + b^i b^j gamma_{ij} + 2 b^t (gamma_{ij} beta^j) b^i\n",
    "  //     = - (alpha b^t)^2 + gamma_{ij} ((b^t)^2 beta^i beta^j + b^i b^j + 2 b^t beta^j b^i)\n",
    "  //     = - (alpha b^t)^2 + gamma_{ij} ((b^t)^2 beta^i beta^j + 2 b^t beta^j b^i + b^i b^j)\n",
    "  //     = - (alpha b^t)^2 + gamma_{ij} (b^i + b^t beta^i) (b^j + b^t beta^j)\n",
    "  CCTK_REAL bx_plus_shiftx_bt = smallb[SMALLBX]+METRIC[SHIFTX]*smallb[SMALLBT];\n",
    "  CCTK_REAL by_plus_shifty_bt = smallb[SMALLBY]+METRIC[SHIFTY]*smallb[SMALLBT];\n",
    "  CCTK_REAL bz_plus_shiftz_bt = smallb[SMALLBZ]+METRIC[SHIFTZ]*smallb[SMALLBT];\n",
    "  smallb[SMALLB2] = -SQR(METRIC_LAP_PSI4[LAPSE]*smallb[SMALLBT]) +\n",
    "    (  METRIC[GXX]*SQR(bx_plus_shiftx_bt) + METRIC[GYY]*SQR(by_plus_shifty_bt) + METRIC[GZZ]*SQR(bz_plus_shiftz_bt) +\n",
    "       2.0*( METRIC[GXY]*(bx_plus_shiftx_bt)*(by_plus_shifty_bt) +\n",
    "             METRIC[GXZ]*(bx_plus_shiftx_bt)*(bz_plus_shiftz_bt) +\n",
    "             METRIC[GYZ]*(by_plus_shifty_bt)*(bz_plus_shiftz_bt) )  ) * METRIC_LAP_PSI4[PSI4]; // mult by psi4 because METRIC[GIJ] is the conformal metric.\n",
    "  /***********************************************************/\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 9: Code validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "First, we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Validation test for inlined_functions.C: FAILED!\n",
      "Diff:\n",
      "1,4c1\n",
      "< static inline CCTK_REAL fasterpow_ppm_reconstruct(CCTK_REAL inputvar,CCTK_REAL inputpow) {\n",
      "<   if(inputpow==2.0) return SQR(inputvar);\n",
      "<   return pow(inputvar,inputpow);\n",
      "< }\n",
      "---\n",
      "> \n",
      "10c7\n",
      "<   //     one, which overestimates the max. speeds by a factor of ~2. \n",
      "---\n",
      ">   //     one, which overestimates the max. speeds by a factor of ~2.\n",
      "15c12\n",
      "<   // kcm^2 = K_{\\mu} K^{\\mu}, \n",
      "---\n",
      ">   // kcm^2 = K_{\\mu} K^{\\mu},\n",
      "38a36\n",
      "> \n",
      "43a42\n",
      "> \n",
      "49c48\n",
      "<   \n",
      "---\n",
      "> \n",
      "59c58,59\n",
      "< static inline void compute_v02(CCTK_REAL dPcold_drho,CCTK_REAL gamma_th,CCTK_REAL eps_th,CCTK_REAL h,CCTK_REAL *smallb,CCTK_REAL *U,  CCTK_REAL &v02L) {\n",
      "---\n",
      "> \n",
      "> static inline void compute_v02(CCTK_REAL dPcold_drho,CCTK_REAL Gamma_th,CCTK_REAL eps_th,CCTK_REAL h,CCTK_REAL *smallb,CCTK_REAL *U,  CCTK_REAL &v02L) {\n",
      "64c64,65\n",
      "<   CCTK_REAL c_s_squared  = (dPcold_drho + gamma_th*(gamma_th-1.0)*eps_th)/(h);\n",
      "---\n",
      ">   CCTK_REAL c_s_squared  = (dPcold_drho + Gamma_th*(Gamma_th-1.0)*eps_th)/(h);\n",
      "> \n",
      "66a68\n",
      "> \n",
      "70,84c72,180\n",
      "< static inline void compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(CCTK_REAL *U, eos_struct &eos,\n",
      "<                                                                          CCTK_REAL &P_cold,CCTK_REAL &eps_cold,CCTK_REAL &dPcold_drho,CCTK_REAL &eps_th,CCTK_REAL &h,\n",
      "<                                                                          CCTK_REAL &gamma_cold) {\n",
      "<   // This code handles equations of state of the form defined\n",
      "<   // in Eqs 13-16 in http://arxiv.org/pdf/0802.0200.pdf\n",
      "< \n",
      "<   if(U[RHOB]==0) {\n",
      "<     P_cold = 0.0;\n",
      "<     eps_cold = 0.0;\n",
      "<     dPcold_drho = 0.0;\n",
      "<     eps_th = 0.0;\n",
      "<     h = 0.0;\n",
      "<     gamma_cold = eos.gamma_tab[0];\n",
      "<     return;\n",
      "<   }\n",
      "---\n",
      "> /* Function    : font_fix__rhob_loop()\n",
      ">  * Authors     : Leo Werneck\n",
      ">  * Description : Determines rhob using the font fix prescription\n",
      ">  * Dependencies: find_polytropic_K_and_Gamma_index()\n",
      ">  *             : compute_P_cold__eps_cold()\n",
      ">  * Reference   : Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf]\n",
      ">  *\n",
      ">  * Inputs      : maxits          - maximum number of iterations allowed\n",
      ">  *             : tol             - font fix tolerance\n",
      ">  *             : W               - See eq. (A26)\n",
      ">  *             : Sf2             - S_{fluid}^{2}, see eq. (A24)\n",
      ">  *             : Psim6           - This is equal to sqrt(\\gamma)\n",
      ">  *             : sdots           - \\tilde{S}_{\\mu}\\tilde{S}^{\\mu}\n",
      ">  *             : BbardotS2       - (\\bar{B}^{\\mu}S_{\\mu})^{2},\n",
      ">  *             : B2bar           - \\bar{B}^{2}, see eq. (A28)\n",
      ">  *             : CONSERVS        - Array of conservative variables\n",
      ">  *             : eos             - Struct of EOS parameters\n",
      ">  *             : rhob_in         - Initial value of rhob\n",
      ">  *             : rhob_out        - Output variable\n",
      ">  *\n",
      ">  * Outputs     : rhob_out        - Updated value of rhob\n",
      ">  *             : return value: 0 - Font fix worked\n",
      ">  *             : return value: 1 - Font fix failed\n",
      ">  */\n",
      "> inline int font_fix__rhob_loop( int maxits, CCTK_REAL tol,\n",
      ">                                 CCTK_REAL W, CCTK_REAL Sf2, CCTK_REAL Psim6, CCTK_REAL sdots, CCTK_REAL BbardotS2, CCTK_REAL B2bar,\n",
      ">                                 CCTK_REAL *CONSERVS,\n",
      ">                                 eos_struct eos, CCTK_REAL rhob_in, CCTK_REAL &rhob_out ) {\n",
      "> \n",
      ">   /* Declare basic variables */\n",
      ">   bool fontcheck=true;\n",
      ">   int itcount = 0, j0, j1;\n",
      ">   CCTK_REAL W0, Sf20, rhob0, rhob1, h, P_cold, eps_cold;\n",
      "> \n",
      ">   //////////////////////\n",
      ">   // OUTER LOOP START //\n",
      ">   //////////////////////\n",
      ">   while(fontcheck && itcount < maxits) {\n",
      "> \n",
      ">     /* Set variables to their input values */\n",
      ">     itcount++;\n",
      ">     W0    = W;\n",
      ">     Sf20  = Sf2;\n",
      ">     rhob1 = rhob_in;\n",
      "> \n",
      ">     /* Based on rhob_in (i.e. rhob1), determine the\n",
      ">      * polytropic index j1\n",
      ">      */\n",
      ">     j1 = find_polytropic_K_and_Gamma_index(eos,rhob1);\n",
      "> \n",
      ">     //////////////////////\n",
      ">     // INNER LOOP START //\n",
      ">     //////////////////////\n",
      ">     do {\n",
      "> \n",
      ">       /* Set rhob0/j0 to be equal to the rhob/j used\n",
      ">        * in the previous iteration, i.e. rhob1/j1.\n",
      ">        */\n",
      ">       rhob0 = rhob1;\n",
      ">       j0    = j1;\n",
      "> \n",
      ">       /* Compute h using h_cold and our polytropic EOS\n",
      ">        * .------------------------------------------.\n",
      ">        * | h = h_cold = 1 + eps_cold + P_cold/rhob. |\n",
      ">        * .------------------------------------------.\n",
      ">        */\n",
      ">       compute_P_cold__eps_cold(eos,rhob0, P_cold, eps_cold);\n",
      ">       h = 1.0 + eps_cold + P_cold/rhob0;\n",
      "> \n",
      ">       /* Update rhob using eq. (A62) in Etienne et al. (2011)\n",
      ">        *          https://arxiv.org/pdf/1112.0568.pdf\n",
      ">        * .---------------------------------------------------------------------------.\n",
      ">        * | rhob = rho_star * Psi^{-6} / sqrt( 1 + S_fluid^{2}/( (rho_star*h)^{2} ) ) |\n",
      ">        * .---------------------------------------------------------------------------.\n",
      ">        */\n",
      ">       rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n",
      "> \n",
      ">       /* Update j1 */\n",
      ">       j1 = find_polytropic_K_and_Gamma_index(eos,rhob1);\n",
      "> \n",
      ">     }  while( fabs(rhob1-rhob0) > rhob1*tol || j1 != j0);\n",
      ">     //////////////////////\n",
      ">     //  INNER LOOP END  //\n",
      ">     //////////////////////\n",
      "> \n",
      ">     /* Output the last value of rhob */\n",
      ">     rhob_out = rhob1;\n",
      "> \n",
      ">     /* Perform physical checks on the variables\n",
      ">      * and output the last value of h obtained\n",
      ">      */\n",
      ">     compute_P_cold__eps_cold(eos,rhob_out, P_cold, eps_cold);\n",
      ">     h = 1.0 + eps_cold + P_cold/rhob_out;\n",
      "> \n",
      ">     /* Set W based on eq. (A60) in Etienne et al. (2011)\n",
      ">      *       https://arxiv.org/pdf/1112.0568.pdf\n",
      ">      * .-------------------------------------------------------.\n",
      ">      * | W = psi^{-6} * sqrt( S_fluid^{2} + (rho_star*h)^{2} ) |\n",
      ">      * .-------------------------------------------------------.\n",
      ">      */\n",
      ">     W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n",
      "> \n",
      ">     /* Then update S_{fluid}^{2} using eq. (A61) in Etienne et al. (2011)\n",
      ">      *           https://arxiv.org/pdf/1112.0568.pdf\n",
      ">      * .---------------------------------------------------------------------------.\n",
      ">      * | S_fluid^{2} = ( W^{2}*S^{2} + (B.S)^2*(B^{2} + 2W) )/( ( W + B^{2} )^{2} )|\n",
      ">      * .---------------------------------------------------------------------------.\n",
      ">      */\n",
      ">     Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n",
      "86c182\n",
      "<   CCTK_REAL U_RHOB_inv = 1.0/U[RHOB];\n",
      "---\n",
      ">     if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n",
      "88,111d183\n",
      "<   if(eos.neos==1) {\n",
      "<     // Eq. 14 of http://arxiv.org/pdf/0802.0200.pdf :\n",
      "<     // P_{cold} = K_i rho_i^{\\Gamma_i}\n",
      "<     P_cold = eos.k_tab[0]*fasterpow_ppm_reconstruct(U[RHOB],eos.gamma_tab[0]);\n",
      "<     // Eq. 16 of http://arxiv.org/pdf/0802.0200.pdf :\n",
      "<     // \\epsilon_{cold} = \\int ( P_{cold}(rho) / rho^2 ) drho \n",
      "<     //                 = \\int ( K_0 \\rho^{\\Gamma_0 - 2} ) drho \n",
      "<     //                 = ( K_0 \\rho^{\\Gamma_0 - 1} ) / (\\Gamma_0 - 1)\n",
      "<     //                 = ( P_{cold} / rho          ) / (\\Gamma_0 - 1)\n",
      "<     eps_cold = P_cold*U_RHOB_inv/(eos.gamma_tab[0]-1.0);\n",
      "<     // dPcold/drho = K_i \\Gamma_i rho_i^{\\Gamma_i-1} = \\Gamma_i P_{cold} / rho\n",
      "<     dPcold_drho = eos.gamma_tab[0]*P_cold*U_RHOB_inv;\n",
      "<     // Eq. 15 of http://arxiv.org/pdf/0802.0200.pdf :\n",
      "<     // P_{th} = (\\Gamma_{th} - 1) \\rho_0 \\epsilon_{th},\n",
      "<     // Eq. 13 of http://arxiv.org/pdf/0802.0200.pdf :\n",
      "<     // P_{th} = P - P_{cold}\n",
      "<     // -> P - P_{cold} = (\\Gamma_{th} - 1) \\rho_0 \\epsilon_{th}\n",
      "<     // -> \\epsilon_{th} = ( P - P_{cold} ) / [ (\\Gamma_{th} - 1) \\rho_0 ]\n",
      "<     eps_th = (U[PRESSURE] - P_cold)/(eos.gamma_th-1.0)*U_RHOB_inv;\n",
      "<     // Just below Eq. 16 in http://arxiv.org/pdf/astro-ph/0503420.pdf :\n",
      "<     // h = 1 + \\epsilon + P/rho\n",
      "<     h = 1.0 + eps_cold + eps_th + U[PRESSURE]*U_RHOB_inv;\n",
      "<     gamma_cold = eos.gamma_tab[0];\n",
      "<     return;\n",
      "113,125c185,194\n",
      "< \n",
      "<   // See comments above for the eos.neos==1 case for relevant \n",
      "<   // equations & references; the extension to arbitrary \"nn\"\n",
      "<   // is straightforward.\n",
      "<   for(int nn=1;nn<eos.neos;nn++) {\n",
      "<     if (U[RHOB] <= eos.rho_tab[nn]  &&  U[RHOB] > eos.rho_tab[nn-1]) {\n",
      "<       P_cold = eos.k_tab[nn]*fasterpow_ppm_reconstruct(U[RHOB],eos.gamma_tab[nn]);\n",
      "<       eps_cold = eos.eps_tab[nn-1] + (P_cold*U_RHOB_inv - eos.P_tab[nn-1]/eos.rho_tab[nn-1])/(eos.gamma_tab[nn]-1.0);\n",
      "<       dPcold_drho = eos.gamma_tab[nn]*P_cold*U_RHOB_inv;\n",
      "<       eps_th = (U[PRESSURE] - P_cold)/(eos.gamma_th-1.0)*U_RHOB_inv;\n",
      "<       h = 1.0 + eps_cold + eps_th + U[PRESSURE]*U_RHOB_inv;\n",
      "<       gamma_cold = eos.gamma_tab[nn];\n",
      "<     }\n",
      "---\n",
      ">   //////////////////////\n",
      ">   //  OUTER LOOP END  //\n",
      ">   //////////////////////\n",
      "> \n",
      ">   /* If the code converged before the max\n",
      ">    * number of iterations were exceeded,\n",
      ">    * return 0, otherwise return 1.\n",
      ">    */\n",
      ">   if(fontcheck || itcount >= maxits) {\n",
      ">     return 1;\n",
      "127,133c196,197\n",
      "<   if (U[RHOB] > eos.rho_tab[eos.neos-1]) {\n",
      "<     P_cold = eos.k_tab[eos.neos]*fasterpow_ppm_reconstruct(U[RHOB],eos.gamma_tab[eos.neos]);\n",
      "<     eps_cold = eos.eps_tab[eos.neos-1] + (P_cold*U_RHOB_inv - eos.P_tab[eos.neos-1]/eos.rho_tab[eos.neos-1])/(eos.gamma_tab[eos.neos]-1.0);\n",
      "<     dPcold_drho = eos.gamma_tab[eos.neos]*P_cold*U_RHOB_inv;\n",
      "<     eps_th = (U[PRESSURE] - P_cold)/(eos.gamma_th-1.0)*U_RHOB_inv;\n",
      "<     h = 1.0 + eps_cold + eps_th + U[PRESSURE]*U_RHOB_inv;\n",
      "<     gamma_cold = eos.gamma_tab[eos.neos];\n",
      "---\n",
      ">   else {\n",
      ">     return 0;\n",
      "136a201\n",
      "> \n",
      "149a215\n",
      "> \n",
      "150a217,218\n",
      "> \n",
      "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
      "151a220,221\n",
      "> #endif\n",
      "> \n",
      "154c224\n",
      "<   // \\gamma_{ij} (v^i + \\beta^i)(v^j + \\beta^j)/(\\alpha)^2 \n",
      "---\n",
      ">   // \\gamma_{ij} (v^i + \\beta^i)(v^j + \\beta^j)/(\\alpha)^2\n",
      "165a236\n",
      "> \n",
      "176a248\n",
      "> \n",
      "187,195c259,267\n",
      "< // one_over_alpha_u0  = sqrt(1.0-one_minus_one_over_alpha_u0_squared);  \n",
      "< /* Proof of following line: */                                          \n",
      "< /*   [        1-1/(alphau0)^2         ] / [ 1/(alphau0) (1 + 1/(alphau0)) ] */ \n",
      "< /* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [ 1/(alphau0) + 1/(alphau0)^2   ] */ \n",
      "< /* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [  (alphau0 + 1)/(alphau0)^2    ] */ \n",
      "< /* = [        (alphau0)^2 - 1)        ] / [        (alphau0 + 1)          ] */  \n",
      "< /*   [   (alphau0 + 1) (alphau0 - 1)  ] / [        (alphau0 + 1)          ] */ \n",
      "< /* =                              alphau0 - 1                               */ \n",
      "< //alpha_u0_minus_one = one_minus_one_over_alpha_u0_squared/one_over_alpha_u0/(1.0+one_over_alpha_u0); \n",
      "---\n",
      "> // one_over_alpha_u0  = sqrt(1.0-one_minus_one_over_alpha_u0_squared);\n",
      "> /* Proof of following line: */\n",
      "> /*   [        1-1/(alphau0)^2         ] / [ 1/(alphau0) (1 + 1/(alphau0)) ] */\n",
      "> /* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [ 1/(alphau0) + 1/(alphau0)^2   ] */\n",
      "> /* = [ (alphau0)^2 - 1)/((alphau0)^2) ] / [  (alphau0 + 1)/(alphau0)^2    ] */\n",
      "> /* = [        (alphau0)^2 - 1)        ] / [        (alphau0 + 1)          ] */\n",
      "> /*   [   (alphau0 + 1) (alphau0 - 1)  ] / [        (alphau0 + 1)          ] */\n",
      "> /* =                              alphau0 - 1                               */\n",
      "> //alpha_u0_minus_one = one_minus_one_over_alpha_u0_squared/one_over_alpha_u0/(1.0+one_over_alpha_u0);\n",
      "197a270\n",
      "> \n",
      "212a286\n",
      "> \n",
      "224c298\n",
      "<     \n",
      "---\n",
      "> \n",
      "235c309,310\n",
      "< static inline void compute_smallba_b2_and_u_i_over_u0_psi4(CCTK_REAL *METRIC,CCTK_REAL *METRIC_LAP_PSI4,CCTK_REAL *U,CCTK_REAL u0L,CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI,  \n",
      "---\n",
      "> \n",
      "> static inline void compute_smallba_b2_and_u_i_over_u0_psi4(CCTK_REAL *METRIC,CCTK_REAL *METRIC_LAP_PSI4,CCTK_REAL *U,CCTK_REAL u0L,CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI,\n",
      "252a328\n",
      "> \n",
      "265a342\n",
      "> \n",
      "276c353\n",
      "<   CCTK_REAL bz_plus_shiftz_bt = smallb[SMALLBZ]+METRIC[SHIFTZ]*smallb[SMALLBT];    \n",
      "---\n",
      ">   CCTK_REAL bz_plus_shiftz_bt = smallb[SMALLBZ]+METRIC[SHIFTZ]*smallb[SMALLBT];\n",
      "283a361\n",
      "> \n"
     ]
    }
   ],
   "source": [
    "# Verify if the code generated by this tutorial module\n",
    "# matches the original IllinoisGRMHD source code\n",
    "\n",
    "# First download the original IllinoisGRMHD source code\n",
    "import urllib\n",
    "from os import path\n",
    "\n",
    "original_IGM_file_url  = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/inlined_functions.C\"\n",
    "original_IGM_file_name = \"inlined_functions-original.C\"\n",
    "original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
    "\n",
    "# Then download the original IllinoisGRMHD source code\n",
    "# We try it here in a couple of ways in an attempt to keep\n",
    "# the code more portable\n",
    "try:\n",
    "    original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "    # Write down the file the original IllinoisGRMHD source code\n",
    "    with open(original_IGM_file_path,\"w\") as file:\n",
    "        file.write(original_IGM_file_code)\n",
    "except:\n",
    "    try:\n",
    "        original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "        # Write down the file the original IllinoisGRMHD source code\n",
    "        with open(original_IGM_file_path,\"w\") as file:\n",
    "            file.write(original_IGM_file_code)\n",
    "    except:\n",
    "        # If all else fails, hope wget does the job\n",
    "        !wget -O $original_IGM_file_path $original_IGM_file_url\n",
    "\n",
    "# Perform validation\n",
    "Validation__inlined_functions__C  = !diff $original_IGM_file_path $outfile_path__inlined_functions__C\n",
    "\n",
    "if Validation__inlined_functions__C == []:\n",
    "    # If the validation passes, we do not need to store the original IGM source code file\n",
    "    !rm $original_IGM_file_path\n",
    "    print(\"Validation test for inlined_functions.C: PASSED!\")\n",
    "else:\n",
    "    # If the validation fails, we keep the original IGM source code file\n",
    "    print(\"Validation test for inlined_functions.C: FAILED!\")\n",
    "    # We also print out the difference between the code generated\n",
    "    # in this tutorial module and the original IGM source code\n",
    "    print(\"Diff:\")\n",
    "    for diff_line in Validation__inlined_functions__C:\n",
    "        print(diff_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 10: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-IllinoisGRMHD__inlined_functions.pdf](Tutorial-IllinoisGRMHD__inlined_functions.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
    "#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__inlined_functions.ipynb\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__inlined_functions.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__inlined_functions.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__inlined_functions.tex\n",
    "!rm -f Tut*.out Tut*.aux Tut*.log"
   ]
  }
 ],
 "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}