Tutorial-ETK_thorn-FishboneMoncriefID.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",
    "# `FishboneMoncriefID`: An Einstein Toolkit Initial Data Thorn for Fishbone-Moncrief initial data\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "[comment]: <> (Abstract: TODO)\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** Agrees with trusted Fishbone-Moncrief initial data module in HARM3D. Also generates results in agreement with trusted version sent to Event Horizon Telescope (EHT) GRMHD code comparison project collaborators. This thorn was used for the [IllinoisGRMHD](http://illinoisgrmhd.net) contribution to the [EHT GRMHD code comparison project](https://arxiv.org/abs/1904.04923).\n",
    "\n",
    "### NRPy+ Source Code for this module: [FishboneMoncriefID/FishboneMoncriefID.py](../edit/FishboneMoncriefID/FishboneMoncriefID.py) [\\[tutorial\\]](Tutorial-FishboneMoncriefID.ipynb) Constructs SymPy expressions for [Fishbone-Moncrief initial data](Tutorial-FishboneMoncriefID.ipynb)\n",
    "\n",
    "## Introduction:\n",
    "In this part of the tutorial, we will construct an Einstein Toolkit (ETK) thorn (module) that will set up Fishbone-Moncrief initial data. In the [Tutorial-FishboneMoncriefID](Tutorial-FishboneMoncriefID.ipynb) tutorial notebook, we used NRPy+ to construct the SymPy expressions for Fishbone-Moncrief 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 convert the SymPy expression for the Fishbone-Moncrief initial data into a C-code kernel\n",
    "1. [Step 2](#einstein): Interfacing with the Einstein Toolkit\n",
    "    1. [Step 2.a](#einstein_c): Constructing the Einstein Toolkit C-code calling functions that include the C code kernels\n",
    "    1. [Step 2.b](#einstein_ccl): CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure\n",
    "    1. [Step 2.c](#einstein_list): Add the C code to the Einstein Toolkit compilation list\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Call on NRPy+ to convert the SymPy expression for the Fishbone-Moncrief initial data into a C-code kernel \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "After importing the core modules, we will set `GridFuncMemAccess` to `ETK`. SymPy expressions for Fishbone-Moncrief initial data are written inside [FishboneMoncriefID/FishboneMoncriefID.py](../edit/FishboneMoncriefID/FishboneMoncriefID.py), and we simply import them for use here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Call on NRPy+ to convert the SymPy expression for the\n",
    "#         Fishbone-Moncrief initial data into a C-code kernel\n",
    "\n",
    "# Step 1a: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh,outputC # NRPy+: Core C code output module\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 loop as lp                # NRPy+: Generate C code loops\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions\n",
    "import FishboneMoncriefID.FishboneMoncriefID as fmid # Stores closed-form SymPy expressions for F-M initial data.\n",
    "\n",
    "# Step 1b: 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 1c: Within the ETK, the 3D gridfunctions x, y, and z store the\n",
    "#         Cartesian grid coordinates. Setting the gri.xx[] arrays\n",
    "#         to point to these gridfunctions forces NRPy+ to treat\n",
    "#         the Cartesian coordinate gridfunctions properly --\n",
    "#         reading them from memory as needed.\n",
    "xcoord,ycoord,zcoord = gri.register_gridfunctions(\"AUX\",[\"xcoord\",\"ycoord\",\"zcoord\"])\n",
    "gri.xx[0] = xcoord\n",
    "gri.xx[1] = ycoord\n",
    "gri.xx[2] = zcoord\n",
    "\n",
    "# Step 1d: Call the FishboneMoncriefID() function from within the\n",
    "#          FishboneMoncriefID/FishboneMoncriefID.py module. This\n",
    "#          sets all the ID gridfunctions.\n",
    "fmid.FishboneMoncriefID()\n",
    "Valencia3velocityU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"Valencia3velocityU\")\n",
    "\n",
    "# -={ Spacetime quantities: Generate C code from expressions and output to file }=-\n",
    "KerrSchild_to_print = [\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"alpha\"),rhs=fmid.IDalpha),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU0\"),rhs=fmid.IDbetaU[0]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU1\"),rhs=fmid.IDbetaU[1]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"betaU2\"),rhs=fmid.IDbetaU[2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD00\"),rhs=fmid.IDgammaDD[0][0]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD01\"),rhs=fmid.IDgammaDD[0][1]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD02\"),rhs=fmid.IDgammaDD[0][2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD11\"),rhs=fmid.IDgammaDD[1][1]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD12\"),rhs=fmid.IDgammaDD[1][2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"gammaDD22\"),rhs=fmid.IDgammaDD[2][2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD00\"),rhs=fmid.IDKDD[0][0]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD01\"),rhs=fmid.IDKDD[0][1]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD02\"),rhs=fmid.IDKDD[0][2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD11\"),rhs=fmid.IDKDD[1][1]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD12\"),rhs=fmid.IDKDD[1][2]),\\\n",
    "                     lhrh(lhs=gri.gfaccess(\"out_gfs\",\"KDD22\"),rhs=fmid.IDKDD[2][2]),\\\n",
    "                     ]\n",
    "# Force outCverbose=False for this module to avoid gigantic C files\n",
    "# filled with the non-CSE expressions.\n",
    "KerrSchild_CcodeKernel = fin.FD_outputC(\"returnstring\",KerrSchild_to_print,params=\"outCverbose=False\")\n",
    "\n",
    "# -={ GRMHD quantities: Generate C code from expressions and output to file }=-\n",
    "FMdisk_GRHD_rho_initial_to_print = [lhrh(lhs=gri.gfaccess(\"out_gfs\",\"rho_initial\"),rhs=fmid.rho_initial)]\n",
    "FMdisk_GRHD_rho_initial_CcodeKernel = fin.FD_outputC(\"returnstring\",FMdisk_GRHD_rho_initial_to_print)\n",
    "\n",
    "FMdisk_GRHD_velocities_to_print = [\\\n",
    "                                 lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU0\"),rhs=fmid.IDValencia3velocityU[0]),\\\n",
    "                                 lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU1\"),rhs=fmid.IDValencia3velocityU[1]),\\\n",
    "                                 lhrh(lhs=gri.gfaccess(\"out_gfs\",\"Valencia3velocityU2\"),rhs=fmid.IDValencia3velocityU[2]),\\\n",
    "                                 ]\n",
    "FMdisk_GRHD_velocities_CcodeKernel = fin.FD_outputC(\"returnstring\",FMdisk_GRHD_velocities_to_print)\n",
    "\n",
    "# Step 1f: Create directories for the thorn if they don't exist.\n",
    "Ccodesdir = \"FishboneMoncriefID\"\n",
    "cmd.mkdir(Ccodesdir)\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"src\"))\n",
    "\n",
    "# Step 1g: Write the C code kernel to file.\n",
    "with open(os.path.join(Ccodesdir,\"src\",\"KerrSchild.h\"), \"w\") as file:\n",
    "    file.write(str(KerrSchild_CcodeKernel.replace(\"time\",\"cctk_time\")))\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"src\",\"FMdisk_GRHD_velocities.h\"), \"w\") as file:\n",
    "    file.write(str(FMdisk_GRHD_velocities_CcodeKernel.replace(\"time\",\"cctk_time\")))\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"src\",\"FMdisk_GRHD_rho_initial.h\"), \"w\") as file:\n",
    "    file.write(str(FMdisk_GRHD_rho_initial_CcodeKernel.replace(\"time\",\"cctk_time\")))\n",
    "\n",
    "hm1string = outputC(fmid.hm1,\"hm1\",filename=\"returnstring\")\n",
    "with open(os.path.join(Ccodesdir,\"src\",\"FMdisk_GRHD_hm1.h\"), \"w\") as file:\n",
    "    file.write(str(hm1string))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein'></a>\n",
    "\n",
    "# Step 2: 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 2.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",
    "Here we construct `InitialData.c`, which contains C driver functions that pull in the necessary NRPy+ C-code kernels.\n",
    "\n",
    "First we set up driver routines to specify the Kerr-Schild metric and the Fishbone-Moncrief disk velocity at a given gridpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting FishboneMoncriefID/src/InitialData.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/InitialData.c\n",
    "#include <math.h>\n",
    "#include <stdio.h>\n",
    "#include <stdbool.h>\n",
    "#include <stdlib.h> // Needed for rand()\n",
    "\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "#include \"cctk_Arguments.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 FishboneMoncrief_KerrSchild(const cGH* restrict const cctkGH,const CCTK_INT *cctk_lsh,\n",
    "                                 const CCTK_INT i0,const CCTK_INT i1,const CCTK_INT i2,\n",
    "                                 const CCTK_REAL *xcoordGF,const CCTK_REAL *ycoordGF,const CCTK_REAL *zcoordGF,\n",
    "                                 CCTK_REAL *alphaGF,CCTK_REAL *betaU0GF,CCTK_REAL *betaU1GF,CCTK_REAL *betaU2GF,\n",
    "                                 CCTK_REAL *gammaDD00GF,CCTK_REAL *gammaDD01GF,CCTK_REAL *gammaDD02GF,CCTK_REAL *gammaDD11GF,CCTK_REAL *gammaDD12GF,CCTK_REAL *gammaDD22GF,\n",
    "                                 CCTK_REAL     *KDD00GF,CCTK_REAL     *KDD01GF,CCTK_REAL     *KDD02GF,CCTK_REAL     *KDD11GF,CCTK_REAL     *KDD12GF,CCTK_REAL     *KDD22GF)\n",
    "{\n",
    "\n",
    "  DECLARE_CCTK_PARAMETERS\n",
    "\n",
    "#include \"KerrSchild.h\"\n",
    "\n",
    "}\n",
    "\n",
    "void FishboneMoncrief_FMdisk_GRHD_velocities(const cGH* restrict const cctkGH,const CCTK_INT *cctk_lsh,\n",
    "                                             const CCTK_INT i0,const CCTK_INT i1,const CCTK_INT i2,\n",
    "                                             const CCTK_REAL *xcoordGF,const CCTK_REAL *ycoordGF,const CCTK_REAL *zcoordGF,\n",
    "                                             CCTK_REAL *Valencia3velocityU0GF, CCTK_REAL *Valencia3velocityU1GF, CCTK_REAL *Valencia3velocityU2GF)\n",
    "{\n",
    "\n",
    "  DECLARE_CCTK_PARAMETERS\n",
    "\n",
    "#include \"FMdisk_GRHD_velocities.h\"\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we set up the driver function for setting all metric and hydrodynamical fields $\\rho,P,\\epsilon,v^i$.\n",
    "\n",
    "**Important**: Suppose the Fishbone-Moncrief initial data yield a density $\\rho(r,\\theta)$ (which is valid for all Fishbone-Moncrief disks centered at the origin, $r=0$, as F-M disks are axisymmetric). Then the disk will have pressure\n",
    "$$\n",
    "P = \\kappa \\rho^\\Gamma.\n",
    "$$\n",
    "\n",
    "Since the disk is not self-gravitating, we are allowed to rescale the maximum density in the disk to be one in code units; i.e., $\\rho_{\\rm max}=1$. This may be incompatible with the initial choice of polytropic constant $\\kappa$, as rescaling the density results in a rescaling of pressure $P$, as follows.\n",
    "\n",
    "When we rescale $\\rho$ so that the maximum density in the disk is one, we make the following transformation:\n",
    "$$\n",
    "\\rho \\to \\rho' = \\frac{\\rho}{\\rho_{\\rm max}}.\n",
    "$$\n",
    "Since pressure has units of $\\rho c^2$, and we use $G=c=1$ units, pressure must therefore be rescaled by the same factor:\n",
    "\\begin{align}\n",
    "P \\to P' &= \\frac{P}{\\rho_{\\rm max}} \\\\\n",
    "&= \\frac{\\kappa \\rho^\\Gamma}{\\rho_{\\rm max}} \\\\\n",
    "&= \\kappa \\frac{\\rho^\\Gamma}{\\rho_{\\rm max}} \\\\\n",
    "&= \\kappa \\frac{(\\rho' \\rho_{\\rm max})^\\Gamma}{\\rho_{\\rm max}} \\\\\n",
    "&= \\kappa \\rho_{\\rm max}^{\\Gamma-1} (\\rho')^\\Gamma \\\\\n",
    "&= \\kappa' (\\rho')^\\Gamma\n",
    "\\end{align}\n",
    "\n",
    "Thus the polytropic equation of state is still valid, but only if \n",
    "$$\n",
    "\\kappa' = \\kappa \\rho_{\\rm max}^{\\Gamma-1} = \\frac{P_{\\rm max}}{\\rho_{\\rm max}}.\n",
    "$$\n",
    "As e.g., `IllinoisGRMHD` requires that the initial $P'$ be given as a polytropic equation of state, with $P'_{\\rm cold} = \\kappa' (\\rho')^\\Gamma$, $\\kappa'$ must be input into the `FishboneMoncriefID` (and `IllinoisGRMHD`) thorns instead of $\\kappa$. If this does not happen, the code will error out, providing the correct value for $\\kappa'$ that must be set in the parameter file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to FishboneMoncriefID/src/InitialData.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $Ccodesdir/src/InitialData.c\n",
    "\n",
    "void FishboneMoncrief_ET_GRHD_initial(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  CCTK_VINFO(\"Fishbone-Moncrief Disk Initial data.\");\n",
    "  CCTK_VINFO(\"Using input parameters of\\n a = %e,\\n M = %e,\\nr_in = %e,\\nr_at_max_density = %e\\nkappa = %e\\ngamma = %e\",a,M,r_in,r_at_max_density,kappa,gamma);\n",
    "\n",
    "  // First compute maximum pressure and density\n",
    "  CCTK_REAL P_max, rho_max;\n",
    "  {\n",
    "    CCTK_REAL hm1;\n",
    "    CCTK_REAL xcoord = r_at_max_density;\n",
    "    CCTK_REAL ycoord = 0.0;\n",
    "    CCTK_REAL zcoord = 0.0;\n",
    "    {\n",
    "#include \"FMdisk_GRHD_hm1.h\"\n",
    "    }\n",
    "    rho_max = pow( hm1 * (gamma-1.0) / (kappa*gamma), 1.0/(gamma-1.0) );\n",
    "    P_max   = kappa * pow(rho_max, gamma);\n",
    "  }\n",
    "\n",
    "  // We enforce units such that rho_max = 1.0; if these units are not obeyed, then\n",
    "  //    we error out. If we did not error out, then the value of kappa used in all\n",
    "  //    EOS routines would need to be changed, and generally these appear as\n",
    "  //    read-only parameters.\n",
    "  if(fabs(P_max/rho_max - kappa) > 1e-8) {\n",
    "    printf(\"Error: To ensure that P = kappa*rho^Gamma, where rho_max = 1.0,\\n\");\n",
    "    printf(\"       you must set (in your parfile) the polytropic constant kappa = P_max/rho_max = %.15e\\n\\n\",P_max/rho_max);\n",
    "    printf(\" Needed values for kappa, for common values of Gamma:\\n\");\n",
    "    printf(\" For Gamma =4/3, use kappa=K_initial=K_poly = 4.249572342020724e-03 to ensure rho_max = 1.0\\n\");\n",
    "    printf(\" For Gamma =5/3, use kappa=K_initial=K_poly = 6.799315747233158e-03 to ensure rho_max = 1.0\\n\");\n",
    "    printf(\" For Gamma = 2,  use kappa=K_initial=K_poly = 8.499144684041449e-03 to ensure rho_max = 1.0\\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "#pragma omp parallel for\n",
    "  for(CCTK_INT k=0;k<cctk_lsh[2];k++) for(CCTK_INT j=0;j<cctk_lsh[1];j++) for(CCTK_INT i=0;i<cctk_lsh[0];i++) {\n",
    "        CCTK_INT idx = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        CCTK_REAL xcoord = x[idx];\n",
    "        CCTK_REAL ycoord = y[idx];\n",
    "        CCTK_REAL zcoord = z[idx];\n",
    "        CCTK_REAL rr = r[idx];\n",
    "\n",
    "        FishboneMoncrief_KerrSchild(cctkGH,cctk_lsh,\n",
    "                                    i,j,k,\n",
    "                                    x,y,z,\n",
    "                                    alp,betax,betay,betaz,\n",
    "                                    gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "                                    kxx,kxy,kxz,kyy,kyz,kzz);\n",
    "\n",
    "        CCTK_REAL hm1;\n",
    "        bool set_to_atmosphere=false;\n",
    "        if(rr > r_in) {\n",
    "          {\n",
    "#include \"FMdisk_GRHD_hm1.h\"\n",
    "          }\n",
    "          if(hm1 > 0) {\n",
    "            rho[idx] = pow( hm1 * (gamma-1.0) / (kappa*gamma), 1.0/(gamma-1.0) ) / rho_max;\n",
    "            press[idx] = kappa*pow(rho[idx], gamma);\n",
    "            // P = (\\Gamma - 1) rho epsilon\n",
    "            eps[idx] = press[idx] / (rho[idx] * (gamma - 1.0));\n",
    "            FishboneMoncrief_FMdisk_GRHD_velocities(cctkGH,cctk_lsh,\n",
    "                                                    i,j,k,\n",
    "                                                    x,y,z,\n",
    "                                                    velx,vely,velz);\n",
    "          } else {\n",
    "            set_to_atmosphere=true;\n",
    "          }\n",
    "        } else {\n",
    "          set_to_atmosphere=true;\n",
    "        }\n",
    "        // Outside the disk? Set to atmosphere all hydrodynamic variables!\n",
    "        if(set_to_atmosphere) {\n",
    "          // Choose an atmosphere such that\n",
    "          //   rho =       1e-5 * r^(-3/2), and\n",
    "          //   P   = k rho^gamma\n",
    "          // Add 1e-100 or 1e-300 to rr or rho to avoid divisions by zero.\n",
    "          rho[idx] = 1e-5 * pow(rr + 1e-100,-3.0/2.0);\n",
    "          press[idx] = kappa*pow(rho[idx], gamma);\n",
    "          eps[idx] = press[idx] / ((rho[idx] + 1e-300) * (gamma - 1.0));\n",
    "          w_lorentz[idx] = 1.0;\n",
    "          velx[idx] = 0.0;\n",
    "          vely[idx] = 0.0;\n",
    "          velz[idx] = 0.0;\n",
    "        }\n",
    "      }\n",
    "\n",
    "  CCTK_INT final_idx = CCTK_GFINDEX3D(cctkGH,cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1);\n",
    "  CCTK_VINFO(\"=====   OUTPUTS   =====\");\n",
    "  CCTK_VINFO(\"betai: %e %e %e \\ngij: %e %e %e %e %e %e \\nKij: %e %e %e %e %e %e\\nalp: %e\\n\",betax[final_idx],betay[final_idx],betaz[final_idx],gxx[final_idx],gxy[final_idx],gxz[final_idx],gyy[final_idx],gyz[final_idx],gzz[final_idx],kxx[final_idx],kxy[final_idx],kxz[final_idx],kyy[final_idx],kyz[final_idx],kzz[final_idx],alp[final_idx]);\n",
    "  CCTK_VINFO(\"rho: %.15e\\nPressure: %.15e\\nvx: %.15e\\nvy: %.15e\\nvz: %.15e\",rho[final_idx],press[final_idx],velx[final_idx],vely[final_idx],velz[final_idx]);\n",
    "}\n",
    "\n",
    "void FishboneMoncrief_ET_GRHD_initial__perturb_pressure(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  for(CCTK_INT k=0;k<cctk_lsh[2];k++) for(CCTK_INT j=0;j<cctk_lsh[1];j++) for(CCTK_INT i=0;i<cctk_lsh[0];i++) {\n",
    "        CCTK_INT idx = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        // Generate random number in range [0,1),\n",
    "        // snippet courtesy http://daviddeley.com/random/crandom.htm\n",
    "        CCTK_REAL random_number_between_0_and_1 = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );\n",
    "\n",
    "        CCTK_REAL random_number_between_min_and_max = random_min + (random_max - random_min)*random_number_between_0_and_1;\n",
    "        press[idx] = press[idx]*(1.0 + random_number_between_min_and_max);\n",
    "        // Add 1e-300 to rho to avoid division by zero when density is zero.\n",
    "        eps[idx] = press[idx] / ((rho[idx] + 1e-300) * (gamma - 1.0));\n",
    "      }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_ccl'></a>\n",
    "\n",
    "## Step 2.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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting FishboneMoncriefID/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/interface.ccl\n",
    "implements: FishboneMoncriefID\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": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting FishboneMoncriefID/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/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",
    "  \"FishboneMoncriefID\" :: \"Initial data from FishboneMoncriefID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_lapse\n",
    "{\n",
    "  \"FishboneMoncriefID\" :: \"Initial lapse from FishboneMoncriefID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_shift\n",
    "{\n",
    "  \"FishboneMoncriefID\" :: \"Initial shift from FishboneMoncriefID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_dtlapse\n",
    "{\n",
    "  \"FishboneMoncriefID\" :: \"Initial dtlapse from FishboneMoncriefID solution\"\n",
    "}\n",
    "EXTENDS KEYWORD initial_dtshift\n",
    "{\n",
    "  \"FishboneMoncriefID\" :: \"Initial dtshift from FishboneMoncriefID solution\"\n",
    "}\n",
    "\n",
    "shares: HydroBase\n",
    "EXTENDS KEYWORD initial_hydro\n",
    "{\n",
    "  \"FishboneMoncriefID\" :: \"Initial GRHD data from FishboneMoncriefID solution\"\n",
    "}\n",
    "\n",
    "#[\"r_in\",\"r_at_max_density\",\"a\",\"M\"] A_b, kappa, gamma\n",
    "restricted:\n",
    "CCTK_REAL r_in \"Fixes the inner edge of the disk\"\n",
    "{\n",
    " 0.0:* :: \"Must be positive\"\n",
    "} 6.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL r_at_max_density \"Radius at maximum disk density. Needs to be > r_in\"\n",
    "{\n",
    " 0.0:* :: \"Must be positive\"\n",
    "} 12.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL a \"The spin parameter of the black hole\"\n",
    "{\n",
    " 0:1.0 :: \"Positive values, up to 1. Negative disallowed, as certain roots are chosen in the hydro fields setup. Check those before enabling negative spins!\"\n",
    "} 0.9375\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL M \"Kerr-Schild BH mass. Probably should always set M=1.\"\n",
    "{\n",
    " 0.0:* :: \"Must be positive\"\n",
    "} 1.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL A_b \"Scaling factor for the vector potential\"\n",
    "{\n",
    " *:* :: \"\"\n",
    "} 1.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL kappa \"Equation of state: P = kappa * rho^gamma\"\n",
    "{\n",
    " 0.0:* :: \"Positive values\"\n",
    "} 1.0e-3\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL gamma \"Equation of state: P = kappa * rho^gamma\"\n",
    "{\n",
    " 0.0:* :: \"Positive values\"\n",
    "} 1.3333333333333333333333333333\n",
    "\n",
    "##################################\n",
    "# PRESSURE PERTURBATION PARAMETERS\n",
    "private:\n",
    "CCTK_REAL random_min \"Floor value of random perturbation to initial pressure, where perturbed pressure = pressure*(1.0 + (random_min + (random_max-random_min)*RAND[0,1)))\"\n",
    "{\n",
    "  *:*      :: \"Any value\"\n",
    "} -0.02\n",
    "\n",
    "private:\n",
    "CCTK_REAL random_max \"Ceiling value of random perturbation to initial pressure, where perturbed pressure = pressure*(1.0 + (random_min + (random_max-random_min)*RAND[0,1)))\"\n",
    "{\n",
    "  *:*      :: \"Any value\"\n",
    "} 0.02\n",
    "\n"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting FishboneMoncriefID/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/schedule.ccl\n",
    "STORAGE: ADMBase::metric[metric_timelevels], ADMBase::curv[metric_timelevels], ADMBase::lapse[lapse_timelevels], ADMBase::shift[shift_timelevels]\n",
    "\n",
    "schedule FishboneMoncrief_ET_GRHD_initial IN HydroBase_Initial\n",
    "{\n",
    "  LANG: C\n",
    "  READS: grid::x(Everywhere)\n",
    "  READS: grid::y(Everywhere)\n",
    "  READS: grid::z(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(Everywhere) # Note that vel is a vector gridfunction.\n",
    "  WRITES: hydrobase::rho(Everywhere)\n",
    "  WRITES: hydrobase::eps(Everywhere)\n",
    "  WRITES: hydrobase::press(Everywhere)\n",
    "} \"Set up general relativistic hydrodynamic (GRHD) fields for Fishbone-Moncrief disk\"\n",
    "\n",
    "schedule FishboneMoncrief_ET_GRHD_initial__perturb_pressure IN CCTK_INITIAL AFTER Seed_Magnetic_Fields BEFORE IllinoisGRMHD_ID_Converter\n",
    "{\n",
    "    LANG: C\n",
    "} \"Add random perturbation to initial pressure, after seed magnetic fields have been set up (in case we'd like the seed magnetic fields to depend on the pristine pressures)\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_list'></a>\n",
    "\n",
    "## Step 2.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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting FishboneMoncriefID/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/make.code.defn\n",
    "SRCS = InitialData.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: 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-FishboneMoncriefID.pdf](Tutorial-ETK_thorn-FishboneMoncriefID.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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-FishboneMoncriefID.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-ETK_thorn-FishboneMoncriefID.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-FishboneMoncriefID\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}