Notebook Status: Validated
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).
The module is organized as follows
We will begin by considering the Klein-Gordon equation for a massless scalar field, φ,
∇μ(∇μφ)=0 .We then define the auxiliary field
Π≡−1α(∂tφ−βi∂iφ) ,so that the Klein-Gordon equation is decomposed into two first order equations (cf. eqs. (5.252) and (5.253) in B&S)
∂tφ=βi∂iφ−αΠ ,∂tΠ=βi∂iΠ+αKΠ−γij(∂jφ∂iα−αΓk ij∂kφ+α∂i∂jφ) ,where K=γijKij is the trace of the extrinsic curvature Kij and γij the inverse of the physical spatial metric γij. We will choose not to define the auxiliary variables φi≡∂iφ in our code and, instead, leave the equations in terms of second derivatives of φ.
To write the equations in terms of BSSN variables (see this tutorial notebook for a review) we start by considering the conformal metric, ˉγij, related to the physical metric by
γij=e4ϕˉγij ,and its inverse
γij=e−4ϕˉγij .Let us also look at equation (3.7) of B&S (with i↔k and lnψ=ϕ, for convenience)
Γk ij=ˉΓk ij+2(δk iˉDjϕ+δk jˉDiϕ−ˉγijˉγkℓˉDℓϕ) .Then let us consider the term that contains Γk ij on the right-hand side of ∂tΠ:
αγijΓk ij∂kφ=αe−4ϕˉγij[ˉΓk ij+2(δk iˉDjϕ+δk jˉDiϕ−ˉγijˉγkℓˉDℓϕ)]∂kφ .Focusing on the term in parenthesis, we have (ignoring, for now, a few non-essential multiplicative terms and replacing ˉDiϕ=∂iϕ)
2ˉγij(δk iˉDjϕ+δk jˉDiϕ−ˉγijˉγkℓˉDℓϕ)∂kφ=2(ˉγkj∂jϕ+ˉγki∂iϕ−3ˉγkℓ∂ℓϕ)∂kφ=2(ˉγij∂iϕ+ˉγij∂iϕ−3ˉγij∂iϕ)∂jφ=−2ˉγij∂jφ∂iϕ ,so that
αγijΓk ij∂kφ=e−4ϕˉγij(αˉΓk ij∂kφ−2α∂jφ∂iϕ)For the rest of the equation, all we need to do is replace γij→e−4ϕˉγij, so that the Klein-Gordon equation becomes
∂tΠ=βi∂iΠ+αKΠ−e−4ϕˉγij(∂jφ∂iα−αˉΓk ij∂kφ+α∂i∂jφ+2α∂jφ∂iϕ) .Note that the evolution equation for φ is left unchanged
∂tφ=βi∂iφ−αΠ .# Step 1.a: import all needed modules from NRPy+:
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par # NRPy+: Parameter interface
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: BSSN quantities
# Step 1.b: Set the coordinate system for the numerical grid
coordsystem = "Spherical"
par.set_parval_from_str("reference_metric::CoordSystem",coordsystem)
# Step 1.c: Then we set the theta and phi axes to be the symmetry
# axes; i.e., axis "1" and "2", corresponding to the
# i1 and i2 directions.
# This sets all spatial derivatives in the theta and
# phi directions to zero (analytically).
# par.set_parval_from_str("indexedexp::symmetry_axes")
# Step 1.d: 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.e: 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.f: Import all needed BSSN quantities
Bq.BSSN_basic_tensors()
trK = Bq.trK
alpha = Bq.alpha
betaU = Bq.betaU
Bq.gammabar__inverse_and_derivs()
gammabarUU = Bq.gammabarUU
Let us label each of the terms in the RHS of the ∂tφ equation so that it is easier to understand their implementation:
∂tφ=−αΠ⏟Term 1+βi∂iφ⏟Term 2 .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.
A note on notation: We choose to declare the scalar field variable, φ, as sf and the scalar field conjugate Momentum variable, Π, as sfM to avoid possible conflicts with other variables which might be commonly denoted by ψ or Π.
# Step 2.a: Declare grid functions for varphi and Pi
sf, sfM = gri.register_gridfunctions("EVOL",["sf", "sfM"])
# Step 2.a: Add Term 1 to sf_rhs: -alpha*Pi
sf_rhs = - alpha * sfM
# Step 2.b: Add Term 2 to sf_rhs: beta^{i}\partial_{i}\varphi
sf_dupD = ixp.declarerank1("sf_dupD")
for i in range(DIM):
sf_rhs += betaU[i] * sf_dupD[i]
# Step 3a: Adding Term 1 to sfM_rhs: alpha * K * Pi
sfM_rhs = alpha * trK * sfM
# Step 3b: Adding Term 2 to sfM_rhs: beta^{i}\partial_{i}Pi
sfM_dupD = ixp.declarerank1("sfM_dupD")
for i in range(DIM):
sfM_rhs += betaU[i] * sfM_dupD[i]
Now let's focus on Term 3:
e−4ϕ(−ˉγij∂iα∂jφ⏟Term 3a−αˉγij∂i∂jφ⏟Term 3b−2αˉγij∂jφ∂iϕ⏟Term 3c)# Step 3c: Adding Term 3 to sfM_rhs
# Step 3c.i: Term 3a: gammabar^{ij}\partial_{i}\alpha\partial_{j}\varphi
alpha_dD = ixp.declarerank1("alpha_dD")
sf_dD = ixp.declarerank1("sf_dD")
sfMrhsTerm3 = sp.sympify(0)
for i in range(DIM):
for j in range(DIM):
sfMrhsTerm3 += - gammabarUU[i][j] * alpha_dD[i] * sf_dD[j]
# Step 3c.ii: Term 3b: \alpha*gammabar^{ij}\partial_{i}\partial_{j}\varphi
sf_dDD = ixp.declarerank2("sf_dDD","sym01")
for i in range(DIM):
for j in range(DIM):
sfMrhsTerm3 += - alpha * gammabarUU[i][j] * sf_dDD[i][j]
# Step 3c.iii: Term 3c: 2*alpha*gammabar^{ij}\partial_{j}\varphi\partial_{i}\phi
Bq.phi_and_derivs() # sets exp^{-4phi} = exp_m4phi and \partial_{i}phi = phi_dD[i]
for i in range(DIM):
for j in range(DIM):
sfMrhsTerm3 += - 2 * alpha * gammabarUU[i][j] * sf_dD[j] * Bq.phi_dD[i]
# Step 3c.iv: Multiplying Term 3 by e^{-4phi} and adding it to sfM_rhs
sfMrhsTerm3 *= Bq.exp_m4phi
sfM_rhs += sfMrhsTerm3
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))
ΔΓk ij≡ˉΓk ij−ˆΓk ij ,ˉΛk≡ˉγijΔΓk ij .Then we have
Term 4=e−4ϕαˉγijˉΓk ij∂kφ=e−4ϕαˉγijΔΓk ij∂kφ+e−4ϕαˉγijˆΓk ij∂kφ=e−4ϕ(αˉΛi∂iφ⏟Term 4a+αˉγijˆΓk ij∂kφ⏟Term 4b)# Step 3d: Adding Term 4 to sfM_rhs
# Step 3d.i: Term 4a: \alpha \bar\Lambda^{i}\partial_{i}\varphi
LambdabarU = Bq.LambdabarU
sfMrhsTerm4 = sp.sympify(0)
for i in range(DIM):
sfMrhsTerm4 += alpha * LambdabarU[i] * sf_dD[i]
# Step 3d.ii: Evaluating \bar\gamma^{ij}\hat\Gamma^{k}_{ij}
GammahatUDD = rfm.GammahatUDD
gammabarGammahatContractionU = ixp.zerorank1()
for k in range(DIM):
for i in range(DIM):
for j in range(DIM):
gammabarGammahatContractionU[k] += gammabarUU[i][j] * GammahatUDD[k][i][j]
# Step 3d.iii: Term 4b: \alpha \bar\gamma^{ij}\hat\Gamma^{k}_{ij}\partial_{k}\varphi
for i in range(DIM):
sfMrhsTerm4 += alpha * gammabarGammahatContractionU[i] * sf_dD[i]
# Step 3d.iii: Multplying Term 4 by e^{-4phi} and adding it to sfM_rhs
sfMrhsTerm4 *= Bq.exp_m4phi
sfM_rhs += sfMrhsTerm4
Here we perform a code validation. We verify agreement in the SymPy expressions for the RHSs of the scalar field equations between
By default, we analyze the RHSs in Spherical coordinates, though other coordinate systems may be chosen.
# Step 4: Code validation against NRPy+ module
# Step 4.a: Import the ScalarFieldCollapse module and
# run the ScalarFieldCollapse.scalarfield_RHSs()
# function to evaluate the RHSs.
import ScalarField.ScalarField_RHSs as sfrhs # NRPyCritCol: Scalar field right-hand sides
sfrhs.ScalarField_RHSs()
# Step 4.b: Perform the consistency check by subtracting
# the RHSs computed in this tutorial from the
# ones computed in the ScalarFieldCollapse module
print("Consistency check between Scalar Field RHSs tutorial and NRPy+ module")
print(" sf_rhs - sfrhs.sf_rhs = " + str( sf_rhs - sfrhs.sf_rhs )+" Should be zero.")
print("sfM_rhs - sfrhs.sfM_rhs = " + str(sfM_rhs - sfrhs.sfM_rhs)+" Should be zero.")
Consistency check between Scalar Field RHSs tutorial and NRPy+ module sf_rhs - sfrhs.sf_rhs = 0 Should be zero. sfM_rhs - sfrhs.sfM_rhs = 0 Should be zero.
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-ScalarField_RHSs.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-ScalarField_RHSs")
Created Tutorial-ScalarField_RHSs.tex, and compiled LaTeX file to PDF file Tutorial-ScalarField_RHSs.pdf