Notebook

Weyl Scalars and Invariants: An Introduction to Einstein Toolkit Diagnostic Thorns¶

Author: Patrick Nelson & Zach Etienne¶

Formatting improvements courtesy Brandon Clark¶

This notebook outlines the construction of an Einstein Toolkit diagnostic module to compute the Weyl scalars and invariants using ADMBase gridfunctions as input, with expressions generated through SymPy and NRPy+. The process includes converting SymPy expressions into C-code kernels and integrating these within the larger Einstein Toolkit infrastructure, providing a methodological approach to develop and integrate custom diagnostic tools.¶

Notebook Status: Validated

Validation Notes: Numerical results from this module have been confirmed to agree with the trusted WeylScal4 Einstein Toolkit thorn to roundoff error.

NRPy+ Source Code for this module:¶

which are fully documented in the NRPy+ Tutorial-WeylScalars-Cartesian module on using NRPy+ to construct the Weyl scalars and invariants as SymPy expressions.

Introduction:¶

In the previous tutorial notebook, we constructed within SymPy full expressions for the real and imaginary components of all five Weyl scalars $\psi_0$, $\psi_1$, $\psi_2$, $\psi_3$, and $\psi_4$ as well as the Weyl invariants. So that we can easily access these expressions, we have ported the Python code needed to generate the Weyl scalar SymPy expressions to WeylScal4NRPD/WeylScalars_Cartesian.py, and the Weyl invariant SymPy expressions to WeylScal4NRPD/WeylScalarInvariants_Cartesian.py.

Here we will work through the steps necessary to construct an Einstein Toolkit diagnostic thorn (module), starting from these SymPy expressions, which computes these expressions using ADMBase gridfunctions as input. This tutorial is in two steps:

Call on NRPy+ to convert the SymPy expressions for the Weyl Scalars and associated Invariants into one C-code kernel for each.
Write the C code and build up the needed Einstein Toolkit infrastructure (i.e., the .ccl files).

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows

Step 1: Call on NRPy+ to convert the SymPy expressions for the Weyl scalars and associated invariants into one C-code kernel for each
Step 2: Interfacing with the Einstein Toolkit
1. Step 2.a: Constructing the Einstein Toolkit C-code calling functions that include the C code kernels
2. Step 2.b: CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure
3. Step 2.c: Add the C file to Einstein Toolkit compilation list
Step 3: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Call on NRPy+ to convert the SymPy expressions for the Weyl scalars and associated invariants into one C-code kernel for each [Back to top]¶

$$\label{nrpy}$$

WARNING: It takes some time to generate the CSE-optimized C code kernels for these quantities, especially the Weyl scalars... expect 5 minutes on a modern computer.

In [1]:

from outputC import *            # NRPy+: Core C code output module
import finite_difference as fin  # NRPy+: Finite difference C code generation module
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 cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface
import loop as lp                # NRPy+: loop infrasructure
import shutil, os, sys, time     # Standard Python modules for multiplatform OS-level functions, benchmarking

# Step 1: Set the coordinate system for the numerical grid to Cartesian.
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.

# Step 2: Set the finite differencing order to FD_order to 4
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", 4)

# Step 3: Create output directories
!mkdir WeylScal4NRPD     2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.
!mkdir WeylScal4NRPD/src 2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.

# Step 4: Generate symbolic expressions
# Since we are writing an Einstein Toolkit thorn, we must set our memory access style to "ETK".
par.set_parval_from_str("grid::GridFuncMemAccess","ETK")
import BSSN.Psi4_tetrads as BP4t
par.set_parval_from_str("BSSN.Psi4_tetrads::TetradChoice","QuasiKinnersley")
#par.set_parval_from_str("BSSN.Psi4_tetrads::UseCorrectUnitNormal","True")
import BSSN.Psi4 as BP4
print("Generating symbolic expressions for psi4...")
start = time.time()
BP4.Psi4()
end = time.time()
print("(BENCH) Finished psi4 symbolic expressions in "+str(end-start)+" seconds.")

psi4r = gri.register_gridfunctions("AUX","psi4r")
psi4r0pt = gri.register_gridfunctions("AUX","psi4r0pt")
psi4r1pt = gri.register_gridfunctions("AUX","psi4r1pt")
psi4r2pt = gri.register_gridfunctions("AUX","psi4r2pt")

# Construct RHSs:
psi4r_lhrh = [lhrh(lhs=gri.gfaccess("out_gfs","psi4r"),rhs=BP4.psi4_re_pt[0]+BP4.psi4_re_pt[1]+BP4.psi4_re_pt[2]),
              lhrh(lhs=gri.gfaccess("out_gfs","psi4r0pt"),rhs=BP4.psi4_re_pt[0]),
              lhrh(lhs=gri.gfaccess("out_gfs","psi4r1pt"),rhs=BP4.psi4_re_pt[1]),
              lhrh(lhs=gri.gfaccess("out_gfs","psi4r2pt"),rhs=BP4.psi4_re_pt[2])]

# Generating the CSE is the slowest
# operation in this notebook, and much of the CSE
# time is spent sorting CSE expressions. Disabling
# this sorting makes the C codegen 3-4x faster,
# but the tradeoff is that every time this is
# run, the CSE patterns will be different 
# (though they should result in mathematically 
# *identical* expressions). You can expect 
# roundoff-level differences as a result.
start = time.time()
print("Generating C code kernel for psi4r...")
psi4r_CcodeKernel = fin.FD_outputC("returnstring",psi4r_lhrh,params="outCverbose=False,CSE_sorting=none")
end = time.time()
print("(BENCH) Finished psi4r C code kernel generation in "+str(end-start)+" seconds.")
psi4r_looped = lp.loop(["i2","i1","i0"],["2","2","2"],["cctk_lsh[2]-2","cctk_lsh[1]-2","cctk_lsh[0]-2"],\
                       ["1","1","1"],["#pragma omp parallel for","",""],"","""
               const CCTK_REAL xx0 = xGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];
               const CCTK_REAL xx1 = yGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];
               const CCTK_REAL xx2 = zGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];
"""+psi4r_CcodeKernel)
with open("WeylScal4NRPD/src/WeylScal4NRPD_psi4r.h", "w") as file:
    file.write(str(psi4r_looped))

Generating symbolic expressions for psi4...
(BENCH) Finished psi4 symbolic expressions in 1.7516753673553467 seconds.
Generating C code kernel for psi4r...
(BENCH) Finished psi4r C code kernel generation in 40.00490069389343 seconds.

Step 2: Interfacing with the Einstein Toolkit [Back to top]¶

$$\label{etk}$$

Step 2.a: Constructing the Einstein Toolkit calling functions that include the C code kernels [Back to top]¶

$$\label{etkc}$$

Now that we have generated the C code kernels (WeylScal4NRPD_psis.h and WeylScal4NRPD_invars.h) express the Weyl scalars and invariants as CSE-optimized finite-difference expressions, we next need to write the C code functions that incorporate these kernels and are called by the Einstein Toolkit scheduler.

In [2]:

%%writefile WeylScal4NRPD/src/WeylScal4NRPD.c

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

void WeylScal4NRPD_calc_psi4r(const cGH* restrict const cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,
const CCTK_REAL invdx0,const CCTK_REAL invdx1,const CCTK_REAL invdx2,
const CCTK_REAL *xGF,const CCTK_REAL *yGF,const CCTK_REAL *zGF,
const CCTK_REAL *hDD00GF,const CCTK_REAL *hDD01GF,const CCTK_REAL *hDD02GF,const CCTK_REAL *hDD11GF,const CCTK_REAL *hDD12GF,const CCTK_REAL *hDD22GF,
const CCTK_REAL *aDD00GF,const CCTK_REAL *aDD01GF,const CCTK_REAL *aDD02GF,const CCTK_REAL *aDD11GF,const CCTK_REAL *aDD12GF,const CCTK_REAL *aDD22GF,
const CCTK_REAL *trKGF,const CCTK_REAL *cfGF,
CCTK_REAL *psi4rGF,
CCTK_REAL *psi4r0ptGF,
CCTK_REAL *psi4r1ptGF,
CCTK_REAL *psi4r2ptGF) {

    DECLARE_CCTK_PARAMETERS;

#include "WeylScal4NRPD_psi4r.h"

}

extern void WeylScal4NRPD_mainfunction(CCTK_ARGUMENTS) {

  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;

  if(cctk_iteration % WeylScal4NRPD_calc_every != 0) { return; }
  
  const CCTK_REAL invdx0 = 1.0 / (CCTK_DELTA_SPACE(0));
  const CCTK_REAL invdx1 = 1.0 / (CCTK_DELTA_SPACE(1));
  const CCTK_REAL invdx2 = 1.0 / (CCTK_DELTA_SPACE(2));

  /* Now, to calculate psi4: */
  WeylScal4NRPD_calc_psi4r(cctkGH,cctk_lsh,cctk_nghostzones,
                             invdx0,invdx1,invdx2,
                             x,y,z,
                             hDD00GF,hDD01GF,hDD02GF,hDD11GF,hDD12GF,hDD22GF,
                             aDD00GF,aDD01GF,aDD02GF,aDD11GF,aDD12GF,aDD22GF,
                             trKGF,cfGF,
                             psi4rGF,
                             psi4r0ptGF,psi4r1ptGF,psi4r2ptGF);
}

Overwriting WeylScal4NRPD/src/WeylScal4NRPD.c

In [3]:

# First we convert from ADM to BSSN, as is required to convert initial data 
#    (given using) ADM quantities, to the BSSN evolved variables
import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as atob
IDhDD,IDaDD,IDtrK,IDvetU,IDbetU,IDalpha,IDcf,IDlambdaU = \
    atob.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear("Cartesian","DoNotOutputADMInputFunction",os.path.join("WeylScal4NRPD","src"))

# Store the original list of registered gridfunctions; we'll want to unregister
#   all the *SphorCart* gridfunctions after we're finished with them below.
orig_glb_gridfcs_list = []
for gf in gri.glb_gridfcs_list:
    orig_glb_gridfcs_list.append(gf)

alphaSphorCart   = gri.register_gridfunctions(                 "AUXEVOL", "alphaSphorCart")
betaSphorCartU   = ixp.register_gridfunctions_for_single_rank1("AUXEVOL", "betaSphorCartU")
BSphorCartU      = ixp.register_gridfunctions_for_single_rank1("AUXEVOL", "BSphorCartU")
gammaSphorCartDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL", "gammaSphorCartDD", "sym01")
KSphorCartDD     = ixp.register_gridfunctions_for_single_rank2("AUXEVOL", "KSphorCartDD", "sym01")

# ADM to BSSN conversion, used for converting ADM initial data into a form readable by this thorn.
# ADM to BSSN, Part 1: Set up function call and pointers to ADM gridfunctions
outstr = """
#include <math.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

void WeylScal4NRPD_ADM_to_BSSN(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    CCTK_REAL *alphaSphorCartGF = alp;
"""
# It's ugly if we output code in the following ordering, so we'll first
#   output to a string and then sort the string to beautify the code a bit.
outstrtmp = []
for i in range(3):
    outstrtmp.append("    CCTK_REAL *betaSphorCartU"+str(i)+"GF = beta"+chr(ord('x')+i)+";\n")
#     outstrtmp.append("    CCTK_REAL *BSphorCartU"+str(i)+"GF = dtbeta"+chr(ord('x')+i)+";\n")
    for j in range(i,3):
        outstrtmp.append("    CCTK_REAL *gammaSphorCartDD"+str(i)+str(j)+"GF = g"+chr(ord('x')+i)+chr(ord('x')+j)+";\n")
        outstrtmp.append("    CCTK_REAL *KSphorCartDD"+str(i)+str(j)+"GF = k"+chr(ord('x')+i)+chr(ord('x')+j)+";\n")
outstrtmp.sort()
for line in outstrtmp:
    outstr += line

# ADM to BSSN, Part 2: Set up ADM to BSSN conversions for BSSN gridfunctions that do not require
#                      finite-difference derivatives (i.e., all gridfunctions except lambda^i (=Gamma^i 
#                      in non-covariant BSSN)):
#                      h_{ij}, a_{ij}, trK, vet^i=beta^i,bet^i=B^i, cf (conformal factor), and alpha
all_but_lambdaU_expressions = [
    lhrh(lhs=gri.gfaccess("in_gfs","hDD00"),rhs=IDhDD[0][0]),
    lhrh(lhs=gri.gfaccess("in_gfs","hDD01"),rhs=IDhDD[0][1]),
    lhrh(lhs=gri.gfaccess("in_gfs","hDD02"),rhs=IDhDD[0][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","hDD11"),rhs=IDhDD[1][1]),
    lhrh(lhs=gri.gfaccess("in_gfs","hDD12"),rhs=IDhDD[1][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","hDD22"),rhs=IDhDD[2][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD00"),rhs=IDaDD[0][0]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD01"),rhs=IDaDD[0][1]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD02"),rhs=IDaDD[0][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD11"),rhs=IDaDD[1][1]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD12"),rhs=IDaDD[1][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","aDD22"),rhs=IDaDD[2][2]),
    lhrh(lhs=gri.gfaccess("in_gfs","trK"),rhs=IDtrK),
    lhrh(lhs=gri.gfaccess("in_gfs","vetU0"),rhs=IDvetU[0]),
    lhrh(lhs=gri.gfaccess("in_gfs","vetU1"),rhs=IDvetU[1]),
    lhrh(lhs=gri.gfaccess("in_gfs","vetU2"),rhs=IDvetU[2]),
    lhrh(lhs=gri.gfaccess("in_gfs","alpha"),rhs=IDalpha),
    lhrh(lhs=gri.gfaccess("in_gfs","cf"),rhs=IDcf)]

outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
all_but_lambdaU_outC = fin.FD_outputC("returnstring",all_but_lambdaU_expressions, outCparams)
outstr += lp.loop(["i2","i1","i0"],["0","0","0"],["cctk_lsh[2]","cctk_lsh[1]","cctk_lsh[0]"],
                   ["1","1","1"],["#pragma omp parallel for","",""],"    ",all_but_lambdaU_outC)

outstr += "} // END void WeylScal4NRPD_ADM_to_BSSN(CCTK_ARGUMENTS)\n"

with open("WeylScal4NRPD/src/ADM_to_BSSN.c", "w") as file:
    file.write(str(outstr))

Step 2.b: CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure [Back to top]¶

$$\label{cclfiles}$$

Writing a module ("thorn") within the Einstein Toolkit requires that three "ccl" files be constructed, all in the root directory of the thorn:

1.interface.ccl: defines the gridfunction groups needed, and provides keywords denoting what this thorn provides and what it should inherit from other thorns.

param.ccl: specifies free parameters within the thorn.
schedule.ccl: allocates storage for gridfunctions, defines how the thorn's functions should be scheduled in a broader simulation, and specifies the regions of memory written to or read from gridfunctions.

Let's start with interface.ccl. The official Einstein Toolkit (Cactus) documentation defines what must/should be included in an interface.ccl file here.

In [4]:

%%writefile WeylScal4NRPD/interface.ccl

# With "implements", we give our thorn its unique name.
implements: WeylScal4NRPD

# By "inheriting" other thorns, we tell the Toolkit that we 
#   will rely on variables/function that exist within those
#   functions. 
inherits:   admbase Boundary Grid methodoflines

# Tell the Toolkit that we want the various Weyl scalars 
#    and invariants to be visible to other thorns by using 
#    the keyword "public". Note that declaring these 
#    gridfunctions *does not* allocate memory for them;
#    that is done by the schedule.ccl file.
public:
CCTK_REAL NRPyPsi4_group type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=0 tensorparity=1'
{
  psi4rGF,psi4r0ptGF,psi4r1ptGF,psi4r2ptGF, psi4iGF
} "Psi4_group"

CCTK_REAL evol_variables type = GF Timelevels=3
{
    aDD00GF,aDD01GF,aDD02GF,aDD11GF,aDD12GF,aDD22GF,alphaGF,cfGF,hDD00GF,hDD01GF,hDD02GF,hDD11GF,hDD12GF,hDD22GF,trKGF,vetU0GF,vetU1GF,vetU2GF
} "BSSN evolved gridfunctions, sans lambdaU and partial t beta"

Overwriting WeylScal4NRPD/interface.ccl

We will now write the file param.ccl. This file allows the listed parameters to be set at runtime. We also give allowed ranges and default values for each parameter. More information on this file's syntax can be found in the official Einstein Toolkit documentation.

The first parameter specifies how many time levels need to be stored. Generally when using the ETK's adaptive-mesh refinement (AMR) driver Carpet, three timelevels are needed so that the diagnostic quantities can be properly interpolated and defined across refinement boundaries.

The second parameter determines how often we will calculate $\psi_4$, and the third parameter indicates whether just $\psi_4$, all Weyl scalars, or all Weyl scalars and invariants are going to be output. The third parameter is currently specified entirely within NRPy+, so by this point, it is not a free parameter. Thus it is not quite correct to include it in this list of free parameters (FIXME).

In [5]:

%%writefile WeylScal4NRPD/param.ccl

restricted:
CCTK_INT timelevels "Number of active timelevels" STEERABLE=RECOVER
{
  0:3 :: ""
} 3

restricted:
CCTK_INT WeylScal4NRPD_calc_every "WeylScal4_psi4_calc_Nth_calc_every" STEERABLE=ALWAYS
{
  *:* :: ""
} 1

Overwriting WeylScal4NRPD/param.ccl

Finally, we will write the file schedule.ccl; its official documentation is found here. This file dictates when the various parts of the thorn will be run. We first assign storage for both the real and imaginary components of $\psi_4$, and then specify that we want our code run in the MoL_PseudoEvolution schedule group (consistent with the original WeylScal4 Einstein Toolkit thorn), after the ADM variables are set. At this step, we declare that we will be writing code in C. We also specify the gridfunctions that we wish to read in from memory--in our case, we need all the components of $K_{ij}$ (the spatial extrinsic curvature) and $\gamma_{ij}$ (the physical [as opposed to conformal] 3-metric), in addition to the coordinate values. Note that the ETK adopts the widely-used convention that components of $\gamma_{ij}$ are prefixed in the code with $\text{g}$ and not $\gamma$.

In [6]:

%%writefile WeylScal4NRPD/schedule.ccl

STORAGE: NRPyPsi4_group[3], evol_variables[3]
STORAGE: ADMBase::metric[3], ADMBase::curv[3], ADMBase::lapse[3], ADMBase::shift[3]

schedule group WeylScal4NRPD_group in MoL_PseudoEvolution after ADMBase_SetADMVars
{
} "Schedule WeylScal4NRPD group"

schedule WeylScal4NRPD_ADM_to_BSSN in WeylScal4NRPD_group before weylscal4_mainfunction
{
  LANG: C
} "Convert ADM into BSSN variables"


schedule WeylScal4NRPD_mainfunction in WeylScal4NRPD_group after WeylScal4NRPD_ADM_to_BSSN
{
  LANG: C
} "Call WeylScal4NRPD main function"

Overwriting WeylScal4NRPD/schedule.ccl

Step 2.c: Tell the Einstein Toolkit to compile the C code [Back to top]¶

$$\label{etk_list}$$

The make.code.defn lists the source files that need to be compiled. Naturally, this thorn has only the one C file $-$ written above $-$ to compile:

In [7]:

%%writefile WeylScal4NRPD/src/make.code.defn

SRCS = WeylScal4NRPD.c ADM_to_BSSN.c

Overwriting WeylScal4NRPD/src/make.code.defn

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

$$\label{latex_pdf_output}$$

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-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [8]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-ETK_thorn-WeylScal4NRPD")

Created Tutorial-ETK_thorn-WeylScal4NRPD.tex, and compiled LaTeX file to PDF file Tutorial-ETK_thorn-WeylScal4NRPD.pdf