Notebook Status: Validated
Validation Notes: This module has been validated as follows:
This module constructs ψ4, a quantity that is immensely useful when extracting gravitational wave content from a numerical relativity simulation. ψ4 is related to the gravitational wave strain via
ψ4=¨h+−i¨h×.We construct ψ4 from the standard ADM spatial metric γij and extrinsic curvature Kij, and their derivatives. The full expression is given by Eq. 5.1 in Baker, Campanelli, Lousto (2001):
ψ4=[Rijkl+2Ki[kKl]j]niˉmjnkˉml−8[Kj[k,l]+Γpj[kKl]p]n[0ˉmj]nkˉml+4[Rjl−KjpKpl+KKjl]n[0ˉmj]n[0ˉml],Note that ψ4 is complex, with the imaginary components originating from the tetrad vector mμ. This module does not specify a tetrad; instead, it only constructs the above expression leaving mμ and nμ unspecified. The next module on tetrads defines these tetrad quantities (currently only a quasi-Kinnersley tetrad is supported).
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 notebook).
This tutorial notebook is organized as follows
BSSN.Psi4
NRPy+ module# Step 1.a: import all needed modules from NRPy+:
import sympy as sp
import NRPy_param_funcs as par
import indexedexp as ixp
import reference_metric as rfm
# Step 1.b: Set the coordinate system for the numerical grid
# Note that this parameter is assumed to be set
# prior to calling the Python Psi4.py module,
# so this Step will not appear there.
par.set_parval_from_str("reference_metric::CoordSystem","Spherical")
# 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()
# Step 1.f: Initialize tetrad vectors.
# mre4U = $\text{Re}(m^\mu)$
# mim4U = $\text{Im}(m^\mu)$, and
# n4U = $n^\mu$
# Note that in the separate Python Psi4.py
# module, these will be set to the tetrad
# chosen within the Psi4_tetrads.py module.
# We choose the most general form for the
# tetrad vectors here instead, to ensure complete
# code validation.
mre4U = ixp.declarerank1("mre4U",DIM=4)
mim4U = ixp.declarerank1("mim4U",DIM=4)
n4U = ixp.declarerank1("n4U" ,DIM=4)
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.,
Rijkl=∂kΓijl−∂lΓijk+ΓimkΓmjl−ΓimlΓmjk,where Γijk is the Christoffel symbol associated with the 3-metric γij:
Γlij=12γlk(γki,j+γkj,i−γij,k)Notice that this equation for the Riemann tensor is equivalent to the equation given in the Wikipedia article on Formulas in Riemannian geometry:
Rℓijk=∂jΓℓik−∂kΓℓij+ΓℓjsΓsik−ΓℓksΓsij,with the replacements i→ℓ, j→i, k→j, l→k, and s→m. Wikipedia also provides a simpler form in terms of second-derivatives of three-metric itself (using the definition of the Christoffel symbol), so that we need not define derivatives of the Christoffel symbol:
Rikℓm=12(γim,kℓ+γkℓ,im−γiℓ,km−γkm,iℓ)+γnp(ΓnkℓΓpim−ΓnkmΓpiℓ).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()
gammaDDdDD = AB.gammaDDdDD
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) * \
(gammaDDdDD[i][m][k][l] + gammaDDdDD[k][l][i][m] - gammaDDdDD[i][l][k][m] - gammaDDdDD[k][m][i][l])
... then we add the term on the right:
# ... then we add the term on the right:
gammaDD = AB.gammaDD
GammaUDD = AB.GammaUDD
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 ψ4 is given by
Rijkl+2Ki[kKl]j=Rijkl+KikKlj−KilKkj# 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
rank4term1DDDD = 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):
rank4term1DDDD[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 ψ4 is given by
−8(Kj[k,l]+Γpj[kKl]p)First let's construct the first term in this sum: Kj[k,l]=12(Kjk,l−Kjl,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
rank3term2DDD = ixp.zerorank3()
KDDdD = AB.KDDdD
for j in range(DIM):
for k in range(DIM):
for l in range(DIM):
rank3term2DDD[j][k][l] = sp.Rational(1,2)*(KDDdD[j][k][l] - KDDdD[j][l][k])
... then we construct the second term in this sum: Γpj[kKl]p=12(ΓpjkKlp−ΓpjlKkp):
# ... 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):
rank3term2DDD[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):
rank3term2DDD[j][k][l] *= sp.sympify(-8)
Following Eq. 5.1 in Baker, Campanelli, Lousto (2001), the rank-2 tensor in the third term of ψ4 is given by
+4(Rjl−KjpKpl+KKjl),where Rjl=Rijil=γimRijmlK=Kii=γimKim
Let's build the components of this term: Rjl, Kpl, 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()
gammaUU = AB.gammaUU
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(Rjl−KjpKpl+KKjl),
# 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),
rank2term3DD = ixp.zerorank2()
for j in range(DIM):
for l in range(DIM):
rank2term3DD[j][l] = RDD[j][l] + trK*KDD[j][l]
for p in range(DIM):
rank2term3DD[j][l] += - KDD[j][p]*KUD[p][l]
# Finally we multiply by +4:
for j in range(DIM):
for l in range(DIM):
rank2term3DD[j][l] *= sp.sympify(4)
Eq. 5.1 in Baker, Campanelli, Lousto (2001) writes ψ4 (which is complex) as the contraction of each of the above terms with products of tetrad vectors:
ψ4=[Rijkl+2Ki[kKl]j]niˉmjnkˉml−8[Kj[k,l]+Γpj[kKl]p]n[0ˉmj]nkˉml+4[Rjl−KjpKpl+KKjl]n[0ˉmj]n[0ˉml],where ˉmμ is the complex conjugate of mμ, and nμ is real. The third term is given by n[0ˉmj]n[0ˉml]=12(n0ˉmj−njˉm0)12(n0ˉml−nlˉm0)=14(n0ˉmj−njˉm0)(n0ˉml−nlˉm0)=14(n0ˉmjn0ˉml−njˉm0n0ˉml−n0ˉmjnlˉm0+njˉm0nlˉm0)
Only mμ is complex, so we can separate the real and imaginary parts of ψ4 by hand, defining Mμ to now be the real part of ˉmμ and Mμ to be the imaginary part. All of the above products are of the form nμˉmνnηˉmδ, so let's evaluate the real and imaginary parts of this product once, for all such terms:
nμˉmνnηˉmδ=nμ(Mν−iMν)nη(Mδ−iMδ)=(nμMνnηMδ−nμMνnηMδ)+i(−nμMνnηMδ−nμMνnηMδ)# Step 6: Construct real & imaginary parts of psi_4
# by contracting constituent rank 2, 3, and 4
# tensors with input tetrads mre4U, mim4U, & 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]
# We split psi_4 into three pieces, to expedite & possibly parallelize C code generation.
psi4_re_pt = [sp.sympify(0),sp.sympify(0),sp.sympify(0)]
psi4_im_pt = [sp.sympify(0),sp.sympify(0),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_pt[0] += rank4term1DDDD[i][j][k][l] * \
tetrad_product__Real_psi4(n4U,mre4U,mim4U, i+1,j+1,k+1,l+1)
psi4_im_pt[0] += rank4term1DDDD[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_pt[1] += rank3term2DDD[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_pt[1] += rank3term2DDD[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_pt[2] += rank2term3DD[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_pt[2] += rank2term3DD[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)))
BSSN.Psi4
NRPy+ module [Back to top]¶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.
# 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.
import BSSN.Psi4 as BP4
BP4.Psi4(specify_tetrad=False)
print("Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO.")
for part in range(3):
print("psi4_im_pt["+str(part)+"] - BP4.psi4_im_pt["+str(part)+"] = " + str(psi4_im_pt[part] - BP4.psi4_im_pt[part]))
print("psi4_re_pt["+str(part)+"] - BP4.psi4_re_pt["+str(part)+"] = " + str(psi4_re_pt[part] - BP4.psi4_re_pt[part]))
Consistency check between this tutorial and BSSN.Psi4 NRPy+ module: ALL SHOULD BE ZERO. psi4_im_pt[0] - BP4.psi4_im_pt[0] = 0 psi4_re_pt[0] - BP4.psi4_re_pt[0] = 0 psi4_im_pt[1] - BP4.psi4_im_pt[1] = 0 psi4_re_pt[1] - BP4.psi4_re_pt[1] = 0 psi4_im_pt[2] - BP4.psi4_im_pt[2] = 0 psi4_re_pt[2] - BP4.psi4_re_pt[2] = 0
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 Tutorial-Psi4.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Psi4")
Created Tutorial-Psi4.tex, and compiled LaTeX file to PDF file Tutorial- Psi4.pdf