Tutorial-BSSN_stress_energy_source_terms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [BSSN](http://www2.yukawa.kyoto-u.ac.jp/~yuichiro.sekiguchi/3+1.pdf) Stress-Energy Source Terms\n",
    "\n",
    "## Author: Zach Etienne\n",
    "\n",
    "## This module constructs the BSSN stress-energy source terms in terms of $T^{\\mu\\nu}$\n",
    "\n",
    "### This module implements the BSSN source terms as prescribed in the reference metric approach of [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632), which builds upon the covariant \"Lagrangian\" BSSN formalism of [Brown (2009)](https://arxiv.org/abs/0902.3652).\n",
    "\n",
    "**Notebook Status:** <font color='orange'><b> Self-validated </b></font>\n",
    "\n",
    "**Validation Notes:** None yet.\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/BSSN_stress_energy_source_terms.py](../edit/BSSN/BSSN_stress_energy_source_terms.py)\n",
    "\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In [the NRPy+ tutorial on the BSSN formulation](Tutorial-BSSN_formulation.ipynb) we outlined the BSSN formulation of Einstein's equations *in the absence of stress-energy* (i.e., in vacuum where Einstein's equations reduce to $G^{\\mu\\nu}=0$). When $T^{\\mu\\nu}$ is nonzero, stress-energy source terms must appear on the right-hand sides of the BSSN equations in order to ensure Einstein's equations of general relativity are satisfied.\n",
    "\n",
    "Analyzing Eqs. 9 of [Baumgarte, Montero, Cordero-Carrión, and Müller](https://arxiv.org/pdf/1211.6632.pdf), we see that adding stress-energy source terms $T_{\\mu\\nu}$ to Einstein's equations in vacuum simply adjusts the right-hand sides of the $\\partial_t \\bar{A}_{ij}$, $\\partial_t K$, and $\\partial_t \\bar{\\Lambda}^i$ equations as follows:\n",
    "\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\ \\partial_t \\bar{A}_{ij} &=& \\left[\\partial_t \\bar{A}_{ij}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 8\\pi \\alpha e^{-4\\phi} \\left(S_{ij}\\right)^{\\rm TF}}} \\\\\n",
    "\\partial_t K &=& \\left[\\partial_t K\\right]_{\\rm vacuum}\\ {\\color{blue}{+\\ 4\\pi \\alpha (\\rho + S)}} \\\\\n",
    "\\partial_t \\bar{\\Lambda}^i &=& \\left[\\partial_t \\bar{\\Lambda}^{i}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 16\\pi \\alpha \\bar{\\gamma}^{ij} S_j}},\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $\\rho$, $S$, $S_i$, and $S_{ij}$ are related to the stress-energy tensor $T^{\\mu\\nu}$ as follows (Eq. 10 of [Baumgarte, Montero, Cordero-Carrión, and Müller](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\ S_{ij} &=& \\gamma_{i \\mu} \\gamma_{j \\nu} T^{\\mu \\nu} \\\\\n",
    "S_{i} &=& -\\gamma_{i\\mu} n_\\nu T^{\\mu\\nu} \\\\\n",
    "S &=& \\gamma^{ij} S_{ij} \\\\\n",
    "\\rho &=& n_\\mu n_\\nu T^{\\mu\\nu},\n",
    "\\end{eqnarray}\n",
    "\n",
    "the unit normal one-form on each spatial slice $n_{\\mu}$ is given by Eq. 10 of [Baumgarte, Montero, Cordero-Carrión, and Müller](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "\n",
    "$$\n",
    "n_\\mu = (-\\alpha,0,0,0),\n",
    "$$\n",
    "\n",
    "and Baumgarte & Shapiro Eq. 2.27 gives $\\gamma_{\\mu\\nu}$:\n",
    "\n",
    "$$\\gamma_{\\mu\\nu} = g_{\\mu\\nu} + n_\\mu n_\\nu.$$\n",
    "\n",
    "Further, analyzing Eqs. 13 & 14 of [Baumgarte, Montero, Cordero-Carrión, and Müller](https://arxiv.org/pdf/1211.6632.pdf) we find that adding stress-energy source terms $T_{\\mu\\nu}$ to Einstein's equations in vacuum adjusts the BSSN constraint equations as follows:\n",
    "\\begin{eqnarray}\n",
    "\\ \\partial_t \\mathcal{H} &=& \\left[\\partial_t \\mathcal{H}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 16\\pi \\rho}} \\\\\n",
    "\\partial_t \\mathcal{M}^i &=& \\left[\\partial_t \\mathcal{M}^i\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 8\\pi S^i}},\n",
    "\\end{eqnarray}\n",
    "\n",
    "This module will construct expressions for $S_{ij}$, $S_i$, $S$, and $\\rho$ in terms of $T^{\\mu\\nu}$, and also add the necessary terms to the BSSN RHSs and constraints.\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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#bssn_sourceterms): BSSN source terms, in terms of $T^{\\mu\\nu}$\n",
    "    1. [Step 2.a](#gamma4dd): Define `gamma4DD[mu][nu]` = $g_{\\mu \\nu} + n_{\\mu} n_{\\nu}$\n",
    "    1. [Step 2.b](#t4uu): Declare `T4UU[mu][nu]`=$T^{\\mu\\nu}$\n",
    "    1. [Step 2.c](#define_bssn_sourceterms): Define BSSN source terms    \n",
    "1. [Step 3](#add_bssn_sourceterms_to_rhss): Add BSSN source terms to BSSN RHSs\n",
    "1. [Step 4](#add_bssn_sourceterms_to_constraints): Add BSSN source terms to BSSN Constraints\n",
    "1. [Step 5](#code_validation): Code Validation against `BSSN.BSSN_stress_energy_source_terms` NRPy+ module\n",
    "1. [Step 6](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize needed Python/NRPy+ modules & set up reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Initialize needed Python/NRPy+ modules\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 BSSN.ADMBSSN_tofrom_4metric as AB4m # NRPy+: ADM/BSSN <-> 4-metric conversions\n",
    "import BSSN.ADM_in_terms_of_BSSN as AitoB  # NRPy+: ADM quantities in terms of BSSN quantities\n",
    "\n",
    "# Step 1.a: Set up reference metric. We'll choose SinhSpherical here, but\n",
    "#           could choose any CoordSystem defined in reference_metric.py:\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "thismodule = \"BSSN.BSSN_stress_energy_source_terms\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn_sourceterms'></a>\n",
    "\n",
    "# Step 2: BSSN source terms, in terms of $T^{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_sourceterms}$$\n",
    "\n",
    "<a id='gamma4dd'></a>\n",
    "\n",
    "## Step 2.a: Define `gamma4DD[mu][nu]` = $g_{\\mu \\nu} + n_{\\mu} n_{\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{gamma4dd}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.a: Define gamma4DD[mu][nu] = g_{mu nu} + n_{mu} n_{nu}\n",
    "alpha = sp.symbols(\"alpha\",real=True)\n",
    "zero  = sp.sympify(0)\n",
    "n4D      = [-alpha, zero, zero ,zero]\n",
    "AB4m.g4DD_ito_BSSN_or_ADM(\"BSSN\")\n",
    "\n",
    "gamma4DD = ixp.zerorank2(DIM=4)\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        gamma4DD[mu][nu] = AB4m.g4DD[mu][nu] + n4D[mu]*n4D[nu]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='t4uu'></a>\n",
    "\n",
    "## Step 2.b: Declare `T4UU[mu][nu]`=$T^{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{t4uu}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.b: Declare T4UU\n",
    "T4UU = ixp.declarerank2(\"T4UU\",\"sym01\",DIM=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='define_bssn_sourceterms'></a>\n",
    "\n",
    "## Step 2.c: Define BSSN source terms \\[Back to [top](#toc)\\]\n",
    "$$\\label{define_bssn_sourceterms}$$\n",
    "\n",
    "Recall from above, we have:\n",
    "\\begin{eqnarray}\n",
    "\\ S_{ij} &=& \\gamma_{i \\mu} \\gamma_{j \\nu} T^{\\mu \\nu} \\\\\n",
    "S_{i} &=& -\\gamma_{i\\mu} n_\\nu T^{\\mu\\nu} \\\\\n",
    "S &=& \\gamma^{ij} S_{ij} \\\\\n",
    "\\rho &=& n_\\mu n_\\nu T^{\\mu\\nu}.\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.c: Define BSSN source terms\n",
    "# Step 2.c.i: S_{ij} = gamma_{i mu} gamma_{j nu} T^{mu nu}\n",
    "SDD = ixp.zerorank2()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        for mu in range(4):\n",
    "            for nu in range(4):\n",
    "                SDD[i][j] += gamma4DD[i+1][mu] * gamma4DD[j+1][nu] * T4UU[mu][nu]\n",
    "# Step 2.c.ii: S_{i} = -gamma_{i mu} n_{nu} T^{mu nu}\n",
    "SD = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            SD[i] += - gamma4DD[i+1][mu] * n4D[nu] * T4UU[mu][nu]\n",
    "# Step 2.c.iii: S = gamma^{ij} S_{ij}\n",
    "AitoB.ADM_in_terms_of_BSSN()\n",
    "S = zero\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        S += AitoB.gammaUU[i][j]*SDD[i][j]\n",
    "# Step 2.c.iv: rho = n_{mu} n_{nu} T^{mu nu}\n",
    "rho = zero\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        rho += n4D[mu]*n4D[nu]*T4UU[mu][nu]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='add_bssn_sourceterms_to_rhss'></a>\n",
    "\n",
    "# Step 3: Add BSSN source terms to BSSN RHSs \\[Back to [top](#toc)\\]\n",
    "$$\\label{add_bssn_sourceterms_to_rhss}$$\n",
    "\n",
    "Recall from above we need to make the following modifications:\n",
    "\\begin{eqnarray}\n",
    "\\ \\partial_t \\bar{A}_{ij} &=& \\left[\\partial_t \\bar{A}_{ij}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 8\\pi \\alpha e^{-4\\phi} \\left(S_{ij}\\right)^{\\rm TF}}} \\\\\n",
    "\\partial_t K &=& \\left[\\partial_t K\\right]_{\\rm vacuum}\\ {\\color{blue}{+\\ 4\\pi \\alpha (\\rho + S)}} \\\\\n",
    "\\partial_t \\bar{\\Lambda}^i &=& \\left[\\partial_t \\bar{\\Lambda}^{i}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 16\\pi \\alpha \\bar{\\gamma}^{ij} S_j}},\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $$\\left(S_{ij}\\right)^{\\rm TF} = S_{ij} - \\frac{1}{3} \\bar{\\gamma}_{ij} \\bar{\\gamma}^{km} S_{km}.$$\n",
    "\n",
    "*Exercise to student:* Prove that replacing the $\\bar{\\gamma}_{ij}$ and $\\bar{\\gamma}^{km}$ with $\\gamma_{ij}$ and $\\gamma^{km}$, respectively, results in exactly the same expression for $\\left(S_{ij}\\right)^{\\rm TF}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3: Add BSSN stress-energy source terms to BSSN RHSs\n",
    "import BSSN.BSSN_quantities as Bq\n",
    "# Can't #declare M_PI here, as it is not SIMD-compatible.\n",
    "PI = par.Cparameters(\"REAL\", thismodule, [\"PI\"], \"3.14159265358979323846264338327950288\")\n",
    "alpha = sp.symbols(\"alpha\",real=True)\n",
    "zero  = sp.sympify(0)\n",
    "\n",
    "# Step 3.a: Initialize RHS source terms to zero.\n",
    "sourceterm_trK_rhs     = zero\n",
    "sourceterm_a_rhsDD     = ixp.zerorank2()\n",
    "sourceterm_lambda_rhsU = ixp.zerorank1()\n",
    "\n",
    "# Step 3.b: trK_rhs\n",
    "sourceterm_trK_rhs = 4*PI*alpha*(rho + S)\n",
    "\n",
    "# Step 3.c: Abar_rhsDD:\n",
    "# Step 3.c.i: Compute trace-free part of S_{ij}:\n",
    "Bq.BSSN_basic_tensors() # Sets gammabarDD\n",
    "gammabarUU, dummydet = ixp.symm_matrix_inverter3x3(Bq.gammabarDD) # Set gammabarUU\n",
    "tracefree_SDD = ixp.zerorank2()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        tracefree_SDD[i][j] = SDD[i][j]\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        for k in range(3):\n",
    "            for m in range(3):\n",
    "                tracefree_SDD[i][j] += -sp.Rational(1,3)*Bq.gammabarDD[i][j]*gammabarUU[k][m]*SDD[k][m]\n",
    "\n",
    "# Step 3.c.ii: Define exp_m4phi = e^{-4 phi}\n",
    "Bq.phi_and_derivs()\n",
    "# Step 3.c.iii: Evaluate RHS\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        Abar_rhsDDij             = -8*PI*alpha*Bq.exp_m4phi*tracefree_SDD[i][j]\n",
    "        sourceterm_a_rhsDD[i][j] =  Abar_rhsDDij / rfm.ReDD[i][j]\n",
    "\n",
    "# Step 3.d: Stress-energy part of Lambdabar_rhsU = stressenergy_Lambdabar_rhsU\n",
    "sourceterm_Lambdabar_rhsU = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        sourceterm_Lambdabar_rhsU[i] += -16*PI*alpha*gammabarUU[i][j]*SD[j]\n",
    "for i in range(3):\n",
    "    sourceterm_lambda_rhsU[i] = sourceterm_Lambdabar_rhsU[i] / rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='add_bssn_sourceterms_to_constraints'></a>\n",
    "\n",
    "# Step 4: Add BSSN source terms to BSSN Constraints \\[Back to [top](#toc)\\]\n",
    "$$\\label{add_bssn_sourceterms_to_constraints}$$\n",
    "\n",
    "Recall from above we need to make the following modifications:\n",
    "\\begin{eqnarray}\n",
    "\\ \\partial_t \\mathcal{H} &=& \\left[\\partial_t \\mathcal{H}\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 16\\pi \\rho}} \\\\\n",
    "\\partial_t \\mathcal{M}^i &=& \\left[\\partial_t \\mathcal{M}^i\\right]_{\\rm vacuum}\\ {\\color{blue}{-\\ 8\\pi S^i}},\n",
    "\\end{eqnarray}\n",
    "\n",
    "where \n",
    "$$\n",
    "S^i = \\gamma^{ij} S_j,\n",
    "$$\n",
    "and $\\gamma^{ij}$ is the inverse ADM 3-metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4: Add BSSN stress-energy source terms to BSSN constraints\n",
    "# Step 4.a: Initialize constraint source terms to zero.\n",
    "sourceterm_H  = sp.sympify(0)\n",
    "sourceterm_MU = ixp.zerorank1()\n",
    "\n",
    "# Step 4.b: Add source term to the Hamiltonian constraint H\n",
    "sourceterm_H = -16*PI*rho\n",
    "\n",
    "# Step 4.c: Add source term to the momentum constraint M^i\n",
    "# Step 4.c.i: Compute gammaUU in terms of BSSN quantities\n",
    "import BSSN.ADM_in_terms_of_BSSN as AitoB\n",
    "AitoB.ADM_in_terms_of_BSSN() # Provides gammaUU\n",
    "# Step 4.c.ii: Raise S_i\n",
    "SU = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        SU[i] += AitoB.gammaUU[i][j]*SD[j]\n",
    "# Step 4.c.iii: Add source term to momentum constraint & rescale:\n",
    "for i in range(3):\n",
    "    sourceterm_MU[i] = -8 * PI * SU[i] / rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 5: Code Validation against `BSSN.BSSN_stress_energy_source_terms` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of the BSSN equations between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.BSSN_stressenergy_source_terms](../edit/BSSN/BSSN_stressenergy_source_terms.py) module.\n",
    "\n",
    "By default, we analyze these expressions in SinhSpherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between BSSN_stress_energy_source_terms tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
      "STRESS-ENERGY SOURCE TERMS:\n",
      "rho - Bsest.rho = 0\n",
      "S   - Bsest.S   = 0\n",
      "SD[0] - Bsest.SD[0] = 0\n",
      "SD[1] - Bsest.SD[1] = 0\n",
      "SD[2] - Bsest.SD[2] = 0\n",
      "SDD[0][0] - Bsest.SDD[0][0] = 0\n",
      "SDD[0][1] - Bsest.SDD[0][1] = 0\n",
      "SDD[0][2] - Bsest.SDD[0][2] = 0\n",
      "SDD[1][0] - Bsest.SDD[1][0] = 0\n",
      "SDD[1][1] - Bsest.SDD[1][1] = 0\n",
      "SDD[1][2] - Bsest.SDD[1][2] = 0\n",
      "SDD[2][0] - Bsest.SDD[2][0] = 0\n",
      "SDD[2][1] - Bsest.SDD[2][1] = 0\n",
      "SDD[2][2] - Bsest.SDD[2][2] = 0\n",
      "\n",
      "BSSN RHSs SOURCE TERMS:\n",
      "sourceterm_trK_rhs - Bsest.sourceterm_trK_rhs = 0\n",
      "sourceterm_a_rhsDD[0][0] - Bsest.sourceterm_a_rhsDD[0][0] = 0\n",
      "sourceterm_a_rhsDD[0][1] - Bsest.sourceterm_a_rhsDD[0][1] = 0\n",
      "sourceterm_a_rhsDD[0][2] - Bsest.sourceterm_a_rhsDD[0][2] = 0\n",
      "sourceterm_a_rhsDD[1][0] - Bsest.sourceterm_a_rhsDD[1][0] = 0\n",
      "sourceterm_a_rhsDD[1][1] - Bsest.sourceterm_a_rhsDD[1][1] = 0\n",
      "sourceterm_a_rhsDD[1][2] - Bsest.sourceterm_a_rhsDD[1][2] = 0\n",
      "sourceterm_a_rhsDD[2][0] - Bsest.sourceterm_a_rhsDD[2][0] = 0\n",
      "sourceterm_a_rhsDD[2][1] - Bsest.sourceterm_a_rhsDD[2][1] = 0\n",
      "sourceterm_a_rhsDD[2][2] - Bsest.sourceterm_a_rhsDD[2][2] = 0\n",
      "sourceterm_lambda_rhsU[0] - Bsest.sourceterm_lambda_rhsU[0] = 0\n",
      "sourceterm_lambda_rhsU[1] - Bsest.sourceterm_lambda_rhsU[1] = 0\n",
      "sourceterm_lambda_rhsU[2] - Bsest.sourceterm_lambda_rhsU[2] = 0\n",
      "\n",
      "BSSN CONSTRAINTS SOURCE TERMS:\n",
      "sourceterm_H - Bsest.sourceterm_H = 0\n",
      "sourceterm_MU[0] - Bsest.sourceterm_MU[0] = 0\n",
      "sourceterm_MU[1] - Bsest.sourceterm_MU[1] = 0\n",
      "sourceterm_MU[2] - Bsest.sourceterm_MU[2] = 0\n"
     ]
    }
   ],
   "source": [
    "# Step 5: Code Validation against BSSN.BSSN_stress_energy_source_terms NRPy+ module\n",
    "\n",
    "# We already have SymPy expressions for BSSN source terms\n",
    "#         in terms of other SymPy variables.\n",
    "#\n",
    "#         Here, we will use the above-defined BSSN stress-energy source term expressions\n",
    "#         to validate against the same expressions in the\n",
    "#         BSSN/BSSN_stress_energy_source_terms.py file, to ensure consistency between\n",
    "#         this tutorial and the module itself.\n",
    "import BSSN.BSSN_stress_energy_source_terms as Bsest\n",
    "\n",
    "print(\"Consistency check between BSSN_stress_energy_source_terms tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
    "\n",
    "print(\"STRESS-ENERGY SOURCE TERMS:\")\n",
    "Bsest.stress_energy_source_terms_ito_T4UU_and_ADM_or_BSSN_metricvars(\"BSSN\")\n",
    "print(\"rho - Bsest.rho = \" + str(rho - Bsest.rho))\n",
    "print(\"S   - Bsest.S   = \" + str(S   - Bsest.S))\n",
    "for i in range(3):\n",
    "    print(\"SD[\"+str(i)+\"] - Bsest.SD[\"+str(i)+\"] = \" + str(SD[i] - Bsest.SD[i]))\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        print(\"SDD[\"+str(i)+\"][\"+str(j)+\"] - Bsest.SDD[\"+str(i)+\"][\"+str(j)+\"] = \" + str(SDD[i][j] - Bsest.SDD[i][j]))\n",
    "\n",
    "print(\"\\nBSSN RHSs SOURCE TERMS:\")\n",
    "Bsest.BSSN_source_terms_for_BSSN_RHSs()\n",
    "print(\"sourceterm_trK_rhs - Bsest.sourceterm_trK_rhs = \" + str(sourceterm_trK_rhs - Bsest.sourceterm_trK_rhs))\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        print(\"sourceterm_a_rhsDD[\"+str(i)+\"][\"+str(j)+\"] - Bsest.sourceterm_a_rhsDD[\"+str(i)+\"][\"+str(j)+\"] = \" +\n",
    "              str(sourceterm_a_rhsDD[i][j] - Bsest.sourceterm_a_rhsDD[i][j]))\n",
    "for i in range(3):\n",
    "    print(\"sourceterm_lambda_rhsU[\"+str(i)+\"] - Bsest.sourceterm_lambda_rhsU[\"+str(i)+\"] = \" +\n",
    "          str(sourceterm_lambda_rhsU[i] - Bsest.sourceterm_lambda_rhsU[i]))\n",
    "\n",
    "\n",
    "print(\"\\nBSSN CONSTRAINTS SOURCE TERMS:\")\n",
    "Bsest.BSSN_source_terms_for_BSSN_constraints()\n",
    "print(\"sourceterm_H - Bsest.sourceterm_H = \" + str(sourceterm_H - Bsest.sourceterm_H))\n",
    "for i in range(3):\n",
    "    print(\"sourceterm_MU[\"+str(i)+\"] - Bsest.sourceterm_MU[\"+str(i)+\"] = \" +\n",
    "          str(sourceterm_MU[i] - Bsest.sourceterm_MU[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-BSSN_stress_energy_source_terms.pdf](Tutorial-BSSN_stress_energy_source_terms.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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_stress_energy_source_terms.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-BSSN_stress_energy_source_terms.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_stress_energy_source_terms\")"
   ]
  }
 ],
 "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
}