Tutorial-BSSN_quantities.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 Quantities\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module documents and constructs a number of quantities useful for building 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) (see also [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632)). \n",
    "\n",
    "**Notebook Status:** <font color='orange'><b> Self-Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n",
    "\n",
    "[comment]: <> (Introduction: TODO)\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",
    "Each family of quantities is constructed within a given function (**boldfaced** below). This notebook is organized as follows\n",
    "\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#declare_bssn_gfs): **`declare_BSSN_gridfunctions_if_not_declared_already()`**: Declare all of the core BSSN variables $\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$ and register them as gridfunctions\n",
    "1. [Step 3](#rescaling_tensors) Rescaling tensors to avoid coordinate singularities\n",
    "    1. [Step 3.a](#bssn_basic_tensors) **`BSSN_basic_tensors()`**: Define all basic conformal BSSN tensors $\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\bar{\\Lambda}^{i}, \\beta^i, B^i\\right\\}$ in terms of BSSN gridfunctions\n",
    "1. [Step 4](#bssn_barred_metric__inverse_and_derivs): **`gammabar__inverse_and_derivs()`**: $\\bar{\\gamma}^{ij}$, and spatial derivatives of $\\bar{\\gamma}_{ij}$ including $\\bar{\\Gamma}^{i}_{jk}$\n",
    "    1. [Step 4.a](#bssn_barred_metric__inverse): Inverse conformal 3-metric: $\\bar{\\gamma^{ij}}$\n",
    "    1. [Step 4.b](#bssn_barred_metric__derivs): Derivatives of the conformal 3-metric $\\bar{\\gamma}_{ij,k}$ and $\\bar{\\gamma}_{ij,kl}$, and associated \"barred\" Christoffel symbols $\\bar{\\Gamma}^{i}_{jk}$\n",
    "1. [Step 5](#detgammabar_and_derivs): **`detgammabar_and_derivs()`**: $\\det \\bar{\\gamma}_{ij}$ and its derivatives\n",
    "1. [Step 6](#abar_quantities): **`AbarUU_AbarUD_trAbar()`**: Quantities related to conformal traceless extrinsic curvature $\\bar{A}_{ij}$: $\\bar{A}^{ij}$, $\\bar{A}^i_j$, and $\\bar{A}^k_k$\n",
    "1. [Step 7](#rbar): **`RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()`**: The conformal (\"barred\") Ricci tensor $\\bar{R}_{ij}$ and associated quantities\n",
    "    1. [Step 7.a](#rbar_part1): Conformal Ricci tensor, part 1: The $\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j}$ term\n",
    "    1. [Step 7.b](#rbar_part2): Conformal Ricci tensor, part 2: The $\\bar{\\gamma}_{k(i} \\hat{D}_{j)} \\bar{\\Lambda}^{k}$ term\n",
    "    1. [Step 7.c](#rbar_part3): Conformal Ricci tensor, part 3: The $\\Delta^{k} \\Delta_{(i j) k}  + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right )$ terms\n",
    "    1. [Step 7.d](#summing_rbar_terms): Summing the terms and defining $\\bar{R}_{ij}$\n",
    "1. [Step 8](#beta_derivs): **`betaU_derivs()`**: Unrescaled shift vector $\\beta^i$ and spatial derivatives $\\beta^i_{,j}$ and $\\beta^i_{,jk}$\n",
    "1. [Step 9](#phi_and_derivs): **`phi_and_derivs()`**: Standard BSSN conformal factor $\\phi$, and its derivatives $\\phi_{,i}$, $\\phi_{,ij}$, $\\bar{D}_j \\phi_i$, and $\\bar{D}_j\\bar{D}_k \\phi_i$\n",
    "    1. [Step 9.a](#phi_ito_cf): $\\phi$ in terms of the chosen (possibly non-standard) conformal factor variable `cf` (e.g., `cf`$=W=e^{-4\\phi}$)\n",
    "    1. [Step 9.b](#phi_covariant_derivs): Partial and covariant derivatives of $\\phi$\n",
    "1. [Step 10](#code_validation): Code Validation against `BSSN.BSSN_quantities` NRPy+ module\n",
    "1. [Step 11](#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",
    "\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Import all needed modules from NRPy+:\n",
    "import NRPy_param_funcs as par\n",
    "import sympy as sp\n",
    "import indexedexp as ixp\n",
    "import grid as gri\n",
    "import reference_metric as rfm\n",
    "import sys\n",
    "\n",
    "# Step 1.a: Set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\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()\n",
    "\n",
    "# Step 1.c: 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",
    "par.set_parval_from_str(\"grid::DIM\",DIM)\n",
    "\n",
    "# Step 1.d: Declare/initialize parameters for this module\n",
    "thismodule = \"BSSN_quantities\"\n",
    "par.initialize_param(par.glb_param(\"char\", thismodule, \"EvolvedConformalFactor_cf\", \"W\"))\n",
    "par.initialize_param(par.glb_param(\"bool\", thismodule, \"detgbarOverdetghat_equals_one\", \"True\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='declare_bssn_gfs'></a>\n",
    "\n",
    "# Step 2: `declare_BSSN_gridfunctions_if_not_declared_already()`: Declare all of the core BSSN variables $\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$ and register them as gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{declare_bssn_gfs}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Register all needed BSSN gridfunctions.\n",
    "# Step 2.a: Register indexed quantities, using ixp.register_... functions\n",
    "hDD = ixp.register_gridfunctions_for_single_rank2(\"EVOL\", \"hDD\", \"sym01\")\n",
    "aDD = ixp.register_gridfunctions_for_single_rank2(\"EVOL\", \"aDD\", \"sym01\")\n",
    "lambdaU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"lambdaU\")\n",
    "vetU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"vetU\")\n",
    "betU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"betU\")\n",
    "# Step 2.b: Register scalar quantities, using gri.register_gridfunctions()\n",
    "trK, cf, alpha = gri.register_gridfunctions(\"EVOL\",[\"trK\", \"cf\", \"alpha\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rescaling_tensors'></a>\n",
    "\n",
    "# Step 3: Rescaling tensors to avoid coordinate singularities \\[Back to [top](#toc)\\]\n",
    "$$\\label{rescaling_tensors}$$\n",
    "\n",
    "While the [covariant form of the BSSN evolution equations](Tutorial-BSSNCurvilinear.ipynb) are properly covariant (with the potential exception of the shift evolution equation, since the shift is a [freely specifiable gauge quantity](https://en.wikipedia.org/wiki/Gauge_fixing)), components of the rank-1 and rank-2 tensors $\\varepsilon_{i j}$, $\\bar{A}_{i j}$, and $\\bar{\\Lambda}^{i}$ will drop to zero (destroying information) or diverge (to $\\infty$) at coordinate singularities. \n",
    "\n",
    "The good news is, this singular behavior is well-understood in terms of the scale factors of the reference metric, enabling us to define rescaled version of these quantities that are well behaved (so that, e.g., they can be finite differenced).\n",
    "\n",
    "For example, given a smooth vector *in a 3D Cartesian basis* $\\bar{\\Lambda}^{i}$, all components $\\bar{\\Lambda}^{x}$, $\\bar{\\Lambda}^{y}$, and $\\bar{\\Lambda}^{z}$ will be smooth (by assumption). When changing the basis to spherical coordinates (applying the appropriate Jacobian matrix transformation), we will find that since $\\phi = \\arctan(y/x)$, $\\bar{\\Lambda}^{\\phi}$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{\\Lambda}^{\\phi} &= \\frac{\\partial \\phi}{\\partial x} \\bar{\\Lambda}^{x} + \n",
    "\\frac{\\partial \\phi}{\\partial y} \\bar{\\Lambda}^{y} + \n",
    "\\frac{\\partial \\phi}{\\partial z} \\bar{\\Lambda}^{z} \\\\\n",
    "&= -\\frac{y}{\\sqrt{x^2+y^2}} \\bar{\\Lambda}^{x} + \n",
    "\\frac{x}{\\sqrt{x^2+y^2}} \\bar{\\Lambda}^{y} \\\\\n",
    "&= -\\frac{y}{r \\sin\\theta} \\bar{\\Lambda}^{x} + \n",
    "\\frac{x}{r \\sin\\theta} \\bar{\\Lambda}^{y}.\n",
    "\\end{align}\n",
    "\n",
    "Thus $\\bar{\\Lambda}^{\\phi}$ diverges at all points where $r\\sin\\theta=0$ due to the $\\frac{1}{r\\sin\\theta}$ that appear in the Jacobian transformation. \n",
    "\n",
    "This divergence might pose no problem on cell-centered grids that avoid $r \\sin\\theta=0$, except that the BSSN equations require that *first and second derivatives* of these quantities be taken. Usual strategies for numerical approximation of these derivatives (e.g., finite difference methods) will \"see\" these divergences and errors generally will not drop to zero with increased numerical sampling of the functions at points near where the functions diverge.\n",
    "\n",
    "However, notice that if we define $\\lambda^{\\phi}$ such that\n",
    "\n",
    "$$\\bar{\\Lambda}^{\\phi} = \\frac{1}{r\\sin\\theta} \\lambda^{\\phi},$$\n",
    "\n",
    "then $\\lambda^{\\phi}$ will be smooth as well. \n",
    "\n",
    "Avoiding such singularities can be generalized to other coordinate systems, so long as $\\lambda^i$ is defined as:\n",
    "\n",
    "$$\\bar{\\Lambda}^{i} = \\frac{\\lambda^i}{\\text{scalefactor[i]}} ,$$\n",
    "\n",
    "where scalefactor\\[i\\] is the $i$th scale factor in the given coordinate system. In an identical fashion, we define the smooth versions of $\\beta^i$ and $B^i$ to be $\\mathcal{V}^i$ and $\\mathcal{B}^i$, respectively. We refer to $\\mathcal{V}^i$ and $\\mathcal{B}^i$ as vet\\[i\\] and bet\\[i\\] respectively in the code after the Hebrew letters that bear some resemblance. \n",
    "\n",
    "Similarly, we define the smooth versions of $\\bar{A}_{ij}$ and $\\varepsilon_{ij}$ ($a_{ij}$ and $h_{ij}$, respectively) via\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{A}_{ij} &= \\text{scalefactor[i]}\\ \\text{scalefactor[j]}\\  a_{ij} \\\\\n",
    "\\varepsilon_{ij} &= \\text{scalefactor[i]}\\ \\text{scalefactor[j]}\\  h_{ij},\n",
    "\\end{align}\n",
    "\n",
    "where in this case we *multiply* due to the fact that these tensors are purely covariant (as opposed to contravariant). To slightly simplify the notation, in NRPy+ we define the *rescaling matrices* `ReU[i]` and `ReDD[i][j]`, such that\n",
    "\n",
    "\\begin{align}\n",
    "\\text{ReU[i]} &= 1 / \\text{scalefactor[i]} \\\\\n",
    "\\text{ReDD[i][j]} &= \\text{scalefactor[i] scalefactor[j]}.\n",
    "\\end{align}\n",
    "\n",
    "Thus, for example, $\\bar{A}_{ij}$ and $\\bar{\\Lambda}^i$ can be expressed as the [Hadamard product](https://en.wikipedia.org/w/index.php?title=Hadamard_product_(matrices)&oldid=852272177) of matrices :\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{A}_{ij} &= \\mathbf{ReDD}\\circ\\mathbf{a} = \\text{ReDD[i][j]} a_{ij} \\\\\n",
    "\\bar{\\Lambda}^{i} &= \\mathbf{ReU}\\circ\\mathbf{\\lambda} = \\text{ReU[i]} \\lambda^i,\n",
    "\\end{align}\n",
    "where no sums are implied by the repeated indices.\n",
    "\n",
    "Further, since the scale factors are *time independent*, \n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t \\bar{A}_{ij} &= \\text{ReDD[i][j]}\\  \\partial_t a_{ij} \\\\\n",
    "\\partial_t \\bar{\\gamma}_{ij} &= \\partial_t \\left(\\varepsilon_{ij} + \\hat{\\gamma}_{ij}\\right)\\\\\n",
    "&= \\partial_t \\varepsilon_{ij} \\\\\n",
    "&= \\text{scalefactor[i]}\\ \\text{scalefactor[j]}\\ \\partial_t h_{ij}.\n",
    "\\end{align}\n",
    "\n",
    "Thus instead of taking space or time derivatives of BSSN quantities\n",
    "\n",
    "$$\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}, \\alpha, \\beta^i, B^i\\right\\},$$ \n",
    "\n",
    "across coordinate singularities, we instead factor out the singular scale factors according to this prescription so that space or time derivatives of BSSN quantities are written in terms of finite-difference derivatives of the *rescaled* variables \n",
    "\n",
    "$$\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\},$$ \n",
    "\n",
    "and *exact* expressions for (spatial) derivatives of scale factors. Note that `cf` is the chosen conformal factor (supported choices for `cf` are discussed in [Step 6.a](#phi_ito_cf)). \n",
    "\n",
    "As an example, let's evaluate $\\bar{\\Lambda}^{i}_{\\, ,\\, j}$ according to this prescription:\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{\\Lambda}^{i}_{\\, ,\\, j} &= -\\frac{\\lambda^i}{(\\text{ReU[i]})^2}\\ \\partial_j \\left(\\text{ReU[i]}\\right) + \\frac{\\partial_j \\lambda^i}{\\text{ReU[i]}} \\\\\n",
    "&= -\\frac{\\lambda^i}{(\\text{ReU[i]})^2}\\ \\text{ReUdD[i][j]} + \\frac{\\partial_j \\lambda^i}{\\text{ReU[i]}}.\n",
    "\\end{align}\n",
    "\n",
    "Here, the derivative `ReUdD[i][j]` **is computed symbolically and exactly** using SymPy, and the derivative $\\partial_j \\lambda^i$ represents a derivative of a *smooth* quantity (so long as $\\bar{\\Lambda}^{i}$ is smooth in the Cartesian basis)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn_basic_tensors'></a>\n",
    "\n",
    "## Step 3.a: `BSSN_basic_tensors()`: Define all basic conformal BSSN tensors $\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\bar{\\Lambda}^{i}, \\beta^i, B^i\\right\\}$ in terms of BSSN gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_basic_tensors}$$\n",
    "\n",
    "The `BSSN_vars__tensors()` function defines the tensorial BSSN quantities $\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\bar{\\Lambda}^{i}, \\beta^i, B^i\\right\\}$, in terms of the rescaled \"base\" tensorial quantities $\\left\\{h_{i j},a_{i j}, \\lambda^{i}, \\mathcal{V}^i, \\mathcal{B}^i\\right\\},$ respectively:\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{\\gamma}_{i j} &= \\hat{\\gamma}_{ij} + \\varepsilon_{ij}, \\text{ where } \\varepsilon_{ij} = h_{ij} \\circ \\text{ReDD[i][j]} \\\\\n",
    "\\bar{A}_{i j} &= a_{ij} \\circ \\text{ReDD[i][j]} \\\\\n",
    "\\bar{\\Lambda}^{i} &= \\lambda^i \\circ \\text{ReU[i]} \\\\\n",
    "\\beta^{i} &= \\mathcal{V}^i \\circ \\text{ReU[i]} \\\\\n",
    "B^{i} &= \\mathcal{B}^i \\circ \\text{ReU[i]}\n",
    "\\end{align}\n",
    "\n",
    "Rescaling vectors and tensors are built upon the scale factors for the chosen (in general, singular) coordinate system, which are defined in NRPy+'s [reference_metric.py](../edit/reference_metric.py) ([Tutorial](Tutorial-Reference_Metric.ipynb)), and the rescaled variables are defined in the stub function [BSSN/BSSN_rescaled_vars.py](../edit/BSSN/BSSN_rescaled_vars.py). \n",
    "\n",
    "Here we implement `BSSN_vars__tensors()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3.a: Define all basic conformal BSSN tensors in terms of BSSN gridfunctions\n",
    "\n",
    "# Step 3.a.i: gammabarDD and AbarDD:\n",
    "gammabarDD = ixp.zerorank2()\n",
    "AbarDD     = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # gammabar_{ij}  = h_{ij}*ReDD[i][j] + gammahat_{ij}\n",
    "        gammabarDD[i][j] = hDD[i][j]*rfm.ReDD[i][j] + rfm.ghatDD[i][j]\n",
    "        # Abar_{ij}      = a_{ij}*ReDD[i][j]\n",
    "        AbarDD[i][j]     = aDD[i][j]*rfm.ReDD[i][j]\n",
    "\n",
    "# Step 3.a.ii: LambdabarU, betaU, and BU:\n",
    "LambdabarU = ixp.zerorank1()\n",
    "betaU      = ixp.zerorank1()\n",
    "BU         = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    LambdabarU[i] = lambdaU[i]*rfm.ReU[i]\n",
    "    betaU[i]      = vetU[i]   *rfm.ReU[i]\n",
    "    BU[i]         = betU[i]   *rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn_barred_metric__inverse_and_derivs'></a>\n",
    "\n",
    "# Step 4: `gammabar__inverse_and_derivs()`: $\\bar{\\gamma}^{ij}$, and spatial derivatives of $\\bar{\\gamma}_{ij}$ including $\\bar{\\Gamma}^{i}_{jk}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_barred_metric__inverse_and_derivs}$$\n",
    "\n",
    "<a id='bssn_barred_metric__inverse'></a>\n",
    "\n",
    "## Step 4.a: Inverse conformal 3-metric: $\\bar{\\gamma^{ij}}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_barred_metric__inverse}$$\n",
    "\n",
    "Since $\\bar{\\gamma}^{ij}$ is the inverse of $\\bar{\\gamma}_{ij}$, we apply a $3\\times 3$ symmetric matrix inversion to compute $\\bar{\\gamma}^{ij}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4.a: Inverse conformal 3-metric gammabarUU:\n",
    "# Step 4.a.i: gammabarUU:\n",
    "gammabarUU, dummydet = ixp.symm_matrix_inverter3x3(gammabarDD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn_barred_metric__derivs'></a>\n",
    "\n",
    "## Step 4.b: Derivatives of the conformal 3-metric $\\bar{\\gamma}_{ij,k}$ and $\\bar{\\gamma}_{ij,kl}$, and associated \"barred\" Christoffel symbols $\\bar{\\Gamma}^{i}_{jk}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_barred_metric__derivs}$$\n",
    "\n",
    "In the BSSN-in-curvilinear coordinates formulation, all quantities must be defined in terms of rescaled quantities $h_{ij}$ and their derivatives (evaluated using finite differences), as well as reference-metric quantities and their derivatives (evaluated exactly using SymPy). \n",
    "\n",
    "For example, $\\bar{\\gamma}_{ij,k}$ is given by:\n",
    "\\begin{align}\n",
    "\\bar{\\gamma}_{ij,k} &= \\partial_k \\bar{\\gamma}_{ij} \\\\\n",
    "&= \\partial_k \\left(\\hat{\\gamma}_{ij} + \\varepsilon_{ij}\\right) \\\\\n",
    "&= \\partial_k \\left(\\hat{\\gamma}_{ij} + h_{ij} \\text{ReDD[i][j]}\\right) \\\\\n",
    "&= \\hat{\\gamma}_{ij,k} + h_{ij,k} \\text{ReDD[i][j]} + h_{ij} \\text{ReDDdD[i][j][k]},\n",
    "\\end{align}\n",
    "where `ReDDdD[i][j][k]` is computed within `rfm.reference_metric()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4.b.i gammabarDDdD[i][j][k]\n",
    "#            = \\hat{\\gamma}_{ij,k} + h_{ij,k} \\text{ReDD[i][j]} + h_{ij} \\text{ReDDdD[i][j][k]}.\n",
    "\n",
    "gammabarDD_dD = ixp.zerorank3()\n",
    "hDD_dD = ixp.declarerank3(\"hDD_dD\",\"sym01\")\n",
    "hDD_dupD = ixp.declarerank3(\"hDD_dupD\",\"sym01\")\n",
    "gammabarDD_dupD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            gammabarDD_dD[i][j][k] = rfm.ghatDDdD[i][j][k] + \\\n",
    "                                     hDD_dD[i][j][k]*rfm.ReDD[i][j] + hDD[i][j]*rfm.ReDDdD[i][j][k]\n",
    "\n",
    "            # Compute associated upwinded derivative, needed for the \\bar{\\gamma}_{ij} RHS\n",
    "            gammabarDD_dupD[i][j][k] = rfm.ghatDDdD[i][j][k] + \\\n",
    "                                       hDD_dupD[i][j][k]*rfm.ReDD[i][j] + hDD[i][j]*rfm.ReDDdD[i][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By extension, the second derivative $\\bar{\\gamma}_{ij,kl}$ is given by\n",
    "\\begin{align}\n",
    "\\bar{\\gamma}_{ij,kl} &= \\partial_l \\left(\\hat{\\gamma}_{ij,k} + h_{ij,k} \\text{ReDD[i][j]} + h_{ij} \\text{ReDDdD[i][j][k]}\\right)\\\\\n",
    "&= \\hat{\\gamma}_{ij,kl} + h_{ij,kl} \\text{ReDD[i][j]} + h_{ij,k} \\text{ReDDdD[i][j][l]} + h_{ij,l} \\text{ReDDdD[i][j][k]} + h_{ij} \\text{ReDDdDD[i][j][k][l]}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4.b.ii: Compute gammabarDD_dDD in terms of the rescaled BSSN quantity hDD\n",
    "#      and its derivatives, as well as the reference metric and rescaling\n",
    "#      matrix, and its derivatives (expression given below):\n",
    "hDD_dDD = ixp.declarerank4(\"hDD_dDD\",\"sym01_sym23\")\n",
    "gammabarDD_dDD = ixp.zerorank4()\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",
    "                # gammabar_{ij,kl} = gammahat_{ij,kl}\n",
    "                #                  + h_{ij,kl} ReDD[i][j]\n",
    "                #                  + h_{ij,k} ReDDdD[i][j][l] + h_{ij,l} ReDDdD[i][j][k]\n",
    "                #                  + h_{ij} ReDDdDD[i][j][k][l]\n",
    "                gammabarDD_dDD[i][j][k][l]  = rfm.ghatDDdDD[i][j][k][l]\n",
    "                gammabarDD_dDD[i][j][k][l] += hDD_dDD[i][j][k][l]*rfm.ReDD[i][j]\n",
    "                gammabarDD_dDD[i][j][k][l] += hDD_dD[i][j][k]*rfm.ReDDdD[i][j][l] + \\\n",
    "                                              hDD_dD[i][j][l]*rfm.ReDDdD[i][j][k]\n",
    "                gammabarDD_dDD[i][j][k][l] += hDD[i][j]*rfm.ReDDdDD[i][j][k][l]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we compute the Christoffel symbol associated with the barred 3-metric: $\\bar{\\Gamma}^{i}_{kl}$:\n",
    "$$\n",
    "\\bar{\\Gamma}^{i}_{kl} = \\frac{1}{2} \\bar{\\gamma}^{im} \\left(\\bar{\\gamma}_{mk,l} + \\bar{\\gamma}_{ml,k} - \\bar{\\gamma}_{kl,m} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4.b.iii: Define barred Christoffel symbol \\bar{\\Gamma}^{i}_{kl} = GammabarUDD[i][k][l] (see expression below)\n",
    "GammabarUDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for m in range(DIM):\n",
    "                # Gammabar^i_{kl} = 1/2 * gammabar^{im} ( gammabar_{mk,l} + gammabar_{ml,k} - gammabar_{kl,m}):\n",
    "                GammabarUDD[i][k][l] += sp.Rational(1,2)*gammabarUU[i][m]* \\\n",
    "                                        (gammabarDD_dD[m][k][l] + gammabarDD_dD[m][l][k] - gammabarDD_dD[k][l][m])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='detgammabar_and_derivs'></a>\n",
    "\n",
    "# Step 5: `detgammabar_and_derivs()`: $\\det \\bar{\\gamma}_{ij}$ and its derivatives \\[Back to [top](#toc)\\]\n",
    "$$\\label{detgammabar_and_derivs}$$\n",
    "\n",
    "\n",
    "As described just before Section III of [Baumgarte *et al* (2012)](https://arxiv.org/pdf/1211.6632.pdf), we are free to choose $\\det \\bar{\\gamma}_{ij}$, which should remain fixed in time.\n",
    "\n",
    "As in [Baumgarte *et al* (2012)](https://arxiv.org/pdf/1211.6632.pdf) generally we make the choice $\\det \\bar{\\gamma}_{ij} = \\det \\hat{\\gamma}_{ij}$, but *this need not be the case; we could choose to set $\\det \\bar{\\gamma}_{ij}$ to another expression.*\n",
    "\n",
    "In case we do not choose to set $\\det \\bar{\\gamma}_{ij}/\\det \\hat{\\gamma}_{ij}=1$, below we begin the implementation of a gridfunction, `detgbarOverdetghat`, which defines an alternative expression in its place. \n",
    "\n",
    "$\\det \\bar{\\gamma}_{ij}/\\det \\hat{\\gamma}_{ij}$=`detgbarOverdetghat`$\\ne 1$ is not yet implemented. However, we can define `detgammabar` and its derivatives in terms of a generic `detgbarOverdetghat` and $\\det \\hat{\\gamma}_{ij}$ and their derivatives:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{detgammabar} &= \\det \\bar{\\gamma}_{ij} = \\text{detgbarOverdetghat} \\cdot \\left(\\det \\hat{\\gamma}_{ij}\\right) \\\\\n",
    "\\text{detgammabar}\\_\\text{dD[k]} &= \\left(\\det \\bar{\\gamma}_{ij}\\right)_{,k} = \\text{detgbarOverdetghat}\\_\\text{dD[k]} \\det \\hat{\\gamma}_{ij} +  \\text{detgbarOverdetghat} \\left(\\det \\hat{\\gamma}_{ij}\\right)_{,k} \\\\\n",
    "\\end{align}\n",
    "https://en.wikipedia.org/wiki/Determinant#Properties_of_the_determinant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5: det(gammabarDD) and its derivatives\n",
    "detgbarOverdetghat = sp.sympify(1)\n",
    "detgbarOverdetghat_dD = ixp.zerorank1()\n",
    "detgbarOverdetghat_dDD = ixp.zerorank2()\n",
    "\n",
    "if par.parval_from_str(thismodule+\"::detgbarOverdetghat_equals_one\") == \"False\":\n",
    "    print(\"Error: detgbarOverdetghat_equals_one=\\\"False\\\" is not fully implemented yet.\")\n",
    "    sys.exit(1)\n",
    "## Approach for implementing detgbarOverdetghat_equals_one=False:\n",
    "#     detgbarOverdetghat = gri.register_gridfunctions(\"AUX\", [\"detgbarOverdetghat\"])\n",
    "#     detgbarOverdetghatInitial = gri.register_gridfunctions(\"AUX\", [\"detgbarOverdetghatInitial\"])\n",
    "#     detgbarOverdetghat_dD = ixp.declarerank1(\"detgbarOverdetghat_dD\")\n",
    "#     detgbarOverdetghat_dDD = ixp.declarerank2(\"detgbarOverdetghat_dDD\", \"sym01\")\n",
    "\n",
    "# Step 5.b: Define detgammabar, detgammabar_dD, and detgammabar_dDD (needed for\n",
    "#           \\partial_t \\bar{\\Lambda}^i below)detgammabar = detgbarOverdetghat * rfm.detgammahat\n",
    "detgammabar = detgbarOverdetghat * rfm.detgammahat\n",
    "\n",
    "detgammabar_dD = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    detgammabar_dD[i] = detgbarOverdetghat_dD[i] * rfm.detgammahat + detgbarOverdetghat * rfm.detgammahatdD[i]\n",
    "\n",
    "detgammabar_dDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        detgammabar_dDD[i][j] = detgbarOverdetghat_dDD[i][j] * rfm.detgammahat + \\\n",
    "                                detgbarOverdetghat_dD[i] * rfm.detgammahatdD[j] + \\\n",
    "                                detgbarOverdetghat_dD[j] * rfm.detgammahatdD[i] + \\\n",
    "                                detgbarOverdetghat * rfm.detgammahatdDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='abar_quantities'></a>\n",
    "\n",
    "# Step 6: `AbarUU_AbarUD_trAbar_AbarDD_dD()`: Quantities related to conformal traceless extrinsic curvature $\\bar{A}_{ij}$: $\\bar{A}^{ij}$, $\\bar{A}^i_j$, and $\\bar{A}^k_k$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{abar_quantities}$$\n",
    "\n",
    "$\\bar{A}^{ij}$ is given by application of the raising operators (a.k.a., the inverse 3-metric) $\\bar{\\gamma}^{jk}$ on both of the covariant (\"down\") components:\n",
    "$$\n",
    "\\bar{A}^{ij} = \\bar{\\gamma}^{ik}\\bar{\\gamma}^{jl} \\bar{A}_{kl}.\n",
    "$$\n",
    "\n",
    "$\\bar{A}^i_j$ is given by a single application of the raising operator (a.k.a., the inverse 3-metric) $\\bar{\\gamma}^{ik}$ on $\\bar{A}_{kj}$:\n",
    "$$\n",
    "\\bar{A}^i_j = \\bar{\\gamma}^{ik}\\bar{A}_{kj}.\n",
    "$$\n",
    "\n",
    "The trace of $\\bar{A}_{ij}$, $\\bar{A}^k_k$, is given by a contraction with the barred 3-metric:\n",
    "$$\n",
    "\\text{Tr}(\\bar{A}_{ij}) = \\bar{A}^k_k = \\bar{\\gamma}^{kj}\\bar{A}_{jk}.\n",
    "$$\n",
    "\n",
    "Note that while $\\bar{A}_{ij}$ is defined as the *traceless* conformal extrinsic curvature, it may acquire a nonzero trace (assuming the initial data impose tracelessness) due to numerical error. $\\text{Tr}(\\bar{A}_{ij})$ is included in the BSSN equations to drive $\\text{Tr}(\\bar{A}_{ij})$ to zero.\n",
    "\n",
    "In terms of rescaled BSSN quantities, $\\bar{A}_{ij}$ is given by\n",
    "$$\n",
    "\\bar{A}_{ij} = \\text{ReDD[i][j]} a_{ij},\n",
    "$$\n",
    "so in terms of the same quantities, $\\bar{A}_{ij,k}$ is given by\n",
    "$$\n",
    "\\bar{A}_{ij,k} = \\text{ReDDdD[i][j][k]} a_{ij} + \\text{ReDD[i][j]} a_{ij,k}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6: Quantities related to conformal traceless extrinsic curvature\n",
    "\n",
    "# Step 6.a.i: Compute Abar^{ij} in terms of Abar_{ij} and gammabar^{ij}\n",
    "AbarUU = ixp.zerorank2()\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",
    "                # Abar^{ij} = gammabar^{ik} gammabar^{jl} Abar_{kl}\n",
    "                AbarUU[i][j] += gammabarUU[i][k]*gammabarUU[j][l]*AbarDD[k][l]\n",
    "\n",
    "# Step 6.a.ii: Compute Abar^i_j in terms of Abar_{ij} and gammabar^{ij}\n",
    "AbarUD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            # Abar^i_j = gammabar^{ik} Abar_{kj}\n",
    "            AbarUD[i][j] += gammabarUU[i][k]*AbarDD[k][j]\n",
    "\n",
    "# Step 6.a.iii: Compute Abar^k_k = trace of Abar:\n",
    "trAbar = sp.sympify(0)\n",
    "for k in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # Abar^k_k = gammabar^{kj} Abar_{jk}\n",
    "        trAbar += gammabarUU[k][j]*AbarDD[j][k]\n",
    "\n",
    "# Step 6.a.iv: Compute Abar_{ij,k}\n",
    "AbarDD_dD   = ixp.zerorank3()\n",
    "AbarDD_dupD = ixp.zerorank3()\n",
    "aDD_dD   = ixp.declarerank3(\"aDD_dD\"  ,\"sym01\")\n",
    "aDD_dupD = ixp.declarerank3(\"aDD_dupD\",\"sym01\")\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            AbarDD_dupD[i][j][k] = rfm.ReDDdD[i][j][k]*aDD[i][j] + rfm.ReDD[i][j]*aDD_dupD[i][j][k]\n",
    "            AbarDD_dD[i][j][k]   = rfm.ReDDdD[i][j][k]*aDD[i][j] + rfm.ReDD[i][j]*aDD_dD[  i][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rbar'></a>\n",
    "\n",
    "# Step 7: `RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()`: The conformal (\"barred\") Ricci tensor $\\bar{R}_{ij}$ and associated quantities \\[Back to [top](#toc)\\]\n",
    "$$\\label{rbar}$$\n",
    "\n",
    "Let's compute perhaps the most complicated expression in the BSSN evolution equations, the conformal Ricci tensor:\n",
    "\n",
    "\\begin{align}\n",
    "  \\bar{R}_{i j} {} = {} & - \\frac{1}{2} \\bar{\\gamma}^{k l} \\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} + \\bar{\\gamma}_{k(i} \\hat{D}_{j)} \\bar{\\Lambda}^{k} + \\Delta^{k} \\Delta_{(i j) k} \\nonumber \\\\\n",
    "  & + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right ) \\; .\n",
    "\\end{align}\n",
    "\n",
    "Let's tackle the $\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j}$ term first:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rbar_part1'></a>\n",
    "\n",
    "## Step 7.a: Conformal Ricci tensor, part 1: The $\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j}$ term \\[Back to [top](#toc)\\]\n",
    "$$\\label{rbar_part1}$$\n",
    "\n",
    "First note that the covariant derivative of a metric with respect to itself is zero\n",
    "$$\\hat{D}_{l} \\hat{\\gamma}_{ij} = 0,$$\n",
    "so \n",
    "$$\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} = \\hat{D}_{k} \\hat{D}_{l} \\left(\\hat{\\gamma}_{i j} + \\varepsilon_{ij}\\right) = \\hat{D}_{k} \\hat{D}_{l} \\varepsilon_{ij}.$$\n",
    "\n",
    "Next, the covariant derivative of a tensor is given by (from the [wikipedia article on covariant differentiation](https://en.wikipedia.org/wiki/Covariant_derivative)):\n",
    "\\begin{align}\n",
    "  {(\\nabla_{e_c} T)^{a_1 \\ldots a_r}}_{b_1 \\ldots b_s} = {}\n",
    "    &\\frac{\\partial}{\\partial x^c}{T^{a_1 \\ldots a_r}}_{b_1 \\ldots b_s} \\\\\n",
    "    &+ \\,{\\Gamma ^{a_1}}_{dc} {T^{d a_2 \\ldots a_r}}_{b_1 \\ldots b_s} + \\cdots + {\\Gamma^{a_r}}_{dc} {T^{a_1 \\ldots a_{r-1}d}}_{b_1 \\ldots b_s} \\\\\n",
    "    &-\\,{\\Gamma^d}_{b_1 c} {T^{a_1 \\ldots a_r}}_{d b_2 \\ldots b_s} - \\cdots - {\\Gamma^d}_{b_s c} {T^{a_1 \\ldots a_r}}_{b_1 \\ldots b_{s-1} d}.\n",
    "\\end{align}\n",
    "\n",
    "Therefore, \n",
    "$$\\hat{D}_{l} \\bar{\\gamma}_{i j} = \\hat{D}_{l} \\varepsilon_{i j} = \\varepsilon_{i j,l} - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j} -\\hat{\\Gamma}^m_{j l} \\varepsilon_{i m}.$$\n",
    "\n",
    "Since the covariant first derivative is a tensor, the covariant second derivative is given by (same as [Eq. 27 in Baumgarte et al (2012)](https://arxiv.org/pdf/1211.6632.pdf))\n",
    "\\begin{align}\n",
    "\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} &= \\hat{D}_{k} \\hat{D}_{l} \\varepsilon_{i j}  \\\\\n",
    "&= \\partial_k \\hat{D}_{l} \\varepsilon_{i j}\n",
    " - \\hat{\\Gamma}^m_{lk} \\left(\\hat{D}_{m} \\varepsilon_{i j}\\right) \n",
    " - \\hat{\\Gamma}^m_{ik} \\left(\\hat{D}_{l} \\varepsilon_{m j}\\right)\n",
    " - \\hat{\\Gamma}^m_{jk} \\left(\\hat{D}_{l} \\varepsilon_{i m}\\right),\n",
    "\\end{align}\n",
    "\n",
    "where the first term is the partial derivative of the expression already derived for $\\hat{D}_{l} \\varepsilon_{i j}$:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_k \\hat{D}_{l} \\varepsilon_{i j} &= \\partial_k \\left(\\varepsilon_{ij,l} - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j} -\\hat{\\Gamma}^m_{j l} \\varepsilon_{i m} \\right) \\\\\n",
    "&= \\varepsilon_{ij,lk} - \\hat{\\Gamma}^m_{i l,k} \\varepsilon_{m j} - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j,k} - \\hat{\\Gamma}^m_{j l,k} \\varepsilon_{i m} - \\hat{\\Gamma}^m_{j l} \\varepsilon_{i m,k}.\n",
    "\\end{align}\n",
    "\n",
    "In terms of the evolved quantity $h_{ij}$, the derivatives of $\\varepsilon_{ij}$ are given by:\n",
    "\\begin{align}\n",
    "\\varepsilon_{ij,k} &= \\partial_k \\left(h_{ij} \\text{ReDD[i][j]}\\right) \\\\\n",
    "&= h_{ij,k} \\text{ReDD[i][j]} + h_{ij} \\text{ReDDdD[i][j][k]},\n",
    "\\end{align}\n",
    "and\n",
    "\\begin{align}\n",
    "\\varepsilon_{ij,kl} &= \\partial_l \\left(h_{ij,k} \\text{ReDD[i][j]} + h_{ij} \\text{ReDDdD[i][j][k]} \\right)\\\\\n",
    "&= h_{ij,kl} \\text{ReDD[i][j]} + h_{ij,k} \\text{ReDDdD[i][j][l]} + h_{ij,l} \\text{ReDDdD[i][j][k]} + h_{ij} \\text{ReDDdDD[i][j][k][l]}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7: Conformal Ricci tensor, part 1: The \\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} term\n",
    "\n",
    "# Step 7.a.i: Define \\varepsilon_{ij} = epsDD[i][j]\n",
    "epsDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        epsDD[i][j] = hDD[i][j]*rfm.ReDD[i][j]\n",
    "\n",
    "# Step 7.a.ii: Define epsDD_dD[i][j][k]\n",
    "hDD_dD = ixp.declarerank3(\"hDD_dD\",\"sym01\")\n",
    "epsDD_dD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            epsDD_dD[i][j][k] = hDD_dD[i][j][k]*rfm.ReDD[i][j] + hDD[i][j]*rfm.ReDDdD[i][j][k]\n",
    "\n",
    "# Step 7.a.iii: Define epsDD_dDD[i][j][k][l]\n",
    "hDD_dDD = ixp.declarerank4(\"hDD_dDD\",\"sym01_sym23\")\n",
    "epsDD_dDD = ixp.zerorank4()\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",
    "                epsDD_dDD[i][j][k][l] = hDD_dDD[i][j][k][l]*rfm.ReDD[i][j] + \\\n",
    "                                        hDD_dD[i][j][k]*rfm.ReDDdD[i][j][l] + \\\n",
    "                                        hDD_dD[i][j][l]*rfm.ReDDdD[i][j][k] + \\\n",
    "                                        hDD[i][j]*rfm.ReDDdDD[i][j][k][l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next compute three quantities derived above:\n",
    "\n",
    "* `gammabarDD_DhatD[i][j][l]` = $\\hat{D}_{l} \\bar{\\gamma}_{i j} = \\hat{D}_{l} \\varepsilon_{i j} = \\varepsilon_{i j,l} - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j} -\\hat{\\Gamma}^m_{j l} \\varepsilon_{i m}$,\n",
    "* `gammabarDD_DhatD\\_dD[i][j][l][k]` = $\\partial_k \\hat{D}_{l} \\bar{\\gamma}_{i j} = \\partial_k \\hat{D}_{l} \\varepsilon_{i j} = \\varepsilon_{ij,lk} - \\hat{\\Gamma}^m_{i l,k} \\varepsilon_{m j} - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j,k} - \\hat{\\Gamma}^m_{j l,k} \\varepsilon_{i m} - \\hat{\\Gamma}^m_{j l} \\varepsilon_{i m,k}$, and\n",
    "* `gammabarDD_DhatDD[i][j][l][k]` = $\\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} = \\partial_k \\hat{D}_{l} \\varepsilon_{i j} - \\hat{\\Gamma}^m_{lk} \\left(\\hat{D}_{m} \\varepsilon_{i j}\\right)  - \\hat{\\Gamma}^m_{ik} \\left(\\hat{D}_{l} \\varepsilon_{m j}\\right) - \\hat{\\Gamma}^m_{jk} \\left(\\hat{D}_{l} \\varepsilon_{i m}\\right)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7.a.iv: DhatgammabarDDdD[i][j][l] = \\bar{\\gamma}_{ij;\\hat{l}}\n",
    "# \\bar{\\gamma}_{ij;\\hat{l}} = \\varepsilon_{i j,l}\n",
    "#                           - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j}\n",
    "#                           - \\hat{\\Gamma}^m_{j l} \\varepsilon_{i m}\n",
    "gammabarDD_dHatD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            gammabarDD_dHatD[i][j][l] = epsDD_dD[i][j][l]\n",
    "            for m in range(DIM):\n",
    "                gammabarDD_dHatD[i][j][l] += - rfm.GammahatUDD[m][i][l]*epsDD[m][j] \\\n",
    "                                             - rfm.GammahatUDD[m][j][l]*epsDD[i][m]\n",
    "\n",
    "# Step 7.a.v: \\bar{\\gamma}_{ij;\\hat{l},k} = DhatgammabarDD_dHatD_dD[i][j][l][k]:\n",
    "#        \\bar{\\gamma}_{ij;\\hat{l},k} = \\varepsilon_{ij,lk}\n",
    "#                                      - \\hat{\\Gamma}^m_{i l,k} \\varepsilon_{m j}\n",
    "#                                      - \\hat{\\Gamma}^m_{i l} \\varepsilon_{m j,k}\n",
    "#                                      - \\hat{\\Gamma}^m_{j l,k} \\varepsilon_{i m}\n",
    "#                                      - \\hat{\\Gamma}^m_{j l} \\varepsilon_{i m,k}\n",
    "gammabarDD_dHatD_dD = ixp.zerorank4()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for k in range(DIM):\n",
    "                gammabarDD_dHatD_dD[i][j][l][k] = epsDD_dDD[i][j][l][k]\n",
    "                for m in range(DIM):\n",
    "                    gammabarDD_dHatD_dD[i][j][l][k] += -rfm.GammahatUDDdD[m][i][l][k]*epsDD[m][j]  \\\n",
    "                                                       -rfm.GammahatUDD[m][i][l]*epsDD_dD[m][j][k] \\\n",
    "                                                       -rfm.GammahatUDDdD[m][j][l][k]*epsDD[i][m]  \\\n",
    "                                                       -rfm.GammahatUDD[m][j][l]*epsDD_dD[i][m][k]\n",
    "\n",
    "# Step 7.a.vi: \\bar{\\gamma}_{ij;\\hat{l}\\hat{k}} = DhatgammabarDD_dHatDD[i][j][l][k]\n",
    "#          \\bar{\\gamma}_{ij;\\hat{l}\\hat{k}} = \\partial_k \\hat{D}_{l} \\varepsilon_{i j}\n",
    "#                                           - \\hat{\\Gamma}^m_{lk} \\left(\\hat{D}_{m} \\varepsilon_{i j}\\right)\n",
    "#                                           - \\hat{\\Gamma}^m_{ik} \\left(\\hat{D}_{l} \\varepsilon_{m j}\\right)\n",
    "#                                           - \\hat{\\Gamma}^m_{jk} \\left(\\hat{D}_{l} \\varepsilon_{i m}\\right)\n",
    "gammabarDD_dHatDD = ixp.zerorank4()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for l in range(DIM):\n",
    "            for k in range(DIM):\n",
    "                gammabarDD_dHatDD[i][j][l][k] = gammabarDD_dHatD_dD[i][j][l][k]\n",
    "                for m in range(DIM):\n",
    "                    gammabarDD_dHatDD[i][j][l][k] += - rfm.GammahatUDD[m][l][k]*gammabarDD_dHatD[i][j][m] \\\n",
    "                                                     - rfm.GammahatUDD[m][i][k]*gammabarDD_dHatD[m][j][l] \\\n",
    "                                                     - rfm.GammahatUDD[m][j][k]*gammabarDD_dHatD[i][m][l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rbar_part2'></a>\n",
    "\n",
    "## Step 7.b: Conformal Ricci tensor, part 2: The $\\bar{\\gamma}_{k(i} \\hat{D}_{j)} \\bar{\\Lambda}^{k}$ term \\[Back to [top](#toc)\\]\n",
    "$$\\label{rbar_part2}$$\n",
    "\n",
    "By definition, the index symmetrization operation is given by:\n",
    "$$\\bar{\\gamma}_{k(i} \\hat{D}_{j)} \\bar{\\Lambda}^{k} = \\frac{1}{2} \\left( \\bar{\\gamma}_{ki} \\hat{D}_{j} \\bar{\\Lambda}^{k} + \\bar{\\gamma}_{kj} \\hat{D}_{i} \\bar{\\Lambda}^{k} \\right),$$\n",
    "\n",
    "and $\\bar{\\gamma}_{ij}$ is trivially computed ($=\\varepsilon_{ij} + \\hat{\\gamma}_{ij}$) so the only nontrival part to computing this term is in evaluating $\\hat{D}_{j} \\bar{\\Lambda}^{k}$.\n",
    "\n",
    "The covariant derivative is with respect to the hatted metric (i.e. the reference metric), so\n",
    "$$\\hat{D}_{j} \\bar{\\Lambda}^{k} = \\partial_j \\bar{\\Lambda}^{k} + \\hat{\\Gamma}^{k}_{mj} \\bar{\\Lambda}^m,$$\n",
    "except we cannot take derivatives of $\\bar{\\Lambda}^{k}$ directly due to potential issues with coordinate singularities. Instead we write it in terms of the rescaled quantity $\\lambda^k$ via\n",
    "$$\\bar{\\Lambda}^{k} = \\lambda^k \\text{ReU[k]}.$$\n",
    "\n",
    "Then the expression for $\\hat{D}_{j} \\bar{\\Lambda}^{k}$ becomes\n",
    "$$\n",
    "\\hat{D}_{j} \\bar{\\Lambda}^{k} = \\lambda^{k}_{,j} \\text{ReU[k]} + \\lambda^{k} \\text{ReUdD[k][j]} + \\hat{\\Gamma}^{k}_{mj} \\lambda^{m} \\text{ReU[m]},\n",
    "$$\n",
    "and the NRPy+ code for this expression is written"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7.b: Second term of RhatDD: compute \\hat{D}_{j} \\bar{\\Lambda}^{k} = LambarU_dHatD[k][j]\n",
    "lambdaU_dD = ixp.declarerank2(\"lambdaU_dD\",\"nosym\")\n",
    "LambarU_dHatD = ixp.zerorank2()\n",
    "for j in range(DIM):\n",
    "    for k in range(DIM):\n",
    "        LambarU_dHatD[k][j] = lambdaU_dD[k][j]*rfm.ReU[k] + lambdaU[k]*rfm.ReUdD[k][j]\n",
    "        for m in range(DIM):\n",
    "            LambarU_dHatD[k][j] += rfm.GammahatUDD[k][m][j]*lambdaU[m]*rfm.ReU[m]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rbar_part3'></a>\n",
    "\n",
    "## Step 7.c: Conformal Ricci tensor, part 3: The $\\Delta^{k} \\Delta_{(i j) k}  + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right )$ terms \\[Back to [top](#toc)\\]\n",
    "$$\\label{rbar_part3}$$\n",
    "\n",
    "Our goal here is to compute the quantities appearing as the final terms of the conformal Ricci tensor:\n",
    "$$\n",
    "\\Delta^{k} \\Delta_{(i j) k}  + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right).\n",
    "$$\n",
    "\n",
    "* `DGammaUDD[k][i][j]`$= \\Delta^k_{ij}$ is simply the difference in Christoffel symbols: $\\Delta^{k}_{ij} = \\bar{\\Gamma}^i_{jk} - \\hat{\\Gamma}^i_{jk}$, and \n",
    "* `DGammaU[k]`$= \\Delta^k$ is the contraction: $\\bar{\\gamma}^{ij} \\Delta^{k}_{ij}$\n",
    "\n",
    "Adding these expressions to Ricci is straightforward, since $\\bar{\\Gamma}^i_{jk}$ and $\\bar{\\gamma}^{ij}$ were defined above in [Step 4](#bssn_barred_metric__inverse_and_derivs), and $\\hat{\\Gamma}^i_{jk}$ was computed within NRPy+'s `reference_metric()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7.c: Conformal Ricci tensor, part 3: The \\Delta^{k} \\Delta_{(i j) k}\n",
    "#           + \\bar{\\gamma}^{k l}*(2 \\Delta_{k(i}^{m} \\Delta_{j) m l}\n",
    "#           + \\Delta_{i k}^{m} \\Delta_{m j l}) terms\n",
    "\n",
    "# Step 7.c.i: Define \\Delta^i_{jk} = \\bar{\\Gamma}^i_{jk} - \\hat{\\Gamma}^i_{jk} = DGammaUDD[i][j][k]\n",
    "DGammaUDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            DGammaUDD[i][j][k] = GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k]\n",
    "\n",
    "# Step 7.c.ii: Define \\Delta^i = \\bar{\\gamma}^{jk} \\Delta^i_{jk}\n",
    "DGammaU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            DGammaU[i] += gammabarUU[j][k] * DGammaUDD[i][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define $\\Delta_{ijk}=\\bar{\\gamma}_{im}\\Delta^m_{jk}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7.c.iii: Define \\Delta_{ijk} = \\bar{\\gamma}_{im} \\Delta^m_{jk}\n",
    "DGammaDDD = 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",
    "                DGammaDDD[i][j][k] += gammabarDD[i][m] * DGammaUDD[m][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='summing_rbar_terms'></a>\n",
    "\n",
    "## Step 7.d: Summing the terms and defining $\\bar{R}_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{summing_rbar_terms}$$\n",
    "\n",
    "We have now constructed all of the terms going into $\\bar{R}_{ij}$:\n",
    "\n",
    "\\begin{align}\n",
    "  \\bar{R}_{i j} {} = {} & - \\frac{1}{2} \\bar{\\gamma}^{k l} \\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j} + \\bar{\\gamma}_{k(i} \\hat{D}_{j)} \\bar{\\Lambda}^{k} + \\Delta^{k} \\Delta_{(i j) k} \\nonumber \\\\\n",
    "  & + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right ) \\; .\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 7.d: Summing the terms and defining \\bar{R}_{ij}\n",
    "\n",
    "# Step 7.d.i: Add the first term to RbarDD:\n",
    "#         Rbar_{ij} += - \\frac{1}{2} \\bar{\\gamma}^{k l} \\hat{D}_{k} \\hat{D}_{l} \\bar{\\gamma}_{i j}\n",
    "RbarDD = ixp.zerorank2()\n",
    "RbarDDpiece = ixp.zerorank2()\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",
    "                RbarDD[i][j] += -sp.Rational(1,2) * gammabarUU[k][l]*gammabarDD_dHatDD[i][j][l][k]\n",
    "                RbarDDpiece[i][j] += -sp.Rational(1,2) * gammabarUU[k][l]*gammabarDD_dHatDD[i][j][l][k]\n",
    "\n",
    "# Step 7.d.ii: Add the second term to RbarDD:\n",
    "#         Rbar_{ij} += (1/2) * (gammabar_{ki} Lambar^k_{;\\hat{j}} + gammabar_{kj} Lambar^k_{;\\hat{i}})\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            RbarDD[i][j] += sp.Rational(1,2) * (gammabarDD[k][i]*LambarU_dHatD[k][j] + \\\n",
    "                                                gammabarDD[k][j]*LambarU_dHatD[k][i])\n",
    "\n",
    "# Step 7.d.iii: Add the remaining term to RbarDD:\n",
    "#      Rbar_{ij} += \\Delta^{k} \\Delta_{(i j) k} = 1/2 \\Delta^{k} (\\Delta_{i j k} + \\Delta_{j i k})\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            RbarDD[i][j] += sp.Rational(1,2) * DGammaU[k] * (DGammaDDD[i][j][k] + DGammaDDD[j][i][k])\n",
    "\n",
    "# Step 7.d.iv: Add the final term to RbarDD:\n",
    "#      Rbar_{ij} += \\bar{\\gamma}^{k l} (\\Delta^{m}_{k i} \\Delta_{j m l}\n",
    "#                   + \\Delta^{m}_{k j} \\Delta_{i m l}\n",
    "#                   + \\Delta^{m}_{i k} \\Delta_{m j l})\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",
    "                for m in range(DIM):\n",
    "                    RbarDD[i][j] += gammabarUU[k][l] * (DGammaUDD[m][k][i]*DGammaDDD[j][m][l] +\n",
    "                                                        DGammaUDD[m][k][j]*DGammaDDD[i][m][l] +\n",
    "                                                        DGammaUDD[m][i][k]*DGammaDDD[m][j][l])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='beta_derivs'></a>\n",
    "\n",
    "# Step 8: **`betaU_derivs()`**: The unrescaled shift vector $\\beta^i$ spatial derivatives: $\\beta^i_{,j}$ & $\\beta^i_{,jk}$, written in terms of the rescaled shift vector $\\mathcal{V}^i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{beta_derivs}$$\n",
    "\n",
    "This step, which documents the function `betaUbar_and_derivs()` inside the [BSSN.BSSN_unrescaled_and_barred_vars](../edit/BSSN/BSSN_unrescaled_and_barred_vars) module, defines three quantities:\n",
    "\n",
    "[comment]: <> (Fix Link Above: TODO)\n",
    "\n",
    "* `betaU_dD[i][j]`$=\\beta^i_{,j} = \\left(\\mathcal{V}^i \\circ \\text{ReU[i]}\\right)_{,j} = \\mathcal{V}^i_{,j} \\circ \\text{ReU[i]} + \\mathcal{V}^i \\circ \\text{ReUdD[i][j]}$\n",
    "* `betaU_dupD[i][j]`: the same as above, except using *upwinded* finite-difference derivatives to compute $\\mathcal{V}^i_{,j}$ instead of *centered* finite-difference derivatives.\n",
    "* `betaU_dDD[i][j][k]`$=\\beta^i_{,jk} = \\mathcal{V}^i_{,jk} \\circ \\text{ReU[i]} + \\mathcal{V}^i_{,j} \\circ \\text{ReUdD[i][k]} + \\mathcal{V}^i_{,k} \\circ \\text{ReUdD[i][j]}+\\mathcal{V}^i \\circ \\text{ReUdDD[i][j][k]}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 8: The unrescaled shift vector betaU spatial derivatives:\n",
    "#         betaUdD & betaUdDD, written in terms of the\n",
    "#         rescaled shift vector vetU\n",
    "vetU_dD = ixp.declarerank2(\"vetU_dD\",\"nosym\")\n",
    "vetU_dupD = ixp.declarerank2(\"vetU_dupD\",\"nosym\") # Needed for upwinded \\beta^i_{,j}\n",
    "vetU_dDD = ixp.declarerank3(\"vetU_dDD\",\"sym12\")   # Needed for \\beta^i_{,j}\n",
    "betaU_dD = ixp.zerorank2()\n",
    "betaU_dupD = ixp.zerorank2() # Needed for, e.g., \\beta^i RHS\n",
    "betaU_dDD = ixp.zerorank3() # Needed for, e.g., \\bar{\\Lambda}^i RHS\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        betaU_dD[i][j] = vetU_dD[i][j]*rfm.ReU[i] + vetU[i]*rfm.ReUdD[i][j]\n",
    "        betaU_dupD[i][j] = vetU_dupD[i][j]*rfm.ReU[i] + vetU[i]*rfm.ReUdD[i][j] # Needed for \\beta^i RHS\n",
    "        for k in range(DIM):\n",
    "            # Needed for, e.g., \\bar{\\Lambda}^i RHS:\n",
    "            betaU_dDD[i][j][k] = vetU_dDD[i][j][k]*rfm.ReU[i]  + vetU_dD[i][j]*rfm.ReUdD[i][k] + \\\n",
    "                                 vetU_dD[i][k]*rfm.ReUdD[i][j] + vetU[i]*rfm.ReUdDD[i][j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='phi_and_derivs'></a>\n",
    "\n",
    "# Step 9: **`phi_and_derivs()`**: Standard BSSN conformal factor $\\phi$, and its derivatives $\\phi_{,i}$, $\\phi_{,ij}$, $\\bar{D}_j \\phi_i$, and $\\bar{D}_j\\bar{D}_k \\phi_i$, all written in terms of BSSN gridfunctions like $\\text{cf}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{phi_and_derivs}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='phi_ito_cf'></a>\n",
    "\n",
    "## Step 9.a: $\\phi$ in terms of the chosen (possibly non-standard) conformal factor variable $\\text{cf}$ (e.g., $\\text{cf}=\\chi=e^{-4\\phi}$) \\[Back to [top](#toc)\\]\n",
    "$$\\label{phi_ito_cf}$$\n",
    "\n",
    "When solving the BSSN time evolution equations across the coordinate singularity (i.e., the \"puncture\") inside puncture black holes for example, the standard conformal factor $\\phi$ becomes very sharp, whereas $\\chi=e^{-4\\phi}$ is far smoother (see, e.g., [Campanelli, Lousto, Marronetti, and Zlochower (2006)](https://arxiv.org/abs/gr-qc/0511048) for additional discussion). Thus if we choose to rewrite derivatives of $\\phi$ in the BSSN equations in terms of finite-difference derivatives `cf`$=\\chi$, numerical errors will be far smaller near the puncture.\n",
    "\n",
    "The BSSN modules in NRPy+ support three options for the conformal factor variable `cf`:\n",
    "\n",
    "1. `cf`$=\\phi$,\n",
    "1. `cf`$=\\chi=e^{-4\\phi}$, and\n",
    "1. `cf`$=W = e^{-2\\phi}$.\n",
    "\n",
    "The BSSN equations are written in terms of $\\phi$ (actually only $e^{-4\\phi}$ appears) and derivatives of $\\phi$, we now define $e^{-4\\phi}$ and derivatives of $\\phi$ in terms of the chosen `cf`.\n",
    "\n",
    "First, we define the base variables needed within the BSSN equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 9: Standard BSSN conformal factor phi,\n",
    "#         and its partial and covariant derivatives,\n",
    "#         all in terms of BSSN gridfunctions like cf\n",
    "\n",
    "# Step 9.a.i: Define partial derivatives of \\phi in terms of evolved quantity \"cf\":\n",
    "cf_dD = ixp.declarerank1(\"cf_dD\")\n",
    "cf_dupD = ixp.declarerank1(\"cf_dupD\") # Needed for \\partial_t \\phi next.\n",
    "cf_dDD = ixp.declarerank2(\"cf_dDD\",\"sym01\")\n",
    "phi_dD = ixp.zerorank1()\n",
    "phi_dupD = ixp.zerorank1()\n",
    "phi_dDD = ixp.zerorank2()\n",
    "exp_m4phi = sp.sympify(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define $\\phi_{,i}$, $\\phi_{,ij}$, and $e^{-4\\phi}$ for each of the choices of `cf`.\n",
    "\n",
    "For `cf`$=\\phi$, this is trivial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 9.a.ii: Assuming cf=phi, define exp_m4phi, phi_dD,\n",
    "#              phi_dupD (upwind finite-difference version of phi_dD), and phi_DD\n",
    "if par.parval_from_str(thismodule+\"::EvolvedConformalFactor_cf\") == \"phi\":\n",
    "    for i in range(DIM):\n",
    "        phi_dD[i] = cf_dD[i]\n",
    "        phi_dupD[i] = cf_dupD[i]\n",
    "        for j in range(DIM):\n",
    "            phi_dDD[i][j] = cf_dDD[i][j]\n",
    "    exp_m4phi = sp.exp(-4*cf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `cf`$=W=e^{-2\\phi}$, we have\n",
    "\n",
    "* $\\phi_{,i} = -\\text{cf}_{,i} / (2 \\text{cf})$\n",
    "* $\\phi_{,ij} = (-\\text{cf}_{,ij} + \\text{cf}_{,i}\\text{cf}_{,j}/\\text{cf}) / (2 \\text{cf})$\n",
    "* $e^{-4\\phi} = \\text{cf}^2$\n",
    "\n",
    "***Exercise to student: Prove the above relations***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 9.a.iii: Assuming cf=W=e^{-2 phi}, define exp_m4phi, phi_dD,\n",
    "#               phi_dupD (upwind finite-difference version of phi_dD), and phi_DD\n",
    "if par.parval_from_str(thismodule+\"::EvolvedConformalFactor_cf\") == \"W\":\n",
    "    # \\partial_i W = \\partial_i (e^{-2 phi}) = -2 e^{-2 phi} \\partial_i phi\n",
    "    # -> \\partial_i phi = -\\partial_i cf / (2 cf)\n",
    "    for i in range(DIM):\n",
    "        phi_dD[i] = - cf_dD[i] / (2*cf)\n",
    "        phi_dupD[i] = - cf_dupD[i] / (2*cf)\n",
    "        for j in range(DIM):\n",
    "            # \\partial_j \\partial_i phi = - \\partial_j [\\partial_i cf / (2 cf)]\n",
    "            #                           = - cf_{,ij} / (2 cf) + \\partial_i cf \\partial_j cf / (2 cf^2)\n",
    "            phi_dDD[i][j] = (- cf_dDD[i][j] + cf_dD[i]*cf_dD[j] / cf) / (2*cf)\n",
    "    exp_m4phi = cf*cf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `cf`$=W=e^{-4\\phi}$, we have\n",
    "\n",
    "* $\\phi_{,i} = -\\text{cf}_{,i} / (4 \\text{cf})$\n",
    "* $\\phi_{,ij} = (-\\text{cf}_{,ij} + \\text{cf}_{,i}\\text{cf}_{,j}/\\text{cf}) / (4 \\text{cf})$\n",
    "* $e^{-4\\phi} = \\text{cf}$\n",
    "\n",
    "***Exercise to student: Prove the above relations***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 9.a.iv: Assuming cf=chi=e^{-4 phi}, define exp_m4phi, phi_dD,\n",
    "#              phi_dupD (upwind finite-difference version of phi_dD), and phi_DD\n",
    "if par.parval_from_str(thismodule+\"::EvolvedConformalFactor_cf\") == \"chi\":\n",
    "    # \\partial_i chi = \\partial_i (e^{-4 phi}) = -4 e^{-4 phi} \\partial_i phi\n",
    "    # -> \\partial_i phi = -\\partial_i cf / (4 cf)\n",
    "    for i in range(DIM):\n",
    "        phi_dD[i] = - cf_dD[i] / (4*cf)\n",
    "        phi_dupD[i] = - cf_dupD[i] / (4*cf)\n",
    "        for j in range(DIM):\n",
    "            # \\partial_j \\partial_i phi = - \\partial_j [\\partial_i cf / (4 cf)]\n",
    "            #                           = - cf_{,ij} / (4 cf) + \\partial_i cf \\partial_j cf / (4 cf^2)\n",
    "            phi_dDD[i][j] = (- cf_dDD[i][j] + cf_dD[i]*cf_dD[j] / cf) / (4*cf)\n",
    "    exp_m4phi = cf\n",
    "\n",
    "# Step 9.a.v: Error out if unsupported EvolvedConformalFactor_cf choice is made:\n",
    "cf_choice = par.parval_from_str(thismodule+\"::EvolvedConformalFactor_cf\")\n",
    "if cf_choice not in ('phi', 'W', 'chi'):\n",
    "    print(\"Error: EvolvedConformalFactor_cf == \"+par.parval_from_str(thismodule+\"::EvolvedConformalFactor_cf\")+\" unsupported!\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='phi_covariant_derivs'></a>\n",
    "\n",
    "## Step 9.b: Covariant derivatives of $\\phi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{phi_covariant_derivs}$$\n",
    "\n",
    "Since $\\phi$ is a scalar, $\\bar{D}_i \\phi = \\partial_i \\phi$.\n",
    "\n",
    "Thus the second covariant derivative is given by\n",
    "\\begin{align}\n",
    "\\bar{D}_i \\bar{D}_j \\phi &= \\phi_{;\\bar{i}\\bar{j}} = \\bar{D}_i \\phi_{,j}\\\\\n",
    "                         &= \\phi_{,ij} - \\bar{\\Gamma}^k_{ij} \\phi_{,k}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 9.b: Define phi_dBarD = phi_dD (since phi is a scalar) and phi_dBarDD (covariant derivative)\n",
    "#          \\bar{D}_i \\bar{D}_j \\phi = \\phi_{;\\bar{i}\\bar{j}} = \\bar{D}_i \\phi_{,j}\n",
    "#                                   = \\phi_{,ij} - \\bar{\\Gamma}^k_{ij} \\phi_{,k}\n",
    "phi_dBarD = phi_dD\n",
    "phi_dBarDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        phi_dBarDD[i][j] = phi_dDD[i][j]\n",
    "        for k in range(DIM):\n",
    "            phi_dBarDD[i][j] += - GammabarUDD[k][i][j]*phi_dD[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 10: Code validation against `BSSN.BSSN_quantities` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "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_quantities](../edit/BSSN/BSSN_quantities.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": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "all_passed=True\n",
    "def comp_func(expr1,expr2,basename,prefixname2=\"Bq.\"):\n",
    "    if str(expr1-expr2)!=\"0\":\n",
    "        print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n",
    "        all_passed=False\n",
    "\n",
    "def gfnm(basename,idx1,idx2=None,idx3=None):\n",
    "    if idx2 is None:\n",
    "        return basename+\"[\"+str(idx1)+\"]\"\n",
    "    if idx3 is None:\n",
    "        return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"]\"\n",
    "    return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"][\"+str(idx3)+\"]\"\n",
    "\n",
    "expr_list = []\n",
    "exprcheck_list = []\n",
    "namecheck_list = []\n",
    "\n",
    "# Step 3:\n",
    "import BSSN.BSSN_quantities as Bq\n",
    "Bq.BSSN_basic_tensors()\n",
    "for i in range(DIM):\n",
    "    namecheck_list.extend([gfnm(\"LambdabarU\",i),gfnm(\"betaU\",i),gfnm(\"BU\",i)])\n",
    "    exprcheck_list.extend([Bq.LambdabarU[i],Bq.betaU[i],Bq.BU[i]])\n",
    "    expr_list.extend([LambdabarU[i],betaU[i],BU[i]])\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"gammabarDD\",i,j),gfnm(\"AbarDD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.gammabarDD[i][j],Bq.AbarDD[i][j]])\n",
    "        expr_list.extend([gammabarDD[i][j],AbarDD[i][j]])\n",
    "\n",
    "# Step 4:\n",
    "Bq.gammabar__inverse_and_derivs()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"gammabarUU\",i,j)])\n",
    "        exprcheck_list.extend([Bq.gammabarUU[i][j]])\n",
    "        expr_list.extend([gammabarUU[i][j]])\n",
    "        for k in range(DIM):\n",
    "            namecheck_list.extend([gfnm(\"gammabarDD_dD\",i,j,k),\n",
    "                                   gfnm(\"gammabarDD_dupD\",i,j,k),\n",
    "                                   gfnm(\"GammabarUDD\",i,j,k)])\n",
    "            exprcheck_list.extend([Bq.gammabarDD_dD[i][j][k],Bq.gammabarDD_dupD[i][j][k],Bq.GammabarUDD[i][j][k]])\n",
    "            expr_list.extend(        [gammabarDD_dD[i][j][k],gammabarDD_dupD[i][j][k],GammabarUDD[i][j][k]])\n",
    "\n",
    "# Step 5:\n",
    "Bq.detgammabar_and_derivs()\n",
    "namecheck_list.extend([\"detgammabar\"])\n",
    "exprcheck_list.extend([Bq.detgammabar])\n",
    "expr_list.extend([detgammabar])\n",
    "for i in range(DIM):\n",
    "    namecheck_list.extend([gfnm(\"detgammabar_dD\",i)])\n",
    "    exprcheck_list.extend([Bq.detgammabar_dD[i]])\n",
    "    expr_list.extend([detgammabar_dD[i]])\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"detgammabar_dDD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.detgammabar_dDD[i][j]])\n",
    "        expr_list.extend([detgammabar_dDD[i][j]])\n",
    "\n",
    "# Step 6:\n",
    "Bq.AbarUU_AbarUD_trAbar_AbarDD_dD()\n",
    "namecheck_list.extend([\"trAbar\"])\n",
    "exprcheck_list.extend([Bq.trAbar])\n",
    "expr_list.extend([trAbar])\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"AbarUU\",i,j),gfnm(\"AbarUD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.AbarUU[i][j],Bq.AbarUD[i][j]])\n",
    "        expr_list.extend([AbarUU[i][j],AbarUD[i][j]])\n",
    "        for k in range(DIM):\n",
    "            namecheck_list.extend([gfnm(\"AbarDD_dD\",i,j,k)])\n",
    "            exprcheck_list.extend([Bq.AbarDD_dD[i][j][k]])\n",
    "            expr_list.extend([AbarDD_dD[i][j][k]])\n",
    "\n",
    "# Step 7:\n",
    "Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "for i in range(DIM):\n",
    "    namecheck_list.extend([gfnm(\"DGammaU\",i)])\n",
    "    exprcheck_list.extend([Bq.DGammaU[i]])\n",
    "    expr_list.extend([DGammaU[i]])\n",
    "\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"RbarDD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.RbarDD[i][j]])\n",
    "        expr_list.extend([RbarDD[i][j]])\n",
    "        for k in range(DIM):\n",
    "            namecheck_list.extend([gfnm(\"DGammaUDD\",i,j,k),gfnm(\"gammabarDD_dHatD\",i,j,k)])\n",
    "            exprcheck_list.extend([Bq.DGammaUDD[i][j][k],Bq.gammabarDD_dHatD[i][j][k]])\n",
    "            expr_list.extend([DGammaUDD[i][j][k],gammabarDD_dHatD[i][j][k]])\n",
    "\n",
    "# Step 8:\n",
    "Bq.betaU_derivs()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"betaU_dD\",i,j),gfnm(\"betaU_dupD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.betaU_dD[i][j],Bq.betaU_dupD[i][j]])\n",
    "        expr_list.extend([betaU_dD[i][j],betaU_dupD[i][j]])\n",
    "        for k in range(DIM):\n",
    "            namecheck_list.extend([gfnm(\"betaU_dDD\",i,j,k)])\n",
    "            exprcheck_list.extend([Bq.betaU_dDD[i][j][k]])\n",
    "            expr_list.extend([betaU_dDD[i][j][k]])\n",
    "\n",
    "# Step 9:\n",
    "Bq.phi_and_derivs()\n",
    "#phi_dD,phi_dupD,phi_dDD,exp_m4phi,phi_dBarD,phi_dBarDD\n",
    "namecheck_list.extend([\"exp_m4phi\"])\n",
    "exprcheck_list.extend([Bq.exp_m4phi])\n",
    "expr_list.extend([exp_m4phi])\n",
    "for i in range(DIM):\n",
    "    namecheck_list.extend([gfnm(\"phi_dD\",i),gfnm(\"phi_dupD\",i),gfnm(\"phi_dBarD\",i)])\n",
    "    exprcheck_list.extend([Bq.phi_dD[i],Bq.phi_dupD[i],Bq.phi_dBarD[i]])\n",
    "    expr_list.extend(        [phi_dD[i],phi_dupD[i],phi_dBarD[i]])\n",
    "    for j in range(DIM):\n",
    "        namecheck_list.extend([gfnm(\"phi_dDD\",i,j),gfnm(\"phi_dBarDD\",i,j)])\n",
    "        exprcheck_list.extend([Bq.phi_dDD[i][j],Bq.phi_dBarDD[i][j]])\n",
    "        expr_list.extend([phi_dDD[i][j],phi_dBarDD[i][j]])\n",
    "\n",
    "for i in range(len(expr_list)):\n",
    "    comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])\n",
    "\n",
    "if all_passed:\n",
    "    print(\"ALL TESTS PASSED!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 11: 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_quantities.pdf](Tutorial-BSSN_quantities.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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_quantities.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-BSSN_quantities.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_quantities\")"
   ]
  }
 ],
 "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
}