Tutorial-ADM_Initial_Data-StaticTrumpet.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",
    "# [Static Trumpet Black Hole](https://arxiv.org/abs/1403.5484) Initial data\n",
    "\n",
    "## Authors: Terrence Pierre Jacques, Zach Etienne & Ian Ruchlin\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module sets up Static Trumpet Black Hole initial data ([Dennison and Baumgarte, 2014 Class. Quantum Grav. 31 117001](https://arxiv.org/abs/1403.5484)). \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 constraint violation at the expected order to the exact solution (see plot 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 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/StaticTrumpet.py](../edit/BSSN/StaticTrumpet.py)\n",
    "\n",
    "## Introduction:\n",
    "These initial data are derived from a family of analytical coordinate systems representing the Schwarzschild spacetime. The coordinates extend smoothly through the black hole event horizon, the spatial coordinates are isotropic (so that the spatial metric can be written as a conformal factor to some power times a flat spatial metric), and, for almost all members of the family, the spatial slices take a so-called $\\textit{trumpet geometry}$. Moreover, all expressions are surprisingly simple. This module sets the static trumpet 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](#initialize_nrpy): Set up the needed NRPy+ infrastructure and declare core gridfunctions\n",
    "1. [Step 2](#conformal_factor_psi): The conformal factor $\\psi$\n",
    "    1. [Step 2.a](#define_psi): Define the conformal factor $\\psi$\n",
    "    1. [Step 2.b](#nonzero_gamma): Define and construct nonzero components of $\\gamma_{ij}$\n",
    "1. [Step 3](#extrinsic_curvature): Define and construct nonzero components of the extrinsic curvature $K_{ij}$, at the radius $R_0 = M$\n",
    "1. [Step 4](#lapse_shift): Construct Lapse function $\\alpha$ and components of shift vector $\\beta$\n",
    "1. [Step 5](#code_validation): Code Validation against `BSSN.StaticTrumpet` NRPy+ module\n",
    "1. [Step 6](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initialize_nrpy'></a>\n",
    "\n",
    "# Step 1: Set up the needed NRPy+ infrastructure and declare core gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{initialize_nrpy}$$\n",
    "\n",
    "First, we will import the core modules of Python/NRPy+ and specify the main gridfunctions that we will need. \n",
    "\n",
    "**Input for initial data**:\n",
    "\n",
    "* The black hole mass $M$.\n"
   ]
  },
  {
   "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 = \"StaticTrumpet\"\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",
    "# Input parameters:\n",
    "M = par.Cparameters(\"REAL\", thismodule, [\"M\"], 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='conformal_factor_psi'></a>\n",
    "\n",
    "# Step 2: The conformal factor $\\psi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{conformal_factor_psi}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='define_psi'></a>\n",
    "\n",
    "## Step 2.a: Define the conformal factor $\\psi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{define_psi}$$\n",
    "\n",
    "\n",
    "The conformal factor, defined in equation 13 of [Dennison and Baumgarte (2014)](https://arxiv.org/abs/1403.5484), setting $R_0 = M$,\n",
    "$$ \\psi = \\sqrt{1 + \\frac{M}{r}}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Set psi, the StaticTrumpet conformal factor\n",
    "# Dennison and Baumgarte (2014) Eq. 13\n",
    "# https://arxiv.org/pdf/1403.5484.pdf\n",
    "\n",
    "# psi = sqrt{1 + M/r }\n",
    "psi0 = sp.sqrt(1 + M/r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nonzero_gamma'></a>\n",
    "\n",
    "## Step 2.b: Define and construct nonzero components of $\\gamma_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{nonzero_gamma}$$\n",
    "\n",
    "The spatial metric, defined in equation 15 of [Dennison and Baumgarte (2014)](https://arxiv.org/abs/1403.5484),\n",
    "$$ \\gamma_{ij} = \\psi^4 \\eta_{ij}, $$\n",
    "\n",
    "where $\\eta_{ij}$ is the flat metric in spherical polar coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# *** The physical spatial metric in spherical basis ***\n",
    "# Set the upper-triangle of the matrix...\n",
    "# Eq. 15\n",
    "# eta_00 = 1, eta_11 = r^2, eta_22 = r^2 * sin^2 (theta)\n",
    "gammaSphDD = ixp.zerorank2()\n",
    "gammaSphDD[0][0] = psi0**4\n",
    "gammaSphDD[1][1] = psi0**4 * r**2\n",
    "gammaSphDD[2][2] = psi0**4 * r**2*sp.sin(th)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='extrinsic_curvature'></a>\n",
    "\n",
    "# Step 3: Define and construct nonzero components of the extrinsic curvature $K_{ij}$, at the radius $R_0 = M$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{extrinsic_curvature}$$\n",
    "\n",
    "Components of the extrinsic curvature in spherical basis, defined in equations 19 and 20 of [Dennison and Baumgarte (2014)](https://arxiv.org/abs/1403.5484),\n",
    "\n",
    "$$ K_{rr} = - \\frac{r \\left( M-R_0 \\right) + MR_0}{r^2 f_1}, $$\n",
    "\n",
    "<br>\n",
    "\n",
    "$$ K_{\\theta\\theta} = \\frac{K_{\\phi\\phi}}{\\sin^2 \\theta} = f_1, $$\n",
    "\n",
    "<br>\n",
    "\n",
    "where $f_1 = \\sqrt{2r \\left( M-R_0 \\right) + R_0 \\left( 2M-R_0 \\right).}$ Setting $R_0 = M$, these equations reduce to\n",
    "\n",
    "<br>\n",
    "\n",
    "$$ K_{rr} = -\\frac{M}{r^2}, $$\n",
    "\n",
    "<br>\n",
    "\n",
    "$$ K_{\\theta\\theta} = \\frac{K_{\\phi\\phi}}{\\sin^2 \\theta} = M. $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# *** The physical trace-free extrinsic curvature in spherical basis ***\n",
    "# Set the upper-triangle of the matrix...\n",
    "\n",
    "# Eq.19 and 20\n",
    "KSphDD = ixp.zerorank2()\n",
    "\n",
    "# K_{rr} = - M / r^2\n",
    "KSphDD[0][0] = -M / r**2\n",
    "\n",
    "# K_{theta theta} = K_{phi phi} / sin^2 theta = M\n",
    "KSphDD[1][1] = M\n",
    "\n",
    "KSphDD[2][2] = M * sp.sin(th)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lapse_shift'></a>\n",
    "\n",
    "# Step 4: Construct Lapse function $\\alpha$ and components of shift vector $\\beta$\\[Back to [top](#toc)\\]\n",
    "$$\\label{lapse_shift}$$\n",
    "\n",
    "Laspe function and shift vector components, equation 15 of [Dennison and Baumgarte (2014)](https://arxiv.org/abs/1403.5484), setting $R_0 = M$,\n",
    "$$ \\alpha = \\frac{r}{r+M}, $$\n",
    "<br>\n",
    "$$ \\beta^r = \\frac{Mr}{\\left (r+M \\right)^2}, $$\n",
    "<br>\n",
    "$$\\beta^\\theta = \\beta^\\phi = 0. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Eq. 15\n",
    "# alpha = r / (r+M)\n",
    "alphaSph = r / (r + M)\n",
    "\n",
    "betaSphU = ixp.zerorank1()\n",
    "# beta^r = Mr / (r + M)^2\n",
    "betaSphU[0] = M*r / (r + M)**2\n",
    "\n",
    "BSphU = ixp.zerorank1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 5: Code Validation against `BSSN.StaticTrumpet` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for static trumpet black hole initial data between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.StaticTrumpet](../edit/BSSN/StaticTrumpet.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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 - st.alphaSph = 0\n",
      "betaSphU[0] - st.betaSphU[0] = 0\n",
      "BSphU[0] - st.BaSphU[0] = 0\n",
      "gammaSphDD[0][0] - st.gammaSphDD[0][0] = 0\n",
      "KSphDD[0][0] - st.KSphDD[0][0] = 0\n",
      "gammaSphDD[0][1] - st.gammaSphDD[0][1] = 0\n",
      "KSphDD[0][1] - st.KSphDD[0][1] = 0\n",
      "gammaSphDD[0][2] - st.gammaSphDD[0][2] = 0\n",
      "KSphDD[0][2] - st.KSphDD[0][2] = 0\n",
      "betaSphU[1] - st.betaSphU[1] = 0\n",
      "BSphU[1] - st.BaSphU[1] = 0\n",
      "gammaSphDD[1][0] - st.gammaSphDD[1][0] = 0\n",
      "KSphDD[1][0] - st.KSphDD[1][0] = 0\n",
      "gammaSphDD[1][1] - st.gammaSphDD[1][1] = 0\n",
      "KSphDD[1][1] - st.KSphDD[1][1] = 0\n",
      "gammaSphDD[1][2] - st.gammaSphDD[1][2] = 0\n",
      "KSphDD[1][2] - st.KSphDD[1][2] = 0\n",
      "betaSphU[2] - st.betaSphU[2] = 0\n",
      "BSphU[2] - st.BaSphU[2] = 0\n",
      "gammaSphDD[2][0] - st.gammaSphDD[2][0] = 0\n",
      "KSphDD[2][0] - st.KSphDD[2][0] = 0\n",
      "gammaSphDD[2][1] - st.gammaSphDD[2][1] = 0\n",
      "KSphDD[2][1] - st.KSphDD[2][1] = 0\n",
      "gammaSphDD[2][2] - st.gammaSphDD[2][2] = 0\n",
      "KSphDD[2][2] - st.KSphDD[2][2] = 0\n"
     ]
    }
   ],
   "source": [
    "# First we import reference_metric, which is\n",
    "#   needed since BSSN.BrillLindquist 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.StaticTrumpet as st\n",
    "st.StaticTrumpet()\n",
    "\n",
    "print(\"Consistency check between Brill-Lindquist tutorial and NRPy+ BSSN.BrillLindquist module. ALL SHOULD BE ZERO.\")\n",
    "print(\"alphaSph - st.alphaSph = \"+str(sp.simplify(alphaSph - st.alphaSph)))\n",
    "for i in range(DIM):\n",
    "    print(\"betaSphU[\"+str(i)+\"] - st.betaSphU[\"+str(i)+\"] = \"+\\\n",
    "          str(sp.simplify(betaSphU[i] - st.betaSphU[i])))\n",
    "    print(\"BSphU[\"+str(i)+\"] - st.BaSphU[\"+str(i)+\"] = \"+str(sp.simplify(BSphU[i] - st.BSphU[i])))\n",
    "    for j in range(DIM):\n",
    "        print(\"gammaSphDD[\"+str(i)+\"][\"+str(j)+\"] - st.gammaSphDD[\"+str(i)+\"][\"+str(j)+\"] = \"+\\\n",
    "              str(sp.simplify(gammaSphDD[i][j] - st.gammaSphDD[i][j])))\n",
    "        print(\"KSphDD[\"+str(i)+\"][\"+str(j)+\"] - st.KSphDD[\"+str(i)+\"][\"+str(j)+\"] = \"+\\\n",
    "              str(sp.simplify(KSphDD[i][j] - st.KSphDD[i][j])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: Output this notebook to $\\LaTeX$-formatted PDF \\[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-StaticTrumpet.pdf](Tutorial-ADM_Initial_Data-StaticTrumpet.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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data-StaticTrumpet.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-ADM_Initial_Data-StaticTrumpet.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-StaticTrumpet\")"
   ]
  }
 ],
 "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
}