Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.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 Time-Evolution Equations for the Gauge Fields $\\alpha$ and $\\beta^i$\n",
    "\n",
    "## Authors: Zach Etienne & Terrence Pierre Jacques\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "[comment]: <> (Abstract: TODO, or make the introduction an abstract and addiotnal notes section, and write a new Introduction)\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 (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_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py)\n",
    "\n",
    "\n",
    "## Introduction:\n",
    "This tutorial notebook constructs SymPy expressions for the right-hand sides of the time-evolution equations for the gauge fields $\\alpha$ (the lapse, governing how much proper time elapses at each point between one timestep in a 3+1 solution to Einstein's equations and the next) and $\\beta^i$ (the shift, governing how much proper distance numerical grid points move from one timestep in a 3+1 solution to Einstein's equations and the next).\n",
    "\n",
    "Though we are completely free to choose gauge conditions (i.e., free to choose the form of the right-hand sides of the gauge time evolution equations), very few have been found robust in the presence of (puncture) black holes. So we focus here only on a few of the most stable choices.\n",
    "\n"
   ]
  },
  {
   "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](#lapseconditions): Right-hand side of $\\partial_t \\alpha$\n",
    "    1. [Step 2.a](#onepluslog): $1+\\log$ lapse\n",
    "    1. [Step 2.b](#harmonicslicing): Harmonic slicing\n",
    "    1. [Step 2.c](#frozen): Frozen lapse\n",
    "    1. [Step 2.d](#statictrumpet_onepluslog): Alternative 1+log condition for Static Trumpet initial data\n",
    "1. [Step 3](#shiftconditions): Right-hand side of $\\partial_t \\beta^i$: Second-order Gamma-driving shift conditions\n",
    "    1. [Step 3.a](#origgammadriving): Original, non-covariant Gamma-driving shift condition\n",
    "    1. [Step 3.b](#covgammadriving): [Brown](https://arxiv.org/abs/0902.3652)'s suggested covariant Gamma-driving shift condition\n",
    "        1. [Step 3.b.i](#partial_beta): The right-hand side of the $\\partial_t \\beta^i$ equation\n",
    "        1. [Step 3.b.ii](#partial_upper_b): The right-hand side of the $\\partial_t B^i$ equation\n",
    "    1. [Step 3.c](#statictrumpet_nonadvecgammadriving): Non-advecting Gamma-driving shift condition (used for evolving \"Static Trumpet\" initial data)\n",
    "1. [Step 4](#rescale): Rescale right-hand sides of BSSN gauge equations\n",
    "1. [Step 5](#code_validation): Code Validation against `BSSN.BSSN_gauge_RHSs` NRPy+ module\n",
    "1. [Step 6](#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",
    "Let's start by importing all the needed modules from Python/NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:55.800214Z",
     "iopub.status.busy": "2021-03-07T17:21:55.799384Z",
     "iopub.status.idle": "2021-03-07T17:21:58.416995Z",
     "shell.execute_reply": "2021-03-07T17:21:58.416408Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Import all needed modules from NRPy+:\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 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 BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities\n",
    "import BSSN.BSSN_RHSs as Brhs     # NRPy+: Constructs BSSN right-hand-side expressions\n",
    "import sys                        # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step 1.c: Declare/initialize parameters for this module\n",
    "thismodule = \"BSSN_gauge_RHSs\"\n",
    "par.initialize_param(par.glb_param(\"char\", thismodule, \"LapseEvolutionOption\", \"OnePlusLog\"))\n",
    "par.initialize_param(par.glb_param(\"char\", thismodule, \"ShiftEvolutionOption\", \"GammaDriving2ndOrder_Covariant\"))\n",
    "\n",
    "# Step 1.d: 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.e: 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()\n",
    "\n",
    "# Step 1.f: Define BSSN scalars & tensors (in terms of rescaled BSSN quantities)\n",
    "import BSSN.BSSN_quantities as Bq\n",
    "Bq.BSSN_basic_tensors()\n",
    "Bq.betaU_derivs()\n",
    "\n",
    "import BSSN.BSSN_RHSs as Brhs\n",
    "Brhs.BSSN_RHSs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lapseconditions'></a>\n",
    "\n",
    "# Step 2: Right-hand side of $\\partial_t \\alpha$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{lapseconditions}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='onepluslog'></a>\n",
    "\n",
    "## Step 2.a: $1+\\log$ lapse \\[Back to [top](#toc)\\]\n",
    "$$\\label{onepluslog}$$\n",
    "\n",
    "The [$1=\\log$ lapse condition](https://arxiv.org/abs/gr-qc/0206072) is a member of the [Bona-Masso family of lapse choices](https://arxiv.org/abs/gr-qc/9412071), which has the desirable property of singularity avoidance. As is common (e.g., see [Campanelli *et al* (2005)](https://arxiv.org/abs/gr-qc/0511048)), we make the replacement $\\partial_t \\to \\partial_0 = \\partial_t - \\beta^i \\partial_i$ to ensure lapse characteristics advect with the shift. The bracketed term in the $1+\\log$ lapse condition below encodes the shift advection term:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_0 \\alpha &= -2 \\alpha K \\\\\n",
    "\\implies \\partial_t \\alpha &= \\left[\\beta^i \\partial_i \\alpha\\right] - 2 \\alpha K\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.425515Z",
     "iopub.status.busy": "2021-03-07T17:21:58.424823Z",
     "iopub.status.idle": "2021-03-07T17:21:58.427465Z",
     "shell.execute_reply": "2021-03-07T17:21:58.426931Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a: The 1+log lapse condition:\n",
    "#   \\partial_t \\alpha = \\beta^i \\alpha_{,i} - 2*\\alpha*K\n",
    "# First import expressions from BSSN_quantities\n",
    "cf    = Bq.cf\n",
    "trK   = Bq.trK\n",
    "alpha = Bq.alpha\n",
    "betaU = Bq.betaU\n",
    "\n",
    "# Implement the 1+log lapse condition\n",
    "if par.parval_from_str(thismodule+\"::LapseEvolutionOption\") == \"OnePlusLog\":\n",
    "    alpha_rhs = -2*alpha*trK\n",
    "    alpha_dupD = ixp.declarerank1(\"alpha_dupD\")\n",
    "    for i in range(DIM):\n",
    "        alpha_rhs += betaU[i]*alpha_dupD[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='harmonicslicing'></a>\n",
    "\n",
    "## Step 2.b: Harmonic slicing \\[Back to [top](#toc)\\]\n",
    "$$\\label{harmonicslicing}$$\n",
    "\n",
    "As defined on Pg 2 of https://arxiv.org/pdf/gr-qc/9902024.pdf , this is given by \n",
    "\n",
    "$$\n",
    "\\partial_t \\alpha = \\partial_t e^{6 \\phi} = 6 e^{6 \\phi} \\partial_t \\phi\n",
    "$$\n",
    "\n",
    "If \n",
    "\n",
    "$$\\text{cf} = W = e^{-2 \\phi},$$ \n",
    "\n",
    "then\n",
    "\n",
    "$$\n",
    "6 e^{6 \\phi} \\partial_t \\phi = 6 W^{-3} \\partial_t \\phi.\n",
    "$$\n",
    "\n",
    "However,\n",
    "$$\n",
    "\\partial_t \\phi = -\\partial_t \\text{cf} / (2 \\text{cf})$$\n",
    "\n",
    "(as described above), so if `cf`$=W$, then\n",
    "\\begin{align}\n",
    "\\partial_t \\alpha &= 6 e^{6 \\phi} \\partial_t \\phi \\\\\n",
    "&= 6 W^{-3} \\left(-\\frac{\\partial_t W}{2 W}\\right) \\\\\n",
    "&= -3 \\text{cf}^{-4} \\text{cf}\\_\\text{rhs}\n",
    "\\end{align}\n",
    "\n",
    "**Exercise to students: Implement Harmonic slicing for `cf`$=\\chi$** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.433453Z",
     "iopub.status.busy": "2021-03-07T17:21:58.432662Z",
     "iopub.status.idle": "2021-03-07T17:21:58.434925Z",
     "shell.execute_reply": "2021-03-07T17:21:58.435395Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.b: Implement the harmonic slicing lapse condition\n",
    "if par.parval_from_str(thismodule+\"::LapseEvolutionOption\") == \"HarmonicSlicing\":\n",
    "    if par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\") == \"W\":\n",
    "        alpha_rhs = -3*cf**(-4)*Brhs.cf_rhs\n",
    "    elif par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\") == \"phi\":\n",
    "        alpha_rhs = 6*sp.exp(6*cf)*Brhs.cf_rhs\n",
    "    else:\n",
    "        print(\"Error LapseEvolutionOption==HarmonicSlicing unsupported for EvolvedConformalFactor_cf!=(W or phi)\")\n",
    "        sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='frozen'></a>\n",
    "\n",
    "## Step 2.c: Frozen lapse \\[Back to [top](#toc)\\]\n",
    "$$\\label{frozen}$$\n",
    "\n",
    "This slicing condition is given by\n",
    "$$\\partial_t \\alpha = 0,$$\n",
    "\n",
    "which is rarely a stable lapse condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.440683Z",
     "iopub.status.busy": "2021-03-07T17:21:58.439922Z",
     "iopub.status.idle": "2021-03-07T17:21:58.442993Z",
     "shell.execute_reply": "2021-03-07T17:21:58.443943Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.c: Frozen lapse\n",
    "#    \\partial_t \\alpha = 0\n",
    "if par.parval_from_str(thismodule+\"::LapseEvolutionOption\") == \"Frozen\":\n",
    "    alpha_rhs = sp.sympify(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='statictrumpet_onepluslog'></a>\n",
    "\n",
    "## Step 2.d: Alternative $1+\\log$ lapse for Static Trumpet initial data \\[Back to [top](#toc)\\]\n",
    "$$\\label{statictrumpet_onepluslog}$$\n",
    "\n",
    "An alternative to the standard 1+log condition to be used with Static Trumpet initial data, the lapse is evolved according to a\n",
    "condition consistent with staticity, given by equation 67 in [Ruchlin, Etienne, & Baumgarte (2018)](https://arxiv.org/pdf/1712.07658.pdf)\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_0 \\alpha &= -\\alpha(1 - \\alpha) K \\\\\n",
    "\\implies \\partial_t \\alpha &= \\left[\\beta^i \\partial_i \\alpha\\right] -\\alpha(1 - \\alpha) K\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.452895Z",
     "iopub.status.busy": "2021-03-07T17:21:58.452086Z",
     "iopub.status.idle": "2021-03-07T17:21:58.454568Z",
     "shell.execute_reply": "2021-03-07T17:21:58.455119Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.d: Alternative 1+log lapse condition:\n",
    "#   \\partial_t \\alpha = \\beta^i \\alpha_{,i} -\\alpha*(1 - \\alpha)*K\n",
    "\n",
    "# Implement the alternative 1+log lapse condition\n",
    "if par.parval_from_str(thismodule+\"::LapseEvolutionOption\") == \"OnePlusLogAlt\":\n",
    "    alpha_rhs = -alpha*(1 - alpha)*trK\n",
    "    alpha_dupD = ixp.declarerank1(\"alpha_dupD\")\n",
    "    for i in range(DIM):\n",
    "        alpha_rhs += betaU[i]*alpha_dupD[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='shiftconditions'></a>\n",
    "\n",
    "# Step 3: Right-hand side of $\\partial_t \\beta^i$: Second-order Gamma-driving shift conditions \\[Back to [top](#toc)\\]\n",
    "$$\\label{shiftconditions}$$\n",
    "\n",
    "The motivation behind Gamma-driving shift conditions are well documented in the book [*Numerical Relativity* by Baumgarte & Shapiro](https://www.amazon.com/Numerical-Relativity-Einsteins-Equations-Computer/dp/052151407X/).\n",
    "\n",
    "<a id='origgammadriving'></a>\n",
    "\n",
    "## Step 3.a: Original, non-covariant Gamma-driving shift condition \\[Back to [top](#toc)\\]\n",
    "$$\\label{origgammadriving}$$\n",
    "\n",
    "**Option 1: Non-Covariant, Second-Order Shift**\n",
    "\n",
    "We adopt the [*shifting (i.e., advecting) shift*](https://arxiv.org/abs/gr-qc/0605030) non-covariant, second-order shift condition:\n",
    "\\begin{align}\n",
    "\\partial_0 \\beta^i &= B^{i} \\\\\n",
    "\\partial_0 B^i &= \\frac{3}{4} \\partial_{0} \\bar{\\Lambda}^{i} - \\eta B^{i} \\\\\n",
    "\\implies \\partial_t \\beta^i &= \\left[\\beta^j \\partial_j \\beta^i\\right] + B^{i} \\\\\n",
    "\\partial_t B^i &= \\left[\\beta^j \\partial_j B^i\\right] + \\frac{3}{4} \\partial_{0} \\bar{\\Lambda}^{i} - \\eta B^{i},\n",
    "\\end{align}\n",
    "where $\\eta$ is the shift damping parameter, and $\\partial_{0} \\bar{\\Lambda}^{i}$ in the right-hand side of the $\\partial_{0} B^{i}$ equation is computed by adding $\\beta^j \\partial_j \\bar{\\Lambda}^i$ to the right-hand side expression given for $\\partial_t \\bar{\\Lambda}^i$ in the BSSN time-evolution equations as listed [here](Tutorial-BSSN_formulation.ipynb), so no explicit time dependence occurs in the right-hand sides of the BSSN evolution equations and the Method of Lines can be applied directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.468721Z",
     "iopub.status.busy": "2021-03-07T17:21:58.467870Z",
     "iopub.status.idle": "2021-03-07T17:21:58.470825Z",
     "shell.execute_reply": "2021-03-07T17:21:58.470193Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.a: Set \\partial_t \\beta^i\n",
    "# First import expressions from BSSN_quantities\n",
    "BU         = Bq.BU\n",
    "betU       = Bq.betU\n",
    "betaU_dupD = Bq.betaU_dupD\n",
    "# Define needed quantities\n",
    "beta_rhsU = ixp.zerorank1()\n",
    "B_rhsU = ixp.zerorank1()\n",
    "if par.parval_from_str(thismodule+\"::ShiftEvolutionOption\") == \"GammaDriving2ndOrder_NoCovariant\":\n",
    "    # Step 3.a.i: Compute right-hand side of beta^i\n",
    "    # *  \\partial_t \\beta^i = \\beta^j \\beta^i_{,j} + B^i\n",
    "    for i in range(DIM):\n",
    "        beta_rhsU[i] += BU[i]\n",
    "        for j in range(DIM):\n",
    "            beta_rhsU[i] += betaU[j]*betaU_dupD[i][j]\n",
    "    # Compute right-hand side of B^i:\n",
    "    eta = par.Cparameters(\"REAL\", thismodule, [\"eta\"],2.0)\n",
    "\n",
    "    # Step 3.a.ii: Compute right-hand side of B^i\n",
    "    # *  \\partial_t B^i     = \\beta^j B^i_{,j} + 3/4 * \\partial_0 \\Lambda^i - eta B^i\n",
    "    # Step 3.a.iii: Define BU_dupD, in terms of derivative of rescaled variable \\bet^i\n",
    "    BU_dupD = ixp.zerorank2()\n",
    "    betU_dupD = ixp.declarerank2(\"betU_dupD\",\"nosym\")\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            BU_dupD[i][j] = betU_dupD[i][j]*rfm.ReU[i] + betU[i]*rfm.ReUdD[i][j]\n",
    "\n",
    "    # Step 3.a.iv: Compute \\partial_0 \\bar{\\Lambda}^i = (\\partial_t - \\beta^i \\partial_i) \\bar{\\Lambda}^j\n",
    "    Lambdabar_partial0 = ixp.zerorank1()\n",
    "    for i in range(DIM):\n",
    "        Lambdabar_partial0[i] = Brhs.Lambdabar_rhsU[i]\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            Lambdabar_partial0[j] += -betaU[i]*Brhs.LambdabarU_dupD[j][i]\n",
    "\n",
    "    # Step 3.a.v: Evaluate RHS of B^i:\n",
    "    for i in range(DIM):\n",
    "        B_rhsU[i] += sp.Rational(3,4)*Lambdabar_partial0[i] - eta*BU[i]\n",
    "        for j in range(DIM):\n",
    "            B_rhsU[i] += betaU[j]*BU_dupD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='covgammadriving'></a>\n",
    "\n",
    "## Step 3.b: [Brown](https://arxiv.org/abs/0902.3652)'s suggested covariant Gamma-driving shift condition \\[Back to [top](#toc)\\]\n",
    "$$\\label{covgammadriving}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='partial_beta'></a>\n",
    "\n",
    "### Step 3.b.i: The right-hand side of the $\\partial_t \\beta^i$ equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{partial_beta}$$\n",
    "\n",
    "This is [Brown's](https://arxiv.org/abs/0902.3652) suggested formulation (Eq. 20b; note that Eq. 20a is the same as our lapse condition, as $\\bar{D}_j \\alpha = \\partial_j \\alpha$ for scalar $\\alpha$):\n",
    "$$\\partial_t \\beta^i = \\left[\\beta^j \\bar{D}_j \\beta^i\\right] + B^{i}$$\n",
    "Based on the definition of covariant derivative, we have\n",
    "$$\n",
    "\\bar{D}_{j} \\beta^{i} = \\beta^i_{,j} + \\bar{\\Gamma}^i_{mj} \\beta^m,\n",
    "$$\n",
    "so the above equation becomes\n",
    "\\begin{align}\n",
    "\\partial_t \\beta^i &= \\left[\\beta^j \\left(\\beta^i_{,j} + \\bar{\\Gamma}^i_{mj} \\beta^m\\right)\\right] + B^{i}\\\\\n",
    "&= {\\underbrace {\\textstyle \\beta^j \\beta^i_{,j}}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\beta^j \\bar{\\Gamma}^i_{mj} \\beta^m}_{\\text{Term 2}}} + \n",
    "{\\underbrace {\\textstyle B^i}_{\\text{Term 3}}} \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.545297Z",
     "iopub.status.busy": "2021-03-07T17:21:58.509414Z",
     "iopub.status.idle": "2021-03-07T17:21:58.603517Z",
     "shell.execute_reply": "2021-03-07T17:21:58.602939Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.b: The right-hand side of the \\partial_t \\beta^i equation\n",
    "if par.parval_from_str(thismodule+\"::ShiftEvolutionOption\") == \"GammaDriving2ndOrder_Covariant\":\n",
    "    # Step 3.b Option 2: \\partial_t \\beta^i = \\left[\\beta^j \\bar{D}_j \\beta^i\\right] + B^{i}\n",
    "    # First we need GammabarUDD, defined in Bq.gammabar__inverse_and_derivs()\n",
    "    Bq.gammabar__inverse_and_derivs()\n",
    "    GammabarUDD = Bq.GammabarUDD\n",
    "    # Then compute right-hand side:\n",
    "    # Term 1: \\beta^j \\beta^i_{,j}\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            beta_rhsU[i] += betaU[j]*betaU_dupD[i][j]\n",
    "\n",
    "    # Term 2: \\beta^j \\bar{\\Gamma}^i_{mj} \\beta^m\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                beta_rhsU[i] += betaU[j]*GammabarUDD[i][m][j]*betaU[m]\n",
    "    # Term 3: B^i\n",
    "    for i in range(DIM):\n",
    "        beta_rhsU[i] += BU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='partial_upper_b'></a>\n",
    "\n",
    "### Step 3.b.ii: The right-hand side of the $\\partial_t B^i$ equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{partial_upper_b}$$\n",
    "\n",
    "$$\\partial_t B^i = \\left[\\beta^j \\bar{D}_j B^i\\right] + \\frac{3}{4}\\left( \\partial_t \\bar{\\Lambda}^{i} - \\beta^j \\bar{D}_j \\bar{\\Lambda}^{i} \\right) - \\eta B^{i}$$\n",
    "\n",
    "Based on the definition of covariant derivative, we have for vector $V^i$\n",
    "$$\n",
    "\\bar{D}_{j} V^{i} = V^i_{,j} + \\bar{\\Gamma}^i_{mj} V^m,\n",
    "$$\n",
    "so the above equation becomes\n",
    "\\begin{align}\n",
    "\\partial_t B^i &= \\left[\\beta^j \\left(B^i_{,j} + \\bar{\\Gamma}^i_{mj} B^m\\right)\\right] + \\frac{3}{4}\\left[ \\partial_t \\bar{\\Lambda}^{i} - \\beta^j \\left(\\bar{\\Lambda}^i_{,j} + \\bar{\\Gamma}^i_{mj} \\bar{\\Lambda}^m\\right) \\right] - \\eta B^{i} \\\\\n",
    "&= {\\underbrace {\\textstyle \\beta^j B^i_{,j}}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\beta^j \\bar{\\Gamma}^i_{mj} B^m}_{\\text{Term 2}}} + \n",
    "{\\underbrace {\\textstyle \\frac{3}{4}\\partial_t \\bar{\\Lambda}^{i}}_{\\text{Term 3}}} -\n",
    "{\\underbrace {\\textstyle \\frac{3}{4}\\beta^j \\bar{\\Lambda}^i_{,j}}_{\\text{Term 4}}} -\n",
    "{\\underbrace {\\textstyle \\frac{3}{4}\\beta^j \\bar{\\Gamma}^i_{mj} \\bar{\\Lambda}^m}_{\\text{Term 5}}} -\n",
    "{\\underbrace {\\textstyle \\eta B^i}_{\\text{Term 6}}}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.678281Z",
     "iopub.status.busy": "2021-03-07T17:21:58.642565Z",
     "iopub.status.idle": "2021-03-07T17:21:58.700215Z",
     "shell.execute_reply": "2021-03-07T17:21:58.700720Z"
    }
   },
   "outputs": [],
   "source": [
    "if par.parval_from_str(thismodule+\"::ShiftEvolutionOption\") == \"GammaDriving2ndOrder_Covariant\":\n",
    "    # Step 3.c: Covariant option:\n",
    "    #  \\partial_t B^i = \\beta^j \\bar{D}_j B^i\n",
    "    #               + \\frac{3}{4} ( \\partial_t \\bar{\\Lambda}^{i} - \\beta^j \\bar{D}_j \\bar{\\Lambda}^{i} )\n",
    "    #               - \\eta B^{i}\n",
    "    #                 = \\beta^j B^i_{,j} + \\beta^j \\bar{\\Gamma}^i_{mj} B^m\n",
    "    #               + \\frac{3}{4}[ \\partial_t \\bar{\\Lambda}^{i}\n",
    "    #                            - \\beta^j (\\bar{\\Lambda}^i_{,j} + \\bar{\\Gamma}^i_{mj} \\bar{\\Lambda}^m)]\n",
    "    #               - \\eta B^{i}\n",
    "    # Term 1, part a: First compute B^i_{,j} using upwinded derivative\n",
    "    BU_dupD = ixp.zerorank2()\n",
    "    betU_dupD = ixp.declarerank2(\"betU_dupD\",\"nosym\")\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            BU_dupD[i][j] = betU_dupD[i][j]*rfm.ReU[i] + betU[i]*rfm.ReUdD[i][j]\n",
    "    # Term 1: \\beta^j B^i_{,j}\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            B_rhsU[i] += betaU[j]*BU_dupD[i][j]\n",
    "    # Term 2: \\beta^j \\bar{\\Gamma}^i_{mj} B^m\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                B_rhsU[i] += betaU[j]*GammabarUDD[i][m][j]*BU[m]\n",
    "    # Term 3: \\frac{3}{4}\\partial_t \\bar{\\Lambda}^{i}\n",
    "    for i in range(DIM):\n",
    "        B_rhsU[i] += sp.Rational(3,4)*Brhs.Lambdabar_rhsU[i]\n",
    "    # Term 4: -\\frac{3}{4}\\beta^j \\bar{\\Lambda}^i_{,j}\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            B_rhsU[i] += -sp.Rational(3,4)*betaU[j]*Brhs.LambdabarU_dupD[i][j]\n",
    "    # Term 5: -\\frac{3}{4}\\beta^j \\bar{\\Gamma}^i_{mj} \\bar{\\Lambda}^m\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                B_rhsU[i] += -sp.Rational(3,4)*betaU[j]*GammabarUDD[i][m][j]*Bq.LambdabarU[m]\n",
    "    # Term 6: - \\eta B^i\n",
    "    # eta is a free parameter; we declare it here:\n",
    "    eta = par.Cparameters(\"REAL\", thismodule, [\"eta\"],2.0)\n",
    "    for i in range(DIM):\n",
    "        B_rhsU[i] += -eta*BU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='statictrumpet_nonadvecgammadriving'></a>\n",
    "\n",
    "## Step 3.c: Non-advecting Gamma-driving shift condition (used for evolving \"Static Trumpet\" initial data) \\[Back to [top](#toc)\\]\n",
    "$$\\label{statictrumpet_nonadvecgammadriving}$$\n",
    "\n",
    "\n",
    "For the shift vector evolution equation, we desire only that the right-hand sides vanish analytically (although numerical error is expected to result in specious evolution). To this end, we adopt the nonadvecting Gamma-driver condition, given by equations 68a and 68b in [Ruchlin, Etienne, & Baumgarte (2018)](https://arxiv.org/pdf/1712.07658.pdf)\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t \\beta^i &= B^{i} \\\\\n",
    "\\partial_t B^i &= \\frac{3}{4} \\partial_{t} \\bar{\\Lambda}^{i} - \\eta B^{i},\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.707000Z",
     "iopub.status.busy": "2021-03-07T17:21:58.706332Z",
     "iopub.status.idle": "2021-03-07T17:21:58.708570Z",
     "shell.execute_reply": "2021-03-07T17:21:58.709124Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.c: Set \\partial_t \\beta^i\n",
    "\n",
    "if par.parval_from_str(thismodule+\"::ShiftEvolutionOption\") == \"NonAdvectingGammaDriving\":\n",
    "    # Step 3.c.i: Compute right-hand side of beta^i\n",
    "    # *  \\partial_t \\beta^i = B^i\n",
    "    for i in range(DIM):\n",
    "        beta_rhsU[i] += BU[i]\n",
    "\n",
    "    # Compute right-hand side of B^i:\n",
    "    eta = par.Cparameters(\"REAL\", thismodule, [\"eta\"],2.0)\n",
    "\n",
    "    # Step 3.c.ii: Compute right-hand side of B^i\n",
    "    # *  \\partial_t B^i     = 3/4 * \\partial_t \\Lambda^i - eta B^i\n",
    "    # Step 3.c.iii: Evaluate RHS of B^i:\n",
    "    for i in range(DIM):\n",
    "        B_rhsU[i] += sp.Rational(3,4)*Brhs.Lambdabar_rhsU[i] - eta*BU[i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rescale'></a>\n",
    "\n",
    "# Step 4: Rescale right-hand sides of BSSN gauge equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{rescale}$$\n",
    "\n",
    "Next we rescale the right-hand sides of the BSSN equations so that the evolved variables are $\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.762715Z",
     "iopub.status.busy": "2021-03-07T17:21:58.747634Z",
     "iopub.status.idle": "2021-03-07T17:21:58.765074Z",
     "shell.execute_reply": "2021-03-07T17:21:58.765648Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Step 4: Rescale the BSSN gauge RHS quantities so that the evolved\n",
    "#         variables may remain smooth across coord singularities\n",
    "vet_rhsU    = ixp.zerorank1()\n",
    "bet_rhsU    = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    vet_rhsU[i]    =   beta_rhsU[i] / rfm.ReU[i]\n",
    "    bet_rhsU[i]    =      B_rhsU[i] / rfm.ReU[i]\n",
    "#print(str(Abar_rhsDD[2][2]).replace(\"**\",\"^\").replace(\"_\",\"\").replace(\"xx\",\"x\").replace(\"sin(x2)\",\"Sin[x2]\").replace(\"sin(2*x2)\",\"Sin[2*x2]\").replace(\"cos(x2)\",\"Cos[x2]\").replace(\"detgbaroverdetghat\",\"detg\"))\n",
    "#print(str(Dbarbetacontraction).replace(\"**\",\"^\").replace(\"_\",\"\").replace(\"xx\",\"x\").replace(\"sin(x2)\",\"Sin[x2]\").replace(\"detgbaroverdetghat\",\"detg\"))\n",
    "#print(betaU_dD)\n",
    "#print(str(trK_rhs).replace(\"xx2\",\"xx3\").replace(\"xx1\",\"xx2\").replace(\"xx0\",\"xx1\").replace(\"**\",\"^\").replace(\"_\",\"\").replace(\"sin(xx2)\",\"Sinx2\").replace(\"xx\",\"x\").replace(\"sin(2*x2)\",\"Sin2x2\").replace(\"cos(x2)\",\"Cosx2\").replace(\"detgbaroverdetghat\",\"detg\"))\n",
    "#print(str(bet_rhsU[0]).replace(\"xx2\",\"xx3\").replace(\"xx1\",\"xx2\").replace(\"xx0\",\"xx1\").replace(\"**\",\"^\").replace(\"_\",\"\").replace(\"sin(xx2)\",\"Sinx2\").replace(\"xx\",\"x\").replace(\"sin(2*x2)\",\"Sin2x2\").replace(\"cos(x2)\",\"Cosx2\").replace(\"detgbaroverdetghat\",\"detg\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 5: Code Validation against `BSSN.BSSN_gauge_RHSs` 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 gauge equations between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.BSSN_gauge_RHSs](../edit/BSSN/BSSN_gauge_RHSs.py) module.\n",
    "\n",
    "By default, we analyze the RHSs in Spherical coordinates and with the covariant Gamma-driving second-order shift condition, though other coordinate systems & gauge conditions may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:58.772900Z",
     "iopub.status.busy": "2021-03-07T17:21:58.772178Z",
     "iopub.status.idle": "2021-03-07T17:21:59.150639Z",
     "shell.execute_reply": "2021-03-07T17:21:59.151117Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between BSSN.BSSN_gauge_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
      "alpha_rhs - bssnrhs.alpha_rhs = 0\n",
      "vet_rhsU[0] - bssnrhs.vet_rhsU[0] = 0\n",
      "bet_rhsU[0] - bssnrhs.bet_rhsU[0] = 0\n",
      "vet_rhsU[1] - bssnrhs.vet_rhsU[1] = 0\n",
      "bet_rhsU[1] - bssnrhs.bet_rhsU[1] = 0\n",
      "vet_rhsU[2] - bssnrhs.vet_rhsU[2] = 0\n",
      "bet_rhsU[2] - bssnrhs.bet_rhsU[2] = 0\n"
     ]
    }
   ],
   "source": [
    "# Step 5: We already have SymPy expressions for BSSN gauge RHS expressions\n",
    "#         in terms of other SymPy variables. Even if we reset the\n",
    "#         list of NRPy+ gridfunctions, these *SymPy* expressions for\n",
    "#         BSSN RHS variables *will remain unaffected*.\n",
    "#\n",
    "#         Here, we will use the above-defined BSSN gauge RHS expressions\n",
    "#         to validate against the same expressions in the\n",
    "#         BSSN/BSSN_gauge_RHSs.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 will spawn an error.\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "\n",
    "# Step 5.a: Call the BSSN_gauge_RHSs() function from within the\n",
    "#           BSSN/BSSN_gauge_RHSs.py module,\n",
    "#           which should generate exactly the same expressions as above.\n",
    "import BSSN.BSSN_gauge_RHSs as Bgrhs\n",
    "par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption\",\"GammaDriving2ndOrder_Covariant\")\n",
    "Bgrhs.BSSN_gauge_RHSs()\n",
    "\n",
    "print(\"Consistency check between BSSN.BSSN_gauge_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "print(\"alpha_rhs - bssnrhs.alpha_rhs = \" + str(alpha_rhs - Bgrhs.alpha_rhs))\n",
    "\n",
    "for i in range(DIM):\n",
    "    print(\"vet_rhsU[\"+str(i)+\"] - bssnrhs.vet_rhsU[\"+str(i)+\"] = \" + str(vet_rhsU[i] - Bgrhs.vet_rhsU[i]))\n",
    "    print(\"bet_rhsU[\"+str(i)+\"] - bssnrhs.bet_rhsU[\"+str(i)+\"] = \" + str(bet_rhsU[i] - Bgrhs.bet_rhsU[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: 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_time_evolution-BSSN_gauge_RHSs.pdf](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.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:21:59.155623Z",
     "iopub.status.busy": "2021-03-07T17:21:59.154923Z",
     "iopub.status.idle": "2021-03-07T17:22:03.242672Z",
     "shell.execute_reply": "2021-03-07T17:22:03.243409Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.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_time_evolution-BSSN_gauge_RHSs\")"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}