Tutorial-BSSN_time_evolution-BSSN_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",
    "# The BSSN Time-Evolution Equations\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:** 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_RHS.py](../edit/BSSN/BSSN_RHSs.py)\n",
    "\n",
    "## Introduction:\n",
    "This module documents and constructs the time evolution equations of the BSSN formulation of Einstein's equations, as defined  in [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658) (see also [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632)). \n",
    "\n",
    "**This module is part of the following set of NRPy+ tutorial notebooks on the BSSN formulation of general relativity:**\n",
    "\n",
    "* An overview of the BSSN formulation of Einstein's equations, as well as links for background reading/lectures, are provided in [the NRPy+ tutorial notebook on the BSSN formulation](Tutorial-BSSN_formulation.ipynb).\n",
    "* Basic BSSN quantities are defined in the [BSSN quantities NRPy+ tutorial notebook](Tutorial-BSSN_quantities.ipynb).\n",
    "* Other BSSN equation tutorial notebooks:\n",
    "    * [Time-evolution equations the BSSN gauge quantities $\\alpha$, $\\beta^i$, and $B^i$](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb).\n",
    "    * [BSSN Hamiltonian and momentum constraints](Tutorial-BSSN_constraints.ipynb)\n",
    "    * [Enforcing the $\\bar{\\gamma} = \\hat{\\gamma}$ constraint](Tutorial-BSSN-Enforcing_Determinant_gammabar_equals_gammahat_Constraint.ipynb)\n",
    "\n",
    "### A Note on Notation\n",
    "\n",
    "As is standard in NRPy+, \n",
    "\n",
    "* Greek indices refer to four-dimensional quantities where the zeroth component indicates temporal (time) component.\n",
    "* Latin indices refer to three-dimensional quantities. This is somewhat counterintuitive since Python always indexes its lists starting from 0. As a result, the zeroth component of three-dimensional quantities will necessarily indicate the first *spatial* direction.\n",
    "\n",
    "As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial notebook)."
   ]
  },
  {
   "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",
    "0. [Preliminaries](#bssntimeevolequations): BSSN time-evolution equations, as described in the [BSSN formulation NRPy+ tutorial notebook](Tutorial-BSSN_formulation.ipynb)\n",
    "1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules\n",
    "1. [Step 2](#gammabar): Right-hand side of $\\partial_t \\bar{\\gamma}_{ij}$\n",
    "    1. [Step 2.a](#term1_partial_gamma): Term 1 of $\\partial_t \\bar{\\gamma}_{i j}$\n",
    "    1. [Step 2.b](#term2_partial_gamma): Term 2 of $\\partial_t \\bar{\\gamma}_{i j}$\n",
    "    1. [Step 2.c](#term3_partial_gamma): Term 3 of $\\partial_t \\bar{\\gamma}_{i j}$\n",
    "1. [Step 3](#abar): Right-hand side of $\\partial_t \\bar{A}_{ij}$\n",
    "    1. [Step 3.a](#term1_partial_upper_a): Term 1 of $\\partial_t \\bar{A}_{i j}$\n",
    "    1. [Step 3.c](#term2_partial_upper_a): Term 2 of $\\partial_t \\bar{A}_{i j}$\n",
    "    1. [Step 3.c](#term3_partial_upper_a): Term 3 of $\\partial_t \\bar{A}_{i j}$\n",
    "1. [Step 4](#cf): Right-hand side of $\\partial_t \\phi \\to \\partial_t (\\text{cf})$\n",
    "1. [Step 5](#trk): Right-hand side of $\\partial_t \\text{tr} K$\n",
    "1. [Step 6](#lambdabar): Right-hand side of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.a](#term1_partial_lambda): Term 1 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.b](#term2_partial_lambda): Term 2 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.c](#term3_partial_lambda): Term 3 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.d](#term4_partial_lambda): Term 4 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.e](#term5_partial_lambda): Term 5 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.f](#term6_partial_lambda): Term 6 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "    1. [Step 6.g](#term7_partial_lambda): Term 7 of $\\partial_t \\bar{\\Lambda}^i$\n",
    "1. [Step 7](#rescalingrhss): Rescaling the BSSN right-hand sides; rewriting them in terms of the rescaled quantities $\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$\n",
    "1. [Step 8](#code_validation): Code Validation against `BSSN.BSSN_RHSs` NRPy+ module\n",
    "1. [Step 9](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssntimeevolequations'></a>\n",
    "\n",
    "# Preliminaries: BSSN time-evolution equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssntimeevolequations}$$\n",
    "\n",
    "As described in the [BSSN formulation NRPy+ tutorial notebook](Tutorial-BSSN_formulation.ipynb), the BSSN time-evolution equations are given by\n",
    "\n",
    "\\begin{align}\n",
    "  \\partial_t \\bar{\\gamma}_{i j} {} = {} & \\left[\\beta^k \\partial_k \\bar{\\gamma}_{ij} + \\partial_i \\beta^k \\bar{\\gamma}_{kj} + \\partial_j \\beta^k \\bar{\\gamma}_{ik} \\right] + \\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right ) - 2 \\alpha \\bar{A}_{i j} \\; , \\\\\n",
    "  \\partial_t \\bar{A}_{i j} {} = {} & \\left[\\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik} \\right] - \\frac{2}{3} \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} {\\bar{A}^{k}}_{j} + \\alpha \\bar{A}_{i j} K \\nonumber \\\\\n",
    "  & + e^{-4 \\phi} \\left \\{-2 \\alpha \\bar{D}_{i} \\bar{D}_{j} \\phi + 4 \\alpha \\bar{D}_{i} \\phi \\bar{D}_{j} \\phi  + 4 \\bar{D}_{(i} \\alpha \\bar{D}_{j)} \\phi - \\bar{D}_{i} \\bar{D}_{j} \\alpha + \\alpha \\bar{R}_{i j} \\right \\}^{\\text{TF}} \\; , \\\\\n",
    "  \\partial_t \\phi {} = {} & \\left[\\beta^k \\partial_k \\phi \\right] + \\frac{1}{6} \\left (\\bar{D}_{k} \\beta^{k} - \\alpha K \\right ) \\; , \\\\\n",
    "  \\partial_{t} K {} = {} & \\left[\\beta^k \\partial_k K \\right] + \\frac{1}{3} \\alpha K^{2} + \\alpha \\bar{A}_{i j} \\bar{A}^{i j} - e^{-4 \\phi} \\left (\\bar{D}_{i} \\bar{D}^{i} \\alpha + 2 \\bar{D}^{i} \\alpha \\bar{D}_{i} \\phi \\right ) \\; , \\\\\n",
    "  \\partial_t \\bar{\\Lambda}^{i} {} = {} & \\left[\\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k \\right] + \\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i} + \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j} + \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j} \\nonumber \\\\\n",
    "  & - 2 \\bar{A}^{i j} \\left (\\partial_{j} \\alpha - 6 \\partial_{j} \\phi \\right ) + 2 \\alpha \\bar{A}^{j k} \\Delta_{j k}^{i}  -\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K\n",
    "\\end{align}\n",
    "\n",
    "where the Lie derivative terms (often seen on the left-hand side of these equations) are enclosed in square braces.\n",
    "\n",
    "Notice that the shift advection operator $\\beta^k \\partial_k \\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}\\right\\}$ appears on the right-hand side of *every* expression. As the shift determines how the spatial coordinates $x^i$ move on the next 3D slice of our 4D manifold, we find that representing $\\partial_k$ in these shift advection terms via an *upwinded* finite difference stencil results in far lower numerical errors. This trick is implemented below in all shift advection terms. Upwinded derivatives are indicated in NRPy+ by the `_dupD` variable suffix.\n",
    "\n",
    "\n",
    "As discussed in the [NRPy+ tutorial notebook on BSSN quantities](Tutorial-BSSN_quantities.ipynb), tensorial expressions can diverge at coordinate singularities, so each tensor in the set of BSSN variables\n",
    "\n",
    "$$\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}, \\alpha, \\beta^i, B^i\\right\\},$$\n",
    "\n",
    "is written in terms of the corresponding rescaled quantity in the set\n",
    "\n",
    "$$\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\},$$ \n",
    "\n",
    "respectively, as defined in the [BSSN quantities tutorial](Tutorial-BSSN_quantities.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1.a: import all needed modules from Python/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 sys                        # Standard Python modules for multiplatform OS-level functions\n",
    "import BSSN.BSSN_quantities as Bq # NRPy+: Basic BSSN quantities\n",
    "\n",
    "# Step 1.b: Set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "\n",
    "# Step 1.c: 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.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: Import all basic (unrescaled) BSSN scalars & tensors\n",
    "Bq.BSSN_basic_tensors()\n",
    "gammabarDD = Bq.gammabarDD\n",
    "AbarDD     = Bq.AbarDD\n",
    "LambdabarU = Bq.LambdabarU\n",
    "trK   = Bq.trK\n",
    "alpha = Bq.alpha\n",
    "betaU = Bq.betaU\n",
    "\n",
    "# Step 1.f: Import all needed rescaled BSSN tensors:\n",
    "cf  = Bq.cf\n",
    "lambdaU = Bq.lambdaU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='gammabar'></a>\n",
    "\n",
    "# Step 2: Right-hand side of $\\partial_t \\bar{\\gamma}_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{gammabar}$$\n",
    "\n",
    "Let's start with\n",
    "\n",
    "$$\n",
    "\\partial_t \\bar{\\gamma}_{i j} = \n",
    "{\\underbrace {\\textstyle \\left[\\beta^k \\partial_k \\bar{\\gamma}_{ij} + \\partial_i \\beta^k \\bar{\\gamma}_{kj} + \\partial_j \\beta^k \\bar{\\gamma}_{ik} \\right]}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right )}_{\\text{Term 2}}}\n",
    "{\\underbrace {\\textstyle -2 \\alpha \\bar{A}_{i j}}_{\\text{Term 3}}}.\n",
    "$$\n",
    "\n",
    "\n",
    "First, note that the term containing $\\bar{A}_{k}^{k}$ is there to drive $\\bar{A}_{k}^{k}$ to zero. It was implemented in this form in T. Baumgarte's BSSN-in-spherical-coordinates code, against which NRPy+ and SENR were originally validated. You will find this term in [Brown](https://arxiv.org/abs/0902.3652)'s covariant formulation of the BSSN time-evolution equations (see third term in Eq 21a)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term1_partial_gamma'></a>\n",
    "\n",
    "## Step 2.a: Term 1 of $\\partial_t \\bar{\\gamma}_{i j}$  \\[Back to [top](#toc)\\] \n",
    "$$\\label{term1_partial_gamma}$$\n",
    "\n",
    "Term 1 of $\\partial_t \\bar{\\gamma}_{i j} =$ `gammabar_rhsDD[i][j]`: $\\beta^k \\bar{\\gamma}_{ij,k} + \\beta^k_{,i} \\bar{\\gamma}_{kj} + \\beta^k_{,j} \\bar{\\gamma}_{ik}$\n",
    "\n",
    "\n",
    "First we import derivative expressions for betaU defined in the [NRPy+ BSSN quantities tutorial notebook](Tutorial-BSSN_quantities.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.a.i: Import derivative expressions for betaU defined in the BSSN.BSSN_quantities module:\n",
    "Bq.betaU_derivs()\n",
    "betaU_dD = Bq.betaU_dD\n",
    "betaU_dupD = Bq.betaU_dupD\n",
    "betaU_dDD = Bq.betaU_dDD\n",
    "# Step 2.a.ii: Import derivative expression for gammabarDD\n",
    "Bq.gammabar__inverse_and_derivs()\n",
    "gammabarDD_dupD = Bq.gammabarDD_dupD\n",
    "\n",
    "# Step 2.a.iii: First term of \\partial_t \\bar{\\gamma}_{i j} right-hand side:\n",
    "# \\beta^k \\bar{\\gamma}_{ij,k} + \\beta^k_{,i} \\bar{\\gamma}_{kj} + \\beta^k_{,j} \\bar{\\gamma}_{ik}\n",
    "gammabar_rhsDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            gammabar_rhsDD[i][j] += betaU[k]*gammabarDD_dupD[i][j][k] + betaU_dD[k][i]*gammabarDD[k][j] \\\n",
    "                                                                      + betaU_dD[k][j]*gammabarDD[i][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term2_partial_gamma'></a>\n",
    "\n",
    "## Step 2.b: Term 2 of $\\partial_t \\bar{\\gamma}_{i j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term2_partial_gamma}$$\n",
    "\n",
    "Term 2 of $\\partial_t \\bar{\\gamma}_{i j} =$ `gammabar_rhsDD[i][j]`: $\\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right )$\n",
    "\n",
    "Let's first convert this expression to be in terms of the evolved variables $a_{ij}$ and $\\mathcal{B}^i$, starting with $\\bar{A}_{ij} = a_{ij} \\text{ReDD[i][j]}$. Then $\\bar{A}^k_{k} = \\bar{\\gamma}^{ij} \\bar{A}_{ij}$, and we have already defined $\\bar{\\gamma}^{ij}$ in terms of the evolved quantity $h_{ij}$.\n",
    "\n",
    "Next, we wish to compute \n",
    "\n",
    "$$\\bar{D}_{k} \\beta^{k} = \\beta^k_{,k} + \\frac{\\beta^k \\bar{\\gamma}_{,k}}{2 \\bar{\\gamma}},$$\n",
    "\n",
    "where $\\bar{\\gamma}$ is the determinant of the conformal metric $\\bar{\\gamma}_{ij}$. ***Exercise to student: Prove the above relation.***\n",
    "[Solution.](https://physics.stackexchange.com/questions/81453/general-relativity-christoffel-symbol-identity)\n",
    "\n",
    "Usually (i.e., so long as we make the parameter choice `detgbarOverdetghat_equals_one = False` ) we will choose $\\bar{\\gamma}=\\hat{\\gamma}$, so $\\bar{\\gamma}$ will in general possess coordinate singularities. Thus we would prefer to rewrite derivatives of $\\bar{\\gamma}$ in terms of derivatives of $\\bar{\\gamma}/\\hat{\\gamma} = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.b.i: First import \\bar{A}_{ij} = AbarDD[i][j], and its contraction trAbar = \\bar{A}^k_k\n",
    "#           from BSSN.BSSN_quantities\n",
    "Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()\n",
    "trAbar = Bq.trAbar\n",
    "\n",
    "# Step 2.b.ii: Import detgammabar quantities from BSSN.BSSN_quantities:\n",
    "Bq.detgammabar_and_derivs()\n",
    "detgammabar = Bq.detgammabar\n",
    "detgammabar_dD = Bq.detgammabar_dD\n",
    "\n",
    "# Step 2.b.ii: Compute the contraction \\bar{D}_k \\beta^k = \\beta^k_{,k} + \\frac{\\beta^k \\bar{\\gamma}_{,k}}{2 \\bar{\\gamma}}\n",
    "Dbarbetacontraction = sp.sympify(0)\n",
    "for k in range(DIM):\n",
    "    Dbarbetacontraction += betaU_dD[k][k] + betaU[k]*detgammabar_dD[k]/(2*detgammabar)\n",
    "\n",
    "# Step 2.b.iii: Second term of \\partial_t \\bar{\\gamma}_{i j} right-hand side:\n",
    "# \\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right )\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        gammabar_rhsDD[i][j] += sp.Rational(2,3)*gammabarDD[i][j]*(alpha*trAbar - Dbarbetacontraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term3_partial_gamma'></a>\n",
    "\n",
    "## Step 2.c: Term 3 of $\\partial_t \\bar{\\gamma}_{i j}$ \\[Back to [top](#toc)\\] \n",
    "$$\\label{term3_partial_gamma}$$\n",
    "\n",
    "Term 3 of $\\partial_t \\bar{\\gamma}_{i j}$ = `gammabar_rhsDD[i][j]`: $-2 \\alpha \\bar{A}_{ij}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.c: Third term of \\partial_t \\bar{\\gamma}_{i j} right-hand side:\n",
    "# -2 \\alpha \\bar{A}_{ij}\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        gammabar_rhsDD[i][j] += -2*alpha*AbarDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='abar'></a>\n",
    "\n",
    "# Step 3: Right-hand side of $\\partial_t \\bar{A}_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{abar}$$\n",
    "\n",
    "$$\\partial_t \\bar{A}_{i j} = \n",
    "{\\underbrace {\\textstyle \\left[\\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik} \\right]}_{\\text{Term 1}}}\n",
    "{\\underbrace {\\textstyle - \\frac{2}{3} \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} {\\bar{A}^{k}}_{j} + \\alpha \\bar{A}_{i j} K}_{\\text{Term 2}}} + \n",
    "{\\underbrace {\\textstyle e^{-4 \\phi} \\left \\{-2 \\alpha \\bar{D}_{i} \\bar{D}_{j} \\phi + 4 \\alpha \\bar{D}_{i} \\phi \\bar{D}_{j} \\phi  + 4 \\bar{D}_{(i} \\alpha \\bar{D}_{j)} \\phi - \\bar{D}_{i} \\bar{D}_{j} \\alpha + \\alpha \\bar{R}_{i j} \\right \\}^{\\text{TF}}}_{\\text{Term 3}}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term1_partial_upper_a'></a>\n",
    "\n",
    "## Step 3.a: Term 1 of $\\partial_t \\bar{A}_{i j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term1_partial_upper_a}$$\n",
    "\n",
    "Term 1 of $\\partial_t \\bar{A}_{i j}$ = `Abar_rhsDD[i][j]`: $\\left[\\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik} \\right]$\n",
    "\n",
    "\n",
    "Notice the first subexpression has a $\\beta^k \\partial_k A_{ij}$ advection term, which will be upwinded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3.a: First term of \\partial_t \\bar{A}_{i j}:\n",
    "# \\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik}\n",
    "\n",
    "AbarDD_dupD = Bq.AbarDD_dupD # From Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()\n",
    "\n",
    "Abar_rhsDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            Abar_rhsDD[i][j] += betaU[k]*AbarDD_dupD[i][j][k] + betaU_dD[k][i]*AbarDD[k][j] \\\n",
    "                                                              + betaU_dD[k][j]*AbarDD[i][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term2_partial_upper_a'></a>\n",
    "\n",
    "## Step 3.b: Term 2 of $\\partial_t \\bar{A}_{i j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term2_partial_upper_a}$$\n",
    "\n",
    "Term 2 of $\\partial_t \\bar{A}_{i j}$ = `Abar_rhsDD[i][j]`: $- \\frac{2}{3} \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} \\bar{A}^{k}_{j} + \\alpha \\bar{A}_{i j} K$\n",
    "\n",
    "\n",
    "Note that $\\bar{D}_{k} \\beta^{k}$ was already defined as `Dbarbetacontraction`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3.b: Second term of \\partial_t \\bar{A}_{i j}:\n",
    "# - (2/3) \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} {\\bar{A}^{k}}_{j} + \\alpha \\bar{A}_{i j} K\n",
    "gammabarUU = Bq.gammabarUU # From Bq.gammabar__inverse_and_derivs()\n",
    "AbarUD     = Bq.AbarUD     # From Bq.AbarUU_AbarUD_trAbar()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Abar_rhsDD[i][j] += -sp.Rational(2,3)*AbarDD[i][j]*Dbarbetacontraction + alpha*AbarDD[i][j]*trK\n",
    "        for k in range(DIM):\n",
    "            Abar_rhsDD[i][j] += -2*alpha * AbarDD[i][k]*AbarUD[k][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term3_partial_upper_a'></a>\n",
    "\n",
    "## Step 3.c: Term 3 of $\\partial_t \\bar{A}_{i j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term3_partial_upper_a}$$\n",
    "\n",
    "\n",
    "Term 3 of $\\partial_t \\bar{A}_{i j}$ = `Abar_rhsDD[i][j]`: $e^{-4 \\phi} \\left \\{-2 \\alpha \\bar{D}_{i} \\bar{D}_{j} \\phi + 4 \\alpha \\bar{D}_{i} \\phi \\bar{D}_{j} \\phi  + 4 \\bar{D}_{(i} \\alpha \\bar{D}_{j)} \\phi - \\bar{D}_{i} \\bar{D}_{j} \\alpha + \\alpha \\bar{R}_{i j} \\right \\}^{\\text{TF}}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first covariant derivatives of $\\phi$ and $\\alpha$ are simply partial derivatives. However, $\\phi$ is not a gridfunction; `cf` is. cf = $W$ (default value) denotes that the evolved variable is $W=e^{-2 \\phi}$, which results in smoother spacetime fields around puncture black holes (desirable)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3.c.i: Define partial derivatives of \\phi in terms of evolved quantity \"cf\":\n",
    "Bq.phi_and_derivs()\n",
    "phi_dD   = Bq.phi_dD\n",
    "phi_dupD = Bq.phi_dupD\n",
    "exp_m4phi  = Bq.exp_m4phi\n",
    "phi_dBarD  = Bq.phi_dBarD  # phi_dBarD = Dbar_i phi = phi_dD (since phi is a scalar)\n",
    "phi_dBarDD = Bq.phi_dBarDD # phi_dBarDD = Dbar_i Dbar_j phi (covariant derivative)\n",
    "\n",
    "# Step 3.c.ii: Define RbarDD\n",
    "Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "RbarDD = Bq.RbarDD\n",
    "\n",
    "# Step 3.c.iii: Define first and second derivatives of \\alpha, as well as\n",
    "#         \\bar{D}_i \\bar{D}_j \\alpha, which is defined just like phi\n",
    "alpha_dD = ixp.declarerank1(\"alpha_dD\")\n",
    "alpha_dDD = ixp.declarerank2(\"alpha_dDD\",\"sym01\")\n",
    "alpha_dBarD = alpha_dD\n",
    "alpha_dBarDD = ixp.zerorank2()\n",
    "GammabarUDD = Bq.GammabarUDD # Defined in Bq.gammabar__inverse_and_derivs()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        alpha_dBarDD[i][j] = alpha_dDD[i][j]\n",
    "        for k in range(DIM):\n",
    "            alpha_dBarDD[i][j] += - GammabarUDD[k][i][j]*alpha_dD[k]\n",
    "\n",
    "# Step 3.c.iv: Define the terms in curly braces:\n",
    "curlybrackettermsDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        curlybrackettermsDD[i][j] = -2*alpha*phi_dBarDD[i][j] + 4*alpha*phi_dBarD[i]*phi_dBarD[j] \\\n",
    "                                    +2*alpha_dBarD[i]*phi_dBarD[j] \\\n",
    "                                    +2*alpha_dBarD[j]*phi_dBarD[i] \\\n",
    "                                    -alpha_dBarDD[i][j] + alpha*RbarDD[i][j]\n",
    "\n",
    "# Step 3.c.v: Compute the trace:\n",
    "curlybracketterms_trace = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        curlybracketterms_trace += gammabarUU[i][j]*curlybrackettermsDD[i][j]\n",
    "\n",
    "# Step 3.c.vi: Third and final term of Abar_rhsDD[i][j]:\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Abar_rhsDD[i][j] += exp_m4phi*(curlybrackettermsDD[i][j] - \\\n",
    "                                       sp.Rational(1,3)*gammabarDD[i][j]*curlybracketterms_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cf'></a>\n",
    "\n",
    "# Step 4: Right-hand side of $\\partial_t \\phi \\to \\partial_t (\\text{cf})$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{cf}$$\n",
    "\n",
    "$$\\partial_t \\phi = \n",
    "{\\underbrace {\\textstyle \\left[\\beta^k \\partial_k \\phi \\right]}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\frac{1}{6} \\left (\\bar{D}_{k} \\beta^{k} - \\alpha K \\right)}_{\\text{Term 2}}}$$\n",
    "\n",
    "The right-hand side of $\\partial_t \\phi$ is trivial except for the fact that the actual evolved variable is `cf` (short for conformal factor), which could represent\n",
    "* cf = $\\phi$\n",
    "* cf = $W = e^{-2 \\phi}$ (default)\n",
    "* cf = $\\chi = e^{-4 \\phi}$\n",
    "\n",
    "Thus we are actually computing the right-hand side of the equation $\\partial_t $cf, which is related to $\\partial_t \\phi$ via simple relations:\n",
    "* cf = $\\phi$: $\\partial_t $cf$ = \\partial_t \\phi$ (unchanged)\n",
    "* cf = $W$: $\\partial_t $cf$ = \\partial_t (e^{-2 \\phi}) = -2 e^{-2\\phi}\\partial_t \\phi = -2 W \\partial_t \\phi$. Thus we need to multiply the right-hand side by $-2 W = -2$cf when cf = $W$.\n",
    "* cf = $\\chi$: Same argument as for $W$, except the right-hand side must be multiplied by $-4 \\chi=-4$cf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4: Right-hand side of conformal factor variable \"cf\". Supported\n",
    "#          options include: cf=phi, cf=W=e^(-2*phi) (default), and cf=chi=e^(-4*phi)\n",
    "# \\partial_t phi = \\left[\\beta^k \\partial_k \\phi \\right] <- TERM 1\n",
    "#                  + \\frac{1}{6} \\left (\\bar{D}_{k} \\beta^{k} - \\alpha K \\right ) <- TERM 2\n",
    "\n",
    "cf_rhs = sp.Rational(1,6) * (Dbarbetacontraction - alpha*trK) # Term 2\n",
    "for k in range(DIM):\n",
    "    cf_rhs += betaU[k]*phi_dupD[k] # Term 1\n",
    "\n",
    "# Next multiply to convert phi_rhs to cf_rhs.\n",
    "if par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\") == \"phi\":\n",
    "    pass # do nothing; cf_rhs = phi_rhs\n",
    "elif par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\") == \"W\":\n",
    "    cf_rhs *= -2*cf # cf_rhs = -2*cf*phi_rhs\n",
    "elif par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\") == \"chi\":\n",
    "    cf_rhs *= -4*cf # cf_rhs = -4*cf*phi_rhs\n",
    "else:\n",
    "    print(\"Error: EvolvedConformalFactor_cf == \"+\n",
    "          par.parval_from_str(\"BSSN.BSSN_quantities::EvolvedConformalFactor_cf\")+\" unsupported!\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='trk'></a>\n",
    "\n",
    "# Step 5: Right-hand side of $\\partial_t K$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{trk}$$\n",
    "\n",
    "$$\n",
    "\\partial_{t} K = \n",
    "{\\underbrace {\\textstyle \\left[\\beta^i \\partial_i K \\right]}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\frac{1}{3} \\alpha K^{2}}_{\\text{Term 2}}} +\n",
    "{\\underbrace {\\textstyle \\alpha \\bar{A}_{i j} \\bar{A}^{i j}}_{\\text{Term 3}}}\n",
    "{\\underbrace {\\textstyle - e^{-4 \\phi} \\left (\\bar{D}_{i} \\bar{D}^{i} \\alpha + 2 \\bar{D}^{i} \\alpha \\bar{D}_{i} \\phi \\right )}_{\\text{Term 4}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5: right-hand side of trK (trace of extrinsic curvature):\n",
    "# \\partial_t K = \\beta^k \\partial_k K <- TERM 1\n",
    "#           + \\frac{1}{3} \\alpha K^{2} <- TERM 2\n",
    "#           + \\alpha \\bar{A}_{i j} \\bar{A}^{i j} <- TERM 3\n",
    "#           - - e^{-4 \\phi} (\\bar{D}_{i} \\bar{D}^{i} \\alpha + 2 \\bar{D}^{i} \\alpha \\bar{D}_{i} \\phi ) <- TERM 4\n",
    "# TERM 2:\n",
    "trK_rhs = sp.Rational(1,3)*alpha*trK*trK\n",
    "trK_dupD = ixp.declarerank1(\"trK_dupD\")\n",
    "for i in range(DIM):\n",
    "    # TERM 1:\n",
    "    trK_rhs += betaU[i]*trK_dupD[i]\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # TERM 4:\n",
    "        trK_rhs += -exp_m4phi*gammabarUU[i][j]*(alpha_dBarDD[i][j] + 2*alpha_dBarD[j]*phi_dBarD[i])\n",
    "AbarUU = Bq.AbarUU # From Bq.AbarUU_AbarUD_trAbar()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # TERM 3:\n",
    "        trK_rhs += alpha*AbarDD[i][j]*AbarUU[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lambdabar'></a>\n",
    "\n",
    "# Step 6: Right-hand side of $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{lambdabar}$$\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t \\bar{\\Lambda}^{i} &= \n",
    "{\\underbrace {\\textstyle \\left[\\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k \\right]}_{\\text{Term 1}}} + \n",
    "{\\underbrace {\\textstyle \\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i}}_{\\text{Term 2}}} + \n",
    "{\\underbrace {\\textstyle \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j}}_{\\text{Term 3}}} + \n",
    "{\\underbrace {\\textstyle \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j}}_{\\text{Term 4}}} \\nonumber \\\\\n",
    "  & \n",
    "{\\underbrace {\\textstyle - 2 \\bar{A}^{i j} \\left (\\partial_{j} \\alpha - 6 \\alpha \\partial_{j} \\phi \\right )}_{\\text{Term 5}}} + \n",
    "{\\underbrace {\\textstyle 2 \\alpha \\bar{A}^{j k} \\Delta_{j k}^{i}}_{\\text{Term 6}}} \n",
    "{\\underbrace {\\textstyle -\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K}_{\\text{Term 7}}}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term1_partial_lambda'></a>\n",
    "\n",
    "## Step 6.a: Term 1 of  $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term1_partial_lambda}$$\n",
    "\n",
    "Term 1 of  $\\partial_t \\bar{\\Lambda}^{i}$: $\\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k$\n",
    "\n",
    "Computing this term requires that we define $\\bar{\\Lambda}^i$ and $\\bar{\\Lambda}^i_{,j}$ in terms of the rescaled (i.e., actual evolved) variable $\\lambda^i$ and derivatives: \n",
    "\\begin{align}\n",
    "\\bar{\\Lambda}^i &= \\lambda^i \\text{ReU[i]} \\\\\n",
    "\\bar{\\Lambda}^i_{,\\ j} &= \\lambda^i_{,\\ j} \\text{ReU[i]} + \\lambda^i \\text{ReUdD[i][j]} \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6: right-hand side of \\partial_t \\bar{\\Lambda}^i:\n",
    "# \\partial_t \\bar{\\Lambda}^i = \\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k <- TERM 1\n",
    "#                            + \\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i} <- TERM 2\n",
    "#                            + \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j} <- TERM 3\n",
    "#                            + \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j} <- TERM 4\n",
    "#                            - 2 \\bar{A}^{i j} (\\partial_{j} \\alpha - 6 \\partial_{j} \\phi) <- TERM 5\n",
    "#                            + 2 \\alpha \\bar{A}^{j k} \\Delta_{j k}^{i} <- TERM 6\n",
    "#                            - \\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K <- TERM 7\n",
    "\n",
    "# Step 6.a: Term 1 of \\partial_t \\bar{\\Lambda}^i: \\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k\n",
    "# First we declare \\bar{\\Lambda}^i and \\bar{\\Lambda}^i_{,j} in terms of \\lambda^i and \\lambda^i_{,j}\n",
    "LambdabarU_dupD = ixp.zerorank2()\n",
    "lambdaU_dupD = ixp.declarerank2(\"lambdaU_dupD\",\"nosym\")\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        LambdabarU_dupD[i][j] = lambdaU_dupD[i][j]*rfm.ReU[i] + lambdaU[i]*rfm.ReUdD[i][j]\n",
    "\n",
    "Lambdabar_rhsU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        Lambdabar_rhsU[i] += betaU[k]*LambdabarU_dupD[i][k] - betaU_dD[i][k]*LambdabarU[k] # Term 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term2_partial_lambda'></a>\n",
    "\n",
    "## Step 6.b: Term 2 of $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term2_partial_lambda}$$\n",
    "\n",
    "Term 2 of $\\partial_t \\bar{\\Lambda}^{i}$: $\\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i}$\n",
    "\n",
    "This is a relatively difficult term to compute, as it requires we evaluate the second covariant derivative of the shift vector, with respect to the hatted (i.e., reference) metric.\n",
    "\n",
    "Based on the definition of covariant derivative, we have\n",
    "$$\n",
    "\\hat{D}_{k} \\beta^{i} = \\beta^i_{,k} + \\hat{\\Gamma}^i_{mk} \\beta^m\n",
    "$$\n",
    "\n",
    "Since $\\hat{D}_{k} \\beta^{i}$ is a tensor, the covariant derivative of this will have the same indexing as a tensor $T_k^i$:\n",
    "\n",
    "$$\n",
    "\\hat{D}_{j} T^i_k = T^i_{k,j} + \\hat{\\Gamma}^i_{dj} T^d_k - \\hat{\\Gamma}^d_{kj} T^i_d.\n",
    "$$\n",
    "\n",
    "Therefore,\n",
    "\\begin{align}\n",
    "\\hat{D}_{j} \\left(\\hat{D}_{k} \\beta^{i}\\right) &= \\left(\\beta^i_{,k} + \\hat{\\Gamma}^i_{mk} \\beta^m\\right)_{,j} + \\hat{\\Gamma}^i_{dj} \\left(\\beta^d_{,k} + \\hat{\\Gamma}^d_{mk} \\beta^m\\right) - \\hat{\\Gamma}^d_{kj} \\left(\\beta^i_{,d} + \\hat{\\Gamma}^i_{md} \\beta^m\\right) \\\\\n",
    "&= \\beta^i_{,kj} + \\hat{\\Gamma}^i_{mk,j} \\beta^m + \\hat{\\Gamma}^i_{mk} \\beta^m_{,j} + \\hat{\\Gamma}^i_{dj}\\beta^d_{,k} + \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} \\beta^m - \\hat{\\Gamma}^d_{kj} \\beta^i_{,d} - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} \\beta^m \\\\\n",
    "&= {\\underbrace {\\textstyle \\beta^i_{,kj}}_{\\text{Term 2a}}}\n",
    "{\\underbrace {\\textstyle \\hat{\\Gamma}^i_{mk,j} \\beta^m + \\hat{\\Gamma}^i_{mk} \\beta^m_{,j} + \\hat{\\Gamma}^i_{dj}\\beta^d_{,k} - \\hat{\\Gamma}^d_{kj} \\beta^i_{,d}}_{\\text{Term 2b}}} +\n",
    "{\\underbrace {\\textstyle \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} \\beta^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} \\beta^m}_{\\text{Term 2c}}},\n",
    "\\end{align}\n",
    "\n",
    "where \n",
    "$$\n",
    "\\text{Term 2} = \\bar{\\gamma}^{jk} \\left(\\text{Term 2a} + \\text{Term 2b} + \\text{Term 2c}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.b: Term 2 of \\partial_t \\bar{\\Lambda}^i = \\bar{\\gamma}^{jk} (Term 2a + Term 2b + Term 2c)\n",
    "# Term 2a: \\bar{\\gamma}^{jk} \\beta^i_{,kj}\n",
    "Term2aUDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            Term2aUDD[i][j][k] += betaU_dDD[i][k][j]\n",
    "# Term 2b: \\hat{\\Gamma}^i_{mk,j} \\beta^m + \\hat{\\Gamma}^i_{mk} \\beta^m_{,j}\n",
    "#          + \\hat{\\Gamma}^i_{dj}\\beta^d_{,k} - \\hat{\\Gamma}^d_{kj} \\beta^i_{,d}\n",
    "Term2bUDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                Term2bUDD[i][j][k] += rfm.GammahatUDDdD[i][m][k][j]*betaU[m]    \\\n",
    "                                      + rfm.GammahatUDD[i][m][k]*betaU_dD[m][j] \\\n",
    "                                      + rfm.GammahatUDD[i][m][j]*betaU_dD[m][k] \\\n",
    "                                      - rfm.GammahatUDD[m][k][j]*betaU_dD[i][m]\n",
    "# Term 2c: \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} \\beta^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} \\beta^m\n",
    "Term2cUDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                for d in range(DIM):\n",
    "                    Term2cUDD[i][j][k] += ( rfm.GammahatUDD[i][d][j]*rfm.GammahatUDD[d][m][k] \\\n",
    "                                           -rfm.GammahatUDD[d][k][j]*rfm.GammahatUDD[i][m][d])*betaU[m]\n",
    "\n",
    "Lambdabar_rhsUpieceU = ixp.zerorank1()\n",
    "\n",
    "# Put it all together to get Term 2:\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            Lambdabar_rhsU[i] += gammabarUU[j][k] * (Term2aUDD[i][j][k] + Term2bUDD[i][j][k] + Term2cUDD[i][j][k])\n",
    "            Lambdabar_rhsUpieceU[i] += gammabarUU[j][k] * (Term2aUDD[i][j][k] + Term2bUDD[i][j][k] + Term2cUDD[i][j][k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term3_partial_lambda'></a>\n",
    "\n",
    "## Step 6.c: Term 3 of  $\\partial_t \\bar{\\Lambda}^{i}$: $\\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term3_partial_lambda}$$\n",
    "\n",
    "Term 3 of $\\partial_t \\bar{\\Lambda}^{i}$: $\\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j}$\n",
    " \n",
    "This term is the simplest to implement, as $\\bar{D}_{j} \\beta^{j}$ and $\\Delta^i$ have already been defined, as `Dbarbetacontraction` and `DGammaU[i]`, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.c: Term 3 of \\partial_t \\bar{\\Lambda}^i:\n",
    "#    \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j}\n",
    "DGammaU = Bq.DGammaU # From Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "for i in range(DIM):\n",
    "    Lambdabar_rhsU[i] += sp.Rational(2,3)*DGammaU[i]*Dbarbetacontraction # Term 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term4_partial_lambda'></a>\n",
    "\n",
    "## Step 6.d: Term 4 of  $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term4_partial_lambda}$$\n",
    "\n",
    "$\\partial_t \\bar{\\Lambda}^{i}$: $\\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j}$\n",
    "\n",
    "Recall first that\n",
    "\n",
    "$$\\bar{D}_{k} \\beta^{k} = \\beta^k_{,\\ k} + \\frac{\\beta^k \\bar{\\gamma}_{,k}}{2 \\bar{\\gamma}},$$\n",
    "which is a scalar, so\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{D}_m \\bar{D}_{j} \\beta^{j} &= \\left(\\beta^k_{,\\ k} + \\frac{\\beta^k \\bar{\\gamma}_{,k}}{2 \\bar{\\gamma}}\\right)_{,m} \\\\\n",
    "&= \\beta^k_{\\ ,km} + \\frac{\\beta^k_{\\ ,m} \\bar{\\gamma}_{,k} + \\beta^k \\bar{\\gamma}_{\\ ,km}}{2 \\bar{\\gamma}} - \\frac{\\beta^k \\bar{\\gamma}_{,k} \\bar{\\gamma}_{,m}}{2 \\bar{\\gamma}^2}\n",
    "\\end{align}\n",
    "\n",
    "Thus, \n",
    "\\begin{align}\n",
    "\\bar{D}^i \\bar{D}_{j} \\beta^{j} \n",
    "&= \\bar{\\gamma}^{im} \\bar{D}_m \\bar{D}_{j} \\beta^{j} \\\\\n",
    "&= \\bar{\\gamma}^{im} \\left(\\beta^k_{\\ ,km} + \\frac{\\beta^k_{\\ ,m} \\bar{\\gamma}_{,k} + \\beta^k \\bar{\\gamma}_{\\ ,km}}{2 \\bar{\\gamma}} - \\frac{\\beta^k \\bar{\\gamma}_{,k} \\bar{\\gamma}_{,m}}{2 \\bar{\\gamma}^2} \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.d: Term 4 of \\partial_t \\bar{\\Lambda}^i:\n",
    "#           \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j}\n",
    "detgammabar_dDD = Bq.detgammabar_dDD # From Bq.detgammabar_and_derivs()\n",
    "Dbarbetacontraction_dBarD = ixp.zerorank1()\n",
    "for k in range(DIM):\n",
    "    for m in range(DIM):\n",
    "        Dbarbetacontraction_dBarD[m] += betaU_dDD[k][k][m] + \\\n",
    "                                        (betaU_dD[k][m]*detgammabar_dD[k] +\n",
    "                                         betaU[k]*detgammabar_dDD[k][m])/(2*detgammabar) \\\n",
    "                                        -betaU[k]*detgammabar_dD[k]*detgammabar_dD[m]/(2*detgammabar*detgammabar)\n",
    "for i in range(DIM):\n",
    "    for m in range(DIM):\n",
    "        Lambdabar_rhsU[i] += sp.Rational(1,3)*gammabarUU[i][m]*Dbarbetacontraction_dBarD[m]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term5_partial_lambda'></a>\n",
    "\n",
    "## Step 6.e: Term 5 of  $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term5_partial_lambda}$$\n",
    "\n",
    "Term 5 of  $\\partial_t \\bar{\\Lambda}^{i}$: $- 2 \\bar{A}^{i j} \\left (\\partial_{j} \\alpha - 6\\alpha \\partial_{j} \\phi\\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.e: Term 5 of \\partial_t \\bar{\\Lambda}^i:\n",
    "#           - 2 \\bar{A}^{i j} (\\partial_{j} \\alpha - 6 \\alpha \\partial_{j} \\phi)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Lambdabar_rhsU[i] += -2*AbarUU[i][j]*(alpha_dD[j] - 6*alpha*phi_dD[j])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term6_partial_lambda'></a>\n",
    "\n",
    "## Step 6.f: Term 6 of $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term6_partial_lambda}$$\n",
    "\n",
    "Term 6 of $\\partial_t \\bar{\\Lambda}^{i}$: $2\\alpha \\bar{A}^{j k} \\Delta_{j k}^{i}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.f: Term 6 of \\partial_t \\bar{\\Lambda}^i:\n",
    "#           2 \\alpha \\bar{A}^{j k} \\Delta^{i}_{j k}\n",
    "DGammaUDD = Bq.DGammaUDD # From RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            Lambdabar_rhsU[i] += 2*alpha*AbarUU[j][k]*DGammaUDD[i][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='term7_partial_lambda'></a>\n",
    "\n",
    "## Step 6.g: Term 7 of $\\partial_t \\bar{\\Lambda}^{i}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{term7_partial_lambda}$$\n",
    "\n",
    "$\\partial_t \\bar{\\Lambda}^{i}$: $-\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6.g: Term 7 of \\partial_t \\bar{\\Lambda}^i:\n",
    "#           -\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K\n",
    "trK_dD = ixp.declarerank1(\"trK_dD\")\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Lambdabar_rhsU[i] += -sp.Rational(4,3)*alpha*gammabarUU[i][j]*trK_dD[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rescalingrhss'></a>\n",
    "\n",
    "# Step 7: Rescaling the BSSN right-hand sides; rewriting them in terms of the rescaled quantities $\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{rescalingrhss}$$\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}\\right\\}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Step 7: Rescale the RHS quantities so that the evolved\n",
    "#          variables are smooth across coord singularities\n",
    "h_rhsDD     = ixp.zerorank2()\n",
    "a_rhsDD     = ixp.zerorank2()\n",
    "lambda_rhsU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    lambda_rhsU[i] = Lambdabar_rhsU[i] / rfm.ReU[i]\n",
    "    for j in range(DIM):\n",
    "        h_rhsDD[i][j] = gammabar_rhsDD[i][j] / rfm.ReDD[i][j]\n",
    "        a_rhsDD[i][j] =     Abar_rhsDD[i][j] / rfm.ReDD[i][j]\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 8: Code Validation against `BSSN.BSSN_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 equations between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.BSSN_RHSs](../edit/BSSN/BSSN_RHSs.py) module.\n",
    "\n",
    "By default, we analyze the RHSs in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between BSSN_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
      "trK_rhs - bssnrhs.trK_rhs = 0\n",
      "cf_rhs - bssnrhs.cf_rhs = 0\n",
      "lambda_rhsU[0] - bssnrhs.lambda_rhsU[0] = 0\n",
      "h_rhsDD[0][0] - bssnrhs.h_rhsDD[0][0] = 0\n",
      "a_rhsDD[0][0] - bssnrhs.a_rhsDD[0][0] = 0\n",
      "h_rhsDD[0][1] - bssnrhs.h_rhsDD[0][1] = 0\n",
      "a_rhsDD[0][1] - bssnrhs.a_rhsDD[0][1] = 0\n",
      "h_rhsDD[0][2] - bssnrhs.h_rhsDD[0][2] = 0\n",
      "a_rhsDD[0][2] - bssnrhs.a_rhsDD[0][2] = 0\n",
      "lambda_rhsU[1] - bssnrhs.lambda_rhsU[1] = 0\n",
      "h_rhsDD[1][0] - bssnrhs.h_rhsDD[1][0] = 0\n",
      "a_rhsDD[1][0] - bssnrhs.a_rhsDD[1][0] = 0\n",
      "h_rhsDD[1][1] - bssnrhs.h_rhsDD[1][1] = 0\n",
      "a_rhsDD[1][1] - bssnrhs.a_rhsDD[1][1] = 0\n",
      "h_rhsDD[1][2] - bssnrhs.h_rhsDD[1][2] = 0\n",
      "a_rhsDD[1][2] - bssnrhs.a_rhsDD[1][2] = 0\n",
      "lambda_rhsU[2] - bssnrhs.lambda_rhsU[2] = 0\n",
      "h_rhsDD[2][0] - bssnrhs.h_rhsDD[2][0] = 0\n",
      "a_rhsDD[2][0] - bssnrhs.a_rhsDD[2][0] = 0\n",
      "h_rhsDD[2][1] - bssnrhs.h_rhsDD[2][1] = 0\n",
      "a_rhsDD[2][1] - bssnrhs.a_rhsDD[2][1] = 0\n",
      "h_rhsDD[2][2] - bssnrhs.h_rhsDD[2][2] = 0\n",
      "a_rhsDD[2][2] - bssnrhs.a_rhsDD[2][2] = 0\n"
     ]
    }
   ],
   "source": [
    "# Step 8: We already have SymPy expressions for BSSN 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 RHS expressions\n",
    "#         to validate against the same expressions in the\n",
    "#         BSSN/BSSN_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 9.a: 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_RHSs as bssnrhs\n",
    "bssnrhs.BSSN_RHSs()\n",
    "\n",
    "print(\"Consistency check between BSSN_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "print(\"trK_rhs - bssnrhs.trK_rhs = \" + str(trK_rhs - bssnrhs.trK_rhs))\n",
    "print(\"cf_rhs - bssnrhs.cf_rhs = \" + str(cf_rhs - bssnrhs.cf_rhs))\n",
    "\n",
    "for i in range(DIM):\n",
    "    print(\"lambda_rhsU[\"+str(i)+\"] - bssnrhs.lambda_rhsU[\"+str(i)+\"] = \" +\n",
    "          str(lambda_rhsU[i] - bssnrhs.lambda_rhsU[i]))\n",
    "    for j in range(DIM):\n",
    "        print(\"h_rhsDD[\"+str(i)+\"][\"+str(j)+\"] - bssnrhs.h_rhsDD[\"+str(i)+\"][\"+str(j)+\"] = \"\n",
    "              + str(h_rhsDD[i][j] - bssnrhs.h_rhsDD[i][j]))\n",
    "        print(\"a_rhsDD[\"+str(i)+\"][\"+str(j)+\"] - bssnrhs.a_rhsDD[\"+str(i)+\"][\"+str(j)+\"] = \"\n",
    "              + str(a_rhsDD[i][j] - bssnrhs.a_rhsDD[i][j]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 9: 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_RHSs.pdf](Tutorial-BSSN_time_evolution-BSSN_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": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_time_evolution-BSSN_RHSs.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-BSSN_time_evolution-BSSN_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_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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}