# Converting Numerical ADM Initial Data in the Spherical or Cartesian Basis to BSSN Initial Data in the Desired Curvilinear Basis¶

## Author: Zach Etienne¶

### This module is meant for use only with initial data that can be represented numerically in ADM form, either in the Spherical or Cartesian basis. I.e., the ADM variables are given $\left\{\gamma_{ij}, K_{ij}, \alpha, \beta^i\right\}$ numerically as functions of $(r,\theta,\phi)$ or $(x,y,z)$; e.g., through an initial data solver. If instead the ADM initial data are provided as exact (algebraic) functions of $(r,\theta,\phi)$ or $(x,y,z)$, then is is better to use the Exact-ADM-Spherical/Cartesian-to-BSSNCurvilinear module instead.¶

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)

## Introduction:¶

$$\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:

$$\left\{h_{i j},a_{i j},\phi, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\}.$$

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.

$$\label{toc}$$

This notebook is organized as follows

1. Step 1: Initialize core Python/NRPy+ modules
2. Step 2: Desired output BSSN Curvilinear coordinate system set to Cylindrical, as a proof-of-principle
3. Step 3: Make ADM variables functions of ${\rm xx0},{\rm xx1},{\rm xx2}$ instead of functions of Cartesian or Spherical coordinates
4. Step 4: Applying Jacobian transformations to get in the correct xx0,xx1,xx2 basis
5. Step 5: Call functions within BSSN.BSSN_in_terms_of_ADM (tutorial) to perform the ADM-to-BSSN conversion for all BSSN quantities except $\lambda^i$
6. Step 6: Output all ADM-to-BSSN expressions to a C function
1. Step 6.a: Output the driver function for the above C function
7. Step 7: Compute $\bar{\Lambda}^i$ from finite-difference derivatives of rescaled metric quantities
8. Step 8: Code Validation against BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear NRPy+ module
9. Step 9: Output this notebook to $\LaTeX$-formatted PDF file

$$\label{initializenrpy}$$
In [1]:
# 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:
# 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,
# 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 [Back to top]¶

$$\label{cylindrical}$$
In [2]:
# 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.
pointer_to_ID_inputs = False


# Step 3: Make ADM variables functions of ${\rm xx0},{\rm xx1},{\rm xx2}$ instead of functions of Cartesian or Spherical coordinates [Back to top]¶

$$\label{admxx0xx1xx2}$$

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.

In [3]:
# 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.
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)


# Step 4: Applying Jacobian transformations to get in the correct xx0,xx1,xx2 basis [Back to top]¶

$$\label{adm_jacobian}$$

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

$${\rm Jac\_dUSph\_dDrfmUD[i][j]} = \frac{\partial x^i_{\rm Sph}}{\partial x^j_{\rm rfm}},$$

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:

\begin{align} \beta^i_{\rm rfm} &= \frac{\partial x^i_{\rm rfm}}{\partial x^\ell_{\rm Sph}} \beta^\ell_{\rm Sph}\\ B^i_{\rm rfm} &= \frac{\partial x^i_{\rm rfm}}{\partial x^\ell_{\rm Sph}} B^\ell_{\rm Sph}\\ \gamma^{\rm rfm}_{ij} &= \frac{\partial x^\ell_{\rm Sph}}{\partial x^i_{\rm rfm}} \frac{\partial x^m_{\rm Sph}}{\partial x^j_{\rm rfm}} \gamma^{\rm Sph}_{\ell m}\\ K^{\rm rfm}_{ij} &= \frac{\partial x^\ell_{\rm Sph}}{\partial x^i_{\rm rfm}} \frac{\partial x^m_{\rm Sph}}{\partial x^j_{\rm rfm}} K^{\rm Sph}_{\ell m} \end{align}
In [4]:
# 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()
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):
KDD[i][j]     += Jac_dUSphorCart_dDrfmUD[k][i]*Jac_dUSphorCart_dDrfmUD[l][j] *     KSphorCartDD[k][l]


# Step 5: Call functions within BSSN.BSSN_in_terms_of_ADM (tutorial) to perform the ADM-to-BSSN conversion for all BSSN quantities except $\lambda^i$ [Back to top]¶

$$\label{adm2bssn}$$

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.

In [5]:
# Step 5: Now that we have all ADM quantities in the desired
#         basis, we next perform ADM-to-BSSN conversion:

BitoA.betU_vetU(      betaU,BU)
hDD     = BitoA.hDD
trK     = BitoA.trK
cf      = BitoA.cf
vetU    = BitoA.vetU
betU    = BitoA.betU


$$\label{adm2bssn_c}$$

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.

In [6]:
# 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"
params = "const REAL xx0xx1xx2[3],"+ID_inputs_param+"""
REAL *hDD00,REAL *hDD01,REAL *hDD02,REAL *hDD11,REAL *hDD12,REAL *hDD22,
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",
&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],
trK,  vetU[0],vetU[1],vetU[2],  betU[0],betU[1],betU[2],
alpha, cf],
["*hDD00","*hDD01","*hDD02","*hDD11","*hDD12","*hDD22",
"*trK",  "*vetU0","*vetU1","*vetU2",  "*betU0","*betU1","*betU2",
"*alpha","*cf"],"returnstring",params=outCparams),
enableCparameters=False)

Output C function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs() to file numerical_ADM_to_BSSN_Ccodes/ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h


$$\label{driver}$$

We output the driver function for the above C function: ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs()

In [7]:
# Step 6.a: Output the driver function for the above C function

# Next write ID_BSSN__ALL_BUT_LAMBDAs(), the driver
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};
&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(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)]);
""",

Output C function ID_BSSN__ALL_BUT_LAMBDAs() to file numerical_ADM_to_BSSN_Ccodes/ID_BSSN__ALL_BUT_LAMBDAs.h


# Step 7: Compute $\bar{\Lambda}^i$ from finite-difference derivatives of rescaled metric quantities [Back to top]¶

$$\label{lambda}$$

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.

In [8]:
# 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),

Output C function ID_BSSN_lambdas() to file numerical_ADM_to_BSSN_Ccodes/ID_BSSN_lambdas.h


# Step 8: Code Validation against BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear NRPy+ module [Back to top]¶

$$\label{code_validation}$$

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

1. this tutorial and

By default, we analyze these expressions in Cylindrical coordinates, though other coordinate systems may be chosen.

In [9]:
# Reset the gridfunctions list;
#  below, BSSN_RHSs is called
#          tutorial. This line of code enables us to run
#          without resetting the running Python kernel.
gri.glb_gridfcs_list = []

Ccodesdir=PyModCcodesdir)

print("\n\n ### BEGIN VALIDATION TESTS")
import filecmp
for file in ["ID_BSSN_lambdas.h","ID_BSSN__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)

Output C function ID_BSSN_lambdas() to file numerical_ADM_to_BSSN_Ccodes/PyMod/ID_BSSN_lambdas.h
Output C function ID_BSSN__ALL_BUT_LAMBDAs() to file numerical_ADM_to_BSSN_Ccodes/PyMod/ID_BSSN__ALL_BUT_LAMBDAs.h

### BEGIN VALIDATION TESTS
Validation test PASSED on file: ID_BSSN_lambdas.h
Validation test PASSED on file: ID_BSSN__ALL_BUT_LAMBDAs.h
Validation test PASSED on file: ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h


# Step 9: Output this notebook to $\LaTeX$-formatted PDF file [Back to top]¶

$$\label{latex_pdf_output}$$

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

In [10]:
import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface

Created Tutorial-ADM_Initial_Data-