Tutorial-ADM_Initial_Data-UIUC_BlackHole.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",
    "# [UIUC Black Hole](https://arxiv.org/abs/1001.4077) Initial data\n",
    "\n",
    "## Authors: Terrence Pierre Jacques, Zach Etienne, & Ian Ruchlin\n",
    "\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module sets up UIUC Black Hole initial data ([Liu, Etienne, & Shapiro, PRD 80 121503, 2009](https://arxiv.org/abs/1001.4077)). \n",
    "\n",
    "### We can convert from spherical to any coordinate system defined in [reference_metric.py](../edit/reference_metric.py) (e.g., SinhSpherical, Cylindrical, Cartesian, etc.) using the [Exact ADM Spherical-or-Cartesian-to-BSSNCurvilinear converter module](Tutorial-ADM_Initial_Data-Converting_Exact_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.ipynb)\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to exhibit convergence to zero of the Hamiltonian and momentum constraint violation at the expected order to the exact solution (see plots at bottom of [the exact initial data validation start-to-finish tutorial notebook](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_Exact_Initial_Data.ipynb); momentum constraint violation in non-$\\phi$ directions is zero), and all quantities have been validated against the [original SENR code](https://bitbucket.org/zach_etienne/nrpy).\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/UIUCBlackHole.py](../edit/BSSN/UIUCBlackHole.py)\n",
    "\n",
    "## Introduction:\n",
    "UIUC black holes have the advantage of finite coordinate radius in the maximal spin limit. It is therefore excellent for studying very highly spinning black holes. This module sets the UIUC black hole at the origin. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents:  \n",
    "$$\\label{toc}$$\n",
    "\n",
    "1. [Step 1](#initializenrpy): Set up the needed NRPy+ infrastructure and declare core gridfunctions\n",
    "1. [Step 2](#bl_radius): The Boyer-Lindquist Radius\n",
    "    1. [Step 2.a](#define_inner_outer_radii): Define the inner and outer radii\n",
    "    1. [Step 2.b](#define_bl_radius): Define the Boyer-Lindquist radius\n",
    "1. [Step 3](#line_element): Define the line element, and extract components of $\\gamma_{ij}$\n",
    "1. [Step 4](#extrinsic_curvature): Define and construct nonzero components of the extrinsic curvature $K_{ij}$\n",
    "1. [Step 5](#lapse_shift): Construct Lapse function $\\alpha$ and components of shift vector $\\beta$\n",
    "1. [Step 6](#code_validation): Code Validation against `BSSN.UIUCBlackHole` NRPy+ module\n",
    "1. [Step 7](#latex_pdf_output) Output this notebook to $\\LaTeX$-formatted PDF file\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Set up the needed NRPy+ infrastructure and declare core gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "First, we will import the core modules of Python/NRPy+ and specify the main gridfunctions that we will need.\n",
    "Second, we set some basic NRPy+ parameters. E.g., set the spatial dimension parameter to 3.\n",
    "\n",
    "**Inputs for initial data**:\n",
    "\n",
    "* The black hole mass $M$.\n",
    "* The dimensionless spin parameter $\\chi = a/M$\n",
    "\n",
    "**Additional variables needed for spacetime evolution**:\n",
    "\n",
    "* Desired coordinate system Boyer-Lindquist coordinates $(r_{BL}, \\theta, \\phi)$\n",
    "<br>\n",
    "* Desired initial lapse $\\alpha$ and shift $\\beta^i$. We will choose our gauge conditions as $\\alpha=1$ and $\\beta^i=B^i=0$. $\\alpha = \\psi^{-2}$ will yield much better behavior, but the conformal factor $\\psi$ depends on the desired *destination* coordinate system (which may not be spherical coordinates)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step P0: Load needed modules\n",
    "import sympy as sp             # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par # NRPy+: Parameter interface\n",
    "import indexedexp as ixp       # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "\n",
    "# All gridfunctions will be written in terms of spherical coordinates (r, th, ph):\n",
    "r,th,ph = sp.symbols('r th ph', real=True)\n",
    "\n",
    "thismodule = \"UIUCBlackHole\"\n",
    "\n",
    "# Step 0: Set spatial dimension (must be 3 for BSSN)\n",
    "DIM = 3\n",
    "par.set_parval_from_str(\"grid::DIM\",DIM)\n",
    "\n",
    "# Step 1: Set psi, the conformal factor:\n",
    "\n",
    "# The UIUC initial data represent a Kerr black hole with mass M\n",
    "#  and dimensionless spin chi in UIUC quasi-isotropic coordinates,\n",
    "#   see https://arxiv.org/abs/1001.4077\n",
    "# Input parameters:\n",
    "M,chi = par.Cparameters(\"REAL\", thismodule, [\"M\",\"chi\"],[1.0,0.99])\n",
    "\n",
    "# Spin per unit mass\n",
    "a = M*chi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bl_radius'></a>\n",
    "\n",
    "# Step 2: The Boyer-Lindquist Radius \\[Back to [top](#toc)\\]\n",
    "$$\\label{bl_radius}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='define_inner_outer_radii'></a>\n",
    "\n",
    "## Step 2.a: Defining the Inner and Outer Radii \\[Back to [top](#toc)\\]\n",
    "$$\\label{define_inner_outer_radii}$$\n",
    "\n",
    "\n",
    "\n",
    "Boyer-Lindquist radii of the outer (+) and inner (−) horizons of the BH, defined under equation 1 in [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077) as \n",
    "$$ r_{\\pm} = M \\pm \\sqrt{M^2 - a^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defined under equation 1 in Liu, Etienne, & Shapiro (2009)\n",
    "# https://arxiv.org/pdf/1001.4077.pdf\n",
    "\n",
    "# Boyer - Lindquist outer horizon\n",
    "rp = M + sp.sqrt(M**2 - a**2)\n",
    "# Boyer - Lindquist inner horizon\n",
    "rm = M - sp.sqrt(M**2 - a**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='define_bl_radius'></a>\n",
    "\n",
    "## Step 2.b: Define the Boyer-Lindquist Radius  \\[Back to [top](#toc)\\]\n",
    "$$\\label{define_bl_radius}$$\n",
    "\n",
    "Define $r_{BL}$, equation 11 of [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077), using the radial coordinate $r$:\n",
    "\n",
    "$$  r_{BL} = r \\left( 1 + \\frac{r_+}{4r}\\right)^2.  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Boyer - Lindquist radius in terms of UIUC radius\n",
    "# Eq. 11\n",
    "# r_{BL} = r * ( 1 + r_+ / 4r )^2\n",
    "rBL = r*(1 + rp / (4*r))**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantities used to calculate the spatial metric $\\gamma_{ij}$, found under equation 2 of [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077):\n",
    "$$  \\Sigma = r_{BL}^2 + a^2 \\cos^2 \\theta, $$\n",
    "\n",
    "$$ \\Delta = r_{BL}^2 - 2Mr_{BL} + a^2,  $$ \n",
    "\n",
    "$$  A = \\left(r_{BL}^2 + a^2\\right)^2 - \\Delta a^2 \\sin^2 \\theta.  $$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Expressions found below Eq. 2\n",
    "# Sigma = r_{BL}^2 + a^2 cos^2 theta\n",
    "SIG = rBL**2 + a**2*sp.cos(th)**2\n",
    "\n",
    "# Delta = r_{BL}^2 - 2Mr_{BL} + a^2\n",
    "DEL = rBL**2 - 2*M*rBL + a**2\n",
    "\n",
    "# A = (r_{BL}^2 + a^2)^2 - Delta a^2 sin^2 theta\n",
    "AA = (rBL**2 + a**2)**2 - DEL*a**2*sp.sin(th)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='line_element'></a>\n",
    "\n",
    "# Step 3: Define the Line element and extract components of $\\gamma_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{line_element}$$\n",
    "\n",
    "The line element, defined in equation 13 of [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077):\n",
    "\n",
    "$$ ds^2 = \\frac{\\Sigma\\left(r + \\frac{r_+}{4}\\right)^2 } {r^3 \\left(r_{BL} - r_- \\right)} dr^2 + \\Sigma d\\theta^2  +  \\frac{ A \\sin^2 \\theta  }{\\Sigma} d\\phi^2  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# *** The ADM 3-metric in spherical basis ***\n",
    "gammaSphDD = ixp.zerorank2()\n",
    "# Declare the nonzero components of the 3-metric (Eq. 13):\n",
    "\n",
    "# ds^2 = Sigma (r + r_+/4)^2 / ( r^3 (r_{BL} - r_- ) * dr^2 +\n",
    "# Sigma d theta^2  +  (A sin^2 theta) / Sigma  *  d\\phi^2\n",
    "\n",
    "gammaSphDD[0][0] = ((SIG*(r + rp/4)**2)/(r**3*(rBL - rm)))\n",
    "gammaSphDD[1][1] = SIG\n",
    "gammaSphDD[2][2] = AA/SIG*sp.sin(th)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='extrinsic_curvature'></a>\n",
    "\n",
    "# Step 4: Define and construct nonzero components of extrinsic curvature $K_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{extrinsic_curvature}$$\n",
    "\n",
    "\n",
    "\n",
    "Nonzero components of the extrinsic curvature, equation 14 of [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077):\n",
    "\n",
    "$$ K_{r\\phi} = K_{\\phi r} = \\frac{Ma\\sin^2\\theta}{\\Sigma\\sqrt{A\\Sigma}} \\ \n",
    "    \\left[3r^4_{BL} + 2a^2 r^2_{BL} - a^4 - a^2 \\left(r^2_{BL} - a^2\\right) \\sin^2 \\theta\\right] \\\n",
    "    \\left(1 + \\frac{r_+}{4r}\\right) \\frac{1}{\\sqrt{r(r_{BL} - r_-)}}  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# *** The physical trace-free extrinsic curvature in spherical basis ***\n",
    "# Nonzero components of the extrinsic curvature K, given by\n",
    "# Eq. 14 of Liu, Etienne, & Shapiro, https://arxiv.org/pdf/1001.4077.pdf:\n",
    "KSphDD     = ixp.zerorank2() # K_{ij} = 0 for these initial data\n",
    "\n",
    "\n",
    "# K_{r phi} = K_{phi r} = (Ma sin^2 theta) / (Sigma sqrt{A Sigma}) *\n",
    "#     [3r^4_{BL} + 2a^2 r^2_{BL} - a^4 - a^2 (r^2_{BL} - a^2) sin^2 theta] *\n",
    "#     (1 + r_+ / 4r) (1 / sqrt{r(r_{BL} - r_-)})\n",
    "\n",
    "KSphDD[0][2] = KSphDD[2][0] = (M*a*sp.sin(th)**2)/(SIG*sp.sqrt(AA*SIG))*\\\n",
    "                (3*rBL**4 + 2*a**2*rBL**2 - a**4- a**2*(rBL**2 - a**2)*\\\n",
    "                 sp.sin(th)**2)*(1 + rp/(4*r))*1/sp.sqrt(r*(rBL - rm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonzero components of the extrinsic curvature, equation 15 of [Liu, Etienne, & Shapiro (2009)](https://arxiv.org/abs/1001.4077):\n",
    "\n",
    "$$ K_{\\theta\\phi} = K_{\\phi\\theta} = -\\frac{2a^3 Mr_{BL}\\cos\\theta \\sin^3\\theta} {\\Sigma \\sqrt{A\\Sigma} } \\left(r - \\frac{r_+}{4}\\right) \\sqrt {\\frac{r_{BL} - r_-}{r} }  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Components of the extrinsic curvature K, given by\n",
    "# Eq. 15 of Liu, Etienne, & Shapiro, https://arxiv.org/pdf/1001.4077.pdf:\n",
    "\n",
    "# K_{theta phi} = K_{phi theta} = -(2a^3 Mr_{BL} cos theta sin^3 theta) /\n",
    "#         (Sigma sqrt{A Sigma}) x (r - r_+ / 4) sqrt{(r_{BL} - r_-) / r }\n",
    "\n",
    "KSphDD[1][2] = KSphDD[2][1] = -((2*a**3*M*rBL*sp.cos(th)*sp.sin(th)**3)/ \\\n",
    "                (SIG*sp.sqrt(AA*SIG)))*(r - rp/4)*sp.sqrt((rBL - rm)/r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lapse_shift'></a>\n",
    "\n",
    "# Step 5: Construct Lapse function $\\alpha$ and components of shift vector $\\beta$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{lapse_shift}$$\n",
    "\n",
    "$$\\alpha=1$$ \n",
    "<br>\n",
    "$$\\beta^i=B^i=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphaSph = sp.sympify(1)\n",
    "betaSphU = ixp.zerorank1() # We generally choose \\beta^i = 0 for these initial data\n",
    "BSphU    = ixp.zerorank1() # We generally choose B^i = 0 for these initial data\n",
    "\n",
    "# Validated against original SENR: KSphDD[0][2], KSphDD[1][2], gammaSphDD[2][2], gammaSphDD[0][0], gammaSphDD[1][1]\n",
    "# print(sp.mathematica_code(gammaSphDD[1][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 6: Code Validation against `BSSN.UIUCBlackHole` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for UIUC black hole initial data between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.UIUCBlackHole](../edit/BSSN/UIUCBlackHole.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between Brill-Lindquist tutorial and NRPy+ BSSN.BrillLindquist module. ALL SHOULD BE ZERO.\n",
      "alphaSph - uibh.alphaSph = 0\n",
      "betaSphU[0] - uibh.betaSphU[0] = 0\n",
      "BSphU[0] - uibh.BaSphU[0] = 0\n",
      "gammaSphDD[0][0] - uibh.gammaSphDD[0][0] = 0\n",
      "KSphDD[0][0] - uibh.KSphDD[0][0] = 0\n",
      "gammaSphDD[0][1] - uibh.gammaSphDD[0][1] = 0\n",
      "KSphDD[0][1] - uibh.KSphDD[0][1] = 0\n",
      "gammaSphDD[0][2] - uibh.gammaSphDD[0][2] = 0\n",
      "KSphDD[0][2] - uibh.KSphDD[0][2] = 0\n",
      "betaSphU[1] - uibh.betaSphU[1] = 0\n",
      "BSphU[1] - uibh.BaSphU[1] = 0\n",
      "gammaSphDD[1][0] - uibh.gammaSphDD[1][0] = 0\n",
      "KSphDD[1][0] - uibh.KSphDD[1][0] = 0\n",
      "gammaSphDD[1][1] - uibh.gammaSphDD[1][1] = 0\n",
      "KSphDD[1][1] - uibh.KSphDD[1][1] = 0\n",
      "gammaSphDD[1][2] - uibh.gammaSphDD[1][2] = 0\n",
      "KSphDD[1][2] - uibh.KSphDD[1][2] = 0\n",
      "betaSphU[2] - uibh.betaSphU[2] = 0\n",
      "BSphU[2] - uibh.BaSphU[2] = 0\n",
      "gammaSphDD[2][0] - uibh.gammaSphDD[2][0] = 0\n",
      "KSphDD[2][0] - uibh.KSphDD[2][0] = 0\n",
      "gammaSphDD[2][1] - uibh.gammaSphDD[2][1] = 0\n",
      "KSphDD[2][1] - uibh.KSphDD[2][1] = 0\n",
      "gammaSphDD[2][2] - uibh.gammaSphDD[2][2] = 0\n",
      "KSphDD[2][2] - uibh.KSphDD[2][2] = 0\n"
     ]
    }
   ],
   "source": [
    "# First we import reference_metric, which is\n",
    "#   needed since BSSN.UIUCBlackHole calls\n",
    "#   BSSN.ADM_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear, which\n",
    "#   depends on reference_metric:\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "rfm.reference_metric()\n",
    "\n",
    "import BSSN.UIUCBlackHole as uibh\n",
    "uibh.UIUCBlackHole()\n",
    "\n",
    "print(\"Consistency check between Brill-Lindquist tutorial and NRPy+ BSSN.BrillLindquist module. ALL SHOULD BE ZERO.\")\n",
    "print(\"alphaSph - uibh.alphaSph = \"+str(sp.simplify(alphaSph - uibh.alphaSph)))\n",
    "for i in range(DIM):\n",
    "    print(\"betaSphU[\"+str(i)+\"] - uibh.betaSphU[\"+str(i)+\"] = \"+\\\n",
    "          str(sp.simplify(betaSphU[i] - uibh.betaSphU[i])))\n",
    "    print(\"BSphU[\"+str(i)+\"] - uibh.BaSphU[\"+str(i)+\"] = \"+str(sp.simplify(BSphU[i] - uibh.BSphU[i])))\n",
    "    for j in range(DIM):\n",
    "        print(\"gammaSphDD[\"+str(i)+\"][\"+str(j)+\"] - uibh.gammaSphDD[\"+str(i)+\"][\"+str(j)+\"] = \"+\\\n",
    "              str(sp.simplify(gammaSphDD[i][j] - uibh.gammaSphDD[i][j])))\n",
    "        print(\"KSphDD[\"+str(i)+\"][\"+str(j)+\"] - uibh.KSphDD[\"+str(i)+\"][\"+str(j)+\"] = \"+\\\n",
    "              str(sp.simplify(KSphDD[i][j] - uibh.KSphDD[i][j])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 7: 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-ADM_Initial_Data-UIUC_BlackHole.pdf](Tutorial-ADM_Initial_Data-UIUC_BlackHole.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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data-UIUC_BlackHole.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-ADM_Initial_Data-UIUC_BlackHole.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-UIUC_BlackHole\")"
   ]
  }
 ],
 "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
}