Notebook Status: Self-Validated
Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. Additional validation tests may have been performed, but are as yet, undocumented. (TODO)
Given the ADM variables:
$$\left\{\gamma_{ij}, K_{ij}, \alpha, \beta^i, B^i\right\}$$in the Spherical or Cartesian basis, and as functions of $(r,\theta,\phi)$ or $(x,y,z)$, respectively, this module documents their conversion to the BSSN variables
$$\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\phi, K, \bar{\Lambda}^{i}, \alpha, \beta^i, B^i\right\},$$in the desired curvilinear basis (given by reference_metric::CoordSystem). Then it rescales the resulting BSSNCurvilinear variables (as defined in the BSSN Curvilinear tutorial) into the form needed for BSSNCurvilinear evolutions:
We will use as our core example in this module UIUC initial data, which are (as documented in their NRPy+ initial data module) given in terms of ADM variables in Spherical coordinates.
This notebook is organized as follows
xx0,xx1,xx2 basisBSSN.BSSN_in_terms_of_ADM (tutorial) to perform the ADM-to-BSSN conversion for all BSSN quantities except $\lambda^i$BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear NRPy+ module# Step 1: Initialize core Python/NRPy+ modules
from outputC import outCfunction,outputC,lhrh # NRPy+: Core C code output module
import NRPy_param_funcs as par # NRPy+: Parameter interface
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin # NRPy+: Finite difference C code generation module
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+: Computes useful BSSN quantities; e.g., gammabarUU & GammabarUDD needed below
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
import os, shutil, sys # Standard Python modules for multiplatform OS-level functions
# Step 1.a: Create output directory for C codes generated by this tutorial:
Ccodesdir = os.path.join("numerical_ADM_to_BSSN_Ccodes/")
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
# !rm -r ScalarWaveCurvilinear_Playground_Ccodes
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)
# Step 1.b: Create output directory for C codes generated by the corresponding NRPy+ Python module,
# BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear
PyModCcodesdir = os.path.join("numerical_ADM_to_BSSN_Ccodes/PyMod/")
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
# !rm -r ScalarWaveCurvilinear_Playground_Ccodes
shutil.rmtree(PyModCcodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(PyModCcodesdir)
# Step 2: Desired output BSSN Curvilinear coordinate system set to Cylindrical, as a proof-of-principle
# The ADM & BSSN formalisms only work in 3D; they are 3+1 decompositions of Einstein's equations.
# To implement axisymmetry or spherical symmetry, simply set all spatial derivatives in
# the relevant angular directions to zero; DO NOT SET DIM TO ANYTHING BUT 3.
# Set spatial dimension (must be 3 for BSSN)
DIM = 3
# Set the desired *output* coordinate system to Cylindrical:
par.set_parval_from_str("reference_metric::CoordSystem","Cylindrical")
rfm.reference_metric()
# Set function input parameters to consistent defaults.
ADM_input_function_name = "ID_ADM_SphorCart"
pointer_to_ID_inputs = False
ADM variables are given as functions of $(r,\theta,\phi)$ or $(x,y,z)$. We convert them to functions of (xx0,xx1,xx2) using SymPy's subs() function.
# Step 1: All input quantities are in terms of r,th,ph or x,y,z. We want them in terms
# of xx0,xx1,xx2, so here we call sympify_integers__replace_rthph() to replace
# r,th,ph or x,y,z, respectively, with the appropriate functions of xx0,xx1,xx2
# as defined for this particular reference metric in reference_metric.py's
# xxSph[] or xx_to_Cart[], respectively:
# Define the input variables:
gammaSphorCartDD = ixp.declarerank2("gammaSphorCartDD","sym01")
KSphorCartDD = ixp.declarerank2("KSphorCartDD","sym01")
alphaSphorCart = sp.symbols("alphaSphorCart")
betaSphorCartU = ixp.declarerank1("betaSphorCartU")
BSphorCartU = ixp.declarerank1("BSphorCartU")
# UIUC Black Hole initial data are given in Spherical coordinates.
CoordType_in = "Spherical"
# Make sure that rfm.reference_metric() has been called.
# We'll need the variables it defines throughout this module.
if rfm.have_already_called_reference_metric_function == False:
print("Error. Called Convert_Spherical_ADM_to_BSSN_curvilinear() without")
print(" first setting up reference metric, by calling rfm.reference_metric().")
sys.exit(1)
r_th_ph_or_Cart_xyz_oID_xx = []
if CoordType_in == "Spherical":
r_th_ph_or_Cart_xyz_oID_xx = rfm.xxSph
elif CoordType_in == "Cartesian":
r_th_ph_or_Cart_xyz_oID_xx = rfm.xx_to_Cart
else:
print("Error: Can only convert ADM Cartesian or Spherical initial data to BSSN Curvilinear coords.")
sys.exit(1)
xx0,xx1,xx2 basis [Back to top]¶The following discussion holds for either Spherical or Cartesian input data, so for simplicity let's just assume the data are given in Spherical coordinates.
All ADM tensors and vectors are in the Spherical coordinate basis $x^i_{\rm Sph} = (r,\theta,\phi)$, but we need them in the curvilinear coordinate basis $x^i_{\rm rfm}=$(xx0,xx1,xx2) set by the "reference_metric::CoordSystem" variable. Empirically speaking, it is far easier to write (x(xx0,xx1,xx2),y(xx0,xx1,xx2),z(xx0,xx1,xx2)) than the inverse, so we will compute the Jacobian matrix
via exact differentiation (courtesy SymPy), and the inverse Jacobian $$ {\rm Jac\_dUrfm\_dDSphUD[i][j]} = \frac{\partial x^i_{\rm rfm}}{\partial x^j_{\rm Sph}}, $$
using NRPy+'s generic_matrix_inverter3x3() function. In terms of these, the transformation of BSSN tensors from Spherical to "reference_metric::CoordSystem" coordinates may be written:
# Step 4: All ADM initial data quantities are now functions of xx0,xx1,xx2, but
# they are still in the Spherical or Cartesian basis. We can now directly apply
# Jacobian transformations to get them in the correct xx0,xx1,xx2 basis:
# Next apply Jacobian transformations to convert into the (xx0,xx1,xx2) basis
# alpha is a scalar, so no Jacobian transformation is necessary.
alpha = alphaSphorCart
Jac_dUSphorCart_dDrfmUD = ixp.zerorank2()
for i in range(DIM):
for j in range(DIM):
Jac_dUSphorCart_dDrfmUD[i][j] = sp.diff(r_th_ph_or_Cart_xyz_oID_xx[i],rfm.xx[j])
Jac_dUrfm_dDSphorCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSphorCart_dDrfmUD)
betaU = ixp.zerorank1()
BU = ixp.zerorank1()
gammaDD = ixp.zerorank2()
KDD = ixp.zerorank2()
for i in range(DIM):
for j in range(DIM):
betaU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * betaSphorCartU[j]
BU[i] += Jac_dUrfm_dDSphorCartUD[i][j] * BSphorCartU[j]
for k in range(DIM):
for l in range(DIM):
gammaDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] * gammaSphorCartDD[k][l]
KDD[i][j] += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] * KSphorCartDD[k][l]
BSSN.BSSN_in_terms_of_ADM (tutorial) to perform the ADM-to-BSSN conversion for all BSSN quantities except $\lambda^i$ [Back to top]¶All ADM quantities were input into this function in the Spherical or Cartesian basis, as functions of $r,\theta,\phi$ or $x,y,z$, respectively. In Step 3 and Step 4 above, we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2. Here we convert ADM quantities in the xx0,xx1,xx2 (a.k.a. "rfm") basis to their BSSN Curvilinear counterparts, in the same basis. Note that we withold computation of the BSSN $\lambda^i$ quantities until a later section of this notebook, as they must be evaluated using finite differences.
# Step 5: Now that we have all ADM quantities in the desired
# basis, we next perform ADM-to-BSSN conversion:
import BSSN.BSSN_in_terms_of_ADM as BitoA
BitoA.gammabarDD_hDD( gammaDD)
BitoA.trK_AbarDD_aDD( gammaDD,KDD)
BitoA.cf_from_gammaDD(gammaDD)
BitoA.betU_vetU( betaU,BU)
hDD = BitoA.hDD
trK = BitoA.trK
aDD = BitoA.aDD
cf = BitoA.cf
vetU = BitoA.vetU
betU = BitoA.betU
This function must first call the ID_ADM_SphorCart() defined above. Using these Spherical or Cartesian data, it sets up all quantities needed for BSSNCurvilinear initial data, except $\lambda^i$, which must be computed from numerical data using finite-difference derivatives.
# Step 6: Output all ADM-to-BSSN expressions to a C function
ID_inputs_param = "ID_inputs other_inputs,"
if pointer_to_ID_inputs == True:
ID_inputs_param = "ID_inputs *other_inputs,"
desc="Write BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"
name="ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs"
params = "const REAL xx0xx1xx2[3],"+ID_inputs_param+"""
REAL *hDD00,REAL *hDD01,REAL *hDD02,REAL *hDD11,REAL *hDD12,REAL *hDD22,
REAL *aDD00,REAL *aDD01,REAL *aDD02,REAL *aDD11,REAL *aDD12,REAL *aDD22,
REAL *trK,
REAL *vetU0,REAL *vetU1,REAL *vetU2,
REAL *betU0,REAL *betU1,REAL *betU2,
REAL *alpha, REAL *cf"""
outCparams = "preindent=1,outCverbose=False,includebraces=False"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name, params=params,
body ="""
REAL gammaSphorCartDD00,gammaSphorCartDD01,gammaSphorCartDD02,
gammaSphorCartDD11,gammaSphorCartDD12,gammaSphorCartDD22;
REAL KSphorCartDD00,KSphorCartDD01,KSphorCartDD02,
KSphorCartDD11,KSphorCartDD12,KSphorCartDD22;
REAL alphaSphorCart,betaSphorCartU0,betaSphorCartU1,betaSphorCartU2;
REAL BSphorCartU0,BSphorCartU1,BSphorCartU2;
const REAL xx0 = xx0xx1xx2[0];
const REAL xx1 = xx0xx1xx2[1];
const REAL xx2 = xx0xx1xx2[2];
REAL xyz_or_rthph[3];\n"""+
outputC(r_th_ph_or_Cart_xyz_oID_xx[0:3],["xyz_or_rthph[0]","xyz_or_rthph[1]","xyz_or_rthph[2]"],"returnstring",
outCparams+",CSE_enable=False")+" "+ADM_input_function_name+"""(xyz_or_rthph, other_inputs,
&gammaSphorCartDD00,&gammaSphorCartDD01,&gammaSphorCartDD02,
&gammaSphorCartDD11,&gammaSphorCartDD12,&gammaSphorCartDD22,
&KSphorCartDD00,&KSphorCartDD01,&KSphorCartDD02,
&KSphorCartDD11,&KSphorCartDD12,&KSphorCartDD22,
&alphaSphorCart,&betaSphorCartU0,&betaSphorCartU1,&betaSphorCartU2,
&BSphorCartU0,&BSphorCartU1,&BSphorCartU2);
// Next compute all rescaled BSSN curvilinear quantities:\n"""+
outputC([hDD[0][0],hDD[0][1],hDD[0][2],hDD[1][1],hDD[1][2],hDD[2][2],
aDD[0][0],aDD[0][1],aDD[0][2],aDD[1][1],aDD[1][2],aDD[2][2],
trK, vetU[0],vetU[1],vetU[2], betU[0],betU[1],betU[2],
alpha, cf],
["*hDD00","*hDD01","*hDD02","*hDD11","*hDD12","*hDD22",
"*aDD00","*aDD01","*aDD02","*aDD11","*aDD12","*aDD22",
"*trK", "*vetU0","*vetU1","*vetU2", "*betU0","*betU1","*betU2",
"*alpha","*cf"],"returnstring",params=outCparams),
enableCparameters=False)
# Step 6.a: Output the driver function for the above C function
# Next write ID_BSSN__ALL_BUT_LAMBDAs(), the driver
# function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs():
desc="""Driver function for ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(),
which writes BSSN variables in terms of ADM variables at a given point xx0,xx1,xx2"""
name="ID_BSSN__ALL_BUT_LAMBDAs"
params = "const int Nxx_plus_2NGHOSTS[3],REAL *xx[3],"+ID_inputs_param+"REAL *in_gfs"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name, params=params,
body ="""
const int idx = IDX3(i0,i1,i2);
const REAL xx0xx1xx2[3] = {xx0,xx1,xx2};
ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs(xx0xx1xx2,other_inputs,
&in_gfs[IDX4pt(HDD00GF,idx)],&in_gfs[IDX4pt(HDD01GF,idx)],&in_gfs[IDX4pt(HDD02GF,idx)],
&in_gfs[IDX4pt(HDD11GF,idx)],&in_gfs[IDX4pt(HDD12GF,idx)],&in_gfs[IDX4pt(HDD22GF,idx)],
&in_gfs[IDX4pt(ADD00GF,idx)],&in_gfs[IDX4pt(ADD01GF,idx)],&in_gfs[IDX4pt(ADD02GF,idx)],
&in_gfs[IDX4pt(ADD11GF,idx)],&in_gfs[IDX4pt(ADD12GF,idx)],&in_gfs[IDX4pt(ADD22GF,idx)],
&in_gfs[IDX4pt(TRKGF,idx)],
&in_gfs[IDX4pt(VETU0GF,idx)],&in_gfs[IDX4pt(VETU1GF,idx)],&in_gfs[IDX4pt(VETU2GF,idx)],
&in_gfs[IDX4pt(BETU0GF,idx)],&in_gfs[IDX4pt(BETU1GF,idx)],&in_gfs[IDX4pt(BETU2GF,idx)],
&in_gfs[IDX4pt(ALPHAGF,idx)],&in_gfs[IDX4pt(CFGF,idx)]);
""",
loopopts = "AllPoints,Read_xxs,oldloops",enableCparameters=False)
We compute $\bar{\Lambda}^i$ (Eqs. 4 and 5 of Baumgarte et al.), from finite-difference derivatives of rescaled metric quantities $h_{ij}$:
$$ \bar{\Lambda}^i = \bar{\gamma}^{jk}\left(\bar{\Gamma}^i_{jk} - \hat{\Gamma}^i_{jk}\right). $$The reference_metric.py module provides us with analytic expressions for $\hat{\Gamma}^i_{jk}$, so here we need only compute finite-difference expressions for $\bar{\Gamma}^i_{jk}$, based on the values for $h_{ij}$ provided in the initial data. Once $\bar{\Lambda}^i$ has been computed, we apply the usual rescaling procedure: $$ \lambda^i = \bar{\Lambda}^i/\text{ReU[i]}, $$ and then output the result to a C file using the NRPy+ finite-difference C output routine.
# Step 7: Compute $\bar{\Lambda}^i$ from finite-difference derivatives of rescaled metric quantities
# We will need all BSSN gridfunctions to be defined, as well as
# expressions for gammabarDD_dD in terms of exact derivatives of
# the rescaling matrix and finite-difference derivatives of
# hDD's. This functionality is provided by BSSN.BSSN_unrescaled_and_barred_vars,
# which we call here to overwrite above definitions of gammabarDD,gammabarUU, etc.
Bq.gammabar__inverse_and_derivs() # Provides gammabarUU and GammabarUDD
gammabarUU = Bq.gammabarUU
GammabarUDD = Bq.GammabarUDD
# Next evaluate \bar{\Lambda}^i, based on GammabarUDD above and GammahatUDD
# (from the reference metric):
LambdabarU = ixp.zerorank1()
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])
# Finally apply rescaling:
# lambda^i = Lambdabar^i/\text{ReU[i]}
lambdaU = ixp.zerorank1()
for i in range(DIM):
lambdaU[i] = LambdabarU[i] / rfm.ReU[i]
outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
lambdaU_expressions = [lhrh(lhs=gri.gfaccess("in_gfs","lambdaU0"),rhs=lambdaU[0]),
lhrh(lhs=gri.gfaccess("in_gfs","lambdaU1"),rhs=lambdaU[1]),
lhrh(lhs=gri.gfaccess("in_gfs","lambdaU2"),rhs=lambdaU[2])]
desc="Output lambdaU[i] for BSSN, built using finite-difference derivatives."
name="ID_BSSN_lambdas"
params = "const int Nxx[3],const int Nxx_plus_2NGHOSTS[3],REAL *xx[3],const REAL dxx[3],REAL *in_gfs"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name, params=params,
preloop = """
const REAL invdx0 = 1.0/dxx[0];
const REAL invdx1 = 1.0/dxx[1];
const REAL invdx2 = 1.0/dxx[2];
""",
body = fin.FD_outputC("returnstring",lambdaU_expressions, outCparams),
loopopts = "InteriorPoints,Read_xxs,oldloops",enableCparameters=False)
BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear NRPy+ module [Back to top]¶Here, as a code validation check, we verify agreement in the C codes for converting "numerical" UIUCBlackHole initial data (in Spherical coordinates/basis) to BSSN Curvilinear data in Cylindrical coordinates/basis between
By default, we analyze these expressions in Cylindrical coordinates, though other coordinate systems may be chosen.
# Reset the gridfunctions list;
# in Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear()
# below, BSSN_RHSs is called
# tutorial. This line of code enables us to run
# Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear()
# without resetting the running Python kernel.
gri.glb_gridfcs_list = []
import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as AtoBnum
AtoBnum.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear("Spherical",ADM_input_function_name,
Ccodesdir=PyModCcodesdir)
print("\n\n ### BEGIN VALIDATION TESTS")
import filecmp
for file in ["ID_BSSN_lambdas.h","ID_BSSN__ALL_BUT_LAMBDAs.h",
"ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h"]:
if filecmp.cmp(os.path.join(Ccodesdir,file),
os.path.join(PyModCcodesdir,file)) == False:
print("VALIDATION TEST FAILED ON file: "+file+".")
sys.exit(1)
else:
print("Validation test PASSED on file: "+file)
Once the following code finishes running, the generated PDF may be found at the following location within the directory you have the NRPy+ tutorial saved: Tutorial-ADM_Initial_Data-Converting_Numerical_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.pdf
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-ADM_Initial_Data-Converting_Numerical_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear")