Tutorial-Psi4_tetrads.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",
    "# Tetrads for Evaluating the Outgoing Gravitational Wave Weyl scalar $\\psi_4$\n",
    "\n",
    "## Authors: Patrick Nelson & Zach Etienne\n",
    "\n",
    "## This tutorial demonstrates the construction of quasi-Kinnersley tetrads for the Weyl scalar $\\Psi_4$. It details the process of initializing core NRPy+ modules, setting the numerical grid's coordinate system, defining vectors, and orthogonalizing them using the Levi-Civita symbol.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** See [$\\psi_4$ notebook](Tutorial-Psi4.ipynb), whose Python module depends on the one presented here, for all validation notes.\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/Psi4_tetrads.py](../edit/BSSN/Psi4_tetrads.py)\n",
    "\n",
    "## Introduction: \n",
    "This module constructs tetrad vectors $l^\\mu$, $m^\\mu$, and $n^\\mu$ for the $\\psi_4$ Weyl scalar, a quantity that is immensely useful when extracting gravitational wave content from a numerical relativity simulation. $\\psi_4$ is related to the gravitational wave strain via\n",
    "\n",
    "$$\n",
    "\\psi_4 = \\ddot{h}_+ - i \\ddot{h}_\\times.\n",
    "$$\n",
    "\n",
    "We construct $\\psi_4$ from the standard ADM spatial metric $\\gamma_{ij}$ and extrinsic curvature $K_{ij}$, and their derivatives. The full expression is given by Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf):\n",
    "\n",
    "\\begin{align}\n",
    "\\psi_4 &= \\left[ {R}_{ijkl}+2K_{i[k}K_{l]j}\\right]\n",
    "{n}^i\\bar{m}^j{n}^k\\bar{m}^l  \\\\\n",
    "& -8\\left[ K_{j[k,l]}+{\\Gamma }_{j[k}^pK_{l]p}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^k\\bar{m}^l \\\\\n",
    "& +4\\left[ {R}_{jl}-K_{jp}K_l^p+KK_{jl}\\right]\n",
    "{n}^{[0}\\bar{m}^{j]}{n}^{[0}\\bar{m}^{l]},\n",
    "\\end{align}\n",
    "\n",
    "Note that $\\psi_4$ is complex, with the imaginary components originating from the tetrad vector $m^\\mu$. This module does not specify a tetrad; instead, it only constructs the above expression leaving $m^\\mu$ and $n^\\mu$ unspecified. This module defines these tetrad quantities, implementing the quasi-Kinnersley tetrad of [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), also referred to as \"***the BCL paper***\".\n",
    "\n",
    "### A Note on Notation:\n",
    "\n",
    "As is standard in NRPy+, \n",
    "\n",
    "* Greek indices range from 0 to 3, inclusive, with the zeroth component denoting the temporal (time) component.\n",
    "* Latin indices range from 0 to 2, inclusive, with the zeroth component denoting the first spatial component.\n",
    "\n",
    "As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial notebook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This tutorial notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize needed NRPy+ modules\n",
    "1. [Step 2](#quasikinnersley): The quasi-Kinnersley tetrad\n",
    "1. [Step 3](#code_validation): Code Validation against `BSSN.Psi4_tetrads` NRPy+ module\n",
    "1. [Step 4](#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 core NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:52.011631Z",
     "iopub.status.busy": "2021-03-07T17:15:52.010360Z",
     "iopub.status.idle": "2021-03-07T17:15:53.140162Z",
     "shell.execute_reply": "2021-03-07T17:15:53.140763Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1.a: import all needed modules from NRPy+:\n",
    "import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par    # NRPy+: Parameter interface\n",
    "import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm    # NRPy+: Reference metric support\n",
    "import sys                        # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step 1.b: Set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "\n",
    "# Step 1.c: Given the chosen coordinate system, set up\n",
    "#           corresponding reference metric and needed\n",
    "#           reference metric quantities\n",
    "# The following function call sets up the reference metric\n",
    "#    and related quantities, including rescaling matrices ReDD,\n",
    "#    ReU, and hatted quantities.\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is\n",
    "#           a 3+1-dimensional decomposition of the general\n",
    "#           relativistic field equations)\n",
    "DIM = 3\n",
    "\n",
    "# Step 1.e: Import all ADM quantities as written in terms of BSSN quantities\n",
    "import BSSN.ADM_in_terms_of_BSSN as AB\n",
    "AB.ADM_in_terms_of_BSSN()\n",
    "\n",
    "# Step 1.f: Initialize TetradChoice parameter\n",
    "thismodule = __name__\n",
    "# Current option: QuasiKinnersley = choice made in Baker, Campanelli, and Lousto. PRD 65, 044001 (2002)\n",
    "par.initialize_param(par.glb_param(\"char\", thismodule, \"TetradChoice\", \"QuasiKinnersley\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='quasikinnersley'></a>\n",
    "\n",
    "# Step 2: The quasi-Kinnersley tetrad of [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf) \\[Back to [top](#toc)\\]\n",
    "$$\\label{quasikinnersley}$$\n",
    "\n",
    "To define the Weyl scalars, first, a tetrad must be chosen. Below, for compatibility with the [WeylScal4 diagnostic module](https://bitbucket.org/einsteintoolkit/einsteinanalysis/src/master/WeylScal4/), we implement the quasi-Kinnersley tetrad of [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf).\n",
    "\n",
    "We begin with the vectors given in eqs. 5.6 and 5.7 of the BCL paper, which are orthogonal to each other in flat spacetime; one is in the $\\phi$ direction, one is in $r$, and the third is the cross product of the first two:\n",
    "\\begin{align}\n",
    "    v_1^a &= [-y,x,0] \\\\\n",
    "    v_2^a &= [x,y,z] \\\\\n",
    "    v_3^a &= {\\rm det}(\\gamma)^{1/2} \\gamma^{ad} \\epsilon_{dbc} v_1^b v_2^c,\n",
    "\\end{align}\n",
    "\n",
    "Notice that $v_1^a$ and $v_2^a$ assume the Cartesian basis, but $\\gamma^{ad}$ will be in the $xx^i$ basis given by the chosen `reference_metric::CoordSystem`. Thus to construct $v_3^a$, we must first perform a change of basis on $v_1^a$ and $v_2^a$:\n",
    "\n",
    "$$\n",
    "v_{1,{\\rm xx}}^a = \\frac{\\partial xx^a}{\\partial x_{\\rm Cart}^b} v_{1,{\\rm Cart}}^b.\n",
    "$$\n",
    "This equation is problematic because we generally do not have a closed-form expression for components of the $xx^a$ vector as functions of the Cartesian coordinate vector components $x_{\\rm Cart}^a$. However, we do have closed-form expressions for components of $x_{\\rm Cart}^a$ as functions of $xx^a$. Thus we can construct the needed Jacobian matrix $\\frac{\\partial xx^a}{\\partial x_{\\rm Cart}^b}$ by evaluating the derivative $\\frac{\\partial x_{\\rm Cart}^b}{\\partial xx^a}$ and performing a simple matrix inversion:\n",
    "$$\n",
    "\\frac{\\partial xx^a}{\\partial x_{\\rm Cart}^b} = \\left(\\frac{\\partial x_{\\rm Cart}^b}{\\partial xx^a} \\right)^{-1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:53.215312Z",
     "iopub.status.busy": "2021-03-07T17:15:53.179543Z",
     "iopub.status.idle": "2021-03-07T17:15:55.668569Z",
     "shell.execute_reply": "2021-03-07T17:15:55.669121Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a: Declare the Cartesian x,y,z in terms of\n",
    "#           xx0,xx1,xx2.\n",
    "x = rfm.xx_to_Cart[0]\n",
    "y = rfm.xx_to_Cart[1]\n",
    "z = rfm.xx_to_Cart[2]\n",
    "\n",
    "# Step 2.b: Declare detgamma and gammaUU from\n",
    "#           BSSN.ADM_in_terms_of_BSSN;\n",
    "#           simplify detgamma & gammaUU expressions,\n",
    "#           which expedites Psi4 codegen.\n",
    "detgamma = sp.simplify(AB.detgamma)\n",
    "gammaUU = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        gammaUU[i][j] = sp.simplify(AB.gammaUU[i][j])\n",
    "\n",
    "# Step 2.c: Define v1U and v2U\n",
    "v1UCart = [-y, x, sp.sympify(0)]\n",
    "v2UCart = [x, y, z]\n",
    "\n",
    "# Step 2.d: Construct the Jacobian d x_Cart^i / d xx^j\n",
    "Jac_dUCart_dDrfmUD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Jac_dUCart_dDrfmUD[i][j] = sp.simplify(sp.diff(rfm.xx_to_Cart[i], rfm.xx[j]))\n",
    "\n",
    "# Step 2.e: Invert above Jacobian to get needed d xx^j / d x_Cart^i\n",
    "Jac_dUrfm_dDCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUCart_dDrfmUD)\n",
    "\n",
    "# Step 2.e.i: Simplify expressions for d xx^j / d x_Cart^i:\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Jac_dUrfm_dDCartUD[i][j] = sp.simplify(Jac_dUrfm_dDCartUD[i][j])\n",
    "\n",
    "# Step 2.f: Transform v1U and v2U from the Cartesian to the xx^i basis\n",
    "v1U = ixp.zerorank1()\n",
    "v2U = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        v1U[i] += Jac_dUrfm_dDCartUD[i][j] * v1UCart[j]\n",
    "        v2U[i] += Jac_dUrfm_dDCartUD[i][j] * v2UCart[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... next we construct the third tetrad vector $v_3^a={\\rm det}(\\gamma)^{1/2} \\gamma^{ad} \\epsilon_{dbc} v_1^b v_2^c$, using the Levi-Civita symbol $\\epsilon_{dbc}$ as defined in [indexedexp.py](../edit/indexedexp.py):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:55.708381Z",
     "iopub.status.busy": "2021-03-07T17:15:55.702757Z",
     "iopub.status.idle": "2021-03-07T17:15:55.711174Z",
     "shell.execute_reply": "2021-03-07T17:15:55.710501Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.g: Define v3U\n",
    "v3U = ixp.zerorank1()\n",
    "LeviCivitaSymbolDDD = ixp.LeviCivitaSymbol_dim3_rank3()\n",
    "for a in range(DIM):\n",
    "    for b in range(DIM):\n",
    "        for c in range(DIM):\n",
    "            for d in range(DIM):\n",
    "                v3U[a] += sp.sqrt(detgamma)*gammaUU[a][d]*LeviCivitaSymbolDDD[d][b][c]*v1U[b]*v2U[c]\n",
    "\n",
    "# Step 2.g.i: Simplify expressions for v1U,v2U,v3U. This greatly expedites the C code generation (~10x faster)\n",
    "#             Drat. Simplification with certain versions of SymPy & coord systems results in a hang. Let's just\n",
    "#             evaluate the expressions so the most trivial optimizations can be performed.\n",
    "for a in range(DIM):\n",
    "    v1U[a] = v1U[a].doit() #sp.simplify(v1U[a])\n",
    "    v2U[a] = v2U[a].doit() #sp.simplify(v2U[a])\n",
    "    v3U[a] = v3U[a].doit() #sp.simplify(v3U[a])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As our next step, we carry out the Gram-Schmidt orthonormalization process. The vectors $v_i^a$ are placeholders in the code; the final product of the orthonormalization is the vectors $e_i^a$. So,\n",
    "\\begin{align}\n",
    "e_1^a &= \\frac{v_1^a}{\\sqrt{\\omega_{11}}} \\\\\n",
    "e_2^a &= \\frac{v_2^a - \\omega_{12} e_1^a}{\\sqrt{\\omega_{22}}} \\\\\n",
    "e_3^a &= \\frac{v_3^a - \\omega_{13} e_1^a - \\omega_{23} e_2^a}{\\sqrt{\\omega_{33}}}, \\text{ where}\\\\\n",
    "\\omega_{ij} &= v_i^a v_j^b \\gamma_{ab}\n",
    "\\end{align}\n",
    "\n",
    "Note that the above expressions must be evaluated with the numerators first so that the denominators generate the proper normalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:55.765227Z",
     "iopub.status.busy": "2021-03-07T17:15:55.729307Z",
     "iopub.status.idle": "2021-03-07T17:15:55.808247Z",
     "shell.execute_reply": "2021-03-07T17:15:55.807516Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.h: Define omega_{ij}\n",
    "omegaDD = ixp.zerorank2()\n",
    "gammaDD = AB.gammaDD\n",
    "def v_vectorDU(v1U,v2U,v3U,  i,a):\n",
    "    if i==0:\n",
    "        return v1U[a]\n",
    "    if i==1:\n",
    "        return v2U[a]\n",
    "    if i==2:\n",
    "        return v3U[a]\n",
    "    print(\"ERROR: unknown vector!\")\n",
    "    sys.exit(1)\n",
    "\n",
    "def update_omega(omegaDD, i,j, v1U,v2U,v3U,gammaDD):\n",
    "    omegaDD[i][j] = sp.sympify(0)\n",
    "    for a in range(DIM):\n",
    "        for b in range(DIM):\n",
    "            omegaDD[i][j] += v_vectorDU(v1U,v2U,v3U, i,a)*v_vectorDU(v1U,v2U,v3U, j,b)*gammaDD[a][b]\n",
    "\n",
    "# Step 2.i: Define e^a_i. Note that:\n",
    "#           omegaDD[0][0] = \\omega_{11} above;\n",
    "#           omegaDD[1][1] = \\omega_{22} above, etc.\n",
    "e1U = ixp.zerorank1()\n",
    "e2U = ixp.zerorank1()\n",
    "e3U = ixp.zerorank1()\n",
    "# First e_1^a: Orthogonalize & normalize:\n",
    "update_omega(omegaDD, 0,0, v1U,v2U,v3U,gammaDD)\n",
    "for a in range(DIM):\n",
    "    e1U[a] = v1U[a]/sp.sqrt(omegaDD[0][0])\n",
    "\n",
    "# Next e_2^a: First orthogonalize:\n",
    "update_omega(omegaDD, 0,1, e1U,v2U,v3U,gammaDD)\n",
    "for a in range(DIM):\n",
    "    e2U[a] = (v2U[a] - omegaDD[0][1]*e1U[a])\n",
    "# Then normalize:\n",
    "update_omega(omegaDD, 1,1, e1U,e2U,v3U,gammaDD)\n",
    "for a in range(DIM):\n",
    "    e2U[a] /= sp.sqrt(omegaDD[1][1])\n",
    "\n",
    "# Next e_3^a: First orthogonalize:\n",
    "update_omega(omegaDD, 0,2, e1U,e2U,v3U,gammaDD)\n",
    "update_omega(omegaDD, 1,2, e1U,e2U,v3U,gammaDD)\n",
    "for a in range(DIM):\n",
    "    e3U[a] = (v3U[a] - omegaDD[0][2]*e1U[a] - omegaDD[1][2]*e2U[a])\n",
    "# Then normalize:\n",
    "update_omega(omegaDD, 2,2, e1U,e2U,e3U,gammaDD)\n",
    "for a in range(DIM):\n",
    "    e3U[a] /= sp.sqrt(omegaDD[2][2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have orthogonal, normalized vectors, we can construct the tetrad itself, again drawing on eqs. 5.6. We can draw on SymPy's built-in tools for complex numbers to build the complex vector $m^a$:\n",
    "\\begin{align}\n",
    "    l^\\mu &= \\frac{1}{\\sqrt{2}} \\left(u^\\mu + r^\\mu\\right) \\\\\n",
    "    n^\\mu &= \\frac{1}{\\sqrt{2}} \\left(u^\\mu - r^\\mu\\right) \\\\\n",
    "    \\Re(m^\\mu) &= \\frac{1}{\\sqrt{2}} \\theta^\\mu \\\\\n",
    "    \\Im(m^\\mu) &= \\frac{1}{\\sqrt{2}} \\phi^\\mu,\n",
    "\\end{align}\n",
    "where $r^\\mu=\\{0,e_2^i\\}$, $\\theta^\\mu=\\{0,e_3^i\\}$, $\\phi^\\mu=\\{0,e_1^i\\}$, and $u^\\mu$ is the time-like unit normal to the hypersurface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:55.820257Z",
     "iopub.status.busy": "2021-03-07T17:15:55.819468Z",
     "iopub.status.idle": "2021-03-07T17:15:55.821831Z",
     "shell.execute_reply": "2021-03-07T17:15:55.822345Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Step 2.j: Construct l^mu, n^mu, and m^mu, based on r^mu, theta^mu, phi^mu, and u^mu:\n",
    "r4U     = ixp.zerorank1(DIM=4)\n",
    "u4U     = ixp.zerorank1(DIM=4)\n",
    "theta4U = ixp.zerorank1(DIM=4)\n",
    "phi4U   = ixp.zerorank1(DIM=4)\n",
    "\n",
    "for a in range(DIM):\n",
    "    r4U[    a+1] = e2U[a]\n",
    "    theta4U[a+1] = e3U[a]\n",
    "    phi4U[  a+1] = e1U[a]\n",
    "\n",
    "# FIXME? assumes alpha=1, beta^i = 0\n",
    "u4U[0] = 1\n",
    "\n",
    "l4U = ixp.zerorank1(DIM=4)\n",
    "n4U = ixp.zerorank1(DIM=4)\n",
    "mre4U  = ixp.zerorank1(DIM=4)\n",
    "mim4U  = ixp.zerorank1(DIM=4)\n",
    "\n",
    "# M_SQRT1_2 = 1 / sqrt(2) (defined in math.h on Linux)\n",
    "M_SQRT1_2 = par.Cparameters(\"#define\",thismodule,\"M_SQRT1_2\",\"\")\n",
    "isqrt2 = M_SQRT1_2 #1/sp.sqrt(2) <- SymPy drops precision to 15 sig. digits in unit tests\n",
    "for mu in range(4):\n",
    "    l4U[mu]   = isqrt2*(u4U[mu] + r4U[mu])\n",
    "    n4U[mu]   = isqrt2*(u4U[mu] - r4U[mu])\n",
    "    mre4U[mu] = isqrt2*theta4U[mu]\n",
    "    mim4U[mu] = isqrt2*  phi4U[mu]\n",
    "\n",
    "# ltetU,ntetU,remtetU,immtetU,e1U,e2U,e3U\n",
    "for mu in range(4):\n",
    "    l4U[mu]   = isqrt2*(u4U[mu] + r4U[mu])\n",
    "    n4U[mu]   = isqrt2*(u4U[mu] - r4U[mu])\n",
    "    mre4U[mu] = isqrt2*theta4U[mu]\n",
    "    mim4U[mu] = isqrt2*  phi4U[mu]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code validation against `BSSN.Psi4_tetrads` 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.Psi4_tetrads](../edit/BSSN/Psi4_tetrads.py) module.\n",
    "\n",
    "By default, we compare all quantities in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:55.833390Z",
     "iopub.status.busy": "2021-03-07T17:15:55.832620Z",
     "iopub.status.idle": "2021-03-07T17:15:58.329594Z",
     "shell.execute_reply": "2021-03-07T17:15:58.330100Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALL TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "def comp_func(expr1,expr2,basename,prefixname2=\"BP4T.\"):\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",
    "\n",
    "import BSSN.Psi4_tetrads as BP4T\n",
    "BP4T.Psi4_tetrads()\n",
    "\n",
    "for mu in range(4):\n",
    "    namecheck_list.extend([gfnm(\"l4U\",mu),gfnm(\"n4U\",mu),gfnm(\"mre4U\",mu),gfnm(\"mim4U\",mu)])\n",
    "    exprcheck_list.extend([BP4T.l4U[mu],BP4T.n4U[mu],BP4T.mre4U[mu],BP4T.mim4U[mu]])\n",
    "    expr_list.extend([l4U[mu],n4U[mu],mre4U[mu],mim4U[mu]])\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='latex_pdf_output'></a>\n",
    "\n",
    "# Step 4: 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-Psi4_tetrads.pdf](Tutorial-Psi4_tetrads.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:58.335326Z",
     "iopub.status.busy": "2021-03-07T17:15:58.334227Z",
     "iopub.status.idle": "2021-03-07T17:16:01.828438Z",
     "shell.execute_reply": "2021-03-07T17:16:01.829186Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Psi4_tetrads.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-Psi4_tetrads.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Psi4_tetrads\")"
   ]
  }
 ],
 "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
}