Tutorial-GRFFE_Equations-Cartesian.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",
    "# Equations of General Relativistic Force-Free Electrodynamics (GRFFE)\n",
    "\n",
    "## Authors: Patrick Nelson, Zach Etienne, and Leo Werneck\n",
    "\n",
    "## This notebook documents and constructs a number of quantities useful for building symbolic (SymPy) expressions for the equations of general relativistic force-free electrodynamics (GRFFE), using the same (Valencia) formalism as `IllinoisGRMHD`. The formulation further incorporates the computation of flux and source terms for the induction equation, with detailed code validation against the existing GRFFE equations NRPy+ module.\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",
    "## Introduction\n",
    "\n",
    "We write the equations of general relativistic hydrodynamics in conservative form as follows (adapted from Eqs. 41-44 of [Duez et al](https://arxiv.org/pdf/astro-ph/0503420.pdf)):\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\partial_t \\tilde{S}_i &+& \\partial_j \\left(\\alpha \\sqrt{\\gamma} T^j_{{\\rm EM}i} \\right) = \\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu}_{\\rm EM} g_{\\mu\\nu,i},\n",
    "\\end{eqnarray}\n",
    "where we assume $T^{\\mu\\nu}_{\\rm EM}$ is the electromagnetic stress-energy tensor:\n",
    "$$\n",
    "T^{\\mu\\nu}_{\\rm EM} = b^2 u^{\\mu} u^{\\nu} + \\frac{b^2}{2} g^{\\mu\\nu} - b^\\mu b^\\nu,\n",
    "$$\n",
    "and \n",
    "$$\n",
    "v^j = \\frac{u^j}{u^0} \\\\\n",
    "$$\n",
    "\n",
    "Some quantities can be defined in precisely the same way they are defined in the GRHD equations. ***Therefore, we will not define special functions for generating these quantities, and instead refer the user to the appropriate functions in the [GRHD module](../edit/GRHD/equations.py)*** Namely,\n",
    "\n",
    "* The GRFFE conservative variables:\n",
    "    * $\\tilde{S}_i  = \\alpha \\sqrt{\\gamma} T^0{}_i$, via `GRHD.compute_S_tildeD(alpha, sqrtgammaDET, T4UD)`\n",
    "* The GRFFE fluxes:\n",
    "    * $\\tilde{S}_i$ flux: $\\left(\\alpha \\sqrt{\\gamma} T^j{}_i \\right)$, via `GRHD.compute_S_tilde_fluxUD(alpha, sqrtgammaDET, T4UD)`\n",
    "* GRFFE source terms:\n",
    "    * $\\tilde{S}_i$ source term: $\\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu} g_{\\mu\\nu,i}$, via `GRHD.compute_S_tilde_source_termD(alpha, sqrtgammaDET,g4DD_zerotimederiv_dD, T4UU)`\n",
    "\n",
    "Also, we will write the 4-metric in terms of the ADM 3-metric, lapse, and shift using standard equations.\n",
    "\n",
    "Thus the full set of input variables includes:\n",
    "* Spacetime quantities:\n",
    "    * ADM quantities $\\alpha$, $\\beta^i$, $\\gamma_{ij}$\n",
    "* Hydrodynamical quantities:\n",
    "    * 4-velocity $u^\\mu$\n",
    "* Electrodynamical quantities\n",
    "    * Magnetic field $B^i$\n",
    "\n",
    "Additionally, we need to evolve the vector $A_i$ according to\n",
    "$$\n",
    "\\partial_t A_i = \\epsilon_{ijk} v^j B^k - \\partial_i (\\alpha \\Phi - \\beta^j A_j)\n",
    "$$\n",
    "and the scalar potential $\\Phi$ according to\n",
    "$$\n",
    "\\partial_t [\\sqrt{\\gamma} \\Phi] + \\partial_j (\\alpha\\sqrt{\\gamma}A^j - \\beta^j [\\sqrt{\\gamma} \\Phi]) = - \\xi \\alpha [\\sqrt{\\gamma} \\Phi],\n",
    "$$\n",
    "where $\\xi$ determines the strength of the damping term. This is typically set to $0.1$.\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",
    "For instance, in calculating the first term of $b^2 u^\\mu u^\\nu$, we use Greek indices:\n",
    "\n",
    "```python\n",
    "T4EMUU = ixp.zerorank2(DIM=4)\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        # Term 1: b^2 u^{\\mu} u^{\\nu}\n",
    "        T4EMUU[mu][nu] = smallb2*u4U[mu]*u4U[nu]\n",
    "```\n",
    "\n",
    "When we calculate $\\beta_i = \\gamma_{ij} \\beta^j$, we use Latin indices:\n",
    "```python\n",
    "betaD = ixp.zerorank1(DIM=3)\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        betaD[i] += gammaDD[i][j] * betaU[j]\n",
    "```\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). This can be seen when we handle $\\frac{1}{2} \\alpha \\sqrt{\\gamma} T^{\\mu \\nu}_{\\rm EM} \\partial_i g_{\\mu \\nu}$:\n",
    "```python\n",
    "# \\alpha \\sqrt{\\gamma} T^{\\mu \\nu}_{\\rm EM} \\partial_i g_{\\mu \\nu} / 2\n",
    "for i in range(3):\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            S_tilde_rhsD[i] += alpsqrtgam * T4EMUU[mu][nu] * g4DD_zerotimederiv_dD[mu][nu][i+1] / 2\n",
    "```"
   ]
  },
  {
   "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](#importmodules): Import needed NRPy+ & Python modules\n",
    "1. [Step 2](#u0bu): Define needed quantities for $T^{\\mu\\nu}_{\\rm EM}$, the EM part of the stress-energy tensor\n",
    "1. [Step 3](#stressenergy): **compute_TEM4UU()**, **compute_TEM4UD()**: Define the stress-energy tensor $T^{\\mu\\nu}_{\\rm EM}$ and $T^\\mu_{{\\rm EM}\\nu}$\n",
    "1. [Step 4](#inductioneq): Vector potential induction equation, assuming generalized Lorenz gauge\n",
    "    1. [Step 4.a](#inductionterms) Compute the flux term and the source term for the induction equation\n",
    "        1. **compute_AD_flux_term()**,**compute_AD_source_term_operand_for_FD()**\n",
    "    1. [Step 4.b](#gaugeeq) Compute the damping term and flux term for the gauge equation\n",
    "        1. **compute_psi6Phi_rhs_flux_term_operand()**,**compute_psi6Phi_rhs_damping_term()**\n",
    "1. [Step 5](#declarevarsconstructgrffeeqs): Declare ADM and hydrodynamical input variables, and construct GRFFE equations\n",
    "1. [Step 6](#code_validation): Code Validation against `GRFFE.equations` NRPy+ module for GRMHD\n",
    "1. [Step 7](#code_validation_2): Code Validation against `GRFFE.equations` NRPy+ module for GRFFE\n",
    "1. [Step 8](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='importmodules'></a>\n",
    "\n",
    "# Step 1: Import needed NRPy+ & Python modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{importmodules}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:15.803911Z",
     "iopub.status.busy": "2021-03-07T17:11:15.802519Z",
     "iopub.status.idle": "2021-03-07T17:11:16.139194Z",
     "shell.execute_reply": "2021-03-07T17:11:16.138509Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Import needed core NRPy+ modules\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='u0bu'></a>\n",
    "\n",
    "# Step 2: Define needed quantities for $T^{\\mu\\nu}_{\\rm EM}$, the EM part of the stress-energy tensor \\[Back to [top](#toc)\\]\n",
    "$$\\label{u0bu}$$\n",
    "\n",
    "We are given $B^i$, the magnetic field as measured by a *normal* observer, yet $T^{\\mu\\nu}_{\\rm EM}$ depends on $b^{\\mu}$, the magnetic field as measured by an observer comoving with the plasma $B^{\\mu}_{\\rm (u)}$, divided by $\\sqrt{4\\pi}$.\n",
    "\n",
    "In the ideal MHD limit, $B^{\\mu}_{\\rm (u)}$ is orthogonal to the plasma 4-velocity $u^\\mu$, which sets the $\\mu=0$ component. \n",
    "\n",
    "$B^{\\mu}_{\\rm (u)}$ is related to the magnetic field as measured by a *normal* observer $B^i$ via a simple projection (Eq 21 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)), which results in the expressions (Eqs 23 and 24 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf)):\n",
    "\n",
    "\\begin{align}\n",
    "\\sqrt{4\\pi} b^0 = B^0_{\\rm (u)} &= \\frac{u_j B^j}{\\alpha} \\\\\n",
    "\\sqrt{4\\pi} b^i = B^i_{\\rm (u)} &= \\frac{B^i + (u_j B^j) u^i}{\\alpha u^0}\\\\\n",
    "\\end{align}\n",
    "\n",
    "Further, $B^i$ is related to the actual magnetic field evaluated in IllinoisGRMHD, $\\tilde{B}^i$ via\n",
    "\n",
    "$$B^i = \\frac{\\tilde{B}^i}{\\gamma},$$\n",
    "\n",
    "where $\\gamma$ is the determinant of the spatial 3-metric:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.144888Z",
     "iopub.status.busy": "2021-03-07T17:11:16.143936Z",
     "iopub.status.idle": "2021-03-07T17:11:16.146316Z",
     "shell.execute_reply": "2021-03-07T17:11:16.146814Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a: Define B^i = Btilde^i / sqrt(gamma)\n",
    "def compute_B_notildeU(sqrtgammaDET, B_tildeU):\n",
    "    global B_notildeU\n",
    "    B_notildeU = ixp.zerorank1(DIM=3)\n",
    "    for i in range(3):\n",
    "        B_notildeU[i] = B_tildeU[i]/sqrtgammaDET"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we compute Eqs 23 and 24 in [Duez *et al* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf):\n",
    "\n",
    "\\begin{align}\n",
    "\\sqrt{4\\pi} b^0 = B^0_{\\rm (u)} &= \\frac{u_j B^j}{\\alpha} \\\\\n",
    "\\sqrt{4\\pi} b^i = B^i_{\\rm (u)} &= \\frac{B^i + (u_j B^j) u^i}{\\alpha u^0}.\n",
    "\\end{align}\n",
    "\n",
    "In doing so, we will store the scalar $u_j B^j$ to `u4_dot_B_notilde`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.156012Z",
     "iopub.status.busy": "2021-03-07T17:11:16.155152Z",
     "iopub.status.idle": "2021-03-07T17:11:16.158780Z",
     "shell.execute_reply": "2021-03-07T17:11:16.158130Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.b.i: Define b^mu.\n",
    "def compute_smallb4U(gammaDD,betaU,alpha, u4U,B_notildeU, sqrt4pi):\n",
    "    global smallb4U\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "\n",
    "    u4D = ixp.zerorank1(DIM=4)\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            u4D[mu] += AB4m.g4DD[mu][nu]*u4U[nu]\n",
    "    smallb4U = ixp.zerorank1(DIM=4)\n",
    "    u4_dot_B_notilde = sp.sympify(0)\n",
    "    for i in range(3):\n",
    "        u4_dot_B_notilde += u4D[i+1]*B_notildeU[i]\n",
    "\n",
    "    # b^0 = (u_j B^j)/[alpha * sqrt(4 pi)]\n",
    "    smallb4U[0] = u4_dot_B_notilde / (alpha*sqrt4pi)\n",
    "    # b^i = [B^i + (u_j B^j) u^i]/[alpha * u^0 * sqrt(4 pi)]\n",
    "    for i in range(3):\n",
    "        smallb4U[i+1] = (B_notildeU[i] + u4_dot_B_notilde*u4U[i+1]) / (alpha*u4U[0]*sqrt4pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, in the case of pure GRFFE (that is, if we assume the pressure and density are zero), we can make an additional simplifying assumption. When deriving the equations of GRFFE, one has a choice as to how to define the velocity. Following the example set by [McKinney (2006)](https://arxiv.org/pdf/astro-ph/0601410.pdf), section 2.1, and [Paschalidis (2013)](https://arxiv.org/pdf/1310.3274v2.pdf), Appendix A, we choose a referencing frame comoving with the plasma that fulfills the ideal MHD condition (i.e., the electric field is zero). However, in GRFFE, this choice yields a one-parameter family of timelike vectors. By choosing the drift velocity $v^i = u^i/u^0$ that minimizes the Lorentz factor, we find that the four-velocity is orthogonal to the magnetic field, or $u_j B^j = 0$. With this assumption, \n",
    "\n",
    "\\begin{align}\n",
    "\\sqrt{4\\pi} b^0 &= 0 \\\\\n",
    "\\sqrt{4\\pi} b^i &= \\frac{B^i}{\\alpha u^0}.\n",
    "\\end{align}\n",
    "\n",
    "This simplification also gives the inversion from $\\tilde{S}_i$ to the drift velocity $v^i$ a closed form, greatly simplifying the conservative-to-primitive solver and removing the need to "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.168528Z",
     "iopub.status.busy": "2021-03-07T17:11:16.167373Z",
     "iopub.status.idle": "2021-03-07T17:11:16.170337Z",
     "shell.execute_reply": "2021-03-07T17:11:16.170953Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.b.ii: Define b^mu when u4 and B are orthogonal\n",
    "def compute_smallb4U_with_driftvU_for_FFE(gammaDD,betaU,alpha, u4U,B_notildeU, sqrt4pi):\n",
    "    global smallb4_with_driftv_for_FFE_U\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "\n",
    "    u4D = ixp.zerorank1(DIM=4)\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            u4D[mu] += AB4m.g4DD[mu][nu]*u4U[nu]\n",
    "    smallb4_with_driftv_for_FFE_U = ixp.zerorank1(DIM=4)\n",
    "\n",
    "    # b^0 = 0\n",
    "    smallb4_with_driftv_for_FFE_U[0] = 0\n",
    "    # b^i = B^i / [alpha * u^0 * sqrt(4 pi)]\n",
    "    for i in range(3):\n",
    "        smallb4_with_driftv_for_FFE_U[i+1] = B_notildeU[i] / (alpha*u4U[0]*sqrt4pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we compute `smallbsquared`=$b^2 = b_{\\mu} b^{\\mu} = g_{\\mu \\nu} b^{\\nu}b^{\\mu}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.178212Z",
     "iopub.status.busy": "2021-03-07T17:11:16.177422Z",
     "iopub.status.idle": "2021-03-07T17:11:16.179425Z",
     "shell.execute_reply": "2021-03-07T17:11:16.179911Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.c: Define b^2.\n",
    "def compute_smallbsquared(gammaDD,betaU,alpha, smallb4U):\n",
    "    global smallbsquared\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "\n",
    "    smallbsquared = sp.sympify(0)\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            smallbsquared += AB4m.g4DD[mu][nu]*smallb4U[mu]*smallb4U[nu]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='stressenergy'></a>\n",
    "\n",
    "# Step 3: Define the electromagnetic stress-energy tensor $T^{\\mu\\nu}_{\\rm EM}$ and $T^\\mu_{{\\rm EM}\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{stressenergy}$$\n",
    "\n",
    "Recall from above that\n",
    "\n",
    "$$\n",
    "T^{\\mu\\nu}_{\\rm EM} = b^2 u^{\\mu} u^{\\nu} + \\frac{1}{2} b^2 g^{\\mu\\nu} - b^\\mu b^\\nu.\n",
    "$$\n",
    "Also \n",
    "\n",
    "$$\n",
    "T^\\mu_{{\\rm EM}\\nu} = T^{\\mu\\delta}_{\\rm EM} g_{\\delta \\nu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.189425Z",
     "iopub.status.busy": "2021-03-07T17:11:16.188381Z",
     "iopub.status.idle": "2021-03-07T17:11:16.190879Z",
     "shell.execute_reply": "2021-03-07T17:11:16.191391Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.a: Define T_{EM}^{mu nu} (a 4-dimensional tensor)\n",
    "def compute_TEM4UU(gammaDD,betaU,alpha, smallb4U, smallbsquared,u4U):\n",
    "    global TEM4UU\n",
    "\n",
    "    # Then define g^{mu nu} in terms of the ADM quantities:\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4UU_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "\n",
    "    # Finally compute T^{mu nu}\n",
    "    TEM4UU = ixp.zerorank2(DIM=4)\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            TEM4UU[mu][nu] = smallbsquared*u4U[mu]*u4U[nu] \\\n",
    "                             + sp.Rational(1,2)*smallbsquared*AB4m.g4UU[mu][nu] \\\n",
    "                             - smallb4U[mu]*smallb4U[nu]\n",
    "\n",
    "# Step 3.b: Define T^{mu}_{nu} (a 4-dimensional tensor)\n",
    "def compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU):\n",
    "    global TEM4UD\n",
    "    # Next compute T^mu_nu = T^{mu delta} g_{delta nu}, needed for S_tilde flux.\n",
    "    # First we'll need g_{alpha nu} in terms of ADM quantities:\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",gammaDD,betaU,alpha)\n",
    "    TEM4UD = ixp.zerorank2(DIM=4)\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            for delta in range(4):\n",
    "                TEM4UD[mu][nu] += TEM4UU[mu][delta]*AB4m.g4DD[delta][nu]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='inductioneq'></a>\n",
    "\n",
    "# Step 4: Vector potential induction equation, assuming generalized Lorenz gauge \\[Back to [top](#toc)\\]\n",
    "$$\\label{inductioneq}$$\n",
    "\n",
    "Now, we will turn our attention to the induction equation, which evolves $A_i$, as well as an additional electromagnetic gauge evolution equation, which evolves $\\Phi$. For a cell-centered formulation, they are as follows:\n",
    "$$\n",
    "\\partial_t A_i = \\epsilon_{ijk} v^j B^k - \\partial_i (\\alpha \\Phi - \\beta^j A_j),\n",
    "$$\n",
    "from Eq. 17 of the [original `GiRaFFE` paper](https://arxiv.org/pdf/1704.00599.pdf), and \n",
    "$$\n",
    "\\partial_t [\\sqrt{\\gamma} \\Phi] + \\partial_j (\\alpha\\sqrt{\\gamma}A^j - \\beta^j [\\sqrt{\\gamma} \\Phi]) = - \\xi \\alpha [\\sqrt{\\gamma} \\Phi],\n",
    "$$\n",
    "from Eq. 19. When it comes to taking the derivatives in these equations, we will write these functions with the intention of computing the operand first, storing it as a gridfunction, and then finite differencing that in a later step. \n",
    "\n",
    "<a id='inductionterms'></a>\n",
    "\n",
    "## Step 4.a: Compute the flux term and the source term for the induction equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{inductionterms}$$\n",
    "\n",
    "We'll now take a closer look at the induction equation, the right-hand side of which has two terms:\n",
    "$$\n",
    "\\partial_t A_i = \\underbrace{\\epsilon_{ijk} v^j B^k}_{\\rm \"Flux\"\\ term} - \\partial_i (\\underbrace{\\alpha \\Phi - \\beta^j A_j}_{\\rm \"Source\"\\ term}),\n",
    "$$\n",
    "The flux term here is a simple cross-product between the drift velocity and the magnetic field. The source term is the gradient of a gauge-dependent combination of the lapse $\\alpha$, the shift $\\beta^i$, the vector potential $A_i$, and the scalar potential $\\Phi$; we will write a function to compute that operand, and save the finite-difference derivative until later. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.199410Z",
     "iopub.status.busy": "2021-03-07T17:11:16.198186Z",
     "iopub.status.idle": "2021-03-07T17:11:16.200914Z",
     "shell.execute_reply": "2021-03-07T17:11:16.201415Z"
    }
   },
   "outputs": [],
   "source": [
    "def compute_AD_flux_term(sqrtgammaDET,driftvU,BU):\n",
    "    # Levi-Civita tensor for cross products\n",
    "    LeviCivitaDDD = ixp.LeviCivitaTensorDDD_dim3_rank3(sqrtgammaDET)\n",
    "    global A_fluxD\n",
    "    A_fluxD = ixp.zerorank1()\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                # \\epsilon_{ijk} v^j B^k\n",
    "                A_fluxD[i] += LeviCivitaDDD[i][j][k]*driftvU[j]*BU[k]\n",
    "\n",
    "def compute_AD_source_term_operand_for_FD(sqrtgammaDET,betaU,alpha,psi6Phi,AD):\n",
    "    Phi = psi6Phi/sqrtgammaDET\n",
    "    global AevolParen\n",
    "    # \\alpha \\Phi - \\beta^j A_j\n",
    "    AevolParen = alpha * Phi\n",
    "    for j in range(3):\n",
    "        AevolParen += -betaU[j] * AD[j]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='gaugeeq'></a>\n",
    "\n",
    "## Step 4.b: Compute the damping term and flux term for the gauge equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{gaugeeq}$$\n",
    "\n",
    "Next, we will build the expressions for the RHS of the evolution equation \n",
    "$$\n",
    "\\partial_t [\\sqrt{\\gamma} \\Phi] = -\\partial_j (\\underbrace{\\alpha\\sqrt{\\gamma}A^j - \\beta^j [\\sqrt{\\gamma} \\Phi]}_{\\rm Gauge\\ evolution\\ term}) - \\underbrace{\\xi \\alpha [\\sqrt{\\gamma} \\Phi]}_{\\rm Damping\\ term}.\n",
    "$$\n",
    "Once again, we will do this in two parts. First, we will compute the operand of the divergence in the flux term, leaving the finite-difference derivative for later; then, we will compute the damping term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.209160Z",
     "iopub.status.busy": "2021-03-07T17:11:16.208181Z",
     "iopub.status.idle": "2021-03-07T17:11:16.210948Z",
     "shell.execute_reply": "2021-03-07T17:11:16.210409Z"
    }
   },
   "outputs": [],
   "source": [
    "def compute_psi6Phi_rhs_flux_term_operand(gammaDD,sqrtgammaDET,betaU,alpha,AD,psi6Phi):\n",
    "    gammaUU,_gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "    AU = ixp.zerorank1()\n",
    "    # Raise the index on A in the usual way:\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            AU[i] += gammaUU[i][j] * AD[j]\n",
    "\n",
    "    global PhievolParenU\n",
    "    PhievolParenU = ixp.zerorank1(DIM=3)\n",
    "\n",
    "    for j in range(3):\n",
    "        # \\alpha\\sqrt{\\gamma}A^j - \\beta^j [\\sqrt{\\gamma} \\Phi]\n",
    "        PhievolParenU[j] += alpha*sqrtgammaDET*AU[j] - betaU[j]*psi6Phi\n",
    "\n",
    "def compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping):\n",
    "    # - \\xi \\alpha [\\sqrt{\\gamma} \\Phi]\n",
    "    # Combine the divergence and the damping term\n",
    "    global psi6Phi_damping\n",
    "    psi6Phi_damping = - xi_damping * alpha * psi6Phi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='declarevarsconstructgrffeeqs'></a>\n",
    "\n",
    "# Step 5: Declare ADM and hydrodynamical input variables, and construct GRFFE equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{declarevarsconstructgrffeeqs}$$\n",
    "\n",
    "The GRFFE equations are given by the induction equation (handled later) and the evolution equation for $\\tilde{S}_i$, the Poynting one-form:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\partial_t \\tilde{S}_i &+& \\partial_j \\left(\\alpha \\sqrt{\\gamma} T^j_{{\\rm EM}i} \\right) = \\frac{1}{2} \\alpha\\sqrt{\\gamma} T^{\\mu\\nu}_{\\rm EM} g_{\\mu\\nu,i}.\n",
    "\\end{eqnarray}\n",
    "\n",
    "Notice that all terms in this equation ($\\tilde{S}_i$, $\\left(\\alpha \\sqrt{\\gamma} T^j_{{\\rm EM}i} \\right)$, and the source term) are *identical* to those in the evolution equation for $\\tilde{S}_i$ in [general relativistic hydrodynamics (GRHD)](Tutorial-GRHD_Equations-Cartesian.ipynb); one need only replace $T^{\\mu\\nu}$ of GRHD with the $T^{\\mu\\nu}_{\\rm EM}$ defined above. \n",
    "\n",
    "Thus we will reuse expressions from the [general relativistic hydrodynamics (GRHD)](Tutorial-GRHD_Equations-Cartesian.ipynb) module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.224614Z",
     "iopub.status.busy": "2021-03-07T17:11:16.223898Z",
     "iopub.status.idle": "2021-03-07T17:11:16.477500Z",
     "shell.execute_reply": "2021-03-07T17:11:16.476944Z"
    }
   },
   "outputs": [],
   "source": [
    "# First define hydrodynamical quantities\n",
    "u4U = ixp.declarerank1(\"u4U\", DIM=4)\n",
    "B_tildeU = ixp.declarerank1(\"B_tildeU\", DIM=3)\n",
    "AD       = ixp.declarerank1(\"D\",        DIM=3)\n",
    "B_tildeU = ixp.declarerank1(\"B_tildeU\", DIM=3)\n",
    "psi6Phi = sp.symbols('psi6Phi', real=True)\n",
    "\n",
    "# Then ADM quantities\n",
    "gammaDD = ixp.declarerank2(\"gammaDD\",\"sym01\",DIM=3)\n",
    "betaU   = ixp.declarerank1(\"betaU\", DIM=3)\n",
    "alpha   = sp.symbols('alpha', real=True)\n",
    "\n",
    "# Then numerical constant\n",
    "sqrt4pi = sp.symbols('sqrt4pi', real=True)\n",
    "xi_damping = sp.symbols('xi_damping', real=True)\n",
    "\n",
    "# First compute stress-energy tensor T4UU and T4UD:\n",
    "import GRHD.equations as GRHD\n",
    "GRHD.compute_sqrtgammaDET(gammaDD)\n",
    "compute_B_notildeU(GRHD.sqrtgammaDET, B_tildeU)\n",
    "compute_smallb4U(gammaDD,betaU,alpha, u4U, B_notildeU, sqrt4pi)\n",
    "compute_smallbsquared(gammaDD,betaU,alpha, smallb4U)\n",
    "\n",
    "compute_TEM4UU(gammaDD,betaU,alpha, smallb4U, smallbsquared,u4U)\n",
    "compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU)\n",
    "\n",
    "# Compute conservative variables in terms of primitive variables\n",
    "GRHD.compute_S_tildeD( alpha, GRHD.sqrtgammaDET, TEM4UD)\n",
    "S_tildeD = GRHD.S_tildeD\n",
    "\n",
    "# Next compute fluxes of conservative variables\n",
    "GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET,    TEM4UD)\n",
    "S_tilde_fluxUD = GRHD.S_tilde_fluxUD\n",
    "\n",
    "# Then declare derivatives & compute g4DDdD\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
    "betaU_dD   = ixp.declarerank2(\"betaU_dD\"  ,\"nosym\",DIM=3)\n",
    "alpha_dD   = ixp.declarerank1(\"alpha_dD\"          ,DIM=3)\n",
    "GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
    "\n",
    "# Finally compute source terms on tau_tilde and S_tilde equations\n",
    "GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET, GRHD.g4DD_zerotimederiv_dD, TEM4UU)\n",
    "S_tilde_source_termD = GRHD.S_tilde_source_termD\n",
    "\n",
    "# We must also compute the terms for the induction equation and gauge evolution equation\n",
    "GRHD.compute_vU_from_u4U__no_speed_limit(u4U) # We need the drift velocity\n",
    "compute_AD_flux_term(GRHD.sqrtgammaDET,GRHD.vU,B_notildeU)\n",
    "compute_AD_source_term_operand_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)\n",
    "compute_psi6Phi_rhs_flux_term_operand(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)\n",
    "compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 6: Code Validation against `GRFFE.equations` NRPy+ module for GRMHD \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "As a code validation check, we verify agreement in the SymPy expressions for the GRFFE equations generated in\n",
    "1. this tutorial notebook versus\n",
    "2. the NRPy+ [GRFFE.equations](../edit/GRFFE/equations.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.486832Z",
     "iopub.status.busy": "2021-03-07T17:11:16.486142Z",
     "iopub.status.idle": "2021-03-07T17:11:16.494167Z",
     "shell.execute_reply": "2021-03-07T17:11:16.494692Z"
    }
   },
   "outputs": [],
   "source": [
    "import GRFFE.equations as GRFFE\n",
    "\n",
    "# First compute B^i from Btilde^i:\n",
    "GRFFE.compute_B_notildeU(GRHD.sqrtgammaDET, B_tildeU)\n",
    "\n",
    "# Then compute b^mu and b^2:\n",
    "GRFFE.compute_smallb4U(gammaDD, betaU, alpha, u4U, GRFFE.B_notildeU, sqrt4pi)\n",
    "GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)\n",
    "\n",
    "# Next construct stress-energy tensor T4UU and T4UD:\n",
    "GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,u4U)\n",
    "GRFFE.compute_TEM4UD(gammaDD,betaU,alpha, GRFFE.TEM4UU)\n",
    "\n",
    "# Compute conservative variables in terms of primitive variables\n",
    "GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)\n",
    "Ge_S_tildeD = GRHD.S_tildeD\n",
    "\n",
    "# Next compute fluxes of conservative variables\n",
    "GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)\n",
    "Ge_S_tilde_fluxUD = GRHD.S_tilde_fluxUD\n",
    "\n",
    "# Then declare derivatives & compute g4DDdD\n",
    "# gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
    "# betaU_dD   = ixp.declarerank2(\"betaU_dD\"  ,\"nosym\",DIM=3)\n",
    "# alpha_dD   = ixp.declarerank1(\"alpha_dD\"          ,DIM=3)\n",
    "GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
    "\n",
    "# Finally compute source terms on S_tilde equations\n",
    "GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD,GRFFE.TEM4UU)\n",
    "Ge_S_tilde_source_termD = GRHD.S_tilde_source_termD\n",
    "\n",
    "# We must also compute the terms for the induction equation and gauge evolution equation\n",
    "# GRHD.compute_vU_from_u4U__no_speed_limit(u4U) # We need the drift velocity\n",
    "GRFFE.compute_AD_flux_term(GRHD.sqrtgammaDET,GRHD.vU,B_notildeU)\n",
    "GRFFE.compute_AD_source_term_operand_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)\n",
    "GRFFE.compute_psi6Phi_rhs_flux_term_operand(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)\n",
    "GRFFE.compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.513062Z",
     "iopub.status.busy": "2021-03-07T17:11:16.502671Z",
     "iopub.status.idle": "2021-03-07T17:11:16.527678Z",
     "shell.execute_reply": "2021-03-07T17:11:16.528190Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "def comp_func(expr1,expr2,basename,prefixname2=\"GRFFE.\"):\n",
    "    if str(expr1-expr2)!=\"0\":\n",
    "        print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n",
    "        return 1\n",
    "    return 0\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",
    "# PhievolParenU\n",
    "# psi6Phi_damping\n",
    "# A_fluxD\n",
    "# AevolParen\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        namecheck_list.extend([gfnm(\"TEM4UU\",mu,nu),gfnm(\"TEM4UD\",mu,nu)])\n",
    "        exprcheck_list.extend([GRFFE.TEM4UU[mu][nu],GRFFE.TEM4UD[mu][nu]])\n",
    "        expr_list.extend([TEM4UU[mu][nu],TEM4UD[mu][nu]])\n",
    "\n",
    "for i in range(3):\n",
    "    namecheck_list.extend([gfnm(\"S_tildeD\",i),gfnm(\"S_tilde_source_termD\",i),gfnm(\"A_fluxD\",i),gfnm(\"PhievolParenU\",i)])\n",
    "    exprcheck_list.extend([Ge_S_tildeD[i],Ge_S_tilde_source_termD[i],GRFFE.A_fluxD[i],GRFFE.PhievolParenU[i]])\n",
    "    expr_list.extend([S_tildeD[i],S_tilde_source_termD[i],A_fluxD[i],PhievolParenU[i]])\n",
    "    for j in range(3):\n",
    "        namecheck_list.extend([gfnm(\"S_tilde_fluxUD\",i,j)])\n",
    "        exprcheck_list.extend([Ge_S_tilde_fluxUD[i][j]])\n",
    "        expr_list.extend([S_tilde_fluxUD[i][j]])\n",
    "\n",
    "namecheck_list.extend([\"AevolParen\",\"psi6Phi_damping\"])\n",
    "exprcheck_list.extend([GRFFE.AevolParen,GRFFE.psi6Phi_damping])\n",
    "expr_list.extend([AevolParen,psi6Phi_damping])\n",
    "\n",
    "num_failures = 0\n",
    "for i in range(len(expr_list)):\n",
    "    num_failures += comp_func(expr_list[i],exprcheck_list[i],namecheck_list[i])\n",
    "\n",
    "import sys\n",
    "if num_failures == 0:\n",
    "    print(\"ALL TESTS PASSED!\")\n",
    "else:\n",
    "    print(\"ERROR: \"+str(num_failures)+\" TESTS DID NOT PASS\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation_2'></a>\n",
    "\n",
    "# Step 7: Code Validation against `GRFFE.equations` NRPy+ module for GRFFE \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation_2}$$\n",
    "\n",
    "Additionally, we verify agreement in the SymPy expressions for the GRFFE equations generated in\n",
    "1. this tutorial notebook versus\n",
    "2. the NRPy+ [GRFFE.equations](../edit/GRFFE/equations.py) module\n",
    "for the case $u_j B^j = 0$. \n",
    "\n",
    "We will only recompute $b^\\mu$ and the expressions that depend on it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.604669Z",
     "iopub.status.busy": "2021-03-07T17:11:16.568713Z",
     "iopub.status.idle": "2021-03-07T17:11:16.618441Z",
     "shell.execute_reply": "2021-03-07T17:11:16.618924Z"
    }
   },
   "outputs": [],
   "source": [
    "# Generate the expressions within the tutorial, starting with:\n",
    "compute_smallb4U_with_driftvU_for_FFE(gammaDD,betaU,alpha, u4U, B_notildeU, sqrt4pi)\n",
    "compute_smallbsquared(gammaDD,betaU,alpha, smallb4_with_driftv_for_FFE_U)\n",
    "\n",
    "compute_TEM4UU(gammaDD,betaU,alpha, smallb4_with_driftv_for_FFE_U, smallbsquared,u4U)\n",
    "compute_TEM4UD(gammaDD,betaU,alpha, TEM4UU)\n",
    "\n",
    "# Compute conservative variables in terms of primitive variables\n",
    "GRHD.compute_S_tildeD( alpha, GRHD.sqrtgammaDET, TEM4UD)\n",
    "S_tildeD = GRHD.S_tildeD\n",
    "\n",
    "# Next compute fluxes of conservative variables\n",
    "GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET,    TEM4UD)\n",
    "S_tilde_fluxUD = GRHD.S_tilde_fluxUD\n",
    "\n",
    "# Then declare derivatives & compute g4DDdD\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
    "betaU_dD   = ixp.declarerank2(\"betaU_dD\"  ,\"nosym\",DIM=3)\n",
    "alpha_dD   = ixp.declarerank1(\"alpha_dD\"          ,DIM=3)\n",
    "GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
    "\n",
    "# Finally compute source terms on tau_tilde and S_tilde equations\n",
    "GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET, GRHD.g4DD_zerotimederiv_dD, TEM4UU)\n",
    "S_tilde_source_termD = GRHD.S_tilde_source_termD\n",
    "\n",
    "# Now compute the expressions from the module\n",
    "# Then compute b^mu and b^2:\n",
    "GRFFE.compute_smallb4U_with_driftvU_for_FFE(gammaDD, betaU, alpha, u4U, GRFFE.B_notildeU, sqrt4pi)\n",
    "GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4_with_driftv_for_FFE_U)\n",
    "\n",
    "# Next construct stress-energy tensor T4UU and T4UD:\n",
    "GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4_with_driftv_for_FFE_U, GRFFE.smallbsquared,u4U)\n",
    "GRFFE.compute_TEM4UD(gammaDD,betaU,alpha, GRFFE.TEM4UU)\n",
    "\n",
    "# Compute conservative variables in terms of primitive variables\n",
    "GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)\n",
    "Ge_S_tildeD = GRHD.S_tildeD\n",
    "\n",
    "# Next compute fluxes of conservative variables\n",
    "GRHD.compute_S_tilde_fluxUD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)\n",
    "Ge_S_tilde_fluxUD = GRHD.S_tilde_fluxUD\n",
    "\n",
    "# Then declare derivatives & compute g4DDdD\n",
    "# gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\",DIM=3)\n",
    "# betaU_dD   = ixp.declarerank2(\"betaU_dD\"  ,\"nosym\",DIM=3)\n",
    "# alpha_dD   = ixp.declarerank1(\"alpha_dD\"          ,DIM=3)\n",
    "GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)\n",
    "\n",
    "# Finally compute source terms on S_tilde equations\n",
    "GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD,GRFFE.TEM4UU)\n",
    "Ge_S_tilde_source_termD = GRHD.S_tilde_source_termD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.638877Z",
     "iopub.status.busy": "2021-03-07T17:11:16.638140Z",
     "iopub.status.idle": "2021-03-07T17:11:16.641863Z",
     "shell.execute_reply": "2021-03-07T17:11:16.641318Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "all_passed=True\n",
    "\n",
    "expr_list = []\n",
    "exprcheck_list = []\n",
    "namecheck_list = []\n",
    "\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        namecheck_list.extend([gfnm(\"TEM4UU\",mu,nu),gfnm(\"TEM4UD\",mu,nu)])\n",
    "        exprcheck_list.extend([GRFFE.TEM4UU[mu][nu],GRFFE.TEM4UD[mu][nu]])\n",
    "        expr_list.extend([TEM4UU[mu][nu],TEM4UD[mu][nu]])\n",
    "\n",
    "for i in range(3):\n",
    "    namecheck_list.extend([gfnm(\"S_tildeD\",i),gfnm(\"S_tilde_source_termD\",i)])\n",
    "    exprcheck_list.extend([Ge_S_tildeD[i],Ge_S_tilde_source_termD[i]])\n",
    "    expr_list.extend([S_tildeD[i],S_tilde_source_termD[i]])\n",
    "    for j in range(3):\n",
    "        namecheck_list.extend([gfnm(\"S_tilde_fluxUD\",i,j)])\n",
    "        exprcheck_list.extend([Ge_S_tilde_fluxUD[i][j]])\n",
    "        expr_list.extend([S_tilde_fluxUD[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",
    "import sys\n",
    "if all_passed:\n",
    "    print(\"ALL TESTS PASSED!\")\n",
    "else:\n",
    "    print(\"ERROR: AT LEAST ONE TEST DID NOT PASS\")\n",
    "    sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 8: 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-GRFFE_Equations-Cartesian.pdf](Tutorial-GRFFE_Equations-Cartesian.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": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:11:16.647125Z",
     "iopub.status.busy": "2021-03-07T17:11:16.645938Z",
     "iopub.status.idle": "2021-03-07T17:11:20.521178Z",
     "shell.execute_reply": "2021-03-07T17:11:20.520388Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-GRFFE_Equations-Cartesian.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-GRFFE_Equations-Cartesian.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GRFFE_Equations-Cartesian\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}