Tutorial-ETK_thorn-NRPyPlusTOVID.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",
    "# `NRPyPlusTOVID`: An Einstein Toolkit Thorn for Piecewise-Polytrope TOV neutron star initial data\n",
    "\n",
    "## Author: Zach Etienne and Leo Werneck\n",
    "\n",
    "## This notebook introduces the ETK module `NRPyPlusTOVID`, used for generating initial data for equilibrium neutron stars. The module employs NRPy+ for conversion from SymPy expressions to C-code kernel. It outlines methods for transforming this initial data to conform to the Toolkit's standards, including conversion to spherical and Cartesian ADMBase and HydroBase variables.\n",
    "\n",
    "**Notebook Status:** <font color='orange'><b> Partially Validated </b></font>\n",
    "\n",
    "**Validation Notes:** NRPy+ TOV initial data generation module validated against [Josh Faber's TOV initial data solver](https://ccrg.rit.edu/~jfaber/BNSID/TOV/), as described in the [NRPy+ implementation notes of the TOV solution for piecewise-polytrope neutron stars](Tutorial-TOV-Piecewise_Polytrope_EOSs.ipynb).\n",
    "\n",
    "### NRPy+ Source Code for this module: [TOV/TOV_Solver.py](../edit/TOV/TOV_Solver.py) [\\[tutorial\\]](Tutorial-Tutorial-ADM_Initial_Data-TOV.ipynb) Constructs numerical solution to TOV equations for neutron stars with piecewise polytrope equations of state\n",
    "\n",
    "## Introduction:\n",
    "In this part of the tutorial, we will construct an Einstein Toolkit (ETK) thorn (module) that will set up [TOV initial data](https://en.wikipedia.org/wiki/Tolman–Oppenheimer–Volkoff_equation) for an equilibrium neutron star. As documented in the [Piecewise Polytrope NRPy+ tutorial](Tutorial-TOV-Piecewise_Polytrope_EOSs.ipynb), piecewise-polytrope equations of state are supported, which closely approximate realistic nuclear equations of state appropriate for neutron star matter. In the [Tutorial-Tutorial-ADM_Initial_Data-TOV](Tutorial-Tutorial-ADM_Initial_Data-TOV.ipynb) tutorial notebook, we used NRPy+ to construct the SymPy expressions for these initial data. \n",
    "\n",
    "We will construct this thorn in two steps.\n",
    "\n",
    "1. Call on NRPy+ to convert the SymPy expressions for the initial data into one C-code kernel.\n",
    "1. Write the C code and linkages to the Einstein Toolkit infrastructure (i.e., the .ccl files) to complete this Einstein Toolkit module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$ \n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): **Call on NRPy+ to generate the TOV solution given a piecewise-polytrope equation of state; output the data to a text file**\n",
    "1. [Step 2](#initial_data): **Converting TOV initial data so that it can be used by the Einstein Toolkit**\n",
    "    1. [Step 2.a](#initial_data__interpolation): Interpolate the TOV data file as needed\n",
    "    1. [Step 2.b](#initial_data__tov_to_adm_sph): Converting the TOV variables to ADM variables in Spherical coordinates\n",
    "    1. [Step 2.c](#initial_data__convert_adm_sph_to_admbase): Convert Spherical ADM quantities to `ADMBase` (Cartesian) variables $\\left\\{\\alpha,\\beta^i,\\gamma_{ij},K_{ij}\\right\\}$\n",
    "    1. [Step 2.d](#initial_data__convert_to_hydrobase): Convert TOV solution quantities to `HydroBase` variables $\\left\\{P,\\rho_{\\rm baryonic},\\epsilon,v_{\\rm (n)}^i\\right\\}$\n",
    "1. [Step 3](#einstein): **Interfacing with the Einstein Toolkit**\n",
    "    1. [Step 3.a](#einstein_c): Constructing the Einstein Toolkit C-code calling functions that include the C code kernels\n",
    "    1. [Step 3.b](#einstein_ccl): CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure\n",
    "    1. [Step 3.c](#einstein_list): Add the C code to the Einstein Toolkit compilation list\n",
    "1. [Step 4](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Call on NRPy+ to generate the TOV solution given a piecewise-polytrope equation of state; output the data to a text file \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:50.679110Z",
     "iopub.status.busy": "2021-03-07T17:08:50.677769Z",
     "iopub.status.idle": "2021-03-07T17:08:51.497714Z",
     "shell.execute_reply": "2021-03-07T17:08:51.498165Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1256 1256 1256 1256 1256 1256\n",
      "Just generated a TOV star with\n",
      "* M        = 1.405031497682765e-01 ,\n",
      "* R_Schw   = 9.566039886703138e-01 ,\n",
      "* R_iso    = 8.100079558158075e-01 ,\n",
      "* M/R_Schw = 1.468770268913230e-01 \n",
      "\n",
      "Single Polytrope TOV solution generated in: 0.19237208366394043 s\n",
      "Initial data file: outputTOVpolytrope.txt\n"
     ]
    }
   ],
   "source": [
    "# Step 1: Import needed core NRPy+ modules\n",
    "from outputC import lhrh         # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os                # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Create directory for NRPyPlusTOVID thorn & subdirectories in case they don't exist.\n",
    "outrootdir = \"NRPyPlusTOVID/\"\n",
    "cmd.mkdir(os.path.join(outrootdir))\n",
    "outdir = os.path.join(outrootdir,\"src\") # Main C code output directory\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step 1.a: This is an Einstein Toolkit (ETK) thorn. Here we\n",
    "#           tell NRPy+ that gridfunction memory access will\n",
    "#           therefore be in the \"ETK\" style.\n",
    "par.set_parval_from_str(\"grid::GridFuncMemAccess\",\"ETK\")\n",
    "par.set_parval_from_str(\"grid::DIM\", 3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 1.b: NRPyPlusTOVID uses Cartesian coordinates, so\n",
    "#           we specify the reference metric to be Cartesian here:\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "# ADJUST THIS PARAMETER IF DESIRED.\n",
    "# \"Single\" = Single Polytrope\n",
    "# \"APR4\"   = APR4 Piecewise Polytrope\n",
    "# \"SLY\"    = SLy  Piecewise Polytrope\n",
    "# .---------------------------------------------.\n",
    "# | For all available names please look in the  |\n",
    "# | TOV/Piecewise_Polytrope__dict.py NRPy+ file |\n",
    "# .---------------------------------------------.\n",
    "# vvvvvvvvvvvvvvvv\n",
    "EOSname = \"Single\"\n",
    "# EOSname = \"SLy\"\n",
    "# EOSname = \"APR4\"\n",
    "# ^^^^^^^^^^^^^^^^\n",
    "\n",
    "# Import our TOV solver, which supports both single\n",
    "# and piecewise polytropic EOSs\n",
    "import TOV.TOV_Solver as TOV\n",
    "import TOV.Polytropic_EOSs as poly\n",
    "\n",
    "if EOSname==\"Single\":\n",
    "    # Set neos = 1 (single polytrope)\n",
    "    neos = 1\n",
    "\n",
    "    # Set rho_poly_tab (not needed for a single polytrope)\n",
    "    rho_poly_tab = []\n",
    "\n",
    "    # Set Gamma_poly_tab\n",
    "    Gamma_poly_tab = [2.0]\n",
    "\n",
    "    # Set K_poly_tab0\n",
    "    K_poly_tab0 = 1. # ZACH NOTES: CHANGED FROM 100.\n",
    "\n",
    "    rhob_central = 0.129285309 # M/R_Schw = 1.468770268913230e-01\n",
    "\n",
    "    # Set the eos quantities\n",
    "    eos = poly.set_up_EOS_parameters__complete_set_of_input_variables(neos,rho_poly_tab,Gamma_poly_tab,K_poly_tab0)\n",
    "    import time\n",
    "    start = time.time()\n",
    "    TOV.TOV_Solver(eos,\n",
    "                   outfile=\"outputTOVpolytrope.txt\",\n",
    "                   rho_baryon_central=rhob_central,\n",
    "                   verbose = True)\n",
    "    print(\"Single Polytrope TOV solution generated in: \"+str(time.time()-start)+\" s\")\n",
    "    print(\"Initial data file: outputTOVpolytrope.txt\")\n",
    "else:\n",
    "    # Set up the EOS parameters\n",
    "    eos = poly.set_up_EOS_parameters__Read_et_al_input_variables(EOSname)\n",
    "\n",
    "    # Set up the initial condition for the pressure by\n",
    "    # selecting a central baryon density\n",
    "#     rhob_central = 2.0 # M/R_Schw = 3.303692404611947e-01\n",
    "#     rhob_central = 1.0 # M/R_Schw = 2.051637558540178e-01\n",
    "    rhob_central = 0.8 # M/R_Schw = 1.470662481999595e-01\n",
    "\n",
    "    # Solve the TOV equations given our EOS and central density\n",
    "    import time\n",
    "    start = time.time()\n",
    "    outfilename = \"outputTOVpolytrope-\"+EOSname+\".txt\"\n",
    "    TOV.TOV_Solver(eos,outfile=outfilename,rho_baryon_central=rhob_central,verbose=True)\n",
    "    print(\"PPEOS \"+EOSname+\" TOV solution generated in: \"+str(time.time()-start)+\" s\")\n",
    "    print(\"Initial data file: \"+outfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data'></a>\n",
    "\n",
    "# Step 2: Converting TOV initial data so that it can be used by the Einstein Toolkit \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data}$$\n",
    "\n",
    "Main driver function:\n",
    "\n",
    "* Looping over all gridpoints:\n",
    "    * Read in `const CCTK_REAL rr = r[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)];`\n",
    "    * **Given this radius call interpolation driver to get all the base TOV quantities**\n",
    "    * **Convert TOV spacetime quantities to ADM quantities in *spherical* basis**\n",
    "    * Call the Cartesian ADMBase converter\n",
    "    * Call the HydroBase converter\n",
    "\n",
    "<a id='initial_data__interpolation'></a>\n",
    "\n",
    "## Step 2.a: Interpolate the TOV data file as needed \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data__interpolation}$$\n",
    "\n",
    "We start by interpolating the TOV data file to the gridpoints used by ETK, using the [tov_interp.h](../edit/TOV/tov_interp.h) file, which using Lagrange polynomial interpolation (for more details on the usage of this interpolation file, please look at the [start-to-finish TOV initial data tutorial notebook](Tutorial-Start_to_Finish-BSSNCurvilinear-TOV_initial_data.ipynb)).\n",
    "\n",
    "Keep in mind that the TOV data file just written stored $\\left(r,\\rho(r),\\rho_{\\text{baryonic}}(r),P(r),M(r),e^{\\nu(r)}\\right)$, where $\\rho(r)$ is the total mass-energy density (cf. $\\rho_{\\text{baryonic}}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:51.504730Z",
     "iopub.status.busy": "2021-03-07T17:08:51.503515Z",
     "iopub.status.idle": "2021-03-07T17:08:51.508175Z",
     "shell.execute_reply": "2021-03-07T17:08:51.507439Z"
    }
   },
   "outputs": [],
   "source": [
    "shutil.copy(os.path.join(\"TOV\",\"tov_interp.h\"),outdir)\n",
    "\n",
    "with open(os.path.join(outdir,\"interpolate_TOV_solution_to_point.h\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "/* Load the TOV_interpolate_1D() function */\n",
    "#include \"tov_interp.h\"\n",
    "\n",
    "/* This function returns the TOV quantities at point rr\n",
    " * by interpolating the data in the TOV initial data file.\n",
    " */\n",
    "void interpolate_TOV_solution_to_point(const CCTK_REAL rr, ID_inputs other_inputs,\n",
    "                           CCTK_REAL *exp_4phi, CCTK_REAL *expnu,\n",
    "                           CCTK_REAL *Pressure, CCTK_REAL *rho_baryon, CCTK_REAL *rho__total_energy_density) {\n",
    "\n",
    "  /* The mass valus is not used, but we have to\n",
    "   * store it in this dummy variable because the\n",
    "   * initial data file contains it.\n",
    "   */\n",
    "  CCTK_REAL M;\n",
    "\n",
    "  /* Perform the interpolation, returning:\n",
    "   *  - rho__total_energy_density\n",
    "   *  - rho_baryon\n",
    "   *  - Pressure\n",
    "   *  - Mass (dummy variable, unused)\n",
    "   *  - exp(nu)\n",
    "   *  - exp(4phi)\n",
    "   */\n",
    "  TOV_interpolate_1D(rr,other_inputs.Rbar,other_inputs.Rbar_idx,other_inputs.interp_stencil_size,\n",
    "                     other_inputs.numlines_in_file,\n",
    "                     other_inputs.r_Schw_arr,other_inputs.rho_arr,other_inputs.rho_baryon_arr,other_inputs.P_arr,other_inputs.M_arr,\n",
    "                     other_inputs.expnu_arr,other_inputs.exp4phi_arr,other_inputs.rbar_arr,\n",
    "                     rho__total_energy_density,rho_baryon,Pressure,&M,expnu,exp_4phi);\n",
    "\n",
    "}\\n\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data__tov_to_adm_sph'></a>\n",
    "\n",
    "## Step 2.b: Converting the TOV variables to ADM variables in Spherical coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data__tov_to_adm_sph}$$\n",
    "\n",
    "Now we perform the interpolation of the TOV quantities to ADM quantities in spherical coordinates, using (see [the TOV initial data tutorial notebook](Tutorial-ADM_Initial_Data-TOV.ipynb) for more details):\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\boxed{\n",
    "\\begin{aligned}\n",
    "\\alpha &= e^{\\nu(\\bar{r})/2} \\\\\n",
    "\\beta^k &= 0 \\\\\n",
    "\\gamma_{\\bar{r}\\bar{r}} &= e^{4\\phi}\\\\\n",
    "\\gamma_{\\theta\\theta} &= e^{4\\phi} \\bar{r}^2 \\\\\n",
    "\\gamma_{\\phi\\phi} &= e^{4\\phi} \\bar{r}^2 \\sin^2 \\theta \\\\\n",
    "\\end{aligned}\n",
    "}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:51.515653Z",
     "iopub.status.busy": "2021-03-07T17:08:51.514623Z",
     "iopub.status.idle": "2021-03-07T17:08:51.517412Z",
     "shell.execute_reply": "2021-03-07T17:08:51.517995Z"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(outdir,\"convert_TOV_spacetime_vars_to_ADM_vars.h\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "/* This function converts TOV quantities into\n",
    " * ADM quantities in Spherical coordinates.\n",
    " */\n",
    "void convert_TOV_spacetime_vars_to_ADM_vars( const CCTK_REAL rr, const CCTK_REAL th,\n",
    "                           const CCTK_REAL IDexp_4phi, const CCTK_REAL IDexpnu,\n",
    "                           CCTK_REAL *IDalpha,\n",
    "                           CCTK_REAL *IDgammaDD00, CCTK_REAL *IDgammaDD01, CCTK_REAL *IDgammaDD02,\n",
    "                           CCTK_REAL *IDgammaDD11, CCTK_REAL *IDgammaDD12, CCTK_REAL *IDgammaDD22) {\n",
    "\n",
    "  /***************************************************************\n",
    "   * Convert TOV quantities to ADM quantities in Spherical basis *\n",
    "   ***************************************************************\n",
    "   *\n",
    "   * First we convert the lapse function:\n",
    "   * .------------------.\n",
    "   * | alpha = e^(nu/2) |\n",
    "   * .------------------.\n",
    "   */\n",
    "   *IDalpha = sqrt(IDexpnu);\n",
    "\n",
    "  /* Next we convert the metric function:\n",
    "   * .----------------------------------------.\n",
    "   * | gamma_{00} = e^{4phi}                  |\n",
    "   * .----------------------------------------.\n",
    "   * | gamma_{11} = e^{4phi} r^2              |\n",
    "   * .----------------------------------------.\n",
    "   * | gamma_{22} = e^{4phi} r^2 sin^2(theta) |\n",
    "   * .----------------------------------------.\n",
    "   * | All other components are zero.         |\n",
    "   * .----------------------------------------.\n",
    "   */\n",
    "   *IDgammaDD00 = IDexp_4phi;\n",
    "   *IDgammaDD11 = IDexp_4phi * rr * rr;\n",
    "   *IDgammaDD22 = IDexp_4phi * rr * rr * sin(th) * sin(th);\n",
    "   *IDgammaDD01 = 0.0;\n",
    "   *IDgammaDD02 = 0.0;\n",
    "   *IDgammaDD12 = 0.0;\n",
    "\n",
    "}\\n\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data__convert_adm_sph_to_admbase'></a>\n",
    "\n",
    "## Step 2.c: Convert Spherical ADM quantities to `ADMBase` (Cartesian) variables $\\left\\{\\alpha,\\beta^i,\\gamma_{ij},K_{ij}\\right\\}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data__convert_adm_sph_to_admbase}$$\n",
    "\n",
    "The [TOV line element](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) in *Schwarzschild coordinates* is written (in the $-+++$ form):\n",
    "$$\n",
    "ds^2 = - c^2 e^\\nu dt^2 + \\left(1 - \\frac{2GM}{rc^2}\\right)^{-1} dr^2 + r^2 d\\Omega^2.\n",
    "$$\n",
    "\n",
    "In *isotropic coordinates* with $G=c=1$ (i.e., the initial coordinate slicing and units we prefer to use), the ($-+++$ form) line element is written:\n",
    "$$\n",
    "ds^2 = - e^{\\nu} dt^2 + e^{4\\phi} \\left(d\\bar{r}^2 + \\bar{r}^2 d\\Omega^2\\right),\n",
    "$$\n",
    "where $\\phi$ here is the *conformal factor*.\n",
    "\n",
    "The ADM 3+1 line element for this diagonal metric in isotropic spherical coordinates is given by:\n",
    "$$\n",
    "ds^2 = (-\\alpha^2 + \\beta_k \\beta^k) dt^2 + \\gamma_{\\bar{r}\\bar{r}} d\\bar{r}^2 + \\gamma_{\\theta\\theta} d\\theta^2+ \\gamma_{\\phi\\phi} d\\phi^2,\n",
    "$$\n",
    "\n",
    "from which we can immediately read off the ADM quantities:\n",
    "\\begin{align}\n",
    "\\alpha &= e^{\\nu(\\bar{r})/2} \\\\\n",
    "\\beta^k &= 0 \\\\\n",
    "\\gamma_{\\bar{r}\\bar{r}} &= e^{4\\phi}\\\\\n",
    "\\gamma_{\\theta\\theta} &= e^{4\\phi} \\bar{r}^2 \\\\\n",
    "\\gamma_{\\phi\\phi} &= e^{4\\phi} \\bar{r}^2 \\sin^2 \\theta \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:51.526646Z",
     "iopub.status.busy": "2021-03-07T17:08:51.525935Z",
     "iopub.status.idle": "2021-03-07T17:08:51.528591Z",
     "shell.execute_reply": "2021-03-07T17:08:51.528089Z"
    }
   },
   "outputs": [],
   "source": [
    "thismodule   = __name__\n",
    "\n",
    "IDalpha   = par.Cparameters(\"REAL\", thismodule, \"IDalpha\", 1e300) # IDalpha must be set in C\n",
    "IDbetaU   = ixp.zerorank1() # beta^i is zero\n",
    "IDgammaDD = ixp.zerorank2()\n",
    "for i in range(3):\n",
    "    for j in range(i,3):\n",
    "        IDgammaDD[i][j] = par.Cparameters(\"REAL\", thismodule, \"IDgammaDD\"+str(i)+str(j), 1e300) # IDgammaDD must be set in C\n",
    "        IDgammaDD[j][i] = IDgammaDD[i][j]\n",
    "IDKDD     = ixp.zerorank2() # K_{ij} is zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this ETK module expects Cartesian coordinates, and the TOV solution above is in the spherical basis, we next perform the Jacobian transformations necessary to convert into the Cartesian basis:\n",
    "\n",
    "All ADM tensors and vectors are in the Spherical coordinate basis $x^i_{\\rm Sph} = (r,\\theta,\\phi)$, but we need them in the Cartesian coordinate basis $x^i_{\\rm Cart}=$`(xx0,xx1,xx2)` set by the `\"reference_metric::CoordSystem\"` variable. Empirically speaking, it is far easier to write `(x(xx0,xx1,xx2),y(xx0,xx1, xx2),z(xx0,xx1,xx2))` than the inverse, so we will compute the Jacobian matrix\n",
    "\n",
    "$$\n",
    "{\\rm Jac\\_dUSph\\_dDrfmUD[i][j]} = \\frac{\\partial x^i_{\\rm Sph}}{\\partial x^j_{\\rm Cart}},\n",
    "$$\n",
    "\n",
    "via exact differentiation (courtesy SymPy), and the inverse Jacobian\n",
    "$$\n",
    "{\\rm Jac\\_dUrfm\\_dDSphUD[i][j]} = \\frac{\\partial x^i_{\\rm Cart}}{\\partial x^j_{\\rm Sph}},\n",
    "$$\n",
    "\n",
    "using NRPy+'s `generic_matrix_inverter3x3()` function. \n",
    "\n",
    "In terms of these, the transformation of ADM tensors from Spherical to `\"reference_metric::CoordSystem==Cartesian\"` coordinates may be written:\n",
    "\n",
    "\\begin{align}\n",
    "\\gamma^{\\rm Cart}_{ij} &= \n",
    "\\frac{\\partial x^\\ell_{\\rm Cart}}{\\partial x^i_{\\rm Sph}}\n",
    "\\frac{\\partial x^m_{\\rm Cart}}{\\partial x^j_{\\rm Sph}} \\gamma^{\\rm Sph}_{\\ell m}\n",
    "\\end{align}\n",
    "\n",
    "Since $\\beta^i=K_{ij}=0$ in this case, and $\\alpha$ is not a tensor, only the above Jacobian transformation need be performed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:51.584176Z",
     "iopub.status.busy": "2021-03-07T17:08:51.546783Z",
     "iopub.status.idle": "2021-03-07T17:08:52.914454Z",
     "shell.execute_reply": "2021-03-07T17:08:52.913827Z"
    }
   },
   "outputs": [],
   "source": [
    "# Transform initial data to our coordinate system:\n",
    "# First compute Jacobian and its inverse\n",
    "drrefmetric__dx_0UDmatrix = sp.Matrix([[sp.diff(rfm.xxSph[0],rfm.xx[0]), sp.diff(rfm.xxSph[0],rfm.xx[1]), sp.diff(rfm.xxSph[0],rfm.xx[2])],\n",
    "                                       [sp.diff(rfm.xxSph[1],rfm.xx[0]), sp.diff(rfm.xxSph[1],rfm.xx[1]), sp.diff(rfm.xxSph[1],rfm.xx[2])],\n",
    "                                       [sp.diff(rfm.xxSph[2],rfm.xx[0]), sp.diff(rfm.xxSph[2],rfm.xx[1]), sp.diff(rfm.xxSph[2],rfm.xx[2])]])\n",
    "dx__drrefmetric_0UDmatrix = drrefmetric__dx_0UDmatrix.inv()\n",
    "\n",
    "# Declare as gridfunctions the final quantities we will output for the initial data\n",
    "alpha   = gri.register_gridfunctions(\"EVOL\",\"alpha\")\n",
    "betaU   = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"betaU\")\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"EVOL\",\"gammaDD\",\"sym01\")\n",
    "KDD     = ixp.register_gridfunctions_for_single_rank2(\"EVOL\",\"KDD\",\"sym01\")\n",
    "\n",
    "alpha = IDalpha # No Jacobian necessary!\n",
    "betaU = IDbetaU # Because beta^i = 0\n",
    "KDD   = IDKDD   # Because K_{ij} = 0\n",
    "\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        # Matrices are stored in row, column format, so (i,j) <-> (row,column)\n",
    "        gammaDD[i][j] = 0\n",
    "        for k in range(3):\n",
    "            for l in range(3):\n",
    "                gammaDD[i][j] += drrefmetric__dx_0UDmatrix[(k,i)]*drrefmetric__dx_0UDmatrix[(l,j)]*IDgammaDD[k][l]\n",
    "\n",
    "# -={ Spacetime quantities: Generate C code from expressions and output to file }=-\n",
    "ADMQuantities_to_print = [\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"alpha\"),rhs=alpha),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU0\"),rhs=betaU[0]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU1\"),rhs=betaU[1]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU2\"),rhs=betaU[2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD00\"),rhs=gammaDD[0][0]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD01\"),rhs=gammaDD[0][1]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD02\"),rhs=gammaDD[0][2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD11\"),rhs=gammaDD[1][1]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD12\"),rhs=gammaDD[1][2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD22\"),rhs=gammaDD[2][2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD00\"),rhs=KDD[0][0]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD01\"),rhs=KDD[0][1]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD02\"),rhs=KDD[0][2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD11\"),rhs=KDD[1][1]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD12\"),rhs=KDD[1][2]),\\\n",
    "                          lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD22\"),rhs=KDD[2][2])\n",
    "                          ]\n",
    "\n",
    "with open(os.path.join(outdir,\"ADMQuantities.h\"),\"w\") as file:\n",
    "    ADMQuantities_CcodeKernel = fin.FD_outputC(\"returnstring\",ADMQuantities_to_print,\n",
    "                                               params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "    file.write(\"\"\"\n",
    "static inline\n",
    "void ADMQuantities(const cGH* restrict const cctkGH, const CCTK_INT i0,const CCTK_INT i1,const CCTK_INT i2,\n",
    "\n",
    "                   const CCTK_REAL *restrict xx0GF,const CCTK_REAL *restrict xx1GF,const CCTK_REAL *restrict xx2GF,\n",
    "\n",
    "                   const CCTK_REAL IDalpha,\n",
    "                   const CCTK_REAL IDgammaDD00,const CCTK_REAL IDgammaDD01, const CCTK_REAL IDgammaDD02,\n",
    "                   const CCTK_REAL IDgammaDD11,const CCTK_REAL IDgammaDD12, const CCTK_REAL IDgammaDD22,\n",
    "\n",
    "                   CCTK_REAL *alphaGF,CCTK_REAL *betaU0GF,CCTK_REAL *betaU1GF,CCTK_REAL *betaU2GF,\n",
    "                   CCTK_REAL *gammaDD00GF, CCTK_REAL *gammaDD01GF, CCTK_REAL *gammaDD02GF,\n",
    "                   CCTK_REAL *gammaDD11GF, CCTK_REAL *gammaDD12GF, CCTK_REAL *gammaDD22GF,\n",
    "                   CCTK_REAL *KDD00GF, CCTK_REAL *KDD01GF, CCTK_REAL *KDD02GF,\n",
    "                   CCTK_REAL *KDD11GF, CCTK_REAL *KDD12GF, CCTK_REAL *KDD22GF) {\n",
    "    const CCTK_REAL xx0 = xx0GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)];\n",
    "    const CCTK_REAL xx1 = xx1GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)];\n",
    "    const CCTK_REAL xx2 = xx2GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)];\n",
    "\n",
    "\"\"\"+ADMQuantities_CcodeKernel+\"\"\"\n",
    "}\n",
    "    \"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data__convert_to_hydrobase'></a>\n",
    "\n",
    "## Step 2.d: Convert TOV solution quantities to `HydroBase` variables $\\left\\{P,\\rho_{\\rm baryonic},\\epsilon,v_{\\rm (n)}^i\\right\\}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data__convert_to_hydrobase}$$\n",
    "\n",
    "\n",
    "The TOV solver outputs pressure $P$, the *total* energy density $\\rho$, and the baryonic density $\\rho_{\\rm baryonic}$ as a function of the stellar radius (in isotropic coordinates by default). \n",
    "\n",
    "Then, the `HydroBase` quantities $\\rho^{\\rm HB}_{\\rm baryonic}$, internal energy $\\epsilon^{\\rm HB}$, and pressure $P^{\\rm HB}$ are given in terms of these variables via\n",
    "\n",
    "\\begin{align}\n",
    "P^{\\rm HB} &= P; \\\\\n",
    "\\rho^{\\rm HB}_{\\rm baryonic} &= \\rho_{\\rm baryonic}, \\\\\n",
    "\\rho &= \\rho_{\\rm baryonic} \\left(1 + \\epsilon_{\\rm cold}\\right) \\\\\n",
    "\\implies \\epsilon_{\\rm cold} &= \\frac{\\rho}{\\rho_{\\rm baryonic}} - 1\\\\\n",
    "\\epsilon^{\\rm HB} &= \\epsilon_{\\rm cold}, \\\\\n",
    "\\end{align}\n",
    "[the NRPy+ piecewise polytrope tutorial notebook](Tutorial-TOV-Piecewise_Polytrope_EOSs.ipynb#rhob_from_pcold). Note that $\\rho_{\\rm baryonic}$ will be floored to a nonzero atmosphere value, so that computing $\\epsilon$ will never involve a division by zero.\n",
    "\n",
    "The TOV star is motionless, with all spatial components of the 4-velocity $u^i=0$ and (as seen above) zero shift $\\beta^i$. Thus the Valencia 3-velocity (i.e., the 3-velocity normal to the spatial slice) $v_{\\rm (n)}^i$ is given by\n",
    "\n",
    "$$\n",
    "v_{\\rm (n)}^{i,{\\rm HB}} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.934126Z",
     "iopub.status.busy": "2021-03-07T17:08:52.925556Z",
     "iopub.status.idle": "2021-03-07T17:08:52.936582Z",
     "shell.execute_reply": "2021-03-07T17:08:52.937095Z"
    }
   },
   "outputs": [],
   "source": [
    "IDValencia3velocityU = ixp.zerorank1() # Valencia 3-velocity is zero\n",
    "IDPressure     = par.Cparameters(\"REAL\", thismodule, \"IDPressure\", 1e300) # IDPressure must be set in C\n",
    "IDrho_baryonic = par.Cparameters(\"REAL\", thismodule, \"IDrho_baryonic\", 1e300) # IDrho_baryonic must be set in C\n",
    "IDrho__total_energy_density = par.Cparameters(\"REAL\", thismodule, \"IDrho__total_energy_density\", 1e300) # IDrho__total_energy_density must be set in C\n",
    "\n",
    "# Declare as gridfunctions the final quantities we will output for the initial data\n",
    "Valencia3velocityU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"Valencia3velocityU\")\n",
    "Pressure, rho_baryonic, epsilon = gri.register_gridfunctions(\"EVOL\",[\"Pressure\", \"rho_baryonic\", \"epsilon\"])\n",
    "\n",
    "Valencia3velocityU = IDValencia3velocityU # Because all components of Valencia3velocityU are *zero*\n",
    "Pressure     = IDPressure\n",
    "rho_baryonic = IDrho_baryonic\n",
    "epsilon      = IDrho__total_energy_density / IDrho_baryonic - sp.sympify(1)\n",
    "\n",
    "# -={ Spacetime quantities: Generate C code from expressions and output to file }=-\n",
    "HydroQuantities_to_print = [\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Pressure\"),rhs=Pressure),\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"rho_baryonic\"),rhs=rho_baryonic),\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"epsilon\"),rhs=epsilon),\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU0\"),rhs=Valencia3velocityU[0]),\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU1\"),rhs=Valencia3velocityU[1]),\\\n",
    "                            lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU2\"),rhs=Valencia3velocityU[2])\n",
    "                           ]\n",
    "\n",
    "with open(os.path.join(outdir,\"HydroQuantities.h\"),\"w\") as file:\n",
    "    HydroQuantities_CcodeKernel = fin.FD_outputC(\"returnstring\",HydroQuantities_to_print,\n",
    "                                               params=\"outCverbose=False,includebraces=False,preindent=2\")\n",
    "    file.write(\"\"\"\n",
    "static inline\n",
    "void HydroQuantities(const cGH* restrict const cctkGH, const CCTK_INT i0,const CCTK_INT i1,const CCTK_INT i2,\n",
    "\n",
    "                   const CCTK_REAL IDPressure, const CCTK_REAL IDrho_baryonic,\n",
    "                   const CCTK_REAL IDrho__total_energy_density,\n",
    "\n",
    "                   CCTK_REAL *PressureGF,CCTK_REAL *rho_baryonicGF,\n",
    "                   CCTK_REAL *epsilonGF,\n",
    "                   CCTK_REAL *Valencia3velocityU0GF,\n",
    "                   CCTK_REAL *Valencia3velocityU1GF,\n",
    "                   CCTK_REAL *Valencia3velocityU2GF) {\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "    if(IDrho__total_energy_density <= 0 || IDrho_baryonic <= 0 || IDPressure <= 0) {\n",
    "        rho_baryonicGF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)] = rho_atmosphere;\n",
    "        PressureGF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)]     = K_atmosphere*pow(rho_atmosphere,Gamma_atmosphere);\n",
    "        epsilonGF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)]      = 0;\n",
    "        Valencia3velocityU0GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)] = 0;\n",
    "        Valencia3velocityU1GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)] = 0;\n",
    "        Valencia3velocityU2GF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)] = 0;\n",
    "    } else {\n",
    "\"\"\"+HydroQuantities_CcodeKernel+\"\"\"\n",
    "        // Apply pressure depletion.\n",
    "        PressureGF[CCTK_GFINDEX3D(cctkGH, i0, i1, i2)] *= (1.0 - Pressure_depletion_factor);\n",
    "    }\n",
    "}\n",
    "    \"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein'></a>\n",
    "\n",
    "# Step 3: Interfacing with the Einstein Toolkit \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_c'></a>\n",
    "\n",
    "## Step 3.a: Constructing the Einstein Toolkit C-code calling functions that include the C code kernels \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein_c}$$\n",
    "\n",
    "We will write another C file with the functions we need here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.945231Z",
     "iopub.status.busy": "2021-03-07T17:08:52.944081Z",
     "iopub.status.idle": "2021-03-07T17:08:52.948471Z",
     "shell.execute_reply": "2021-03-07T17:08:52.947628Z"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(outdir,\"InitialData.c\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "#include <math.h>\n",
    "#include <stdio.h>\n",
    "#include <stdbool.h>\n",
    "\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "\n",
    "// Declare initial data input struct:\n",
    "//          stores data from initial data solver,\n",
    "//          so they can be put on the numerical grid.\n",
    "typedef struct __ID_inputs {\n",
    "    CCTK_REAL Rbar;\n",
    "    int Rbar_idx;\n",
    "    int interp_stencil_size;\n",
    "    int numlines_in_file;\n",
    "    CCTK_REAL *r_Schw_arr,*rho_arr,*rho_baryon_arr,*P_arr,*M_arr,*expnu_arr,*exp4phi_arr,*rbar_arr;\n",
    "} ID_inputs;\n",
    "\n",
    "\n",
    "#include \"ADMQuantities.h\"\n",
    "#include \"HydroQuantities.h\"\n",
    "#include \"interpolate_TOV_solution_to_point.h\"\n",
    "#include \"convert_TOV_spacetime_vars_to_ADM_vars.h\"\n",
    "\n",
    "// Alias for \"vel\" vector gridfunction:\n",
    "#define velx (&vel[0*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]])\n",
    "#define vely (&vel[1*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]])\n",
    "#define velz (&vel[2*cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]])\n",
    "\n",
    "void read_TOV_input_data_from_file(ID_inputs *TOV_in) {\n",
    "\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "    // Step 1: Set up TOV initial data\n",
    "    // Step 1.a: Read TOV initial data from data file\n",
    "    // Open the data file:\n",
    "    char filename[100];\n",
    "    sprintf(filename,\"%s\",TOV_filename); // TOV_filename is a CCTK_PARAMETER\n",
    "    FILE *in1Dpolytrope = fopen(filename, \"r\");\n",
    "    if (in1Dpolytrope == NULL) {\n",
    "        fprintf(stderr,\"ERROR: could not open file %s\\\\n\",filename);\n",
    "        exit(1);\n",
    "    }\n",
    "    // Count the number of lines in the data file:\n",
    "    int numlines_in_file = count_num_lines_in_file(in1Dpolytrope);\n",
    "    // Allocate space for all data arrays:\n",
    "    CCTK_REAL *r_Schw_arr     = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *rho_arr        = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *rho_baryon_arr = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *P_arr          = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *M_arr          = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *expnu_arr      = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *exp4phi_arr    = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "    CCTK_REAL *rbar_arr       = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*numlines_in_file);\n",
    "\n",
    "    // Read from the data file, filling in arrays.\n",
    "    // read_datafile__set_arrays() may be found in TOV/tov_interp.h\n",
    "    if(read_datafile__set_arrays(in1Dpolytrope, r_Schw_arr,rho_arr,rho_baryon_arr,P_arr,M_arr,expnu_arr,exp4phi_arr,rbar_arr) == 1) {\n",
    "        fprintf(stderr,\"ERROR WHEN READING FILE %s!\\\\n\",filename);\n",
    "        exit(1);\n",
    "    }\n",
    "    fclose(in1Dpolytrope);\n",
    "    REAL Rbar = -100;\n",
    "    int Rbar_idx = -100;\n",
    "    for(int i=1;i<numlines_in_file;i++) {\n",
    "        if(rho_arr[i-1]>0 && rho_arr[i]==0) { Rbar = rbar_arr[i-1]; Rbar_idx = i-1; }\n",
    "    }\n",
    "    if(Rbar<0) {\n",
    "        fprintf(stderr,\"Error: could not find rbar=Rbar from data file.\\\\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    TOV_in->Rbar = Rbar;\n",
    "    TOV_in->Rbar_idx = Rbar_idx;\n",
    "\n",
    "    const int interp_stencil_size = 12;\n",
    "    TOV_in->interp_stencil_size = interp_stencil_size;\n",
    "    TOV_in->numlines_in_file = numlines_in_file;\n",
    "\n",
    "    TOV_in->r_Schw_arr     = r_Schw_arr;\n",
    "    TOV_in->rho_arr        = rho_arr;\n",
    "    TOV_in->rho_baryon_arr = rho_baryon_arr;\n",
    "    TOV_in->P_arr          = P_arr;\n",
    "    TOV_in->M_arr          = M_arr;\n",
    "    TOV_in->expnu_arr      = expnu_arr;\n",
    "    TOV_in->exp4phi_arr    = exp4phi_arr;\n",
    "    TOV_in->rbar_arr       = rbar_arr;\n",
    "    /* END TOV INPUT ROUTINE */\n",
    "}\n",
    "\n",
    "void NRPyPlusTOVID_ET_InitialData(CCTK_ARGUMENTS) {\n",
    "    DECLARE_CCTK_ARGUMENTS;\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "    ID_inputs TOV_in;\n",
    "    read_TOV_input_data_from_file(&TOV_in);\n",
    "\n",
    "#pragma omp parallel for\n",
    "  for(CCTK_INT i2=0;i2<cctk_lsh[2];i2++) for(CCTK_INT i1=0;i1<cctk_lsh[1];i1++) for(CCTK_INT i0=0;i0<cctk_lsh[0];i0++) {\n",
    "        CCTK_INT idx = CCTK_GFINDEX3D(cctkGH,i0,i1,i2);\n",
    "        CCTK_REAL rr = r[idx];\n",
    "        CCTK_REAL th = acos(z[idx]/rr);\n",
    "\n",
    "        CCTK_REAL IDexp_4phi,IDnu,IDPressure,IDrho_baryonic,IDrho__total_energy_density;\n",
    "        interpolate_TOV_solution_to_point(rr, TOV_in, &IDexp_4phi,&IDnu,\n",
    "                                          &IDPressure,&IDrho_baryonic,&IDrho__total_energy_density);\n",
    "\n",
    "        CCTK_REAL IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22;\n",
    "        convert_TOV_spacetime_vars_to_ADM_vars(rr, th, IDexp_4phi,IDnu,\n",
    "            &IDalpha,&IDgammaDD00,&IDgammaDD01,&IDgammaDD02,&IDgammaDD11,&IDgammaDD12,&IDgammaDD22);\n",
    "\n",
    "        HydroQuantities(cctkGH, i0,i1,i2,\n",
    "                        IDPressure,IDrho_baryonic,IDrho__total_energy_density,\n",
    "                        press,rho,eps,velx,vely,velz);\n",
    "\n",
    "        ADMQuantities(cctkGH,i0,i1,i2,\n",
    "                      x,y,z,\n",
    "                      IDalpha,IDgammaDD00,IDgammaDD01,IDgammaDD02,IDgammaDD11,IDgammaDD12,IDgammaDD22,\n",
    "                      alp,betax,betay,betaz,\n",
    "                      gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "                      kxx,kxy,kxz,kyy,kyz,kzz);\n",
    "  }\n",
    "\n",
    "  free(TOV_in.r_Schw_arr);\n",
    "  free(TOV_in.rho_arr);\n",
    "  free(TOV_in.rho_baryon_arr);\n",
    "  free(TOV_in.P_arr);\n",
    "  free(TOV_in.M_arr);\n",
    "  free(TOV_in.expnu_arr);\n",
    "  free(TOV_in.exp4phi_arr);\n",
    "  free(TOV_in.rbar_arr);\n",
    "}\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_ccl'></a>\n",
    "\n",
    "## Step 3.b: CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein_ccl}$$\n",
    "\n",
    "Writing a module (\"thorn\") within the Einstein Toolkit requires that three \"ccl\" files be constructed, all in the root directory of the thorn:\n",
    "\n",
    "1. `interface.ccl}`: defines the gridfunction groups needed, and provides keywords denoting what this thorn provides and what it should inherit from other thorns. Specifically, this file governs the interaction between this thorn and others; more information can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-178000D2.2). \n",
    "With \"implements\", we give our thorn its unique name. By \"inheriting\" other thorns, we tell the Toolkit that we will rely on variables that exist and are declared \"public\" within those functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.953787Z",
     "iopub.status.busy": "2021-03-07T17:08:52.953057Z",
     "iopub.status.idle": "2021-03-07T17:08:52.956796Z",
     "shell.execute_reply": "2021-03-07T17:08:52.956267Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting NRPyPlusTOVID//interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outrootdir/interface.ccl\n",
    "implements: NRPyPlusTOVID\n",
    "inherits: admbase grid hydrobase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. `param.ccl`: specifies free parameters within the thorn, enabling them to be set at runtime. It is required to provide allowed ranges and default values for each parameter. More information on this file's syntax can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-183000D2.3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.961943Z",
     "iopub.status.busy": "2021-03-07T17:08:52.961162Z",
     "iopub.status.idle": "2021-03-07T17:08:52.964137Z",
     "shell.execute_reply": "2021-03-07T17:08:52.964583Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting NRPyPlusTOVID//param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outrootdir/param.ccl\n",
    "shares: grid\n",
    "shares: ADMBase\n",
    "USES CCTK_INT lapse_timelevels\n",
    "USES CCTK_INT shift_timelevels\n",
    "USES CCTK_INT metric_timelevels\n",
    "\n",
    "USES KEYWORD metric_type\n",
    "\n",
    "EXTENDS KEYWORD initial_data\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial data from NRPyPlusTOVID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_lapse\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial lapse from NRPyPlusTOVID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_shift\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial shift from NRPyPlusTOVID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_dtlapse\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial dtlapse from NRPyPlusTOVID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_dtshift\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial dtshift from NRPyPlusTOVID solution\"\n",
    "}\n",
    "\n",
    "shares: HydroBase\n",
    "EXTENDS KEYWORD initial_hydro\n",
    "{\n",
    "  \"NRPyPlusTOVID\" :: \"Initial GRHD data from NRPyPlusTOVID solution\"\n",
    "}\n",
    "\n",
    "#[\"r_in\",\"r_at_max_density\",\"a\",\"M\"] A_b, kappa, gamma\n",
    "restricted:\n",
    "CCTK_STRING TOV_filename \"Which interpolator should I use\"\n",
    "{\n",
    " \".+\" :: \"Any nonempty string\"\n",
    "} \"outputTOVpolytrope.txt\"\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL rho_atmosphere \"Atmosphere baryonic density\"\n",
    "{\n",
    " 0:* :: \"Physical values\"\n",
    "-1   :: \"forbidden value to make sure it is explicitly set in the parfile\"\n",
    "} -1\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL K_atmosphere \"Polytropic K to be used with the EOS corresponding to rho_atmosphere\"\n",
    "{\n",
    " 0:* :: \"Physical values\"\n",
    "-1   :: \"forbidden value to make sure it is explicitly set in the parfile\"\n",
    "} -1\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL Gamma_atmosphere \"Polytropic Gamma to be used with the EOS corresponding to rho_atmosphere\"\n",
    "{\n",
    " 0:* :: \"Physical values\"\n",
    "-1   :: \"forbidden value to make sure it is explicitly set in the parfile\"\n",
    "} -1\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL Pressure_depletion_factor \"Pressure depletion factor = Pdf: P => (1-Pdf)*P\"\n",
    "{\n",
    " 0:* :: \"Greater than or equal to zero, where zero is no depletion and default.\"\n",
    "} 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. `schedule.ccl`: allocates storage for gridfunctions, defines how the thorn's functions should be scheduled in a broader simulation, and specifies the regions of memory written to or read from gridfunctions. $\\text{schedule.ccl}$'s official documentation may be found [here](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-186000D2.4). \n",
    "\n",
    "We specify here the standardized ETK \"scheduling bins\" in which we want each of our thorn's functions to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.969633Z",
     "iopub.status.busy": "2021-03-07T17:08:52.968900Z",
     "iopub.status.idle": "2021-03-07T17:08:52.971962Z",
     "shell.execute_reply": "2021-03-07T17:08:52.972443Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting NRPyPlusTOVID//schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outrootdir/schedule.ccl\n",
    "STORAGE: ADMBase::metric[metric_timelevels], ADMBase::curv[metric_timelevels], ADMBase::lapse[lapse_timelevels], ADMBase::shift[shift_timelevels]\n",
    "\n",
    "schedule NRPyPlusTOVID_ET_InitialData IN HydroBase_Initial\n",
    "{\n",
    "  LANG: C\n",
    "  READS: grid::x(Everywhere)\n",
    "  READS: grid::y(Everywhere)\n",
    "  READS: grid::y(Everywhere)\n",
    "  WRITES: admbase::alp(Everywhere)\n",
    "  WRITES: admbase::betax(Everywhere)\n",
    "  WRITES: admbase::betay(Everywhere)\n",
    "  WRITES: admbase::betaz(Everywhere)\n",
    "  WRITES: admbase::kxx(Everywhere)\n",
    "  WRITES: admbase::kxy(Everywhere)\n",
    "  WRITES: admbase::kxz(Everywhere)\n",
    "  WRITES: admbase::kyy(Everywhere)\n",
    "  WRITES: admbase::kyz(Everywhere)\n",
    "  WRITES: admbase::kzz(Everywhere)\n",
    "  WRITES: admbase::gxx(Everywhere)\n",
    "  WRITES: admbase::gxy(Everywhere)\n",
    "  WRITES: admbase::gxz(Everywhere)\n",
    "  WRITES: admbase::gyy(Everywhere)\n",
    "  WRITES: admbase::gyz(Everywhere)\n",
    "  WRITES: admbase::gzz(Everywhere)\n",
    "  WRITES: hydrobase::vel[0](Everywhere)\n",
    "  WRITES: hydrobase::vel[1](Everywhere)\n",
    "  WRITES: hydrobase::vel[2](Everywhere)\n",
    "  WRITES: hydrobase::rho(Everywhere)\n",
    "  WRITES: hydrobase::eps(Everywhere)\n",
    "  WRITES: hydrobase::press(Everywhere)\n",
    "} \"Set up general relativistic hydrodynamic (GRHD) fields for TOV initial data\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_list'></a>\n",
    "\n",
    "## Step 3.c: Add the C code to the Einstein Toolkit compilation list \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein_list}$$\n",
    "\n",
    "We will also need `make.code.defn`, which indicates the list of files that need to be compiled. This thorn only has the one C file to compile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.977643Z",
     "iopub.status.busy": "2021-03-07T17:08:52.976810Z",
     "iopub.status.idle": "2021-03-07T17:08:52.979545Z",
     "shell.execute_reply": "2021-03-07T17:08:52.980036Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting NRPyPlusTOVID/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outdir/make.code.defn\n",
    "SRCS = InitialData.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 4: 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-ETK_thorn-NRPyPlusTOVID.pdf](Tutorial-ETK_thorn-NRPyPlusTOVID.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": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:52.987260Z",
     "iopub.status.busy": "2021-03-07T17:08:52.984043Z",
     "iopub.status.idle": "2021-03-07T17:08:57.285289Z",
     "shell.execute_reply": "2021-03-07T17:08:57.286023Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-NRPyPlusTOVID.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-ETK_thorn-NRPyPlusTOVID.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ETK_thorn-NRPyPlusTOVID\")"
   ]
  }
 ],
 "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
}