Tutorial-ScalarField_RHSs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Tutorial: The right-hand sides of the time-evolution equation for a massless scalar field\n",
    "\n",
    "## Author(s): Leonardo Werneck and Zach Etienne\n",
    "\n",
    "# This tutorial notebook documents and constructs the time evolution equations of the Klein-Gordon equations for a massless scalar field written in terms of the BSSN quantities.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** The expressions generated by the NRPy+ module corresponding to this tutorial notebook were used to generate the results shown in [Werneck *et al.* (in preparation)]().</font>\n",
    "\n",
    "## Python module containing the final expressions constructed here: [ScalarField/ScalarField_RHSs.py](../edit/ScalarField/ScalarField_RHSs.py)\n",
    "\n",
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "The module is organized as follows\n",
    "\n",
    "0. [Step 0](#klein_gordon_eq): The Klein-Gordon equation in curvilinear coordinates\n",
    "1. [Step 1](#initializenrpy): Initialize core NRPy+ modules\n",
    "1. [Step 2](#sf_rhs): Right-hand side of $\\partial_{t} \\varphi$\n",
    "1. [Step 3](#pi_rhs): Right-hand side of $\\partial_{t} \\Pi$\n",
    "1. [Step 4](#code_validation): Code Validation against [ScalarField.ScalarField_RHSs.py](../edit/ScalarField/ScalarField_RHSs.py) NRPy+ module\n",
    "1. [Step 5](#latex_pdf_output): Output this module to $\\LaTeX$-formatted PDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='klein_gordon_eq'></a>\n",
    "\n",
    "# Step 0: The Klein-Gordon equation in curvilinear coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{klein_gordon_eq}$$\n",
    "\n",
    "We will begin by considering the [Klein-Gordon equation](https://en.wikipedia.org/wiki/Klein%E2%80%93Gordon_equation#Gravitational_interaction) for a massless scalar field, $\\varphi$,\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mu}\\left(\\nabla^{\\mu}\\varphi\\right) = 0\\ .\n",
    "$$\n",
    "\n",
    "We then define the auxiliary field\n",
    "\n",
    "$$\n",
    "\\Pi \\equiv -\\frac{1}{\\alpha}\\left(\\partial_{t}\\varphi - \\beta^{i}\\partial_{i}\\varphi\\right)\\ ,\n",
    "$$\n",
    "\n",
    "so that the Klein-Gordon equation is decomposed into two first order equations (cf. eqs. (5.252) and (5.253) in B&S)\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_{t}\\varphi &= \\beta^{i}\\partial_{i}\\varphi - \\alpha\\,\\Pi\\ ,\\nonumber\\\\\n",
    "\\partial_{t}\\Pi &= \\beta^{i}\\partial_{i}\\Pi + \\alpha K\\,\\Pi -\\gamma^{ij}\\left(\\partial_{j}\\varphi\\,\\partial_{i}\\alpha -\\alpha\\,\\Gamma^{k}_{\\ ij}\\, \\partial_{k}\\varphi+\\alpha\\,\\partial_{i}\\partial_{j}\\varphi\\right)\\ ,\n",
    "\\end{align}\n",
    "\n",
    "where $K=\\gamma^{ij}K_{ij}$ is the trace of the extrinsic curvature $K_{ij}$ and $\\gamma^{ij}$ the inverse of the physical spatial metric $\\gamma_{ij}$. We will choose *not* to define the auxiliary variables $\\varphi_{i}\\equiv\\partial_{i}\\varphi$ in our code and, instead, leave the equations in terms of second derivatives of $\\varphi$.\n",
    "\n",
    "To write the equations in terms of BSSN variables (see [this tutorial notebook](Tutorial-BSSN_formulation.ipynb) for a review) we start by considering the conformal metric, $\\bar{\\gamma}_{ij}$, related to the physical metric by\n",
    "\n",
    "$$\n",
    "\\gamma_{ij} = e^{4\\phi}\\bar{\\gamma}_{ij}\\ ,\n",
    "$$\n",
    "\n",
    "and its inverse\n",
    "\n",
    "$$\n",
    "\\gamma^{ij} = e^{-4\\phi}\\bar{\\gamma}^{ij}\\ .\n",
    "$$\n",
    "\n",
    "Let us also look at equation (3.7) of B&S (with $i\\leftrightarrow k$ and $\\ln\\psi = \\phi$, for convenience)\n",
    "\n",
    "$$\n",
    "\\Gamma^{k}_{\\ ij} = \\bar{\\Gamma}^{k}_{\\ ij} + 2\\left(\\delta^{k}_{\\ i}\\bar{D}_{j}\\phi + \\delta^{k}_{\\ j}\\bar{D}_{i}\\phi - \\bar{\\gamma}_{ij}\\bar{\\gamma}^{k\\ell}\\bar{D}_{\\ell}\\phi\\right)\\ .\n",
    "$$\n",
    "\n",
    "Then let us consider the term that contains $\\Gamma^{k}_{\\ ij}$ on the right-hand side of $\\partial_{t}\\Pi$:\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha\\,\\gamma^{ij}\\Gamma^{k}_{\\ ij}\\, \\partial_{k}\\varphi &= \\alpha e^{-4\\phi}\\bar{\\gamma}^{ij}\\left[\\bar{\\Gamma}^{k}_{\\ ij} + 2\\left(\\delta^{k}_{\\ i}\\bar{D}_{j}\\phi + \\delta^{k}_{\\ j}\\bar{D}_{i}\\phi - \\bar{\\gamma}_{ij}\\bar{\\gamma}^{k\\ell}\\bar{D}_{\\ell}\\phi\\right)\\right]\\partial_{k}\\varphi\\ .\n",
    "\\end{align}\n",
    "\n",
    "Focusing on the term in parenthesis, we have (ignoring, for now, a few non-essential multiplicative terms and replacing $\\bar{D}_{i}\\phi = \\partial_{i}\\phi$)\n",
    "\n",
    "\\begin{align}\n",
    "2\\bar{\\gamma}^{ij}\\left(\\delta^{k}_{\\ i}\\bar{D}_{j}\\phi + \\delta^{k}_{\\ j}\\bar{D}_{i}\\phi - \\bar{\\gamma}_{ij}\\bar{\\gamma}^{k\\ell}\\bar{D}_{\\ell}\\phi\\right)\\partial_{k}\\varphi &= 2\\left(\\bar{\\gamma}^{kj}\\partial_{j}\\phi + \\bar{\\gamma}^{ki}\\partial_{i}\\phi - 3\\bar{\\gamma}^{k\\ell}\\partial_{\\ell}\\phi\\right)\\partial_{k}\\varphi\\nonumber\\\\\n",
    "&=2\\left(\\bar{\\gamma}^{ij}\\partial_{i}\\phi + \\bar{\\gamma}^{ij}\\partial_{i}\\phi - 3 \\bar{\\gamma}^{ij}\\partial_{i}\\phi\\right)\\partial_{j}\\varphi\\nonumber\\\\\n",
    "&= -2\\bar{\\gamma}^{ij}\\partial_{j}\\varphi\\partial_{i}\\phi\\ ,\n",
    "\\end{align}\n",
    "\n",
    "so that\n",
    "\n",
    "$$\n",
    "\\alpha\\,\\gamma^{ij}\\Gamma^{k}_{\\ ij}\\, \\partial_{k}\\varphi = e^{-4\\phi}\\bar{\\gamma}^{ij}\\left(\\alpha \\bar{\\Gamma}^{k}_{\\ ij}\\, \\partial_{k}\\varphi - 2\\alpha\\partial_{j}\\varphi\\partial_{i}\\phi\\right)\n",
    "$$\n",
    "\n",
    "For the rest of the equation, all we need to do is replace $\\gamma^{ij}\\to e^{-4\\phi}\\bar{\\gamma}^{ij}$, so that the Klein-Gordon equation becomes\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_{t}\\Pi = \\beta^{i}\\partial_{i}\\Pi + \\alpha K\\,\\Pi -e^{-4\\phi}\\bar{\\gamma}^{ij}\\left(\\partial_{j}\\varphi\\,\\partial_{i}\\alpha - \\alpha\\,\\bar{\\Gamma}^{k}_{\\ ij}\\,\\partial_{k}\\varphi + \\alpha\\,\\partial_{i}\\partial_{j}\\varphi + 2\\alpha\\,\\partial_{j}\\varphi\\,\\partial_{i}\\phi\\right)\\ .\n",
    "\\end{align}\n",
    "\n",
    "Note that the evolution equation for $\\varphi$ is left unchanged\n",
    "\n",
    "$$\n",
    "\\partial_{t}\\varphi = \\beta^{i}\\partial_{i}\\varphi - \\alpha\\,\\Pi\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:15.650713Z",
     "iopub.status.busy": "2021-06-15T10:05:15.650257Z",
     "iopub.status.idle": "2021-06-15T10:05:16.120482Z",
     "shell.execute_reply": "2021-06-15T10:05:16.119907Z"
    }
   },
   "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 grid as gri                # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm    # NRPy+: Reference metric support\n",
    "import BSSN.BSSN_quantities as Bq # NRPy+: BSSN quantities\n",
    "\n",
    "# Step 1.b: Set the coordinate system for the numerical grid\n",
    "coordsystem = \"Spherical\"\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",coordsystem)\n",
    "\n",
    "# Step 1.c: Then we set the theta and phi axes to be the symmetry\n",
    "#           axes; i.e., axis \"1\" and \"2\", corresponding to the\n",
    "#           i1 and i2 directions.\n",
    "#           This sets all spatial derivatives in the theta and\n",
    "#           phi directions to zero (analytically).\n",
    "# par.set_parval_from_str(\"indexedexp::symmetry_axes\")\n",
    "\n",
    "# Step 1.d: 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.e: 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.f: Import all needed BSSN quantities\n",
    "Bq.BSSN_basic_tensors()\n",
    "trK   = Bq.trK\n",
    "alpha = Bq.alpha\n",
    "betaU = Bq.betaU\n",
    "Bq.gammabar__inverse_and_derivs()\n",
    "gammabarUU = Bq.gammabarUU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sf_rhs'></a>\n",
    "\n",
    "# Step 2: The right-hand side of $\\partial_{t}\\varphi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{sf_rhs}$$\n",
    "\n",
    "Let us label each of the terms in the RHS of the $\\partial_{t}\\varphi$ equation so that it is easier to understand their implementation:\n",
    "\n",
    "$$\n",
    "\\partial_{t}\\varphi = - \\underbrace{\\alpha \\Pi}_{\\text{Term 1}} + \\underbrace{\\beta^{i}\\partial_{i}\\varphi}_{\\text{Term 2}} \\ .\n",
    "$$\n",
    "\n",
    "The first term is an advection term and therefore we need to set up the appropriate derivate. We begin by declaring the grid functions that we need to implement this equation.\n",
    "\n",
    "**A note on notation:** We choose to declare the $\\color{red}{\\textbf{s}}$calar $\\color{red}{\\textbf{f}}$ield variable, $\\varphi$, as $\\color{red}{\\rm sf}$ and the $\\color{blue}{\\textbf{s}}$calar $\\color{blue}{\\textbf{f}}$ield conjugate $\\color{blue}{\\textbf{M}}$omentum variable, $\\Pi$, as $\\color{blue}{\\rm sfM}$ to avoid possible conflicts with other variables which might be commonly denoted by $\\psi$ or $\\Pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.124877Z",
     "iopub.status.busy": "2021-06-15T10:05:16.124362Z",
     "iopub.status.idle": "2021-06-15T10:05:16.126468Z",
     "shell.execute_reply": "2021-06-15T10:05:16.125952Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a: Declare grid functions for varphi and Pi\n",
    "sf, sfM = gri.register_gridfunctions(\"EVOL\",[\"sf\", \"sfM\"])\n",
    "\n",
    "# Step 2.a: Add Term 1 to sf_rhs: -alpha*Pi\n",
    "sf_rhs = - alpha * sfM\n",
    "\n",
    "# Step 2.b: Add Term 2 to sf_rhs: beta^{i}\\partial_{i}\\varphi\n",
    "sf_dupD = ixp.declarerank1(\"sf_dupD\")\n",
    "for i in range(DIM):\n",
    "    sf_rhs += betaU[i] * sf_dupD[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='pi_rhs'></a>\n",
    "\n",
    "# Step 3: The right-hand side of $\\partial_{t}\\Pi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{pi_rhs}$$\n",
    "\n",
    "Now let us (slightly) rearrange the RHS of $\\partial_{t}\\Pi$ so that we are able to group relevant terms together as we label them\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_{t}\\Pi &= \\underbrace{\\alpha K\\,\\Pi}_{\\text{Term 1}} + \\underbrace{\\beta^{i}\\partial_{i}\\Pi}_{\\text{Term 2}} - \\underbrace{e^{-4\\phi}\\bar{\\gamma}^{ij}\\left(\\partial_{i}\\alpha\\partial_{j}\\varphi+\\alpha\\partial_{i}\\partial_{j}\\varphi + 2\\alpha\\partial_{j}\\varphi\\partial_{i}\\phi\\right)}_{\\text{Term 3}} + \\underbrace{e^{-4\\phi}\\alpha\\bar{\\gamma}^{ij}\\bar{\\Gamma}^{k}_{ij}\\partial_{k}\\varphi}_{\\text{Term 4}}\\ .\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3a: Term 1 of $\\partial_{t}\\Pi$ = sfM_rhs : $\\alpha K\\, \\Pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.129647Z",
     "iopub.status.busy": "2021-06-15T10:05:16.129144Z",
     "iopub.status.idle": "2021-06-15T10:05:16.131067Z",
     "shell.execute_reply": "2021-06-15T10:05:16.130658Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3a: Adding Term 1 to sfM_rhs: alpha * K * Pi\n",
    "sfM_rhs = alpha * trK * sfM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3b: Term 2 of $\\partial_{t}\\Pi$ = sfM_rhs : $\\beta^{i}\\partial_{i}\\Pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.134298Z",
     "iopub.status.busy": "2021-06-15T10:05:16.133881Z",
     "iopub.status.idle": "2021-06-15T10:05:16.135957Z",
     "shell.execute_reply": "2021-06-15T10:05:16.135605Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3b: Adding Term 2 to sfM_rhs: beta^{i}\\partial_{i}Pi\n",
    "sfM_dupD = ixp.declarerank1(\"sfM_dupD\")\n",
    "for i in range(DIM):\n",
    "    sfM_rhs += betaU[i] * sfM_dupD[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3c: Term 3 of $\\partial_{t}\\Pi$ = sfM_rhs : $-e^{-4\\phi}\\bar{\\gamma}^{ij}\\left(\\partial_{i}\\alpha\\partial_{j}\\varphi+\\alpha\\partial_{i}\\partial_{j}\\varphi + 2\\alpha\\partial_{j}\\varphi\\partial_{i}\\phi\\right)$\n",
    "\n",
    "Now let's focus on Term 3:\n",
    "\n",
    "$$\n",
    "e^{-4\\phi}\\left(-\\underbrace{\\bar{\\gamma}^{ij}\\partial_{i}\\alpha\\partial_{j}\\varphi}_{\\text{Term 3a}}-\\underbrace{\\alpha\\bar{\\gamma}^{ij}\\partial_{i}\\partial_{j}\\varphi}_{\\text{Term 3b}} - \\underbrace{2\\alpha\\bar{\\gamma}^{ij}\\partial_{j}\\varphi\\partial_{i}\\phi}_{\\text{Term 3c}}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.198717Z",
     "iopub.status.busy": "2021-06-15T10:05:16.162880Z",
     "iopub.status.idle": "2021-06-15T10:05:16.221354Z",
     "shell.execute_reply": "2021-06-15T10:05:16.220950Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3c: Adding Term 3 to sfM_rhs\n",
    "# Step 3c.i: Term 3a: gammabar^{ij}\\partial_{i}\\alpha\\partial_{j}\\varphi\n",
    "alpha_dD    = ixp.declarerank1(\"alpha_dD\")\n",
    "sf_dD       = ixp.declarerank1(\"sf_dD\")\n",
    "sfMrhsTerm3 = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        sfMrhsTerm3 += - gammabarUU[i][j] * alpha_dD[i] * sf_dD[j]\n",
    "\n",
    "# Step 3c.ii: Term 3b: \\alpha*gammabar^{ij}\\partial_{i}\\partial_{j}\\varphi\n",
    "sf_dDD = ixp.declarerank2(\"sf_dDD\",\"sym01\")\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        sfMrhsTerm3 += - alpha * gammabarUU[i][j] * sf_dDD[i][j]\n",
    "\n",
    "# Step 3c.iii: Term 3c: 2*alpha*gammabar^{ij}\\partial_{j}\\varphi\\partial_{i}\\phi\n",
    "Bq.phi_and_derivs() # sets exp^{-4phi} = exp_m4phi and \\partial_{i}phi = phi_dD[i]\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        sfMrhsTerm3 += - 2 * alpha * gammabarUU[i][j] * sf_dD[j] * Bq.phi_dD[i]\n",
    "\n",
    "# Step 3c.iv: Multiplying Term 3 by e^{-4phi} and adding it to sfM_rhs\n",
    "sfMrhsTerm3 *= Bq.exp_m4phi\n",
    "sfM_rhs     += sfMrhsTerm3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3d: Term 4 of $\\partial_{t}\\Pi$ = sfM_rhs : $e^{-4\\phi}\\alpha\\bar{\\gamma}^{ij}\\bar{\\Gamma}^{k}_{\\ ij}\\partial_{k}\\varphi$\n",
    "\n",
    "We are now going to rewrite this term a bit before implementation. This rewriting is useful in order to reduce the number of finite differences approximations used to evaluate this term. Let us consider the following definitions (see equations 12a and 12b and the discussion above equation 15 in [Brown (2009)](https://arxiv.org/pdf/0902.3652.pdf))\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta\\Gamma^{k}_{\\ ij} &\\equiv \\bar\\Gamma^{k}_{\\ ij} - \\hat\\Gamma^{k}_{\\ ij}\\ ,\\\\\n",
    "\\bar\\Lambda^{k} &\\equiv \\bar\\gamma^{ij}\\Delta\\Gamma^{k}_{\\ ij}\\ .\n",
    "\\end{align}\n",
    "\n",
    "Then we have\n",
    "\n",
    "\\begin{align}\n",
    "\\text{Term 4} &= e^{-4\\phi}\\alpha\\bar{\\gamma}^{ij}\\bar{\\Gamma}^{k}_{\\ ij}\\partial_{k}\\varphi\\nonumber\\\\\n",
    "&=e^{-4\\phi}\\alpha\\bar\\gamma^{ij}\\Delta\\Gamma^{k}_{\\ ij}\\partial_{k}\\varphi + e^{-4\\phi}\\alpha\\bar{\\gamma}^{ij}\\hat{\\Gamma}^{k}_{\\ ij}\\partial_{k}\\varphi\\nonumber\\\\\n",
    "&=e^{-4\\phi}\\left(\\underbrace{\\alpha\\bar\\Lambda^{i}\\partial_{i}\\varphi}_{\\text{Term 4a}} + \\underbrace{\\alpha\\bar{\\gamma}^{ij}\\hat{\\Gamma}^{k}_{\\ ij}\\partial_{k}\\varphi}_{\\text{Term 4b}}\\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.232834Z",
     "iopub.status.busy": "2021-06-15T10:05:16.232430Z",
     "iopub.status.idle": "2021-06-15T10:05:16.234520Z",
     "shell.execute_reply": "2021-06-15T10:05:16.233994Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3d: Adding Term 4 to sfM_rhs\n",
    "# Step 3d.i: Term 4a: \\alpha \\bar\\Lambda^{i}\\partial_{i}\\varphi\n",
    "LambdabarU  = Bq.LambdabarU\n",
    "sfMrhsTerm4 = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    sfMrhsTerm4 += alpha * LambdabarU[i] * sf_dD[i]\n",
    "\n",
    "# Step 3d.ii: Evaluating \\bar\\gamma^{ij}\\hat\\Gamma^{k}_{ij}\n",
    "GammahatUDD = rfm.GammahatUDD\n",
    "gammabarGammahatContractionU = ixp.zerorank1()\n",
    "for k in range(DIM):\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            gammabarGammahatContractionU[k] += gammabarUU[i][j] * GammahatUDD[k][i][j]\n",
    "\n",
    "# Step 3d.iii: Term 4b: \\alpha \\bar\\gamma^{ij}\\hat\\Gamma^{k}_{ij}\\partial_{k}\\varphi\n",
    "for i in range(DIM):\n",
    "    sfMrhsTerm4 += alpha * gammabarGammahatContractionU[i] * sf_dD[i]\n",
    "\n",
    "# Step 3d.iii: Multplying Term 4 by e^{-4phi} and adding it to sfM_rhs\n",
    "sfMrhsTerm4 *= Bq.exp_m4phi\n",
    "sfM_rhs     += sfMrhsTerm4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against [ScalarField.ScalarField_RHSs.py](../edit/ScalarField/ScalarField_RHSs.py) NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here we perform a code validation. We verify agreement in the SymPy expressions for the RHSs of the scalar field equations between\n",
    "1. this tutorial notebook and \n",
    "2. the [ScalarField.ScalarField_RHSs.py](../edit/ScalarField/ScalarField_RHSs.py) NRPy+ module.\n",
    "\n",
    "By default, we analyze the RHSs in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.237280Z",
     "iopub.status.busy": "2021-06-15T10:05:16.236776Z",
     "iopub.status.idle": "2021-06-15T10:05:16.256084Z",
     "shell.execute_reply": "2021-06-15T10:05:16.255564Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between Scalar Field RHSs tutorial and NRPy+ module\n",
      " sf_rhs - sfrhs.sf_rhs  = 0 Should be zero.\n",
      "sfM_rhs - sfrhs.sfM_rhs = 0 Should be zero.\n"
     ]
    }
   ],
   "source": [
    "# Step 4: Code validation against NRPy+ module\n",
    "# Step 4.a: Import the ScalarFieldCollapse module and\n",
    "#           run the ScalarFieldCollapse.scalarfield_RHSs()\n",
    "#           function to evaluate the RHSs.\n",
    "import ScalarField.ScalarField_RHSs as sfrhs # NRPyCritCol: Scalar field right-hand sides\n",
    "sfrhs.ScalarField_RHSs()\n",
    "\n",
    "# Step 4.b: Perform the consistency check by subtracting\n",
    "#           the RHSs computed in this tutorial from the\n",
    "#           ones computed in the ScalarFieldCollapse module\n",
    "print(\"Consistency check between Scalar Field RHSs tutorial and NRPy+ module\")\n",
    "print(\" sf_rhs - sfrhs.sf_rhs  = \" + str( sf_rhs - sfrhs.sf_rhs )+\" Should be zero.\")\n",
    "print(\"sfM_rhs - sfrhs.sfM_rhs = \" + str(sfM_rhs - sfrhs.sfM_rhs)+\" Should be zero.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this module 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-ScalarField_RHSs.pdf](Tutorial-ScalarField_RHSs.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:05:16.258384Z",
     "iopub.status.busy": "2021-06-15T10:05:16.257883Z",
     "iopub.status.idle": "2021-06-15T10:05:18.497406Z",
     "shell.execute_reply": "2021-06-15T10:05:18.497848Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ScalarField_RHSs.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-ScalarField_RHSs.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ScalarField_RHSs\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}