Tutorial-BSSN_constraints.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",
    "# [BSSN](http://www2.yukawa.kyoto-u.ac.jp/~yuichiro.sekiguchi/3+1.pdf) Hamiltonian and momentum constraint equations, in ***curvilinear*** coordinates, using a covariant reference metric approach: C code generation\n",
    "\n",
    "## Authors: Ian Ruchlin & Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This notebook demonstrates the construction of the BSSN Hamiltonian and momentum constraint equations as symbolic (SymPy) expressions, in terms of the core BSSN quantities $\\left\\{h_{i j},a_{i j},\\phi, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$, as defined in [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658). The module implements a generic curvilinear coordinate reference metric approach which is an extension of the spherical coordinate reference metric approach of [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632), building upon the covariant Lagrangian BSSN formalism of [Brown (2009)](https://arxiv.org/abs/0902.3652). \n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** All expressions generated in this module have been validated against a trusted code where applicable (the original NRPy+/SENR code, which itself was validated against [Baumgarte's code](https://arxiv.org/abs/1211.6632)).\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py)\n",
    "\n",
    "\n",
    "[comment]: <> (Introduction: TODO)"
   ]
  },
  {
   "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): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#hamiltonianconstraint): Construct the Hamiltonian constraint $\\mathcal{H}$.\n",
    "1. [Step 3](#momentumconstraint): Construct the momentum constraint $\\mathcal{M}^i$.\n",
    "1. [Step 4](#code_validation): Code Validation against `BSSN.BSSN_constraints` NRPy+ module\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize needed Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "We start by loading the needed modules. Notably, this module depends on several quantities defined in the [BSSN/BSSN_quantities.py](../edit/BSSN/BSSN_quantities.py) Python code, documented in the NRPy+ [BSSN quantities](Tutorial-BSSN_quantities.ipynb). In [Step 2](#hamiltonianconstraint) we call functions within [BSSN.BSSN_quantities](../edit/BSSN/BSSN_quantities.py) to define quantities needed in this module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:01.823104Z",
     "iopub.status.busy": "2021-03-07T17:21:01.822043Z",
     "iopub.status.idle": "2021-03-07T17:21:02.561250Z",
     "shell.execute_reply": "2021-03-07T17:21:02.560504Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Initialize needed Python/NRPy+ modules\n",
    "import sympy as sp               # SymPy, Python's core symbolic algebra package on which NRPy+ depends\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import BSSN.BSSN_quantities as Bq\n",
    "\n",
    "# Step 1.a: Set spatial dimension (must be 3 for BSSN, as BSSN is\n",
    "#           a 3+1-dimensional decomposition of the general\n",
    "#           relativistic field equations)\n",
    "DIM = 3\n",
    "\n",
    "# Step 1.b: Given the chosen coordinate system, set up\n",
    "#           corresponding reference metric and needed\n",
    "#           reference metric quantities\n",
    "# The following function call sets up the reference metric\n",
    "#    and related quantities, including rescaling matrices ReDD,\n",
    "#    ReU, and hatted quantities.\n",
    "rfm.reference_metric()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='hamiltonianconstraint'></a>\n",
    "\n",
    "# Step 2: $\\mathcal{H}$, the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{hamiltonianconstraint}$$\n",
    "\n",
    "Next we define the Hamiltonian constraint. Eq. 13 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf) yields:\n",
    "$$\n",
    "\\mathcal{H} = {\\underbrace {\\textstyle \\frac{2}{3} K^2}_{\\rm Term\\ 1}} - \n",
    "{\\underbrace {\\textstyle \\bar{A}_{ij} \\bar{A}^{ij}}_{\\rm Term\\ 2}} + \n",
    "{\\underbrace {\\textstyle e^{-4\\phi} \\left(\\bar{R} - 8 \\bar{D}^i \\phi \\bar{D}_i \\phi - 8 \\bar{D}^2 \\phi\\right)}_{\\rm Term\\ 3}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:02.635838Z",
     "iopub.status.busy": "2021-03-07T17:21:02.600086Z",
     "iopub.status.idle": "2021-03-07T17:21:03.413371Z",
     "shell.execute_reply": "2021-03-07T17:21:03.413908Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: The Hamiltonian constraint.\n",
    "# First declare all needed variables\n",
    "Bq.declare_BSSN_gridfunctions_if_not_declared_already() # Sets trK\n",
    "Bq.BSSN_basic_tensors()   # Sets AbarDD\n",
    "Bq.gammabar__inverse_and_derivs() # Sets gammabarUU\n",
    "Bq.AbarUU_AbarUD_trAbar_AbarDD_dD() # Sets AbarUU and AbarDD_dD\n",
    "Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU() # Sets RbarDD\n",
    "Bq.phi_and_derivs() # Sets phi_dBarD & phi_dBarDD\n",
    "\n",
    "# Term 1: 2/3 K^2\n",
    "H = sp.Rational(2,3)*Bq.trK**2\n",
    "\n",
    "# Term 2: -A_{ij} A^{ij}\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        H += -Bq.AbarDD[i][j]*Bq.AbarUU[i][j]\n",
    "\n",
    "# Term 3a: trace(Rbar)\n",
    "Rbartrace = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Rbartrace += Bq.gammabarUU[i][j]*Bq.RbarDD[i][j]\n",
    "\n",
    "# Term 3b: -8 \\bar{\\gamma}^{ij} \\bar{D}_i \\phi \\bar{D}_j \\phi = -8*phi_dBar_times_phi_dBar\n",
    "# Term 3c: -8 \\bar{\\gamma}^{ij} \\bar{D}_i \\bar{D}_j \\phi      = -8*phi_dBarDD_contraction\n",
    "phi_dBar_times_phi_dBar = sp.sympify(0) # Term 3b\n",
    "phi_dBarDD_contraction  = sp.sympify(0) # Term 3c\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        phi_dBar_times_phi_dBar += Bq.gammabarUU[i][j]*Bq.phi_dBarD[i]*Bq.phi_dBarD[j]\n",
    "        phi_dBarDD_contraction  += Bq.gammabarUU[i][j]*Bq.phi_dBarDD[i][j]\n",
    "\n",
    "# Add Term 3:\n",
    "H += Bq.exp_m4phi*(Rbartrace - 8*(phi_dBar_times_phi_dBar + phi_dBarDD_contraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='momentumconstraint'></a>\n",
    "\n",
    "# Step 3: $\\mathcal{M}^i$, the momentum constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{momentumconstraint}$$\n",
    "\n",
    "***Courtesy Ian Ruchlin***\n",
    "\n",
    "The following definition of the momentum constraint is a simplification of Eq. 47 or [Ruchlin, Etienne, & Baumgarte (2018)](https://arxiv.org/pdf/1712.07658.pdf), which itself was a corrected version of the momentum constraint presented in Eq. 14 of [Baumgarte *et al*](https://arxiv.org/pdf/1211.6632.pdf).\n",
    "\n",
    "Start with the physical momentum constraint\n",
    "$$\n",
    "\\mathcal{M}^{i} \\equiv D_{j} \\left ( K^{i j} - \\gamma^{i j} K \\right ) = 0 \\; .\n",
    "$$\n",
    "Expanding and using metric compatibility with the physical covariant derivative $D_{i}$ yields\n",
    "$$\n",
    "\\mathcal{M}^{i} = D_{j} K^{i j} - \\gamma^{i j} \\partial_{j} K \\; .\n",
    "$$\n",
    "The physical extrinsic curvature $K_{i j}$ is related to the trace-free extrinsic curvature $A_{i j}$ by\n",
    "$$\n",
    "K_{i j} = A_{i j} + \\frac{1}{3} \\gamma_{i j} K \\; .\n",
    "$$\n",
    "Thus,\n",
    "$$\n",
    "\\mathcal{M}^{i} = D_{j} A^{i j} - \\frac{2}{3} \\gamma^{i j} \\partial_{j} K \\; .\n",
    "$$\n",
    "The physical metric $\\gamma_{i j}$ is related to the conformal metric $\\bar{\\gamma}_{i j}$ by the conformal rescaling\n",
    "$$\n",
    "\\gamma_{i j} = e^{4 \\phi} \\bar{\\gamma}_{i j} \\; ,\n",
    "$$\n",
    "and similarly for the trace-free extrinsic curvature\n",
    "$$\n",
    "A_{i j} = e^{4 \\phi} \\bar{A}_{i j} \\; .\n",
    "$$\n",
    "It can be shown (Eq. (3.34) in Baumgarte & Shapiro (2010) with $\\alpha = -4$ and $\\psi = e^{\\phi}$) that the physical and conformal covariant derivatives obey\n",
    "$$\n",
    "D_{j} A^{i j} = e^{-10 \\phi} \\bar{D}_{j} \\left (e^{6 \\phi} \\bar{A}^{i j} \\right ) \\; .\n",
    "$$\n",
    "Then, the constraint becomes\n",
    "$$\n",
    "\\mathcal{M}^i = e^{-4\\phi} \\left(\n",
    "{\\underbrace {\\textstyle \\bar{D}_j \\bar{A}^{ij}}_{\\rm Term\\ 1}} + \n",
    "{\\underbrace {\\textstyle 6 \\bar{A}^{ij}\\partial_j \\phi}_{\\rm Term\\ 2}} - \n",
    "{\\underbrace {\\textstyle \\frac{2}{3} \\bar{\\gamma}^{ij}\\partial_j K}_{\\rm Term\\ 3}}\\right) \\; .\n",
    "$$\n",
    "\n",
    "Let's first implement Terms 2 and 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:03.445640Z",
     "iopub.status.busy": "2021-03-07T17:21:03.444942Z",
     "iopub.status.idle": "2021-03-07T17:21:03.447661Z",
     "shell.execute_reply": "2021-03-07T17:21:03.447144Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3: M^i, the momentum constraint\n",
    "\n",
    "MU = ixp.zerorank1()\n",
    "\n",
    "# Term 2: 6 A^{ij} \\partial_j \\phi:\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        MU[i] += 6*Bq.AbarUU[i][j]*Bq.phi_dD[j]\n",
    "\n",
    "# Term 3: -2/3 \\bar{\\gamma}^{ij} K_{,j}\n",
    "trK_dD = ixp.declarerank1(\"trK_dD\") # Not defined in BSSN_RHSs; only trK_dupD is defined there.\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        MU[i] += -sp.Rational(2,3)*Bq.gammabarUU[i][j]*trK_dD[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we turn our attention to Term 1. The covariant divergence involves upper indices in $\\bar{A}^{i j}$, but it would be easier for us to finite difference the rescaled $\\bar{A}_{i j}$. A simple application of the inverse conformal metric yields\n",
    "$$\n",
    "\\bar{D}_{j} \\bar{A}^{i j} = \\bar{\\gamma}^{i k} \\bar{\\gamma}^{j l} \\bar{D}_{j} \\bar{A}_{k l} \\; .\n",
    "$$\n",
    "As usual, the covariant derivative is related to the ordinary derivative using the conformal Christoffel symbols\n",
    "$$\n",
    "\\bar{D}_{k} \\bar{A}_{i j} = \\partial_{k} \\bar{A}_{i j} - \\bar{\\Gamma}^{l}_{k i} \\bar{A}_{l j} - \\bar{\\Gamma}^{l}_{k j} \\bar{A}_{i l} \\; .\n",
    "$$\n",
    "It is the ordinary derivative above that is approximated by finite difference. The BSSN formulation used here does not rely on spatial derivatives $\\partial_{k} \\bar{A}_{i j}$ in any of the right-hand-sides (except for the advection term, which uses the upwinded derivative), and so we must declare a new ordinary, centered stencil derivative field of rank 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:03.522418Z",
     "iopub.status.busy": "2021-03-07T17:21:03.486570Z",
     "iopub.status.idle": "2021-03-07T17:21:03.620282Z",
     "shell.execute_reply": "2021-03-07T17:21:03.619601Z"
    }
   },
   "outputs": [],
   "source": [
    "# First define aDD_dD:\n",
    "aDD_dD = ixp.declarerank3(\"aDD_dD\",\"sym01\")\n",
    "\n",
    "# Then evaluate the conformal covariant derivative \\bar{D}_j \\bar{A}_{lm}\n",
    "AbarDD_dBarD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            AbarDD_dBarD[i][j][k] = Bq.AbarDD_dD[i][j][k]\n",
    "            for l in range(DIM):\n",
    "                AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][i]*Bq.AbarDD[l][j]\n",
    "                AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][j]*Bq.AbarDD[i][l]\n",
    "\n",
    "# Term 1: Contract twice with the metric to make \\bar{D}_{j} \\bar{A}^{ij}\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            for l in range(DIM):\n",
    "                MU[i] += Bq.gammabarUU[i][k]*Bq.gammabarUU[j][l]*AbarDD_dBarD[k][l][j]\n",
    "\n",
    "# Finally, we multiply by e^{-4 phi} and rescale the momentum constraint:\n",
    "for i in range(DIM):\n",
    "    MU[i] *= Bq.exp_m4phi / rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against `BSSN.BSSN_constraints` 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 the RHSs of the BSSN equations between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.BSSN_constraints](../edit/BSSN/BSSN_constraints.py) module.\n",
    "\n",
    "By default, we analyze these expressions in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:03.626554Z",
     "iopub.status.busy": "2021-03-07T17:21:03.625853Z",
     "iopub.status.idle": "2021-03-07T17:21:05.588244Z",
     "shell.execute_reply": "2021-03-07T17:21:05.587701Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between BSSN_constraints tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
      "H - bssncon.H = 0\n",
      "MU[0] - bssncon.MU[0] = 0\n",
      "MU[1] - bssncon.MU[1] = 0\n",
      "MU[2] - bssncon.MU[2] = 0\n"
     ]
    }
   ],
   "source": [
    "# Step 4: Code Validation against BSSN.BSSN_constraints NRPy+ module\n",
    "\n",
    "# We already have SymPy expressions for BSSN constraints\n",
    "#         in terms of other SymPy variables. Even if we reset the\n",
    "#         list of NRPy+ gridfunctions, these *SymPy* expressions for\n",
    "#         BSSN constraint variables *will remain unaffected*.\n",
    "#\n",
    "#         Here, we will use the above-defined BSSN constraint expressions\n",
    "#         to validate against the same expressions in the\n",
    "#         BSSN/BSSN_constraints.py file, to ensure consistency between\n",
    "#         this tutorial and the module itself.\n",
    "#\n",
    "# Reset the list of gridfunctions, as registering a gridfunction\n",
    "#   twice (in the bssnrhs.BSSN_RHSs() call) will spawn an error.\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "# Call the BSSN_RHSs() function from within the\n",
    "#          BSSN/BSSN_RHSs.py module,\n",
    "#          which should do exactly the same as in Steps 1-16 above.\n",
    "import BSSN.BSSN_constraints as bssncon\n",
    "bssncon.BSSN_constraints()\n",
    "\n",
    "print(\"Consistency check between BSSN_constraints tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "print(\"H - bssncon.H = \" + str(H - bssncon.H))\n",
    "for i in range(DIM):\n",
    "    print(\"MU[\"+str(i)+\"] - bssncon.MU[\"+str(i)+\"] = \" + str(MU[i] - bssncon.MU[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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-BSSN_constraints.pdf](Tutorial-BSSN_constraints.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": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:05.593434Z",
     "iopub.status.busy": "2021-03-07T17:21:05.592342Z",
     "iopub.status.idle": "2021-03-07T17:21:08.879753Z",
     "shell.execute_reply": "2021-03-07T17:21:08.880335Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_constraints.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-BSSN_constraints.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-BSSN_constraints\")"
   ]
  }
 ],
 "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
}