This module exists as a modification of the NRPy+ $\psi_4$ in curvilinear coordinates module, writing all spacetime quantities in terms of ADM variables and their derivatives directly.
As is standard in NRPy+,
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 module).
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
$$ \psi_4 = \ddot{h}_+ - i \ddot{h}_\times. $$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):
\begin{align} \psi_4 &= \left[ {R}_{ijkl}+2K_{i[k}K_{l]j}\right] {n}^i\bar{m}^j{n}^k\bar{m}^l \\ & -8\left[ K_{j[k,l]}+{\Gamma }_{j[k}^pK_{l]p}\right] {n}^{[0}\bar{m}^{j]}{n}^k\bar{m}^l \\ & +4\left[ {R}_{jl}-K_{jp}K_l^p+KK_{jl}\right] {n}^{[0}\bar{m}^{j]}{n}^{[0}\bar{m}^{l]}, \end{align}Note that $\psi_4$ is imaginary, 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 (currently only a quasi-Kinnersley tetrad is supported).
This tutorial module is organized as follows
BSSN.Psi4
NRPy+ module# Step 1.a: import all needed modules from NRPy+:
import shutil, os
import sys#TylerK: Add sys to get cmdline_helper from NRPy top directory; remove this line and next when debugged
sys.path.append('../../')
import sympy as sp
from outputC import *
import NRPy_param_funcs as par
import indexedexp as ixp
import grid as gri
import finite_difference as fin
import reference_metric as rfm
# Step 1.b: Set the coordinate system for the numerical grid
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
# Step 1.c: Given the chosen coordinate system, set up
# corresponding reference metric and needed
# reference metric quantities
# The following function call sets up the reference metric
# and related quantities, including rescaling matrices ReDD,
# ReU, and hatted quantities.
rfm.reference_metric()
# Step 1.d: Set spatial dimension (must be 3 for BSSN, as BSSN is
# a 3+1-dimensional decomposition of the general
# relativistic field equations)
DIM = 3
# Step 1.e: Import all ADM quantities as written in terms of BSSN quantities
# import BSSN.ADM_in_terms_of_BSSN as AB
# AB.ADM_in_terms_of_BSSN()
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) 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.
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.,
$$ 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}, $$where $\Gamma^i_{jk}$ is the Christoffel symbol associated with the 3-metric $\gamma_{ij}$:
$$ \Gamma^l_{ij} = \frac{1}{2} \gamma^{lk} \left(\gamma_{ki,j} + \gamma_{kj,i} - \gamma_{ij,k} \right) $$Notice that this equation for the Riemann tensor is equivalent to the equation given in the Wikipedia article on Formulas in Riemannian geometry:
$$ R^\ell{}_{ijk}= \partial_j \Gamma^\ell{}_{ik}-\partial_k\Gamma^\ell{}_{ij} +\Gamma^\ell{}_{js}\Gamma_{ik}^s-\Gamma^\ell{}_{ks}\Gamma^s{}_{ij}, $$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:
$$ R_{ik\ell m}=\frac{1}{2}\left( \gamma_{im,k\ell} + \gamma_{k\ell,im} - \gamma_{i\ell,km} - \gamma_{km,i\ell} \right) +\gamma_{np} \left( \Gamma^n{}_{k\ell} \Gamma^p{}_{im} - \Gamma^n{}_{km} \Gamma^p{}_{i\ell} \right). $$First, we construct the term on the left:
# Step 2: Construct the (rank-4) Riemann curvature tensor associated with the ADM 3-metric:
RDDDD = ixp.zerorank4()
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUX","gammaDD", "sym01") # The AUX or EVOL designation is *not*
# used in diagnostic modules.
kDD = ixp.register_gridfunctions_for_single_rank2("AUX","kDD", "sym01")
gammaDD_dD = ixp.declarerank3("gammaDD_dD","sym01")
gammaDD_dDD = ixp.declarerank4("gammaDD_dDD","sym01_sym23")
# gammaDD_dDD = AB.gammaDD_dDD
for i in range(DIM):
for k in range(DIM):
for l in range(DIM):
for m in range(DIM):
RDDDD[i][k][l][m] = sp.Rational(1,2) * \
(gammaDD_dDD[i][m][k][l] + gammaDD_dDD[k][l][i][m] - gammaDD_dDD[i][l][k][m] - gammaDD_dDD[k][m][i][l])
... then we add the term on the right:
# ... then we add the term on the right:
# Define the Christoffel symbols
GammaUDD = ixp.zerorank3(DIM)
gammaUU,gammadetdummy = ixp.symm_matrix_inverter3x3(gammaDD)
for i in range(DIM):
for k in range(DIM):
for l in range(DIM):
for m in range(DIM):
GammaUDD[i][k][l] += (sp.Rational(1,2))*gammaUU[i][m]*\
(gammaDD_dD[m][k][l] + gammaDD_dD[m][l][k] - gammaDD_dD[k][l][m])
for i in range(DIM):
for k in range(DIM):
for l in range(DIM):
for m in range(DIM):
for n in range(DIM):
for p in range(DIM):
RDDDD[i][k][l][m] += gammaDD[n][p] * \
(GammaUDD[n][k][l]*GammaUDD[p][i][m] - GammaUDD[n][k][m]*GammaUDD[p][i][l])
Following Eq. 5.1 in Baker, Campanelli, Lousto (2001), the rank-4 tensor in the first term of $\psi_4$ is given by
$$ R_{ijkl} + 2 K_{i[k} K_{l]j} = R_{ijkl} + K_{ik} K_{lj} - K_{il} K_{kj} $$# Step 3: Construct the (rank-4) tensor in term 1 of psi_4 (referring to Eq 5.1 in
# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
rank4term1 = ixp.zerorank4()
# kDD = AB.kDD
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
rank4term1[i][j][k][l] = RDDDD[i][j][k][l] + kDD[i][k]*kDD[l][j] - kDD[i][l]*kDD[k][j]
Following Eq. 5.1 in Baker, Campanelli, Lousto (2001), the rank-3 tensor in the second term of $\psi_4$ is given by
$$ -8 \left(K_{j[k,l]} + \Gamma^{p}_{j[k} K_{l]p} \right) $$First let's construct the first term in this sum: $K_{j[k,l]} = \frac{1}{2} (K_{jk,l} - K_{jl,k})$:
# Step 4: Construct the (rank-3) tensor in term 2 of psi_4 (referring to Eq 5.1 in
# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
rank3term2 = ixp.zerorank3()
# kDD_dD = AB.kDD_dD
kDD_dD = ixp.declarerank3("kDD_dD","sym01")
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
rank3term2[j][k][l] = sp.Rational(1,2)*(kDD_dD[j][k][l] - kDD_dD[j][l][k])
... 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})$:
# ... 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}):
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
for p in range(DIM):
rank3term2[j][k][l] += sp.Rational(1,2)*(GammaUDD[p][j][k]*kDD[l][p] - GammaUDD[p][j][l]*kDD[k][p])
Finally, we multiply the term by $-8$:
# Finally, we multiply the term by $-8$:
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
rank3term2[j][k][l] *= sp.sympify(-8)
Following Eq. 5.1 in Baker, Campanelli, Lousto (2001), the rank-2 tensor in the third term of $\psi_4$ is given by
$$ +4 \left(R_{jl} - K_{jp} K^p_l + K K_{jl} \right), $$where \begin{align} R_{jl} &= R^i_{jil} \\ &= \gamma^{im} R_{ijml} \\ K &= K^i_i \\ &= \gamma^{im} K_{im} \end{align}
Let's build the components of this term: $R_{jl}$, $K^p_l$, and $K$, as defined above:
# Step 5: Construct the (rank-2) tensor in term 3 of psi_4 (referring to Eq 5.1 in
# Baker, Campanelli, Lousto (2001); https://arxiv.org/pdf/gr-qc/0104063.pdf
# Step 5.1: Construct 3-Ricci tensor R_{ij} = gamma^{im} R_{ijml}
RDD = ixp.zerorank2()
for j in range(DIM):
for l in range(DIM):
for i in range(DIM):
for m in range(DIM):
RDD[j][l] += gammaUU[i][m]*RDDDD[i][j][m][l]
# Step 5.2: Construct K^p_l = gamma^{pi} K_{il}
KUD = ixp.zerorank2()
for p in range(DIM):
for l in range(DIM):
for i in range(DIM):
KUD[p][l] += gammaUU[p][i]*kDD[i][l]
# Step 5.3: Construct trK = gamma^{ij} K_{ij}
trK = sp.sympify(0)
for i in range(DIM):
for j in range(DIM):
trK += gammaUU[i][j]*kDD[i][j]
Next we put these terms together to construct the entire term: $$ +4 \left(R_{jl} - K_{jp} K^p_l + K K_{jl} \right), $$
# Next we put these terms together to construct the entire term in parentheses:
# +4 \left(R_{jl} - K_{jp} K^p_l + K K_{jl} \right),
rank2term3 = ixp.zerorank2()
for j in range(DIM):
for l in range(DIM):
rank2term3[j][l] = RDD[j][l] + trK*kDD[j][l]
for p in range(DIM):
rank2term3[j][l] += - kDD[j][p]*KUD[p][l]
# Finally we multiply by +4:
for j in range(DIM):
for l in range(DIM):
rank2term3[j][l] *= sp.sympify(4)
Eq. 5.1 in Baker, Campanelli, Lousto (2001) writes $\psi_4$ (which is complex) as the contraction of each of the above terms with products of tetrad vectors:
\begin{align} \psi_4 &= \left[ {R}_{ijkl}+2K_{i[k}K_{l]j}\right] {n}^i\bar{m}^j{n}^k\bar{m}^l \\ & -8\left[ K_{j[k,l]}+{\Gamma }_{j[k}^pK_{l]p}\right] {n}^{[0}\bar{m}^{j]}{n}^k\bar{m}^l \\ & +4\left[ {R}_{jl}-K_{jp}K_l^p+KK_{jl}\right] {n}^{[0}\bar{m}^{j]}{n}^{[0}\bar{m}^{l]}, \end{align}where $\bar{m}^\mu$ is the complex conjugate of $m^\mu$, and $n^\mu$ is real. The third term is given by \begin{align} {n}^{[0}\bar{m}^{j]}{n}^{[0}\bar{m}^{l]} &= \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} )\\ &= \frac{1}{4}({n}^{0}\bar{m}^{j} - {n}^{j}\bar{m}^{0} )({n}^{0}\bar{m}^{l} - {n}^{l}\bar{m}^{0} )\\ &= \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}) \end{align}
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 evaluate the real and imaginary parts of this product once, for all such terms:
\begin{align} {n}^\mu\bar{m}^\nu{n}^\eta\bar{m}^\delta &= {n}^\mu(M^\nu - i \mathcal{M}^\nu){n}^\eta(M^\delta - i \mathcal{M}^\delta) \\ &= \left({n}^\mu M^\nu {n}^\eta M^\delta - {n}^\mu \mathcal{M}^\nu {n}^\eta \mathcal{M}^\delta \right)+ i \left( -{n}^\mu M^\nu {n}^\eta \mathcal{M}^\delta -{n}^\mu \mathcal{M}^\nu {n}^\eta M^\delta \right) \end{align}# mre4U = ixp.declarerank1("mre4U",DIM=4)
# mim4U = ixp.declarerank1("mim4U",DIM=4)
# n4U = ixp.declarerank1("n4U" ,DIM=4)
import BSSN.Psi4_tetrads as P4t
P4t.Psi4_tetrads()
mre4U = P4t.mre4U
mim4U = P4t.mim4U
n4U = P4t.n4U
def tetrad_product__Real_psi4(n,Mre,Mim, mu,nu,eta,delta):
return +n[mu]*Mre[nu]*n[eta]*Mre[delta] - n[mu]*Mim[nu]*n[eta]*Mim[delta]
def tetrad_product__Imag_psi4(n,Mre,Mim, mu,nu,eta,delta):
return -n[mu]*Mre[nu]*n[eta]*Mim[delta] - n[mu]*Mim[nu]*n[eta]*Mre[delta]
psi4_re = sp.sympify(0)
psi4_im = sp.sympify(0)
# First term:
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
psi4_re += rank4term1[i][j][k][l]*tetrad_product__Real_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)
psi4_im += rank4term1[i][j][k][l]*tetrad_product__Imag_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)
# Second term:
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
psi4_re += rank3term2[j][k][l] * \
sp.Rational(1,2)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)
-tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )
psi4_im += rank3term2[j][k][l] * \
sp.Rational(1,2)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,k+1,l+1)
-tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,k+1,l+1) )
# Third term:
for j in range(DIM):
for l in range(DIM):
psi4_re += rank2term3[j][l] * \
(sp.Rational(1,4)*(+tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)
-tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)
-tetrad_product__Real_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)
+tetrad_product__Real_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))
psi4_im += rank2term3[j][l] * \
(sp.Rational(1,4)*(+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,0,l+1)
-tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,0,l+1)
-tetrad_product__Imag_psi4(n4U,mre4U,mim4U, 0,j+1,l+1,0)
+tetrad_product__Imag_psi4(n4U,mre4U,mim4U, j+1,0,l+1,0)))
As a code validation check, we verify agreement in the SymPy expressions for the RHSs of the BSSN equations between
By default, we compare all quantities in Spherical coordinates, though other coordinate systems may be chosen.
outCparams = "preindent=1,outCfileaccess=w,outCverbose=False,includebraces=False"
print("STARTING NEW")
fin.FD_outputC("Psi4_new.h", lhrh(lhs="psi4_real", rhs=psi4_re), outCparams)
print("FINISHED NEW")
gri.glb_gridfcs_list = []
import WeylScal4NRPy.WeylScalars_Cartesian as W4
W4.WeylScalars_Cartesian()
print("STARTING OLD")
fin.FD_outputC("Psi4_old.h", lhrh(lhs="psi4_real", rhs=W4.psi4r), outCparams)
print("FINISHED OLD")
# print("FullSimplify["+str(sp.mathematica_code(psi4_re-W4.psi4r))+"]")
# with open("math.txt","w") as file:
# file.write("FullSimplify["+str(sp.mathematica_code(psi4_re-W4.psi4r))+"]")
# # Call the BSSN_RHSs() function from within the
# # BSSN/BSSN_RHSs.py module,
# # which should do exactly the same as in Steps 1-16 above.
# print("vvv Ignore the minor warnings below. vvv")
# import BSSN.Psi4 as BP4
# BP4.Psi4()
# print("^^^ Ignore the minor warnings above. ^^^\n")
# print("Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO.")
# print("psi4_im - BP4.psi4_im = " + str(psi4_im - BP4.psi4_im))
# print("psi4_re - BP4.psi4_re = " + str(psi4_re - BP4.psi4_re))
STARTING NEW Wrote to file "Psi4_new.h" FINISHED NEW STARTING OLD Wrote to file "Psi4_old.h" FINISHED OLD
!gcc -O2 psi4_tester.c -o psi4_tester -lm
!./psi4_tester 4 4 4
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Psi4-Cartesian_validation")
This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex) restricted \write18 enabled. entering extended mode This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex) restricted \write18 enabled. entering extended mode This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex) restricted \write18 enabled. entering extended mode