{
"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",
"# The Outgoing Gravitational Wave Weyl scalar $\\psi_4$\n",
"\n",
"## Author: Zach Etienne\n",
"\n",
"[comment]: <> (Abstract: TODO)\n",
"\n",
"**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** This module has been validated as follows:\n",
"\n",
"* Agreement (to roundoff error) with the WeylScal4 ETK thorn in Cartesian coordinates (as it agrees to roundoff error with Patrick Nelson's [Cartesian Weyl Scalars & Invariants NRPy+ tutorial notebook](Tutorial-WeylScalarsInvariants-Cartesian.ipynb), which itself was validated against WeylScal4). \n",
"* In SinhSpherical coordinates this module yields amplitude falloff and phase agreement with black hole perturbation theory for the $l=2,m=0$ (dominant, spin-weight -2 spherical harmonic) mode of a ringing Brill-Lindquist black hole remnant to more than 7 decades in amplitude, surpassing the agreement seen in Fig. 6 of [Ruchlin, Etienne, & Baumgarte](https://arxiv.org/pdf/1712.07658.pdf). (as shown in [corresponding start-to-finish notebook](Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4.ipynb); must choose EvolOption = \"high resolution\").\n",
"* Above head-on collision calculation was performed with the [Einstein Toolkit](https://einsteintoolkit.org/), and results were found to match in the case of all spherical harmonic modes.\n",
"* An equal-mass binary black hole calculation (QC-0) was performed using this module for $\\psi_4$ (SinhSpherical coordinates in region where gravitational waves extracted), and excellent agreement was observed for all (spin-weight -2 spherical harmonic) modes up to and including $l=4$, when compared with the same calculation performed with the [Einstein Toolkit](https://einsteintoolkit.org/).\n",
"\n",
"### NRPy+ Source Code for this module: [BSSN/Psi4.py](../edit/BSSN/Psi4.py)\n",
"\n",
"## Introduction:\n",
"This module constructs $\\psi_4$, 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. The [next module on tetrads defines these tetrad quantities](Tutorial-Psi4_tetrads.ipynb) (currently only a quasi-Kinnersley tetrad is supported).\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).\n"
]
},
{
"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](#riemann): Constructing the 3-Riemann tensor $R_{ik\\ell m}$\n",
"1. [Step 3](#rank4termone): Constructing the rank-4 tensor in Term 1 of $\\psi_4$: $R_{ijkl} + 2 K_{i[k} K_{l]j}$\n",
"1. [Step 4](#rank3termtwo): Constructing the rank-3 tensor in Term 2 of $\\psi_4$: $-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)$\n",
"1. [Step 5](#rank2termthree): Constructing the rank-2 tensor in term 3 of $\\psi_4$: $+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right)$\n",
"1. [Step 6](#psifour): Constructing $\\psi_4$ through contractions of the above terms with arbitrary tetrad vectors $n^\\mu$ and $m^\\mu$\n",
"1. [Step 7](#code_validation): Code Validation against `BSSN.Psi4` NRPy+ module\n",
"1. [Step 8](#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}$$\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-03-07T17:15:38.383928Z",
"iopub.status.busy": "2021-03-07T17:15:38.383074Z",
"iopub.status.idle": "2021-03-07T17:15:39.491420Z",
"shell.execute_reply": "2021-03-07T17:15:39.490857Z"
}
},
"outputs": [],
"source": [
"# Step 1.a: import all needed modules from NRPy+:\n",
"import sympy as sp\n",
"import NRPy_param_funcs as par\n",
"import indexedexp as ixp\n",
"import reference_metric as rfm\n",
"\n",
"# Step 1.b: Set the coordinate system for the numerical grid\n",
"# Note that this parameter is assumed to be set\n",
"# prior to calling the Python Psi4.py module,\n",
"# so this Step will not appear there.\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 tetrad vectors.\n",
"# mre4U = $\\text{Re}(m^\\mu)$\n",
"# mim4U = $\\text{Im}(m^\\mu)$, and\n",
"# n4U = $n^\\mu$\n",
"# Note that in the separate Python Psi4.py\n",
"# module, these will be set to the tetrad\n",
"# chosen within the Psi4_tetrads.py module.\n",
"\n",
"# We choose the most general form for the\n",
"# tetrad vectors here instead, to ensure complete\n",
"# code validation.\n",
"mre4U = ixp.declarerank1(\"mre4U\",DIM=4)\n",
"mim4U = ixp.declarerank1(\"mim4U\",DIM=4)\n",
"n4U = ixp.declarerank1(\"n4U\" ,DIM=4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='riemann'></a>\n",
"\n",
"# Step 2: Constructing the 3-Riemann tensor $R_{ik\\ell m}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{riemann}$$\n",
"\n",
"Analogously to Christoffel symbols, the Riemann tensor is a measure of the curvature of an $N$-dimensional manifold. Thus the 3-Riemann tensor is not simply a projection of the 4-Riemann tensor (see e.g., Eq. 2.7 of [Campanelli *et al* (1998)](https://arxiv.org/pdf/gr-qc/9803058.pdf) for the relation between 4-Riemann and 3-Riemann), as $N$-dimensional Riemann tensors are meant to define a notion of curvature given only the associated $N$-dimensional metric. \n",
"\n",
"So, given the ADM 3-metric, the Riemann tensor in arbitrary dimension is given by the 3-dimensional version of Eq. 1.19 in Baumgarte & Shapiro's *Numerical Relativity*. I.e.,\n",
"\n",
"$$\n",
"R^i_{jkl} = \\partial_k \\Gamma^{i}_{jl} - \\partial_l \\Gamma^{i}_{jk} + \\Gamma^i_{mk} \\Gamma^m_{jl} - \\Gamma^{i}_{ml} \\Gamma^{m}_{jk},\n",
"$$\n",
"where $\\Gamma^i_{jk}$ is the Christoffel symbol associated with the 3-metric $\\gamma_{ij}$:\n",
"\n",
"$$\n",
"\\Gamma^l_{ij} = \\frac{1}{2} \\gamma^{lk} \\left(\\gamma_{ki,j} + \\gamma_{kj,i} - \\gamma_{ij,k} \\right) \n",
"$$\n",
"\n",
"Notice that this equation for the Riemann tensor is equivalent to the equation given in the Wikipedia article on [Formulas in Riemannian geometry](https://en.wikipedia.org/w/index.php?title=List_of_formulas_in_Riemannian_geometry&oldid=882667524):\n",
"\n",
"$$\n",
"R^\\ell{}_{ijk}=\n",
"\\partial_j \\Gamma^\\ell{}_{ik}-\\partial_k\\Gamma^\\ell{}_{ij}\n",
"+\\Gamma^\\ell{}_{js}\\Gamma_{ik}^s-\\Gamma^\\ell{}_{ks}\\Gamma^s{}_{ij},\n",
"$$\n",
"with the replacements $i\\to \\ell$, $j\\to i$, $k\\to j$, $l\\to k$, and $s\\to m$. Wikipedia also provides a simpler form in terms of second-derivatives of three-metric itself (using the definition of Christoffel symbol), so that we need not define derivatives of the Christoffel symbol:\n",
"\n",
"$$\n",
"R_{ik\\ell m}=\\frac{1}{2}\\left(\n",
"\\gamma_{im,k\\ell} \n",
"+ \\gamma_{k\\ell,im}\n",
"- \\gamma_{i\\ell,km}\n",
"- \\gamma_{km,i\\ell} \\right)\n",
"+\\gamma_{np} \\left(\n",
"\\Gamma^n{}_{k\\ell} \\Gamma^p{}_{im} - \n",
"\\Gamma^n{}_{km} \\Gamma^p{}_{i\\ell} \\right).\n",
"$$\n",
"\n",
"First we construct the term on the left:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.535770Z",
"iopub.status.busy": "2021-03-07T17:15:39.530217Z",
"iopub.status.idle": "2021-03-07T17:15:39.538210Z",
"shell.execute_reply": "2021-03-07T17:15:39.538758Z"
}
},
"outputs": [],
"source": [
"# Step 2: Construct the (rank-4) Riemann curvature tensor associated with the ADM 3-metric:\n",
"RDDDD = ixp.zerorank4()\n",
"gammaDDdDD = AB.gammaDDdDD\n",
"\n",
"for i in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" for m in range(DIM):\n",
" RDDDD[i][k][l][m] = sp.Rational(1,2) * \\\n",
" (gammaDDdDD[i][m][k][l] + gammaDDdDD[k][l][i][m] - gammaDDdDD[i][l][k][m] - gammaDDdDD[k][m][i][l])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"... then we add the term on the right:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.608168Z",
"iopub.status.busy": "2021-03-07T17:15:39.572029Z",
"iopub.status.idle": "2021-03-07T17:15:39.816676Z",
"shell.execute_reply": "2021-03-07T17:15:39.815987Z"
}
},
"outputs": [],
"source": [
"# ... then we add the term on the right:\n",
"gammaDD = AB.gammaDD\n",
"GammaUDD = AB.GammaUDD\n",
"\n",
"for i in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" for m in range(DIM):\n",
" for n in range(DIM):\n",
" for p in range(DIM):\n",
" RDDDD[i][k][l][m] += gammaDD[n][p] * \\\n",
" (GammaUDD[n][k][l]*GammaUDD[p][i][m] - GammaUDD[n][k][m]*GammaUDD[p][i][l])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='rank4termone'></a>\n",
"\n",
"# Step 3: Constructing the rank-4 tensor in Term 1 of $\\psi_4$: $R_{ijkl} + 2 K_{i[k} K_{l]j}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{rank4termone}$$\n",
"\n",
"Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-4 tensor in the first term of $\\psi_4$ is given by\n",
"\n",
"$$\n",
"R_{ijkl} + 2 K_{i[k} K_{l]j} = R_{ijkl} + K_{ik} K_{lj} - K_{il} K_{kj}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.857250Z",
"iopub.status.busy": "2021-03-07T17:15:39.855280Z",
"iopub.status.idle": "2021-03-07T17:15:39.859439Z",
"shell.execute_reply": "2021-03-07T17:15:39.859932Z"
}
},
"outputs": [],
"source": [
"# Step 3: Construct the (rank-4) tensor in term 1 of psi_4 (referring to Eq 5.1 in\n",
"# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
"rank4term1DDDD = ixp.zerorank4()\n",
"KDD = AB.KDD\n",
"\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" rank4term1DDDD[i][j][k][l] = RDDDD[i][j][k][l] + KDD[i][k]*KDD[l][j] - KDD[i][l]*KDD[k][j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='rank3termtwo'></a>\n",
"\n",
"# Step 4: Constructing the rank-3 tensor in Term 2 of $\\psi_4$: $-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)$ \\[Back to [top](#toc)\\]\n",
"$$\\label{rank3termtwo}$$\n",
"\n",
"Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-3 tensor in the second term of $\\psi_4$ is given by\n",
"\n",
"$$\n",
"-8 \\left(K_{j[k,l]} + \\Gamma^{p}_{j[k} K_{l]p} \\right)\n",
"$$\n",
"First let's construct the first term in this sum: $K_{j[k,l]} = \\frac{1}{2} (K_{jk,l} - K_{jl,k})$:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.876863Z",
"iopub.status.busy": "2021-03-07T17:15:39.876237Z",
"iopub.status.idle": "2021-03-07T17:15:39.878536Z",
"shell.execute_reply": "2021-03-07T17:15:39.879036Z"
}
},
"outputs": [],
"source": [
"# Step 4: Construct the (rank-3) tensor in term 2 of psi_4 (referring to Eq 5.1 in\n",
"# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
"rank3term2DDD = ixp.zerorank3()\n",
"KDDdD = AB.KDDdD\n",
"\n",
"for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" rank3term2DDD[j][k][l] = sp.Rational(1,2)*(KDDdD[j][k][l] - KDDdD[j][l][k])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"... then we construct the second term in this sum: $\\Gamma^{p}_{j[k} K_{l]p} = \\frac{1}{2} (\\Gamma^{p}_{jk} K_{lp}-\\Gamma^{p}_{jl} K_{kp})$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.913905Z",
"iopub.status.busy": "2021-03-07T17:15:39.902213Z",
"iopub.status.idle": "2021-03-07T17:15:39.916163Z",
"shell.execute_reply": "2021-03-07T17:15:39.916634Z"
}
},
"outputs": [],
"source": [
"# ... then we construct the second term in this sum:\n",
"# \\Gamma^{p}_{j[k} K_{l]p} = \\frac{1}{2} (\\Gamma^{p}_{jk} K_{lp}-\\Gamma^{p}_{jl} K_{kp}):\n",
"for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" for p in range(DIM):\n",
" rank3term2DDD[j][k][l] += sp.Rational(1,2)*(GammaUDD[p][j][k]*KDD[l][p] - GammaUDD[p][j][l]*KDD[k][p])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we multiply the term by $-8$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.928985Z",
"iopub.status.busy": "2021-03-07T17:15:39.927994Z",
"iopub.status.idle": "2021-03-07T17:15:39.930690Z",
"shell.execute_reply": "2021-03-07T17:15:39.931248Z"
}
},
"outputs": [],
"source": [
"# Finally, we multiply the term by $-8$:\n",
"for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" rank3term2DDD[j][k][l] *= sp.sympify(-8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='rank2termthree'></a>\n",
"\n",
"# Step 5: Constructing the rank-2 tensor in term 3 of $\\psi_4$: $+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right)$ \\[Back to [top](#toc)\\]\n",
"$$\\label{rank2termthree}$$\n",
"\n",
"Following Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf), the rank-2 tensor in the third term of $\\psi_4$ is given by\n",
"\n",
"$$\n",
"+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
"$$\n",
"where\n",
"\\begin{align}\n",
"R_{jl} &= R^i_{jil} \\\\\n",
"&= \\gamma^{im} R_{ijml} \\\\\n",
"K &= K^i_i \\\\\n",
"&= \\gamma^{im} K_{im}\n",
"\\end{align}\n",
"\n",
"Let's build the components of this term: $R_{jl}$, $K^p_l$, and $K$, as defined above:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.957772Z",
"iopub.status.busy": "2021-03-07T17:15:39.954738Z",
"iopub.status.idle": "2021-03-07T17:15:39.959888Z",
"shell.execute_reply": "2021-03-07T17:15:39.960418Z"
}
},
"outputs": [],
"source": [
"# Step 5: Construct the (rank-2) tensor in term 3 of psi_4 (referring to Eq 5.1 in\n",
"# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf\n",
"\n",
"# Step 5.1: Construct 3-Ricci tensor R_{ij} = gamma^{im} R_{ijml}\n",
"RDD = ixp.zerorank2()\n",
"gammaUU = AB.gammaUU\n",
"for j in range(DIM):\n",
" for l in range(DIM):\n",
" for i in range(DIM):\n",
" for m in range(DIM):\n",
" RDD[j][l] += gammaUU[i][m]*RDDDD[i][j][m][l]\n",
"\n",
"# Step 5.2: Construct K^p_l = gamma^{pi} K_{il}\n",
"KUD = ixp.zerorank2()\n",
"for p in range(DIM):\n",
" for l in range(DIM):\n",
" for i in range(DIM):\n",
" KUD[p][l] += gammaUU[p][i]*KDD[i][l]\n",
"\n",
"# Step 5.3: Construct trK = gamma^{ij} K_{ij}\n",
"trK = sp.sympify(0)\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" trK += gammaUU[i][j]*KDD[i][j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we put these terms together to construct the entire term:\n",
"$$\n",
"+4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:39.991138Z",
"iopub.status.busy": "2021-03-07T17:15:39.990507Z",
"iopub.status.idle": "2021-03-07T17:15:39.993355Z",
"shell.execute_reply": "2021-03-07T17:15:39.992768Z"
}
},
"outputs": [],
"source": [
"# Next we put these terms together to construct the entire term in parentheses:\n",
"# +4 \\left(R_{jl} - K_{jp} K^p_l + K K_{jl} \\right),\n",
"rank2term3DD = ixp.zerorank2()\n",
"for j in range(DIM):\n",
" for l in range(DIM):\n",
" rank2term3DD[j][l] = RDD[j][l] + trK*KDD[j][l]\n",
" for p in range(DIM):\n",
" rank2term3DD[j][l] += - KDD[j][p]*KUD[p][l]\n",
"# Finally we multiply by +4:\n",
"for j in range(DIM):\n",
" for l in range(DIM):\n",
" rank2term3DD[j][l] *= sp.sympify(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='psifour'></a>\n",
"\n",
"# Step 6: Constructing $\\psi_4$ through contractions of the above terms with an arbitrary tetrad vectors $m^\\mu$ and $n^\\mu$ \\[Back to [top](#toc)\\]\n",
"$$\\label{psifour}$$\n",
"\n",
"Eq. 5.1 in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf) writes $\\psi_4$ (which is complex) as the contraction of each of the above terms with products of tetrad vectors:\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",
"where $\\bar{m}^\\mu$ is the complex conjugate of $m^\\mu$, and $n^\\mu$ is real. The third term is given by\n",
"\\begin{align}\n",
"{n}^{[0}\\bar{m}^{j]}{n}^{[0}\\bar{m}^{l]}\n",
"&= \\frac{1}{2}({n}^{0}\\bar{m}^{j} - {n}^{j}\\bar{m}^{0} )\\frac{1}{2}({n}^{0}\\bar{m}^{l} - {n}^{l}\\bar{m}^{0} )\\\\\n",
"&= \\frac{1}{4}({n}^{0}\\bar{m}^{j} - {n}^{j}\\bar{m}^{0} )({n}^{0}\\bar{m}^{l} - {n}^{l}\\bar{m}^{0} )\\\\\n",
"&= \\frac{1}{4}({n}^{0}\\bar{m}^{j}{n}^{0}\\bar{m}^{l} - {n}^{j}\\bar{m}^{0}{n}^{0}\\bar{m}^{l} - {n}^{0}\\bar{m}^{j}{n}^{l}\\bar{m}^{0} + {n}^{j}\\bar{m}^{0}{n}^{l}\\bar{m}^{0})\n",
"\\end{align}\n",
"\n",
"Only $m^\\mu$ is complex, so we can separate the real and imaginary parts of $\\psi_4$ by hand, defining $M^\\mu$ to now be the real part of $\\bar{m}^\\mu$ and $\\mathcal{M}^\\mu$ to be the imaginary part. All of the above products are of the form ${n}^\\mu\\bar{m}^\\nu{n}^\\eta\\bar{m}^\\delta$, so let's evalute the real and imaginary parts of this product once, for all such terms:\n",
"\n",
"\\begin{align}\n",
"{n}^\\mu\\bar{m}^\\nu{n}^\\eta\\bar{m}^\\delta\n",
"&= {n}^\\mu(M^\\nu - i \\mathcal{M}^\\nu){n}^\\eta(M^\\delta - i \\mathcal{M}^\\delta) \\\\\n",
"&= \\left({n}^\\mu M^\\nu {n}^\\eta M^\\delta -\n",
"{n}^\\mu \\mathcal{M}^\\nu {n}^\\eta \\mathcal{M}^\\delta \\right)+\n",
"i \\left(\n",
"-{n}^\\mu M^\\nu {n}^\\eta \\mathcal{M}^\\delta\n",
"-{n}^\\mu \\mathcal{M}^\\nu {n}^\\eta M^\\delta\n",
"\\right)\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:40.018914Z",
"iopub.status.busy": "2021-03-07T17:15:40.017595Z",
"iopub.status.idle": "2021-03-07T17:15:41.585762Z",
"shell.execute_reply": "2021-03-07T17:15:41.585203Z"
}
},
"outputs": [],
"source": [
"# Step 6: Construct real & imaginary parts of psi_4\n",
"# by contracting constituent rank 2, 3, and 4\n",
"# tensors with input tetrads mre4U, mim4U, & n4U.\n",
"\n",
"def tetrad_product__Real_psi4(n,Mre,Mim, mu,nu,eta,delta):\n",
" return +n[mu]*Mre[nu]*n[eta]*Mre[delta] - n[mu]*Mim[nu]*n[eta]*Mim[delta]\n",
"\n",
"def tetrad_product__Imag_psi4(n,Mre,Mim, mu,nu,eta,delta):\n",
" return -n[mu]*Mre[nu]*n[eta]*Mim[delta] - n[mu]*Mim[nu]*n[eta]*Mre[delta]\n",
"\n",
"\n",
"# We split psi_4 into three pieces, to expedite & possibly parallelize C code generation.\n",
"psi4_re_pt = [sp.sympify(0),sp.sympify(0),sp.sympify(0)]\n",
"psi4_im_pt = [sp.sympify(0),sp.sympify(0),sp.sympify(0)]\n",
"\n",
"# First term:\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" psi4_re_pt[0] += rank4term1DDDD[i][j][k][l] * \\\n",
" tetrad_product__Real_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)\n",
" psi4_im_pt[0] += rank4term1DDDD[i][j][k][l] * \\\n",
" tetrad_product__Imag_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)\n",
"\n",
"# Second term:\n",
"for j in range(DIM):\n",
" for k in range(DIM):\n",
" for l in range(DIM):\n",
" psi4_re_pt[1] += rank3term2DDD[j][k][l] * \\\n",
" sp.Rational(1,2)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)\n",
" -tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )\n",
" psi4_im_pt[1] += rank3term2DDD[j][k][l] * \\\n",
" sp.Rational(1,2)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)\n",
" -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )\n",
"# Third term:\n",
"for j in range(DIM):\n",
" for l in range(DIM):\n",
" psi4_re_pt[2] += rank2term3DD[j][l] * \\\n",
" (sp.Rational(1,4)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)\n",
" -tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)\n",
" -tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)\n",
" +tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))\n",
" psi4_im_pt[2] += rank2term3DD[j][l] * \\\n",
" (sp.Rational(1,4)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)\n",
" -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)\n",
" -tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)\n",
" +tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 7: Code validation against `BSSN.Psi4` 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](../edit/BSSN/Psi4.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": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:41.591760Z",
"iopub.status.busy": "2021-03-07T17:15:41.591074Z",
"iopub.status.idle": "2021-03-07T17:15:45.738251Z",
"shell.execute_reply": "2021-03-07T17:15:45.737535Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO.\n",
"psi4_im_pt[0] - BP4.psi4_im_pt[0] = 0\n",
"psi4_re_pt[0] - BP4.psi4_re_pt[0] = 0\n",
"psi4_im_pt[1] - BP4.psi4_im_pt[1] = 0\n",
"psi4_re_pt[1] - BP4.psi4_re_pt[1] = 0\n",
"psi4_im_pt[2] - BP4.psi4_im_pt[2] = 0\n",
"psi4_re_pt[2] - BP4.psi4_re_pt[2] = 0\n"
]
}
],
"source": [
"# Call the BSSN_RHSs() function from within the\n",
"# BSSN/BSSN_RHSs.py module,\n",
"# which should do exactly the same as in Steps 1-16 above.\n",
"import BSSN.Psi4 as BP4\n",
"BP4.Psi4(specify_tetrad=False)\n",
"\n",
"print(\"Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO.\")\n",
"\n",
"for part in range(3):\n",
" print(\"psi4_im_pt[\"+str(part)+\"] - BP4.psi4_im_pt[\"+str(part)+\"] = \" + str(psi4_im_pt[part] - BP4.psi4_im_pt[part]))\n",
" print(\"psi4_re_pt[\"+str(part)+\"] - BP4.psi4_re_pt[\"+str(part)+\"] = \" + str(psi4_re_pt[part] - BP4.psi4_re_pt[part]))"
]
},
{
"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-Psi4.pdf](Tutorial-Psi4.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": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:15:45.742903Z",
"iopub.status.busy": "2021-03-07T17:15:45.742209Z",
"iopub.status.idle": "2021-03-07T17:15:49.435590Z",
"shell.execute_reply": "2021-03-07T17:15:49.436331Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-Psi4.tex, and compiled LaTeX file to PDF file Tutorial-\n",
" Psi4.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\")"
]
}
],
"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.10.0rc2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}