Notebook

`BaikalETK`: NRPy+-Based BSSN Solvers for the Einstein Toolkit¶

Author: Zach Etienne¶

Special thanks to Roland Haas for help in answering many implementation questions¶

This notebook generates `Baikal` and `BaikalVacuum`, Einstein Toolkit thorns for solving Einstein's equations in the BSSN formalism, in Cartesian coordinates. These thorns are highly optimized for modern CPU architectures, featuring SIMD intrinsics and OpenMP support.¶

Notebook Status: Validated against the Einstein Toolkit McLachlan BSSN thorn, both in the context of black hole binary simulations (excellent gravitational waveform agreement) as well as binary neutron star simulations (when parameter enable_stress_energy_source_terms below is set to True). Once plots demonstrating this agreement are added to this tutorial notebook, the validation status will be set to Validated.

Validation Notes: This tutorial notebook has been validated against a trusted Einstein Toolkit thorn, but plots demonstrating its validation have yet to be included in this notebook.

Introduction¶

How often did my soul cry out:
Come back to Baikal once again?
I still do not know this lake:
To see does not mean to know.

Igor Severyanin, [1].

Lake Baikal is home to the nerpa seal, NRPy+'s mascot.

This notebook generates two thorns: Baikal and BaikalVacuum. Baikal contains stress-energy source terms (i.e., the $8\pi T^{\mu\nu}$ part of Einstein's equations) for general relativistic hydrodynamics (GRHD) and magnetohydrodynamics (GRMHD), and the BaikalVacuum contains no such source terms, so can be used for e.g., black hole and binary black hole calculations in which no self-gravitating matter is considered.

Baikal and BaikalVacuum are meant to reproduce the functionality of the McLachlan thorn, generated by the Mathematica-based Kranc code, but using the fully open-source NRPy+/SymPy infrastructure.

Unlike McLachlan, Baikal and BaikalVacuum enable the user at runtime to select only the most popular options, like finite-difference order and Kreiss-Oliger dissipation strength. Further, neither thorn yet supports symmetry options or Jacobians needed for Llama grids. Support for these options might be provided in a future update. As for BSSN gauge choice, Baikal and BaikalVacuum by default support only the BSSN moving-puncture gauge conditions, though this can be changed by selecting one of the other gauge choices implemented within the NRPy+ infrastructure.

Regarding the structure of this notebook:

As generating the optimized C code kernels needed for these thorns can individually take roughly 10 minutes per finite-difference order, there is a strong motivation to parallelize the code generation process.

To this end, this Jupyter notebook does not itself run Python code in most code blocks directly. Instead

code blocks needed to generate the C code kernels (i.e., the C-code representations of the basic equations) are simply output to the $BaikalETKdir/BaikalETK_C_kernels_codegen.py file using the %%writefile IPython magic command;
code blocks needed to generate the Einstein Toolkit (ETK) CCL files are output to $BaikalETKdir/BaikalETK_ETK_ccl_files_codegen.py file using the same IPython magic; and
code blocks needed to generate the C driver routines that call the C code kernels and interface with the rest of the ETK are output to $BaikalETKdir/BaikalETK_C_drivers_codegen.py.

The only non-%%writefile IPython magic codeblocks used in this notebook act as glue so that the above three Python modules can be called in sequence, with particular emphasis on generating the C code kernels in parallel (by calling multiple instances of the Python module $BaikalETKdir/BaikalETK_C_kernels_codegen.py). On a modern, multi-core CPU, this greatly cuts down the time needed to generate the thorns (sometimes by 10x or more), all while ensuring a convenient user interface for adjusting the thorns to suit one's needs.

Associated NRPy+ Source Code & Tutorial Modules for this module:¶

BSSN/ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py; [tutorial]: Spherical/Cartesian ADM$\to$Curvilinear BSSN converter function, for which ADM quantities are assumed given at each gridpoint (i.e., exact, closed-form expressions are not given). This is used to generate BaikalETK's ADM$\to$BSSN function, as in the ETK spacetime evolution modules are to assume that initial data are given as ADM quantities in the Cartesian basis at each gridpoint.
BSSN/ADM_in_terms_of_BSSN.py; [tutorial]: Constructs ADM quantities in terms of BSSN quantities (in arbitrary curvilinear coordinates, though we use Cartesian here). This is used to generate BaikalETK's BSSN$\to$ADM function, which make ADM variables available to diagnostic thorns within the ETK.
BSSN/BSSN_constraints.py; [tutorial]: Hamiltonian constraint in BSSN curvilinear basis/coordinates
BSSN/BSSN_RHSs.py; [tutorial]: Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates
BSSN/BSSN_gauge_RHSs.py; [tutorial]: Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear coordinates

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows

Step 1: Initialize needed Python/NRPy+ modules
Step 2: NRPy+-generated C code kernels for BSSN spacetime solve
1. Step 2.a: BSSN RHS expressions
2. Step 2.b: Hamiltonian & momentum constraints
3. Step 2.c: Enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\det{\hat{\gamma}_{ij}}=1$)
4. Step 2.d: Given WhichPart parameter choice, generate algorithm for calling corresponding function within BaikalETK_C_kernels_codegen_onepart() to generate C code kernel
5. Step 2.e: Generate C code kernels for Baikal and BaikalETK
  1. Step 2.e.i: Set compile-time and runtime parameters for Baikal and BaikalVacuum
  2. Step 2.e.ii: Generate all C-code kernels for Baikal and BaikalVacuum, in parallel if possible
Step 3: CCL files - Define how this module interacts and interfaces with the wider Einstein Toolkit infrastructure
1. Step 3.a: param.ccl: specify free parameters within BaikalETK
2. Step 3.b: interface.ccl: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns
3. Step 3.c: schedule.ccl:schedule all functions used within BaikalETK, specify data dependencies within said functions, and allocate memory for gridfunctions
Step 4: C driver functions for ETK registration & NRPy+-generated kernels
1. Step 4.a: Needed ETK functions: Banner, Symmetry registration, Parameter sanity check, Method of Lines (MoL) registration, Boundary condition
2. Step 4.b: BSSN $\leftrightarrow$ ADM conversions
  1. Step 4.b.i: ADM $\to$ BSSN
  2. Step 4.b.ii: BSSN $\to$ ADM
3. Step 4.c Evaluate BSSN right-hand-sides (RHSs)
  1. Step 4.c.i: Evaluate Ricci tensor
  2. Step 4.c.ii: Evaluate BSSN RHSs, using Ricci tensor as input
4. Step 4.d: Enforcing conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\det{\hat{\gamma}_{ij}}=1$)
5. Step 4.e: Diagnostics: Computing Hamiltonian & momentum constraints
6. Step 4.f: driver_BSSN_T4UU(): Compute $T^{\mu\nu}$ from TmunuBase's $T_{\mu\nu}$
7. Step 4.g: Output all C driver functions needed for Baikal/BaikalVacuum
8. Step 4.h: make.code.defn: List of all C driver functions needed to compile BaikalETK
Step 5: Code Validation
1. Step 5.a: Validation against BaikalETK/ Python modules
Step 6: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Initialize needed Python/NRPy+ modules [Back to top]¶

$$\label{initializenrpy}$$

In [1]:

# Compile-time (i.e., NRPy+-time) parameters for both Baikal & BaikalVacuum:
LapseCondition = "OnePlusLog"
ShiftCondition = "GammaDriving2ndOrder_NoCovariant"

In [2]:

# We output to an BaikalETK_validate so as not to overwrite the "trusted"
#   BaikalETK Python modules directory (generated by an earlier/equivalent
#   version of this notebook, and to validate against these trusted Python
#   modules. To this end, at the bottom of this notebook we perform a diff
#   between this notebook's output against the trusted BaikalETK Python
#   modules and error out if any difference exists. The proper approach
#   for updating BaikalETK Python modules is to first update this notebook
#   then once the updated Baikal/BaikalVacuum ETK thorns have been
#   revalidated, copy the Python modules from BaikalETK_validate to
#   BaikalETK.
BaikalETKdir = "BaikalETK_validate"
import cmdline_helper as cmd  # NRPy+: Multi-platform Python command-line interface
import os, sys, shutil        # Standard Python modules for multiplatform OS-level functions

# Output finite difference stencils as functions instead of inlined expressions.
#   Dramatically speeds up compile times (esp with higher-order finite differences
#   and GCC 9.3+)
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
par.set_parval_from_str("finite_difference::enable_FD_functions", True)

cmd.mkdir(os.path.join(BaikalETKdir))

In [3]:

%%writefile $BaikalETKdir/BaikalETK_C_kernels_codegen.py

# Step 1: Import needed core NRPy+ modules
from outputC import lhrh         # 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 loop as lp                # NRPy+: Generate C code loops
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm   # NRPy+: Reference metric support
import os, sys                   # Standard Python modules for multiplatform OS-level functions
import time                      # Standard Python module; useful for benchmarking
import BSSN.BSSN_RHSs as rhs
import BSSN.BSSN_gauge_RHSs as gaugerhs

def BaikalETK_C_kernels_codegen_onepart(params=
        "WhichPart=BSSN_RHSs,ThornName=Baikal,FD_order=4,enable_stress_energy_source_terms=True"):
    # Set default parameters
    WhichPart = "BSSN_RHSs"
    ThornName = "Baikal"
    FD_order = 4
    enable_stress_energy_source_terms = True
    LapseCondition  = "OnePlusLog" # Set the standard 1+log lapse condition
    # Set the standard, second-order advecting-shift, Gamma-driving shift condition:
    ShiftCondition  = "GammaDriving2ndOrder_NoCovariant"
    # Default Kreiss-Oliger dissipation strength
    default_KO_strength = 0.1

    import re
    if params != "":
        params2 = re.sub("^,","",params)
        params = params2.strip()
        splitstring = re.split("=|,", params)

        if len(splitstring) % 2 != 0:
            print("outputC: Invalid params string: "+params)
            sys.exit(1)

        parnm = []
        value = []
        for i in range(int(len(splitstring)/2)):
            parnm.append(splitstring[2*i])
            value.append(splitstring[2*i+1])

        for i in range(len(parnm)):
            parnm.append(splitstring[2*i])
            value.append(splitstring[2*i+1])

        for i in range(len(parnm)):
            if parnm[i] == "WhichPart":
                WhichPart = value[i]
            elif parnm[i] == "ThornName":
                ThornName = value[i]
            elif parnm[i] == "FD_order":
                FD_order = int(value[i])
            elif parnm[i] == "enable_stress_energy_source_terms":
                enable_stress_energy_source_terms = False
                if value[i] == "True":
                    enable_stress_energy_source_terms = True
            elif parnm[i] == "LapseCondition":
                LapseCondition = value[i]
            elif parnm[i] == "ShiftCondition":
                ShiftCondition = value[i]
            elif parnm[i] == "default_KO_strength":
                default_KO_strength = float(value[i])
            else:
                print("BaikalETK Error: Could not parse input param: "+parnm[i])
                sys.exit(1)

    # Set output directory for C kernels
    outdir = os.path.join(ThornName, "src")  # Main C code output directory

    # Set spatial dimension (must be 3 for BSSN)
    par.set_parval_from_str("grid::DIM", 3)

    # Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,
    #                                 FD order, floating point precision, and CFL factor:
    # Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,
    #              SymTP, SinhSymTP
    # NOTE: Only CoordSystem == Cartesian makes sense here; new
    #       boundary conditions are needed within the ETK for
    #       Spherical, etc. coordinates.
    CoordSystem     = "Cartesian"

    par.set_parval_from_str("reference_metric::CoordSystem", CoordSystem)
    rfm.reference_metric()  # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.

    # Set the gridfunction memory access type to ETK-like, so that finite_difference
    #    knows how to read and write gridfunctions from/to memory.
    par.set_parval_from_str("grid::GridFuncMemAccess", "ETK")

    par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption", ShiftCondition)
    par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::LapseEvolutionOption", LapseCondition)

    T4UU = None
    if enable_stress_energy_source_terms == True:
        registered_already = False
        for i in range(len(gri.glb_gridfcs_list)):
            if gri.glb_gridfcs_list[i].name == "T4UU00":
                registered_already = True
        if not registered_already:
            T4UU = ixp.register_gridfunctions_for_single_rank2("AUXEVOL", "T4UU", "sym01", DIM=4)
        else:
            T4UU = ixp.declarerank2("T4UU", "sym01", DIM=4)

    # Register the BSSN constraints (Hamiltonian & momentum constraints) as gridfunctions.
    registered_already = False
    for i in range(len(gri.glb_gridfcs_list)):
        if gri.glb_gridfcs_list[i].name == "H":
            registered_already = True
    if not registered_already:
        # We ignore return values for register_gridfunctions...() calls below
        #    as they are unused.
        gri.register_gridfunctions("AUX", "H")
        ixp.register_gridfunctions_for_single_rank1("AUX", "MU")

Overwriting BaikalETK_validate/BaikalETK_C_kernels_codegen.py

Step 2: Output C code for BSSN spacetime solve [Back to top]¶

$$\label{bssn}$$

Step 2.a: BSSN RHS expressions [Back to top]¶

$$\label{bssnrhs}$$

BaikalETK implements a fully covariant version of the BSSN 3+1 decomposition of Einstein's equations of general relativity, which is fully documented within NRPy+ (start here). However, especially if you are a newcomer to the field of numerical relativity, you may also find the following lectures and papers useful for understanding the adopted formalism:

Mathematical foundations of BSSN and 3+1 initial value problem decompositions of Einstein's equations:
- Thomas Baumgarte's lectures on mathematical formulation of numerical relativity
- Yuichiro Sekiguchi's introduction to BSSN
Extensions to the standard BSSN approach used in NRPy+

Here, we simply call the BSSN.BSSN_RHSs; [tutorial] and BSSN.BSSN_gauge_RHSs; [tutorial] NRPy+ Python modules to generate the symbolic expressions, add Kreiss-Oliger dissipation, and then output the finite-difference C code form of the equations using NRPy+'s finite_difference (tutorial) C code generation module.

In [4]:

%%writefile -a $BaikalETKdir/BaikalETK_C_kernels_codegen.py

    def BSSN_RHSs__generate_symbolic_expressions():
        ######################################
        # START: GENERATE SYMBOLIC EXPRESSIONS
        print("Generating symbolic expressions for BSSN RHSs...")
        start = time.time()
        # Enable rfm_precompute infrastructure, which results in
        #   BSSN RHSs that are free of transcendental functions,
        #   even in curvilinear coordinates, so long as
        #   ConformalFactor is set to "W" (default).
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")
        par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(outdir,"rfm_files/"))

        # Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:
        import BSSN.BSSN_quantities as Bq
        par.set_parval_from_str("BSSN.BSSN_quantities::LeaveRicciSymbolic","True")

        rhs.BSSN_RHSs()

        if T4UU != None:
            import BSSN.BSSN_stress_energy_source_terms as Bsest
            Bsest.BSSN_source_terms_for_BSSN_RHSs(T4UU)
            rhs.trK_rhs += Bsest.sourceterm_trK_rhs
            for i in range(3):
                # Needed for Gamma-driving shift RHSs:
                rhs.Lambdabar_rhsU[i] += Bsest.sourceterm_Lambdabar_rhsU[i]
                # Needed for BSSN RHSs:
                rhs.lambda_rhsU[i]    += Bsest.sourceterm_lambda_rhsU[i]
                for j in range(3):
                    rhs.a_rhsDD[i][j] += Bsest.sourceterm_a_rhsDD[i][j]

        gaugerhs.BSSN_gauge_RHSs()

        # Add Kreiss-Oliger dissipation to the BSSN RHSs:
        thismodule = "KO_Dissipation"
        diss_strength = par.Cparameters("REAL", thismodule, "diss_strength", default_KO_strength)

        alpha_dKOD = ixp.declarerank1("alpha_dKOD")
        cf_dKOD    = ixp.declarerank1("cf_dKOD")
        trK_dKOD   = ixp.declarerank1("trK_dKOD")
        betU_dKOD    = ixp.declarerank2("betU_dKOD","nosym")
        vetU_dKOD    = ixp.declarerank2("vetU_dKOD","nosym")
        lambdaU_dKOD = ixp.declarerank2("lambdaU_dKOD","nosym")
        aDD_dKOD = ixp.declarerank3("aDD_dKOD","sym01")
        hDD_dKOD = ixp.declarerank3("hDD_dKOD","sym01")
        for k in range(3):
            gaugerhs.alpha_rhs += diss_strength*alpha_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
            rhs.cf_rhs         += diss_strength*   cf_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
            rhs.trK_rhs        += diss_strength*  trK_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
            for i in range(3):
                if "2ndOrder" in ShiftCondition:
                    gaugerhs.bet_rhsU[i] += diss_strength*   betU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
                gaugerhs.vet_rhsU[i]     += diss_strength*   vetU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
                rhs.lambda_rhsU[i]       += diss_strength*lambdaU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
                for j in range(3):
                    rhs.a_rhsDD[i][j] += diss_strength*aDD_dKOD[i][j][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]
                    rhs.h_rhsDD[i][j] += diss_strength*hDD_dKOD[i][j][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]

        # We use betaU as our upwinding control vector:
        Bq.BSSN_basic_tensors()
        betaU = Bq.betaU

        # Now that we are finished with all the rfm hatted
        #           quantities in generic precomputed functional
        #           form, let's restore them to their closed-
        #           form expressions.
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","False") # Reset to False to disable rfm_precompute.
        rfm.ref_metric__hatted_quantities()
        par.set_parval_from_str("BSSN.BSSN_quantities::LeaveRicciSymbolic","False")
        end = time.time()
        print("(BENCH) Finished BSSN RHS symbolic expressions in "+str(end-start)+" seconds.")
        # END: GENERATE SYMBOLIC EXPRESSIONS
        ######################################

        BSSN_RHSs_SymbExpressions = [lhrh(lhs=gri.gfaccess("rhs_gfs","aDD00"),   rhs=rhs.a_rhsDD[0][0]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","aDD01"),   rhs=rhs.a_rhsDD[0][1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","aDD02"),   rhs=rhs.a_rhsDD[0][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","aDD11"),   rhs=rhs.a_rhsDD[1][1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","aDD12"),   rhs=rhs.a_rhsDD[1][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","aDD22"),   rhs=rhs.a_rhsDD[2][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","alpha"),   rhs=gaugerhs.alpha_rhs),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","betU0"),   rhs=gaugerhs.bet_rhsU[0]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","betU1"),   rhs=gaugerhs.bet_rhsU[1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","betU2"),   rhs=gaugerhs.bet_rhsU[2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","cf"),      rhs=rhs.cf_rhs),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD00"),   rhs=rhs.h_rhsDD[0][0]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD01")   ,rhs=rhs.h_rhsDD[0][1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD02"),   rhs=rhs.h_rhsDD[0][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD11"),   rhs=rhs.h_rhsDD[1][1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD12"),   rhs=rhs.h_rhsDD[1][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","hDD22"),   rhs=rhs.h_rhsDD[2][2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","lambdaU0"),rhs=rhs.lambda_rhsU[0]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","lambdaU1"),rhs=rhs.lambda_rhsU[1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","lambdaU2"),rhs=rhs.lambda_rhsU[2]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","trK"),     rhs=rhs.trK_rhs),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","vetU0"),   rhs=gaugerhs.vet_rhsU[0]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","vetU1"),   rhs=gaugerhs.vet_rhsU[1]),
                                     lhrh(lhs=gri.gfaccess("rhs_gfs","vetU2"),   rhs=gaugerhs.vet_rhsU[2]) ]

        return [betaU,BSSN_RHSs_SymbExpressions]

    def Ricci__generate_symbolic_expressions():
        ######################################
        # START: GENERATE SYMBOLIC EXPRESSIONS
        print("Generating symbolic expressions for Ricci tensor...")
        start = time.time()
        # Enable rfm_precompute infrastructure, which results in
        #   BSSN RHSs that are free of transcendental functions,
        #   even in curvilinear coordinates, so long as
        #   ConformalFactor is set to "W" (default).
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")
        par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(outdir,"rfm_files/"))

        # Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:
        import BSSN.BSSN_quantities as Bq

        # Next compute Ricci tensor
        # par.set_parval_from_str("BSSN.BSSN_quantities::LeaveRicciSymbolic","False")
        RbarDD_already_registered = False
        for i in range(len(gri.glb_gridfcs_list)):
            if "RbarDD00" in gri.glb_gridfcs_list[i].name:
                RbarDD_already_registered = True
        if not RbarDD_already_registered:
            # We ignore the return value of ixp.register_gridfunctions_for_single_rank2() below
            #    as it is unused.
            ixp.register_gridfunctions_for_single_rank2("AUXEVOL","RbarDD","sym01")
        rhs.BSSN_RHSs()
        Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()

        # Now that we are finished with all the rfm hatted
        #           quantities in generic precomputed functional
        #           form, let's restore them to their closed-
        #           form expressions.
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","False") # Reset to False to disable rfm_precompute.
        rfm.ref_metric__hatted_quantities()
        end = time.time()
        print("(BENCH) Finished Ricci symbolic expressions in "+str(end-start)+" seconds.")
        # END: GENERATE SYMBOLIC EXPRESSIONS
        ######################################

        Ricci_SymbExpressions = [lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD00"),rhs=Bq.RbarDD[0][0]),
                                 lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD01"),rhs=Bq.RbarDD[0][1]),
                                 lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD02"),rhs=Bq.RbarDD[0][2]),
                                 lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD11"),rhs=Bq.RbarDD[1][1]),
                                 lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD12"),rhs=Bq.RbarDD[1][2]),
                                 lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD22"),rhs=Bq.RbarDD[2][2])]

        return [Ricci_SymbExpressions]

    def BSSN_RHSs__generate_Ccode(all_RHSs_Ricci_exprs_list):
        betaU                     = all_RHSs_Ricci_exprs_list[0]
        BSSN_RHSs_SymbExpressions = all_RHSs_Ricci_exprs_list[1]

        print("Generating C code for BSSN RHSs (FD_order="+str(FD_order)+",Tmunu="+str(enable_stress_energy_source_terms)+") in "+par.parval_from_str("reference_metric::CoordSystem")+" coordinates.")
        start = time.time()

        # Store original finite-differencing order:
        FD_order_orig = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")
        # Set new finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)

        BSSN_RHSs_string = fin.FD_outputC("returnstring",BSSN_RHSs_SymbExpressions,
                                          params="outCverbose=False,enable_SIMD=True,GoldenKernelsEnable=True",
                                          upwindcontrolvec=betaU)

        filename = "BSSN_RHSs_enable_Tmunu_"+str(enable_stress_energy_source_terms)+"_FD_order_"+str(FD_order)+".h"
        with open(os.path.join(outdir,filename), "w") as file:
            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],
           ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],
                               ["1","1","SIMD_width"],
                                ["#pragma omp parallel for",
                                 "#include \"rfm_files/rfm_struct__SIMD_outer_read2.h\"",
                                 r"""    #include "rfm_files/rfm_struct__SIMD_outer_read1.h"
    #if (defined __INTEL_COMPILER && __INTEL_COMPILER_BUILD_DATE >= 20180804)
        #pragma ivdep         // Forces Intel compiler (if Intel compiler used) to ignore certain SIMD vector dependencies
        #pragma vector always // Forces Intel compiler (if Intel compiler used) to vectorize
    #endif"""],"",
                    "#include \"rfm_files/rfm_struct__SIMD_inner_read0.h\"\n"+BSSN_RHSs_string))

        # Restore original finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order_orig)

        end = time.time()
        print("(BENCH) Finished BSSN_RHS C codegen (FD_order="+str(FD_order)+",Tmunu="+str(enable_stress_energy_source_terms)+") in " + str(end - start) + " seconds.")


    def Ricci__generate_Ccode(Ricci_exprs_list):
        Ricci_SymbExpressions     = Ricci_exprs_list[0]

        print("Generating C code for Ricci tensor (FD_order="+str(FD_order)+") in "+par.parval_from_str("reference_metric::CoordSystem")+" coordinates.")
        start = time.time()

        # Store original finite-differencing order:
        FD_order_orig = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")
        # Set new finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)

        Ricci_string = fin.FD_outputC("returnstring", Ricci_SymbExpressions,
                                       params="outCverbose=False,enable_SIMD=True,GoldenKernelsEnable=True")

        filename = "BSSN_Ricci_FD_order_"+str(FD_order)+".h"
        with open(os.path.join(outdir, filename), "w") as file:
            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],
           ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],
                               ["1","1","SIMD_width"],
                                ["#pragma omp parallel for",
                                 "#include \"rfm_files/rfm_struct__SIMD_outer_read2.h\"",
                                 r"""    #include "rfm_files/rfm_struct__SIMD_outer_read1.h"
    #if (defined __INTEL_COMPILER && __INTEL_COMPILER_BUILD_DATE >= 20180804)
        #pragma ivdep         // Forces Intel compiler (if Intel compiler used) to ignore certain SIMD vector dependencies
        #pragma vector always // Forces Intel compiler (if Intel compiler used) to vectorize
    #endif"""],"",
                    "#include \"rfm_files/rfm_struct__SIMD_inner_read0.h\"\n"+Ricci_string))

        # Restore original finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order_orig)

        end = time.time()
        print("(BENCH) Finished Ricci C codegen (FD_order="+str(FD_order)+") in " + str(end - start) + " seconds.")

Appending to BaikalETK_validate/BaikalETK_C_kernels_codegen.py

Step 2.b: Hamiltonian & momentum constraints [Back to top]¶

$$\label{hammomconstraints}$$

Next output the C code for evaluating the Hamiltonian & momentum constraints (Tutorial). In the absence of numerical error, this constraint should evaluate to zero. However it does not due to numerical (typically truncation and roundoff) error. Therefore it is useful to measure the Hamiltonian & momentum constraint violation to gauge the accuracy of our simulation, and, ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected.

In [5]:

%%writefile -a $BaikalETKdir/BaikalETK_C_kernels_codegen.py

    def BSSN_constraints__generate_symbolic_expressions_and_C_code():
        ######################################
        # START: GENERATE SYMBOLIC EXPRESSIONS
        # Define the Hamiltonian constraint and output the optimized C code.
        import BSSN.BSSN_constraints as bssncon

        bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False)
        if T4UU != None:
            import BSSN.BSSN_stress_energy_source_terms as Bsest
            Bsest.BSSN_source_terms_for_BSSN_RHSs(T4UU)
            Bsest.BSSN_source_terms_for_BSSN_constraints(T4UU)
            bssncon.H += Bsest.sourceterm_H
            for i in range(3):
                bssncon.MU[i] += Bsest.sourceterm_MU[i]
        # END: GENERATE SYMBOLIC EXPRESSIONS
        ######################################

        # Store original finite-differencing order:
        FD_order_orig = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")
        # Set new finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)

        start = time.time()
        print("Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.")
        Ham_mom_string = fin.FD_outputC("returnstring",
                                        [lhrh(lhs=gri.gfaccess("aux_gfs", "H"),   rhs=bssncon.H),
                                         lhrh(lhs=gri.gfaccess("aux_gfs", "MU0"), rhs=bssncon.MU[0]),
                                         lhrh(lhs=gri.gfaccess("aux_gfs", "MU1"), rhs=bssncon.MU[1]),
                                         lhrh(lhs=gri.gfaccess("aux_gfs", "MU2"), rhs=bssncon.MU[2])],
                                        params="outCverbose=False")

        with open(os.path.join(outdir,"BSSN_constraints_enable_Tmunu_"+str(enable_stress_energy_source_terms)+"_FD_order_"+str(FD_order)+".h"), "w") as file:
            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],
           ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],
                               ["1","1","1"],["#pragma omp parallel for","",""], "", Ham_mom_string))

        # Restore original finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order_orig)

        end = time.time()
        print("(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order="+str(FD_order)+",Tmunu="+str(enable_stress_energy_source_terms)+") in " + str(end - start) + " seconds.")

Appending to BaikalETK_validate/BaikalETK_C_kernels_codegen.py

Step 2.c: Enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\det{\hat{\gamma}_{ij}}=1$) [Back to top]¶

$$\label{gamconstraint}$$

Then enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (Eq. 53 of Ruchlin, Etienne, and Baumgarte (2018)), as documented in the corresponding NRPy+ tutorial notebook

Applying curvilinear boundary conditions should affect the initial data at the outer boundary, and will in general cause the $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint to be violated there. Thus after we apply these boundary conditions, we must always call the routine for enforcing the $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint:

In [6]:

%%writefile -a $BaikalETKdir/BaikalETK_C_kernels_codegen.py

    def enforce_detgammabar_eq_detgammahat__generate_symbolic_expressions_and_C_code():
        ######################################
        # START: GENERATE SYMBOLIC EXPRESSIONS

        # Enable rfm_precompute infrastructure, which results in
        #   BSSN RHSs that are free of transcendental functions,
        #   even in curvilinear coordinates, so long as
        #   ConformalFactor is set to "W" (default).
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")
        par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(outdir,"rfm_files/"))

        import BSSN.Enforce_Detgammahat_Constraint as EGC
        enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammahat_Constraint_symb_expressions()

        # Now that we are finished with all the rfm hatted
        #           quantities in generic precomputed functional
        #           form, let's restore them to their closed-
        #           form expressions.
        par.set_parval_from_str("reference_metric::enable_rfm_precompute","False") # Reset to False to disable rfm_precompute.
        rfm.ref_metric__hatted_quantities()
        # END: GENERATE SYMBOLIC EXPRESSIONS
        ######################################

        start = time.time()
        print("Generating optimized C code (FD_order="+str(FD_order)+") for gamma constraint. May take a while, depending on CoordSystem.")
        enforce_gammadet_string = fin.FD_outputC("returnstring", enforce_detg_constraint_symb_expressions,
                                                 params="outCverbose=False,preindent=0,includebraces=False")

        with open(os.path.join(outdir,"enforcedetgammabar_constraint.h"), "w") as file:
            file.write(lp.loop(["i2","i1","i0"],["0", "0", "0"],
                               ["cctk_lsh[2]","cctk_lsh[1]","cctk_lsh[0]"],
                               ["1","1","1"],
                                ["#pragma omp parallel for",
                                     "#include \"rfm_files/rfm_struct__read2.h\"",
                                     "#include \"rfm_files/rfm_struct__read1.h\""],"",
                                     "#include \"rfm_files/rfm_struct__read0.h\"\n"+enforce_gammadet_string))
        end = time.time()
        print("(BENCH) Finished gamma constraint C codegen (FD_order="+str(FD_order)+") in " + str(end - start) + " seconds.")

Appending to BaikalETK_validate/BaikalETK_C_kernels_codegen.py

Step 2.d: Given `WhichPart` parameter choice, generate algorithm for calling corresponding function within `BaikalETK_C_kernels_codegen_onepart()` to generate C code kernel [Back to top]¶

$$\label{driver_call_codegen_funcs}$$

In [7]:

%%writefile -a $BaikalETKdir/BaikalETK_C_kernels_codegen.py

    if WhichPart=="BSSN_RHSs":
        BSSN_RHSs__generate_Ccode(BSSN_RHSs__generate_symbolic_expressions())
    elif WhichPart=="Ricci":
        Ricci__generate_Ccode(Ricci__generate_symbolic_expressions())
    elif WhichPart=="BSSN_constraints":
        BSSN_constraints__generate_symbolic_expressions_and_C_code()
    elif WhichPart=="detgammabar_constraint":
        enforce_detgammabar_eq_detgammahat__generate_symbolic_expressions_and_C_code()
    else:
        print("Error: WhichPart = "+WhichPart+" is not recognized.")
        sys.exit(1)

    # Store all NRPy+ environment variables to an output string so NRPy+ environment from within this subprocess can be easily restored
    import pickle  # Standard Python module for converting arbitrary data structures to a uniform format.
    # https://www.pythonforthelab.com/blog/storing-binary-data-and-serializing/
    outstr = []
    outstr.append(pickle.dumps(len(gri.glb_gridfcs_list)))
    for lst in gri.glb_gridfcs_list:
        outstr.append(pickle.dumps(lst.gftype))
        outstr.append(pickle.dumps(lst.name))
        outstr.append(pickle.dumps(lst.rank))
        outstr.append(pickle.dumps(lst.DIM))

    outstr.append(pickle.dumps(len(par.glb_params_list)))
    for lst in par.glb_params_list:
        outstr.append(pickle.dumps(lst.type))
        outstr.append(pickle.dumps(lst.module))
        outstr.append(pickle.dumps(lst.parname))
        outstr.append(pickle.dumps(lst.defaultval))

    outstr.append(pickle.dumps(len(par.glb_Cparams_list)))
    for lst in par.glb_Cparams_list:
        outstr.append(pickle.dumps(lst.type))
        outstr.append(pickle.dumps(lst.module))
        outstr.append(pickle.dumps(lst.parname))
        outstr.append(pickle.dumps(lst.defaultval))

    outstr.append(pickle.dumps(len(outC_function_dict)))
    for Cfuncname, Cfunc in outC_function_dict.items():
        outstr.append(pickle.dumps(Cfuncname))
        outstr.append(pickle.dumps(Cfunc))

    return outstr

Appending to BaikalETK_validate/BaikalETK_C_kernels_codegen.py

Step 2.e: Generate C code kernels for `Baikal` and `BaikalETK` [Back to top]¶

$$\label{kernel_codegen}$$

Here we generate the C code kernels (i.e., the C-code representation of the equations needed) for Baikal and BaikalETK.

Step 2.e.i: Set compile-time and runtime parameters for `Baikal` and `BaikalVacuum` [Back to top]¶

$$\label{feature_choice}$$

NRPy+ is a code generation package that is designed to offer maximum flexibility at the time of C code generation. As a result, although NRPy+ can in principle output an infinite variety of C code kernels for solving Einstein's equations, generally free parameters in each kernel steerable at runtime are restricted to simple scalars. This leads to more optimized kernels (i.e., significantly improved performance as compared to McLachlan), but at the expense of flexibility in choosing e.g., different gauges at runtime (currently only the most commonly used gauge is supported), as well as the need to generate multiple kernels (one per finite-differencing order). Reducing the number of kernels and adding more flexibility at runtime will be a focus of future work.

For now, Baikal and BaikalVacuum support the following runtime options:

Baikal: Enables $T^{\mu\nu}$ source terms; useful for general relativistic hydrodynamics (GRHD) and general relativistic magnetohydrodynamics (GRMHD) simulations.
- Finite differencing of order 2 and 4, via runtime parameter FD_order
  - BSSN RHSs, Ricci tensor, Hamiltonian constraint, and $\hat{\gamma}=\bar{\gamma}$ constraint
- Kreiss-Oliger dissipation strength via runtime parameter diss_strength, which is the exact analogue of the eps_diss parameter of the Dissipation thorn
BaikalVacuum: Disables $T^{\mu\nu}$ source terms; useful for generating dynamical black hole and binary black hole spacetimes, in which matter is not present.
- Finite differencing of orders 4, 6, and 8 via runtime parameter FD_order
  - BSSN RHSs, Ricci tensor, Hamiltonian constraint, and $\hat{\gamma}=\bar{\gamma}$ constraint
- Kreiss-Oliger dissipation strength via runtime parameter diss_strength, which is the exact analogue of the eps_diss parameter of the Dissipation thorn

Both thorns adopt the standard moving-puncture gauge conditions, which include the 1+log slicing condition for the lapse and the non-covariant $\Gamma$-driving shift condition, as documented in this NRPy+ Tutorial notebook, and implemented in the corresponding NRPy+ Python module. In case you'd like to make another gauge choice, you need only choose the desired gauge from the NRPy+ Python module and add it as parameters ShiftCondition and LapseCondition to the main BaikalETK code generation function BaikalETK_C_kernels_codegen_onepart(). You will find that adding user-defined gauge choices is a straightforward process, as the module is easily extended.

Next we set up the default parameter lists for BaikalETK_C_kernels_codegen_onepart() for Baikal and BaikalVacuum thorns. We set these parameter lists as strings to make parallelizing the code generation far easier (easier to pass a list of strings than a list of function arguments to Python's multiprocessing.Pool()).

In [8]:

# Step 2.e.i: Set compile-time and runtime parameters for Baikal and BaikalVacuum
#         Runtime parameters for
#         Baikal:       FD_orders = [2,4]  ; Gamma-driving eta parameter; Kreiss-Oliger dissipation strength
#         BaikalVacuum: FD_orders = [4,6,8]; Gamma-driving eta parameter; Kreiss-Oliger dissipation strength

paramslist = []
FD_orders = [2,4,6,8]
WhichParamSet = 0
for WhichPart in ["BSSN_RHSs","Ricci","BSSN_constraints","detgammabar_constraint"]:
    for FD_order in FD_orders:
        enable_stress_energy_list = [True,False]
        if FD_order == 2:
            enable_stress_energy_list = [True]
        elif FD_order >= 6:
            enable_stress_energy_list = [False]
        for enable_stress_energy in enable_stress_energy_list:
            ThornName = "Baikal"
            if enable_stress_energy == False:
                ThornName = "BaikalVacuum"
            paramstr = "WhichPart="+WhichPart+","
            paramstr+= "ThornName="+ThornName+","
            paramstr+= "FD_order="+str(FD_order)+","
            paramstr+= "LapseCondition="+LapseCondition+","
            paramstr+= "ShiftCondition="+ShiftCondition+","
            paramstr+= "enable_stress_energy_source_terms="+str(enable_stress_energy)
            if (WhichPart != "detgammabar_constraint") \
               or (WhichPart == "detgammabar_constraint" and FD_order==4):
                # Do not output detgammabar_constraint code more than once for each thorn, as
                #    it does not depend on FD_order
                paramslist.append(paramstr)
                WhichParamSet = WhichParamSet + 1

paramslist.sort()  # Sort the list alphabetically.

Step 2.e.ii: Generate all C-code kernels for `Baikal` and `BaikalVacuum`, in parallel if possible [Back to top]¶

$$\label{parallel_codegen}$$

In [9]:

nrpy_dir_path = os.path.join(".")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

# Create all output directories if they do not yet exist
import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
import shutil                   # Standard Python module for multiplatform OS-level functions
for ThornName in ["Baikal","BaikalVacuum"]:
    outrootdir = ThornName
    cmd.mkdir(os.path.join(outrootdir))
    outdir = os.path.join(outrootdir,"src") # Main C code output directory

    # Copy SIMD/SIMD_intrinsics.h to $outdir/SIMD/SIMD_intrinsics.h, replacing
    #   the line "#define REAL_SIMD_ARRAY REAL" with "#define REAL_SIMD_ARRAY CCTK_REAL"
    #   (since REAL is undefined in the ETK, but CCTK_REAL takes its place)
    cmd.mkdir(os.path.join(outdir,"SIMD"))
    import fileinput
    f = fileinput.input(os.path.join(nrpy_dir_path,"SIMD","SIMD_intrinsics.h"))
    with open(os.path.join(outdir,"SIMD","SIMD_intrinsics.h"),"w") as outfile:
        for line in f:
            outfile.write(line.replace("#define REAL_SIMD_ARRAY REAL", "#define REAL_SIMD_ARRAY CCTK_REAL"))

    # Create directory for rfm_files output
    cmd.mkdir(os.path.join(outdir,"rfm_files"))

# Start parallel C code generation (codegen)
# NRPyEnvVars stores the NRPy+ environment from all the subprocesses in the following
#     parallel codegen
NRPyEnvVars = []

import time   # Standard Python module for benchmarking
import logging
start = time.time()
if __name__ == "__main__":
    try:
        if os.name == 'nt':
            # Windows & Jupyter multiprocessing do not mix, so we run in serial on Windows.
            #  Here's why: https://stackoverflow.com/questions/45719956/python-multiprocessing-in-jupyter-on-windows-attributeerror-cant-get-attribut
            raise Exception("Parallel codegen currently not available in Windows")
        # Step 3.d.ii: Import the multiprocessing module.
        import multiprocessing
        print("***************************************")
        print("Starting parallel C kernel codegen...")
        print("***************************************")

        # Step 3.d.iii: Define master function for parallelization.
        #           Note that lambdifying this doesn't work in Python 3
        def master_func(i):
            import BaikalETK.BaikalETK_C_kernels_codegen as BCk  # lgtm [py/repeated-import]
            return BCk.BaikalETK_C_kernels_codegen_onepart(params=paramslist[i])

        # Step 3.d.iv: Evaluate list of functions in parallel if possible;
        #           otherwise fallback to serial evaluation:
        pool = multiprocessing.Pool() #processes=len(paramslist))
        NRPyEnvVars.append(pool.map(master_func,range(len(paramslist))))
        pool.terminate()
        pool.join()
    except:
        logging.exception("Ignore this warning/backtrace if on a system in which serial codegen is necessary:")
        print("***************************************")
        print("Starting serial C kernel codegen...")
        print("(If you were running in parallel before,")
        print(" this means parallel codegen failed)")
        print("***************************************")
        import BaikalETK.BaikalETK_C_kernels_codegen as BCk  # lgtm [py/repeated-import]
        # No need to pickle if doing serial codegen.
        for param in paramslist:
            BCk.BaikalETK_C_kernels_codegen_onepart(params=param)
        NRPyEnvVars = []  # Reset NRPyEnvVars in case multiprocessing wrote to it and failed.

print("(BENCH) Finished C kernel codegen for Baikal and BaikalVacuum in "+str(time.time()-start)+" seconds.")

***************************************
Starting parallel C kernel codegen...
***************************************
Generating symbolic expressions for Ricci tensor...Generating symbolic expressions for BSSN RHSs...Generating symbolic expressions for Ricci tensor...

Generating symbolic expressions for Ricci tensor...Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for BSSN RHSs...Generating symbolic expressions for Ricci tensor...



Generating symbolic expressions for BSSN RHSs...Generating symbolic expressions for Ricci tensor...

Generating optimized C code (FD_order=4) for gamma constraint. May take a while, depending on CoordSystem.Generating optimized C code (FD_order=4) for gamma constraint. May take a while, depending on CoordSystem.

(BENCH) Finished gamma constraint C codegen (FD_order=4) in 0.10375690460205078 seconds.
(BENCH) Finished gamma constraint C codegen (FD_order=4) in 0.12864923477172852 seconds.
(BENCH) Finished BSSN RHS symbolic expressions in 0.9982852935791016 seconds.
Generating C code for BSSN RHSs (FD_order=8,Tmunu=False) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 1.0139648914337158 seconds.
Generating C code for BSSN RHSs (FD_order=6,Tmunu=False) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 1.0944099426269531 seconds.
Generating C code for BSSN RHSs (FD_order=4,Tmunu=False) in Cartesian coordinates.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.

(BENCH) Finished Ricci symbolic expressions in 1.9740955829620361 seconds.
Generating C code for Ricci tensor (FD_order=4) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 1.9901673793792725 seconds.
Generating C code for Ricci tensor (FD_order=6) in Cartesian coordinates.(BENCH) Finished Ricci symbolic expressions in 1.951263427734375 seconds.

Generating C code for Ricci tensor (FD_order=4) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 2.0254504680633545 seconds.
Generating C code for Ricci tensor (FD_order=8) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 2.0379810333251953 seconds.
Generating C code for Ricci tensor (FD_order=2) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 2.175140857696533 seconds.
Generating C code for BSSN RHSs (FD_order=4,Tmunu=True) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 2.1774909496307373 seconds.
Generating C code for BSSN RHSs (FD_order=2,Tmunu=True) in Cartesian coordinates.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
(BENCH) Finished BSSN_RHS C codegen (FD_order=4,Tmunu=False) in 14.319687128067017 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=6,Tmunu=False) in 15.422996997833252 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=2,Tmunu=True) in 15.380635499954224 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=4,Tmunu=True) in 16.61129069328308 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=8,Tmunu=False) in 18.059557676315308 seconds.
(BENCH) Finished Ricci C codegen (FD_order=4) in 21.51375389099121 seconds.
(BENCH) Finished Ricci C codegen (FD_order=4) in 22.40598464012146 seconds.
(BENCH) Finished Ricci C codegen (FD_order=6) in 22.7694993019104 seconds.
(BENCH) Finished Ricci C codegen (FD_order=2) in 23.129387378692627 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=4,Tmunu=False) in 24.090815782546997 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=6,Tmunu=False) in 24.394925117492676 seconds.
(BENCH) Finished Ricci C codegen (FD_order=8) in 23.86538815498352 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=4,Tmunu=True) in 22.553589582443237 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=2,Tmunu=True) in 22.88631319999695 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=8,Tmunu=False) in 26.489030838012695 seconds.
(BENCH) Finished C kernel codegen for Baikal and BaikalVacuum in 27.915515661239624 seconds.

Step 3: ETK `ccl` file generation [Back to top]¶

$$\label{cclfiles}$$

The Einstein Toolkit (ETK) ccl files contain runtime parameters (param.ccl), registered gridfunctions (interface.ccl), and function scheduling (schedule.ccl). As parameters and gridfunctions are registered with NRPy+ when the C-code kernels are generated, and this generation occurs on separate processes in parallel, we store the entire NRPy+ environment for each process. This results in a tremendous amount of duplication, which is sorted out next. Once all duplicated environment variables (e.g., registered gridfunctions) are removed, we replace the current NRPy+ environment with the new one, by setting gri.glb_gridfcs_list[],par.glb_params_list[],par.glb_Cparams_list[].

In [10]:

# Store all NRPy+ environment variables to file so NRPy+ environment from within this subprocess can be easily restored
import pickle                   # Standard Python module for converting arbitrary data structures to a uniform format.
from outputC import outC_function_dict # NRPy: core C code output module
import grid as gri              # NRPy+: Functions having to do with numerical grids
import finite_difference as fin # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par  # NRPy+: Parameter interface

# https://www.pythonforthelab.com/blog/storing-binary-data-and-serializing/
if len(NRPyEnvVars) > 0:
    grfcs_list = []
    param_list = []
    Cparm_list = []

    outCfunc_dict = {}

    for WhichParamSet in NRPyEnvVars[0]:
        # gridfunctions
        i=0
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            grfcs_list.append(gri.glb_gridfc(gftype=pickle.loads(WhichParamSet[i+0]),
                                             name  =pickle.loads(WhichParamSet[i+1]),
                                             rank  =pickle.loads(WhichParamSet[i+2]),
                                             DIM   =pickle.loads(WhichParamSet[i+3]))) ; i+=4
        # parameters
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            param_list.append(par.glb_param(type      =pickle.loads(WhichParamSet[i+0]),
                                            module    =pickle.loads(WhichParamSet[i+1]),
                                            parname   =pickle.loads(WhichParamSet[i+2]),
                                            defaultval=pickle.loads(WhichParamSet[i+3]))) ; i+=4
        # Cparameters
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            Cparm_list.append(par.glb_Cparam(type      =pickle.loads(WhichParamSet[i+0]),
                                             module    =pickle.loads(WhichParamSet[i+1]),
                                             parname   =pickle.loads(WhichParamSet[i+2]),
                                             defaultval=pickle.loads(WhichParamSet[i+3]))) ; i+=4
        # outC_func_dict
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            funcname = pickle.loads(WhichParamSet[i+0])
            funcbody = pickle.loads(WhichParamSet[i+1]) ; i+=2
            outCfunc_dict[funcname] = funcbody

    grfcs_list_uniq = []
    for gf_ntuple_stored in grfcs_list:
        found_gf = False
        for gf_ntuple_new in grfcs_list_uniq:
            if gf_ntuple_new == gf_ntuple_stored:
                found_gf = True
        if found_gf == False:
            grfcs_list_uniq.append(gf_ntuple_stored)

    param_list_uniq = []
    for pr_ntuple_stored in param_list:
        found_pr = False
        for pr_ntuple_new in param_list_uniq:
            if pr_ntuple_new == pr_ntuple_stored:
                found_pr = True
        if found_pr == False:
            param_list_uniq.append(pr_ntuple_stored)

    # Set glb_paramsvals_list:
    # Step 1: Reset all paramsvals to their defaults
    par.glb_paramsvals_list = []
    for parm in param_list_uniq:
        par.glb_paramsvals_list.append(parm.defaultval)

    Cparm_list_uniq = []
    for Cp_ntuple_stored in Cparm_list:
        found_Cp = False
        for Cp_ntuple_new in Cparm_list_uniq:
            if Cp_ntuple_new == Cp_ntuple_stored:
                found_Cp = True
        if found_Cp == False:
            Cparm_list_uniq.append(Cp_ntuple_stored)

    # Dictionary outCfunc_dict (by the nature of Python dictionaries) will not have duplicates!

    gri.glb_gridfcs_list = []
    par.glb_params_list  = []
    par.glb_Cparams_list = []

    gri.glb_gridfcs_list = grfcs_list_uniq
    par.glb_params_list  = param_list_uniq
    par.glb_Cparams_list = Cparm_list_uniq
    for key, item in outCfunc_dict.items():
        outC_function_dict[key] = item

# Set lapse_floor to default to 1e-15
lap_floor = par.Cparameters("REAL", "BaikalETK", "lapse_floor", 1e-15)

# Step 2: Override defaults with values used here.
#         Note that almost no NRPy+ parameters
#         are used after this point (DIM is definitely
#         one exception), so most of these lines have no
#         effect.
# import reference_metric as rfm
# import grid as gri
# import BSSN.BSSN_gauge_RHSs
par.set_parval_from_str("grid::DIM", 3)
import reference_metric as rfm
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.
par.set_parval_from_str("grid::GridFuncMemAccess", "ETK")
par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption", ShiftCondition)
par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::LapseEvolutionOption", LapseCondition)

# Finally, output all functions needed for computing finite-difference stencils
#   to thornname/src/finite_difference_functions.h
for thornname in ["Baikal", "BaikalVacuum"]:
    fin.output_finite_difference_functions_h(os.path.join(thornname,"src"))

Step 3.a: `param.ccl`: specify free parameters within `BaikalETK` [Back to top]¶

$$\label{paramccl}$$

All parameters necessary for the computation of the BSSN right-hand side (RHS) expressions are registered within NRPy+; we use this information to automatically generate param.ccl. NRPy+ also specifies default values for each parameter.

More information on param.ccl syntax can be found in the official Einstein Toolkit documentation.

In [11]:

%%writefile $BaikalETKdir/BaikalETK_ETK_ccl_files_codegen.py

# Step 1: Import needed core NRPy+ modules
import NRPy_param_funcs as par # NRPy+: Parameter interface
import grid as gri             # NRPy+: Functions having to do with numerical grids
import os, sys                 # Standard Python modules for multiplatform OS-level functions

def keep_param__return_type(paramtuple):
    keep_param = True # We'll not set some parameters in param.ccl;
                      #   e.g., those that should be #define'd like M_PI.
    typestring = ""
    # Separate thorns within the ETK take care of grid/coordinate parameters;
    #   thus we ignore NRPy+ grid/coordinate parameters:
    if paramtuple.module == "grid" or paramtuple.module == "reference_metric":
        keep_param = False

    partype = paramtuple.type
    if partype == "bool":
        typestring += "BOOLEAN "
    elif partype == "REAL":
        if paramtuple.defaultval != 1e300: # 1e300 is a magic value indicating that the C parameter should be mutable
            typestring += "CCTK_REAL "
        else:
            keep_param = False
    elif partype == "int":
        typestring += "CCTK_INT "
    elif partype == "#define":
        keep_param = False
    elif partype == "char":
        # FIXME: char/string parameter types should in principle be supported
        print("Error: parameter "+paramtuple.module+"::"+paramtuple.parname+
              " has unsupported type: \""+ paramtuple.type + "\"")
        sys.exit(1)
    else:
        print("Error: parameter "+paramtuple.module+"::"+paramtuple.parname+
              " has unsupported type: \""+ paramtuple.type + "\"")
        sys.exit(1)
    return keep_param, typestring

def output_param_ccl(ThornName="BaikalETK"):
    with open(os.path.join(ThornName,"param.ccl"), "w") as file:
        file.write("""
# This param.ccl file was automatically generated by NRPy+.
#   You are advised against modifying it directly; instead
#   modify the Python code that generates it.

shares: ADMBase   # Extends multiple ADMBase variables:

EXTENDS CCTK_KEYWORD evolution_method "evolution_method"
{
  "BaikalETK" :: ""
}

EXTENDS CCTK_KEYWORD lapse_evolution_method "lapse_evolution_method"
{
  "BaikalETK" :: ""
}

EXTENDS CCTK_KEYWORD shift_evolution_method "shift_evolution_method"
{
  "BaikalETK" :: ""
}

EXTENDS CCTK_KEYWORD dtshift_evolution_method "dtshift_evolution_method"
{
  "BaikalETK" :: ""
}

EXTENDS CCTK_KEYWORD dtlapse_evolution_method "dtlapse_evolution_method"
{
  "BaikalETK" :: ""
}

restricted:

CCTK_INT FD_order "Finite-differencing order"
{\n""".replace("BaikalETK",ThornName))
        FDorders = []
        for _root, _dirs, files in os.walk(os.path.join(ThornName,"src")): # _root,_dirs unused.
            for Cfilename in files:
                if ("BSSN_Ricci_FD_order" in Cfilename) and (".h" in Cfilename):
                    array = Cfilename.replace(".","_").split("_")
                    FDorders.append(int(array[-2]))
        FDorders.sort()
        for order in FDorders:
            file.write(" "+str(order)+":"+str(order)+"   :: \"finite-differencing order = "+str(order)+"\"\n")
        FDorder_default = 4
        if FDorder_default not in FDorders:
            print("WARNING: 4th-order FD kernel was not output!?! Changing default FD order to "+str(FDorders[0]))
            FDorder_default = FDorders[0]
        file.write("} "+str(FDorder_default)+ "\n\n") # choose 4th order by default, consistent with ML_BSSN
        paramccl_str = ""
        for i in range(len(par.glb_Cparams_list)):
            # keep_param is a boolean indicating whether we should accept or reject
            #    the parameter. singleparstring will contain the string indicating
            #    the variable type.
            keep_param, singleparstring = keep_param__return_type(par.glb_Cparams_list[i])

            if keep_param:
                parname = par.glb_Cparams_list[i].parname
                partype = par.glb_Cparams_list[i].type
                singleparstring += parname + " \""+ parname +" (see NRPy+ for parameter definition)\"\n"
                singleparstring += "{\n"
                if partype != "bool":
                    singleparstring += " *:* :: \"All values accepted. NRPy+ does not restrict the allowed ranges of parameters yet.\"\n"
                singleparstring += "} "+str(par.glb_Cparams_list[i].defaultval)+"\n\n"

                paramccl_str += singleparstring
        file.write(paramccl_str)

Overwriting BaikalETK_validate/BaikalETK_ETK_ccl_files_codegen.py

Step 3.b: `interface.ccl`: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns [Back to top]¶

$$\label{interfaceccl}$$

interface.ccl declares all gridfunctions and determines how BaikalETK interacts with other Einstein Toolkit thorns.

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

In [12]:

%%writefile -a $BaikalETKdir/BaikalETK_ETK_ccl_files_codegen.py

# First construct lists of the basic gridfunctions used in NRPy+.
#    Each type will be its own group in BaikalETK.
evol_gfs_list    = []
auxevol_gfs_list = []
aux_gfs_list     = []
for i in range(len(gri.glb_gridfcs_list)):
    if gri.glb_gridfcs_list[i].gftype == "EVOL":
        evol_gfs_list.append(   gri.glb_gridfcs_list[i].name+"GF")

    if gri.glb_gridfcs_list[i].gftype == "AUX":
        aux_gfs_list.append(    gri.glb_gridfcs_list[i].name+"GF")

    if gri.glb_gridfcs_list[i].gftype == "AUXEVOL":
        auxevol_gfs_list.append(gri.glb_gridfcs_list[i].name+"GF")

# NRPy+'s finite-difference code generator assumes gridfunctions
#    are alphabetized; not sorting may result in unnecessary
#    cache misses.
evol_gfs_list.sort()
aux_gfs_list.sort()
auxevol_gfs_list.sort()
rhs_list = []
for gf in evol_gfs_list:
    rhs_list.append(gf.replace("GF","")+"_rhsGF")

def output_interface_ccl(ThornName="BaikalETK",enable_stress_energy_source_terms=False):
    outstr = """
# This interface.ccl file was automatically generated by NRPy+.
#   You are advised against modifying it directly; instead
#   modify the Python code that generates it.

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

# 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\n"""
    if enable_stress_energy_source_terms == True:
        outstr += "inherits:   TmunuBase\n"
    # Need a raw string here due to all the backslashes:
    outstr += r"""
# Needed functions and #include's:
USES INCLUDE: Symmetry.h
USES INCLUDE: Boundary.h

# Needed Method of Lines function
CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, \
                                          CCTK_INT IN RHSIndex)
REQUIRES FUNCTION MoLRegisterEvolvedGroup

# Needed Boundary Conditions function
CCTK_INT FUNCTION GetBoundarySpecification(CCTK_INT IN size, CCTK_INT OUT ARRAY nboundaryzones, CCTK_INT OUT ARRAY is_internal, CCTK_INT OUT ARRAY is_staggered, CCTK_INT OUT ARRAY shiftout)
USES FUNCTION GetBoundarySpecification

CCTK_INT FUNCTION SymmetryTableHandleForGrid(CCTK_POINTER_TO_CONST IN cctkGH)
USES FUNCTION SymmetryTableHandleForGrid

CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, CCTK_STRING IN var_name, CCTK_STRING IN bc_name)
USES FUNCTION Boundary_SelectVarForBC

# Needed for EinsteinEvolve/NewRad outer boundary condition driver:
CCTK_INT FUNCTION                         \
    NewRad_Apply                          \
        (CCTK_POINTER_TO_CONST IN cctkGH, \
         CCTK_REAL ARRAY IN var,          \
         CCTK_REAL ARRAY INOUT rhs,       \
         CCTK_REAL IN var0,               \
         CCTK_REAL IN v0,                 \
         CCTK_INT IN radpower)
REQUIRES FUNCTION NewRad_Apply

# Needed to convert ADM initial data into BSSN initial data (gamma extrapolation)
CCTK_INT FUNCTION                         \
    ExtrapolateGammas                     \
        (CCTK_POINTER_TO_CONST IN cctkGH, \
         CCTK_REAL ARRAY INOUT var)
REQUIRES FUNCTION ExtrapolateGammas

# Tell the Toolkit that we want all gridfunctions
#    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.

# FIXME: add info for symmetry conditions:
#    https://einsteintoolkit.org/thornguide/CactusBase/SymBase/documentation.html
public:
"""

    # Next we declare gridfunctions based on their corresponding gridfunction groups as registered within NRPy+
    def output_list_of_gfs(gfs_list, description="User did not provide description"):
        gfs_list_parsed = []
        for j in range(len(gfs_list)):
            # Add all gridfunctions in the list...
            gfs_list_parsed.append(gfs_list[j])
            # Then remove T4UU gridfunction from list if enable_stress_energy_source_terms==False:
            if "T4UU" in gfs_list_parsed[-1] and enable_stress_energy_source_terms==False:
                del gfs_list_parsed[-1]
        gfsstr = "    "
        for j in range(len(gfs_list_parsed)):
            gfsstr += gfs_list_parsed[j]
            if j != len(gfs_list_parsed)-1:
                gfsstr += ","  # This is a comma-separated list of gridfunctions
            else:
                gfsstr += "\n} \""+description+"\"\n\n"
        return gfsstr
    # First EVOL type:
    outstr += "CCTK_REAL evol_variables type = GF Timelevels=3\n{\n"
    outstr += output_list_of_gfs(evol_gfs_list, "BSSN evolved gridfunctions")
    # Second EVOL right-hand-sides
    outstr += "CCTK_REAL evol_variables_rhs type = GF Timelevels=1 TAGS=\'InterpNumTimelevels=1 prolongation=\"none\"\'\n{\n"
    outstr += output_list_of_gfs(rhs_list, "right-hand-side storage for BSSN evolved gridfunctions")
    # Then AUX type:
    outstr += "CCTK_REAL aux_variables type = GF Timelevels=3\n{\n"
    outstr += output_list_of_gfs(aux_gfs_list, "Auxiliary gridfunctions for BSSN diagnostics")
    # Finally, AUXEVOL type:
    outstr += "CCTK_REAL auxevol_variables type = GF Timelevels=1 TAGS=\'InterpNumTimelevels=1 prolongation=\"none\"\'\n{\n"
    outstr += output_list_of_gfs(auxevol_gfs_list, "Auxiliary gridfunctions needed for evaluating the BSSN RHSs")

    with open(os.path.join(ThornName, "interface.ccl"), "w") as file:
        file.write(outstr.replace("BaikalETK", ThornName))

Appending to BaikalETK_validate/BaikalETK_ETK_ccl_files_codegen.py

Step 3.c: `schedule.ccl`: schedule all functions used within `BaikalETK`, specify data dependencies within said functions, and allocate memory for gridfunctions [Back to top]¶

$$\label{scheduleccl}$$

Official documentation on constructing ETK schedule.ccl files is found here.

In [13]:

%%writefile -a $BaikalETKdir/BaikalETK_ETK_ccl_files_codegen.py

def output_schedule_ccl(ThornName="BaikalETK",enable_stress_energy_source_terms=False):
    outstr = """
# This schedule.ccl file was automatically generated by NRPy+.
#   You are advised against modifying it directly; instead
#   modify the Python code that generates it.

# First allocate storage for one timelevel of ADMBase gridfunctions, which is the
#    bare minimum needed by NRPy+. If another thorn (e.g., ADMBase itself) requests
#    more timelevels of storage, Cactus automatically allocates the maximum requested.
STORAGE: ADMBase::metric[1], ADMBase::curv[1], ADMBase::lapse[1], ADMBase::shift[1]

# Next allocate storage for all 3 gridfunction groups used in BaikalETK
STORAGE: evol_variables[3]     # Evolution variables
STORAGE: evol_variables_rhs[1] # Variables storing right-hand-sides
STORAGE: aux_variables[3]      # Diagnostics variables
STORAGE: auxevol_variables[1]  # Single-timelevel storage of variables needed for evolutions.

# The following scheduler is based on Lean/LeanBSSNMoL/schedule.ccl

schedule BaikalETK_Banner at STARTUP
{
  LANG: C
  OPTIONS: meta
} "Output ASCII art banner"

schedule BaikalETK_RegisterSlicing at STARTUP after BaikalETK_Banner
{
  LANG: C
  OPTIONS: meta
} "Register 3+1 slicing condition"

schedule BaikalETK_Symmetry_registration at BASEGRID
{
  LANG: C
  OPTIONS: Global
} "Register symmetries, the CartGrid3D way."

schedule BaikalETK_zero_rhss at BASEGRID after BaikalETK_Symmetry_registration
{
  LANG: C
} "Idea from Lean: set all rhs functions to zero to prevent spurious nans"

schedule BaikalETK_ADM_to_BSSN at CCTK_INITIAL after ADMBase_PostInitial
{
  LANG: C
  OPTIONS: Local
  SYNC: evol_variables
} "Convert initial data into BSSN variables"

schedule GROUP ApplyBCs as BaikalETK_ApplyBCs at CCTK_INITIAL after BaikalETK_ADM_to_BSSN
{
} "Apply boundary conditions"


# MoL: registration

schedule BaikalETK_MoL_registration in MoL_Register
{
  LANG: C
  OPTIONS: META
} "Register variables for MoL"


# MoL: compute RHSs, etc
"""
    if enable_stress_energy_source_terms == True:
        outstr += """
schedule BaikalETK_driver_BSSN_T4UU in MoL_CalcRHS as BaikalETK_T4UU before BaikalETK_BSSN_to_ADM
{
  LANG: C
} "MoL: Compute T4UU, needed for BSSN RHSs."
"""
    outstr += """
schedule BaikalETK_driver_pt1_BSSN_Ricci in MoL_CalcRHS as BaikalETK_Ricci before BaikalETK_RHS
{
  LANG: C
} "MoL: Compute Ricci tensor"

schedule BaikalETK_driver_pt2_BSSN_RHSs in MoL_CalcRHS as BaikalETK_RHS after BaikalETK_Ricci
{
  LANG: C
} "MoL: Evaluate BSSN RHSs"

schedule BaikalETK_NewRad in MoL_CalcRHS after BaikalETK_RHS
{
  LANG: C
} "NewRad boundary conditions, scheduled right after RHS eval."

schedule BaikalETK_floor_the_lapse in MoL_PostStep before BaikalETK_enforce_detgammahat_constraint before BC_Update
{
  LANG: C
} "Set lapse = max(lapse_floor, lapse)"

schedule BaikalETK_enforce_detgammahat_constraint in MoL_PostStep before BC_Update
{
  LANG: C
} "Enforce detgammabar = detgammahat (= 1 in Cartesian)"

schedule BaikalETK_BoundaryConditions_evolved_gfs in MoL_PostStep
{
  LANG: C
  OPTIONS: LEVEL
  SYNC: evol_variables
} "Apply boundary conditions and perform AMR+interprocessor synchronization"

schedule GROUP ApplyBCs as BaikalETK_ApplyBCs in MoL_PostStep after BaikalETK_BoundaryConditions_evolved_gfs
{
} "Group for applying boundary conditions"


# Next update ADM quantities

schedule BaikalETK_BSSN_to_ADM in MoL_PostStep after BaikalETK_ApplyBCs before ADMBase_SetADMVars
{
  LANG: C
  OPTIONS: Local
} "Perform BSSN-to-ADM conversion. Useful for diagnostics."

# Compute Hamiltonian & momentum constraints
"""
    if enable_stress_energy_source_terms == True:
        outstr += """
schedule BaikalETK_driver_BSSN_T4UU in MoL_PseudoEvolution before BaikalETK_BSSN_constraints
{
  LANG: C
  OPTIONS: Local
} "MoL_PseudoEvolution: Compute T4UU, needed for BSSN constraints"
"""
    outstr += """

schedule BaikalETK_BSSN_constraints in MoL_PseudoEvolution
{
  LANG: C
  OPTIONS: Local
} "Compute BSSN (Hamiltonian and momentum) constraints"

schedule BaikalETK_BoundaryConditions_aux_gfs in MoL_PseudoEvolution after BaikalETK_BSSN_constraints
{
  LANG: C
  OPTIONS: LEVEL
  SYNC: aux_variables
} "Enforce symmetry BCs in constraint computation"

"""
    if enable_stress_energy_source_terms == True:
        outstr += """
schedule BaikalETK_BSSN_to_ADM in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs
{
  LANG: C
  OPTIONS: Local
} "Perform BSSN-to-ADM conversion in MoL_PseudoEvolution. Needed for proper HydroBase integration."
"""
    outstr += """
schedule GROUP ApplyBCs as BaikalETK_auxgfs_ApplyBCs in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs
{
} "Apply boundary conditions"
"""
    with open(os.path.join(ThornName,"schedule.ccl"), "w") as file:
        file.write(outstr.replace("BaikalETK",ThornName))

Appending to BaikalETK_validate/BaikalETK_ETK_ccl_files_codegen.py

In [14]:

import BaikalETK_validate.BaikalETK_ETK_ccl_files_codegen as cclgen

for enable_stress_energy_source_terms in [True,False]:
    ThornName="Baikal"
    if enable_stress_energy_source_terms==False:
        ThornName="BaikalVacuum"
    cclgen.output_param_ccl(ThornName)
    cclgen.output_interface_ccl(ThornName,enable_stress_energy_source_terms)
    cclgen.output_schedule_ccl(ThornName,enable_stress_energy_source_terms)

Step 4: C driver functions for ETK registration & NRPy+-generated kernels [Back to top]¶

$$\label{cdrivers}$$

Now that we have constructed the basic C code kernels and the needed Einstein Toolkit ccl files, we next write the driver functions for registering BaikalETK within the Toolkit and the C code kernels. Each of these driver functions will be called directly from the thorn's schedule.ccl in the ETK.

Step 4.a: Needed ETK functions: Banner, Symmetry registration, Parameter sanity check, Method of Lines (`MoL`) registration, Boundary condition [Back to top]¶

$$\label{etkfunctions}$$

To-do: Parameter sanity check function. E.g., error should be thrown if `cctk_nghostzones[]` is set too small for the chosen finite-differencing order within NRPy+.¶

In [15]:

%%writefile $BaikalETKdir/BaikalETK_C_drivers_codegen.py

# Step 1: Import needed core NRPy+ and Python modules
from outputC import lhrh         # NRPy+: Core C code output module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
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 loop as lp                # NRPy+: Generate C code loops
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import os, sys                   # Standard Python modules for multiplatform OS-level functions
# We need the function keep_param__return_type() from this module:
import BaikalETK_validate.BaikalETK_ETK_ccl_files_codegen as ccl

make_code_defn_list = []
def append_to_make_code_defn_list(filename):
    if filename not in make_code_defn_list:
        make_code_defn_list.append(filename)
    return filename

Overwriting BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [16]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

def driver_C_codes(Csrcdict, ThornName,
                   rhs_list,evol_gfs_list,aux_gfs_list,auxevol_gfs_list,
                   LapseCondition = "OnePlusLog", enable_stress_energy_source_terms=False):
    # First the ETK banner code, proudly showing the NRPy+ banner
    import NRPy_logo as logo
    outstr = """
#include <stdio.h>

void BaikalETK_Banner()
{
    """
    logostr = logo.print_logo(print_to_stdout=False)
    outstr += "printf(\"BaikalETK: another Einstein Toolkit thorn generated by\\n\");\n"
    for line in logostr.splitlines():
        outstr += "    printf(\""+line+"\\n\");\n"
    outstr += "}\n"

    # Finally add C code string to dictionaries (Python dictionaries are immutable)
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("Banner.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [17]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Then the RegisterSlicing() function, needed for other ETK thorns
    outstr = """
#include "cctk.h"

#include "Slicing.h"

int BaikalETK_RegisterSlicing (void)
{
  Einstein_RegisterSlicing ("BaikalETK");
  return 0;
}"""

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("RegisterSlicing.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [18]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next BaikalETK_Symmetry_registration(): Register symmetries

    full_gfs_list = []
    full_gfs_list.extend(evol_gfs_list)
    full_gfs_list.extend(auxevol_gfs_list)
    full_gfs_list.extend(aux_gfs_list)

    outstr = """
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"

void BaikalETK_Symmetry_registration(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  // Stores gridfunction parity across x=0, y=0, and z=0 planes, respectively
  int sym[3];

  // Next register parities for each gridfunction based on its name
  //    (to ensure this algorithm is robust, gridfunctions with integers
  //     in their base names are forbidden in NRPy+).
"""
    outstr += ""
    for gfname in full_gfs_list:
        gfname_without_GFsuffix = gfname[:-2]
        # Do not add T4UU gridfunctions if enable_stress_energy_source_terms==False:
        if not (enable_stress_energy_source_terms == False and "T4UU" in gfname_without_GFsuffix):
            outstr += """
      // Default to scalar symmetry:
      sym[0] = 1; sym[1] = 1; sym[2] = 1;
      // Now modify sym[0], sym[1], and/or sym[2] as needed
      //    to account for gridfunction parity across
      //    x=0, y=0, and/or z=0 planes, respectively
"""
            # If gridfunction name does not end in a digit, by NRPy+ syntax, it must be a scalar
            if gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 1].isdigit() == False:
                outstr += "      // (this gridfunction is a scalar -- no need to change default sym[]'s!)\n"
            elif len(gfname_without_GFsuffix) > 2:
                # Rank-1 indexed expression (e.g., vector)
                if gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 2].isdigit() == False:
                    if int(gfname_without_GFsuffix[-1]) > 2:
                        print("Error: Found invalid gridfunction name: "+gfname)
                        sys.exit(1)
                    symidx = gfname_without_GFsuffix[-1]
                    if int(symidx) < 3:   outstr += "      sym[" + symidx + "] = -1;\n"
                # Rank-2 indexed expression
                elif gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 2].isdigit() == True:
                    if len(gfname_without_GFsuffix) > 3 and gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 3].isdigit() == True:
                        print("Error: Found a Rank-3 or above gridfunction: "+gfname+", which is at the moment unsupported.")
                        print("It should be easy to support this if desired.")
                        sys.exit(1)
                    symidx0 = gfname_without_GFsuffix[-2]
                    if "T4UU" in gfname: symidx0 = str(int(symidx0)-1) # E.g., T4UU23 is T4UUyz, corresponding to directions 1,2
                    if int(symidx0) >= 0: outstr += "      sym[" + symidx0 + "] *= -1;\n"
                    symidx1 = gfname_without_GFsuffix[-1]
                    if "T4UU" in gfname: symidx1 = str(int(symidx1)-1) # E.g., T4UU23 is T4UUyz, corresponding to directions 1,2
                    if int(symidx1) >= 0: outstr += "      sym[" + symidx1 + "] *= -1;\n"
            else:
                print("Don't know how you got this far with a gridfunction named "+gfname+", but I'll take no more of this nonsense.")
                print("   Please follow best-practices and rename your gridfunction to be more descriptive")
                sys.exit(1)
            outstr += "      SetCartSymVN(cctkGH, sym, \"BaikalETK::" + gfname + "\");\n"
    outstr += "}\n"

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("Symmetry_registration_oldCartGrid3D.c")] = \
        outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [19]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next set RHSs to zero
    outstr = """
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"

void BaikalETK_zero_rhss(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
"""
    set_rhss_to_zero = ""
    for gf in rhs_list:
        set_rhss_to_zero += gf+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;\n"

    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","","",],"",set_rhss_to_zero)
    outstr += "}\n"
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("zero_rhss.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [20]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next floor the lapse
    outstr = """
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"

#ifndef MAX
#define MAX(a,b) ( ((a) > (b)) ? (a) : (b) )

void BaikalETK_floor_the_lapse(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
"""
    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","","",],"","""
alphaGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)] = MAX(alphaGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)], lapse_floor);
""")
    outstr += """
}
#undef MAX
#endif\n"""
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("floor_the_lapse.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [21]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next registration with the Method of Lines thorn
    outstr = """
//--------------------------------------------------------------------------
// Register with the Method of Lines time stepper
// (MoL thorn, found in arrangements/CactusBase/MoL)
// MoL documentation located in arrangements/CactusBase/MoL/doc
//--------------------------------------------------------------------------
#include <stdio.h>

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

#include "Symmetry.h"

void BaikalETK_MoL_registration(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  CCTK_INT ierr = 0, group, rhs;

  // Register evolution & RHS gridfunction groups with MoL, so it knows

  group = CCTK_GroupIndex("BaikalETK::evol_variables");
  rhs = CCTK_GroupIndex("BaikalETK::evol_variables_rhs");
  ierr += MoLRegisterEvolvedGroup(group, rhs);

  if (ierr) CCTK_ERROR("Problems registering with MoL");
}
"""
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("MoL_registration.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [22]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next register with the boundary conditions thorns.
    # PART 1: Set BC type to "none" for all variables
    # Since we choose NewRad boundary conditions, we must register all
    #   gridfunctions to have boundary type "none". This is because
    #   NewRad is seen by the rest of the Toolkit as a modification to the
    #   RHSs.

    # This code is based on Kranc's McLachlan/ML_BSSN/src/Boundaries.cc code.
    outstr = """
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Faces.h"
#include "util_Table.h"
#include "Symmetry.h"

// Set `none` boundary conditions on BSSN RHSs, as these are set via NewRad.
void BaikalETK_BoundaryConditions_evolved_gfs(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;
"""
    for gf in evol_gfs_list:
        outstr += """
  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "BaikalETK::"""+gf+"""", "none");
  if (ierr < 0) CCTK_ERROR("Failed to register BC for BaikalETK::"""+gf+"""!");
"""
    outstr += """
}

// Set `flat` boundary conditions on BSSN constraints, similar to what Lean does.
void BaikalETK_BoundaryConditions_aux_gfs(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;

  CCTK_INT bndsize[6];
  CCTK_INT is_internal[6];
  CCTK_INT is_staggered[6];
  CCTK_INT shiftout[6];

  GetBoundarySpecification(6, bndsize, is_internal, is_staggered, shiftout);

"""
    for gf in aux_gfs_list:
        outstr += """
  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, bndsize[0], -1, "BaikalETK::"""+gf+"""", "flat");
  if (ierr < 0) CCTK_ERROR("Failed to register BC for BaikalETK::"""+gf+"""!");
"""
    outstr += "}\n"

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("BoundaryConditions.c")] = outstr.replace("BaikalETK",ThornName)

    # PART 2: Set C code for calling NewRad BCs
    #   As explained in lean_public/LeanBSSNMoL/src/calc_bssn_rhs.F90,
    #   the function NewRad_Apply takes the following arguments:
    #   NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower),
    #     which implement the boundary condition:
    #       var  =  var_at_infinite_r + u(r-var_char_speed*t)/r^var_radpower
    #  Obviously for var_radpower>0, var_at_infinite_r is the value of
    #    the variable at r->infinity. var_char_speed is the propagation
    #    speed at the outer boundary, and var_radpower is the radial
    #    falloff rate.

    outstr = """
#include <math.h>

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

void BaikalETK_NewRad(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

"""
    for gf in evol_gfs_list:
        var_at_infinite_r = "0.0"
        var_char_speed    = "1.0"
        var_radpower      = "1.0"

        if "alpha" in gf:
            var_at_infinite_r = "1.0"
            if LapseCondition == "OnePlusLog":
                var_char_speed = "sqrt(2.0)"
            else:
                pass # 1.0 (default) is fine
        if "aDD" in gf or "trK" in gf: # consistent with Lean code.
            var_radpower = "2.0"

        outstr += "  NewRad_Apply(cctkGH, "+gf+", "+gf.replace("GF","")+"_rhsGF, "+var_at_infinite_r+", "+var_char_speed+", "+var_radpower+");\n"
    outstr += "}\n"

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("BoundaryCondition_NewRad.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.b: BSSN $\leftrightarrow$ ADM conversions [Back to top]¶

$$\label{bssnadmconversions}$$

Step 4.b.i: ADM $\to$ BSSN [Back to top]¶

$$\label{admtobssn}$$

Initial data in the Einstein Toolkit are given in terms of ADM quantities, so a conversion is necessary so the quantities are in terms of BSSN variables used for evolving the initial data forward in time.

In [23]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # 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(ThornName,"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)

    # We ignore the return values for the following register_gridfunctions...() calls,
    #   as they are unused.
    gri.register_gridfunctions(                 "AUXEVOL", "alphaSphorCart")
    ixp.register_gridfunctions_for_single_rank1("AUXEVOL", "betaSphorCartU")
    ixp.register_gridfunctions_for_single_rank1("AUXEVOL", "BSphorCartU")
    ixp.register_gridfunctions_for_single_rank2("AUXEVOL", "gammaSphorCartDD", "sym01")
    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 BaikalETK_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

    # Output finite difference stencils as inlined expressions.
    #   We do this instead of outputting as FD functions, as this function
    #   does not take long to compile, and we have already output all the
    #   FD functions to file, so if this one contains new FD functions,
    #   the compile will fail.
    par.set_parval_from_str("finite_difference::enable_FD_functions", False)

    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", "betU0"), rhs=IDbetU[0]),
        lhrh(lhs=gri.gfaccess("in_gfs", "betU1"), rhs=IDbetU[1]),
        lhrh(lhs=gri.gfaccess("in_gfs", "betU2"), rhs=IDbetU[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)

    # ADM to BSSN, Part 3: Set up ADM to BSSN conversions for BSSN gridfunctions defined from
    #                      finite-difference derivatives: lambda^i, which is Gamma^i in non-covariant BSSN):
    outstr += """
    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);
"""

    path = os.path.join(ThornName,"src")
    BaikalETK_src_filelist = []
    for _root, _dirs, files in os.walk(path):  # _root, _dirs unused.
        for filename in files:
            BaikalETK_src_filelist.append(filename)
    BaikalETK_src_filelist.sort() # Sort the list in place.

    BSSN_FD_orders_output = []
    for filename in BaikalETK_src_filelist:
        if "BSSN_RHSs_" in filename:
            array = filename.replace(".","_").split("_")
            FDorder =  int(array[-2])
            if FDorder not in BSSN_FD_orders_output:
                BSSN_FD_orders_output.append(FDorder)
    BSSN_FD_orders_output.sort()

    for current_FD_order in BSSN_FD_orders_output:
        # Store original finite-differencing order:
        orig_FD_order = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")
        # Set new finite-differencing order:
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", current_FD_order)

        outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
        lambdaU_expressions = [lhrh(lhs=gri.gfaccess("in_gfs","lambdaU0"),rhs=IDlambdaU[0]),
                               lhrh(lhs=gri.gfaccess("in_gfs","lambdaU1"),rhs=IDlambdaU[1]),
                               lhrh(lhs=gri.gfaccess("in_gfs","lambdaU2"),rhs=IDlambdaU[2])]
        lambdaU_expressions_FDout = fin.FD_outputC("returnstring",lambdaU_expressions, outCparams)

        lambdafile = "ADM_to_BSSN__compute_lambdaU_FD_order_"+str(current_FD_order)+".h"
        with open(os.path.join(ThornName,"src",lambdafile),"w") as file:
            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],
                       ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],
                       ["1","1","1"],["#pragma omp parallel for","",""],"",lambdaU_expressions_FDout))

        outstr += "    if(FD_order == "+str(current_FD_order)+") {\n"
        outstr += "        #include \""+lambdafile+"\"\n"
        outstr += "    }\n"
        # Restore original FD order
        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", orig_FD_order)

    outstr += """
    ExtrapolateGammas(cctkGH,lambdaU0GF);
    ExtrapolateGammas(cctkGH,lambdaU1GF);
    ExtrapolateGammas(cctkGH,lambdaU2GF);
}
"""

    # Unregister the *SphorCartGF's.
    gri.glb_gridfcs_list = orig_glb_gridfcs_list

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("ADM_to_BSSN.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.b.ii: BSSN $\to$ ADM [Back to top]¶

$$\label{bssntoadm}$$

All modules (thorns) in the Einstein Toolkit that deal with spacetime quantities do so via the core ADMBase module, which assumes variables are written in ADM form. Therefore, in order for BaikalETK to interface properly with the rest of the Toolkit, its native BSSN variables must be converted to ADM quantities.

In [24]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    import BSSN.ADM_in_terms_of_BSSN as btoa
    import BSSN.BSSN_quantities as Bq

    btoa.ADM_in_terms_of_BSSN()
    Bq.BSSN_basic_tensors() # Gives us betaU & BU

    outstr = """
#include <math.h>

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

void BaikalETK_BSSN_to_ADM(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

"""
    btoa_lhrh = []
    for i in range(3):
        for j in range(i,3):
            btoa_lhrh.append(lhrh(lhs="g"+chr(ord('x')+i)+chr(ord('x')+j)+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]",
                                  rhs=btoa.gammaDD[i][j]))
    for i in range(3):
        for j in range(i,3):
            btoa_lhrh.append(lhrh(lhs="k"+chr(ord('x')+i)+chr(ord('x')+j)+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]",
                                  rhs=btoa.KDD[i][j]))
    btoa_lhrh.append(lhrh(lhs="alp[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]",rhs=Bq.alpha))

    for i in range(3):
        btoa_lhrh.append(lhrh(lhs="beta"+chr(ord('x')+i)+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]",
                                  rhs=Bq.betaU[i]))

    for i in range(3):
        btoa_lhrh.append(lhrh(lhs="dtbeta"+chr(ord('x')+i)+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]",
                                  rhs=Bq.BU[i]))

    outCparams = "preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False"
    bssn_to_adm_Ccode = fin.FD_outputC("returnstring",btoa_lhrh, 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","",""],"",bssn_to_adm_Ccode)

    outstr += "}\n"

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("BSSN_to_ADM.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.c: Evaluate BSSN right-hand-sides (RHSs) [Back to top]¶

$$\label{bssnrhss}$$

Step 4.c.i: Evaluate Ricci tensor [Back to top]¶

$$\label{ricci}$$

To slightly optimize the performance of BaikalETK's BSSN solver, we split the computation of the complicated expressions for the Ricci tensor $\\bar{R}_{ij}$ into its own function, and then use the result when evaluating the BSSN right-hand-side (RHS) expressions.

In [25]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    ###########################
    ###########################
    # BSSN_RHSs and Ricci driver functions
    ###########################
    common_includes = """
#include <math.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "SIMD/SIMD_intrinsics.h"
#include "finite_difference_functions.h"
"""
    common_preloop = """
    DECLARE_CCTK_ARGUMENTS;
    const CCTK_REAL NOSIMDinvdx0 = 1.0/CCTK_DELTA_SPACE(0);
    const REAL_SIMD_ARRAY invdx0 = ConstSIMD(NOSIMDinvdx0);
    const CCTK_REAL NOSIMDinvdx1 = 1.0/CCTK_DELTA_SPACE(1);
    const REAL_SIMD_ARRAY invdx1 = ConstSIMD(NOSIMDinvdx1);
    const CCTK_REAL NOSIMDinvdx2 = 1.0/CCTK_DELTA_SPACE(2);
    const REAL_SIMD_ARRAY invdx2 = ConstSIMD(NOSIMDinvdx2);
"""

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

In [26]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Output the driver code for computing the Ricci tensor:
    outstr = common_includes
    for order in BSSN_FD_orders_output:
        outstr += """extern void """+ThornName+"_"+"BSSN_Ricci_FD_order_"+str(order)+"(CCTK_ARGUMENTS);\n"
    outstr += """
void BaikalETK_driver_pt1_BSSN_Ricci(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    const CCTK_INT *FD_order = CCTK_ParameterGet("FD_order","BaikalETK",NULL);
"""
    for order in BSSN_FD_orders_output:
        outstr += "    if(*FD_order == "+str(order)+") {\n"
        outstr += "        "+ThornName+"_"+"BSSN_Ricci_FD_order_"+str(order)+"(CCTK_PASS_CTOC);\n"
        outstr += "    }\n"
    outstr += "} // END FUNCTION\n"
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("driver_pt1_BSSN_Ricci.c")] = outstr.replace("BaikalETK",ThornName)

    # Create functions for the largest C kernels (BSSN RHSs and Ricci) and output
    #    the .h files to .c files with function wrappers; delete original .h files
    for filename in BaikalETK_src_filelist:
        if ("BSSN_Ricci_FD_order_") in filename and (".h" in filename):
            outstr = common_includes + "void BaikalETK_"+filename.replace(".h","")+"(CCTK_ARGUMENTS) {\n"
            outstr += common_preloop
            with open(os.path.join(path,filename), "r") as currfile:
                outstr += currfile.read()
            # Now that we've inserted the contents of the kernel into this file,
            #     we delete the file containing the kernel
            os.remove(os.path.join(path,filename))
            outstr += "} // END FUNCTION\n"
            # Add C code string to dictionary (Python dictionaries are immutable)
            Csrcdict[append_to_make_code_defn_list(filename.replace(".h",".c"))] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.c.ii: Evaluate BSSN RHSs, using Ricci tensor as input [Back to top]¶

$$\label{bssnrhssricciinput}$$

Next we construct the driver function for evaluating the BSSN RHSs, which make use of the Ricci tensor $\bar{R}_{ij}$, which has just been computed.

In [27]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    ###########################
    # Output BSSN RHSs driver function
    outstr = common_includes
    for filename in BaikalETK_src_filelist:
        if ("BSSN_RHSs_" in filename) and (".h" in filename):
            outstr += """extern void """ + ThornName+"_"+filename.replace(".h", "(CCTK_ARGUMENTS);") + "\n"

    outstr += """
void BaikalETK_driver_pt2_BSSN_RHSs(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;

    const CCTK_INT *FD_order = CCTK_ParameterGet("FD_order","BaikalETK",NULL);

"""
    for filename in BaikalETK_src_filelist:
        if ("BSSN_RHSs_" in filename) and (".h" in filename):
            array = filename.replace(".", "_").split("_")
            outstr += "    if(*FD_order == " + str(array[-2]) + ") {\n"
            outstr += "        " + ThornName+"_"+filename.replace(".h", "(CCTK_PASS_CTOC);") + "\n"
            outstr += "    }\n"
    outstr += "} // END FUNCTION\n"
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("driver_pt2_BSSN_RHSs.c")] = outstr.replace("BaikalETK", ThornName)

    def SIMD_declare_C_params():
        SIMD_declare_C_params_str = ""
        for i in range(len(par.glb_Cparams_list)):
            # keep_param is a boolean indicating whether we should accept or reject
            #    the parameter. singleparstring will contain the string indicating
            #    the variable type.
            keep_param, singleparstring = ccl.keep_param__return_type(par.glb_Cparams_list[i])

            if (keep_param) and ("CCTK_REAL" in singleparstring):
                parname = par.glb_Cparams_list[i].parname
                SIMD_declare_C_params_str += "    const "+singleparstring + "*NOSIMD"+parname+\
                                             " = CCTK_ParameterGet(\""+parname+"\",\"BaikalETK\",NULL);\n"
                SIMD_declare_C_params_str += "    const REAL_SIMD_ARRAY "+parname+" = ConstSIMD(*NOSIMD"+parname+");\n"
        return SIMD_declare_C_params_str

    # Create functions for the largest C kernels (BSSN RHSs and Ricci) and output
    #    the .h files to .c files with function wrappers; delete original .h files
    path = os.path.join(ThornName, "src")
    for filename in BaikalETK_src_filelist:
        if ("BSSN_RHSs_" in filename) and (".h" in filename):
            outstr = common_includes + "void BaikalETK_"+filename.replace(".h","")+"(CCTK_ARGUMENTS) {\n"
            outstr += common_preloop+SIMD_declare_C_params()
            with open(os.path.join(path,filename), "r") as currfile:
                outstr += currfile.read()
            # Now that we've inserted the contents of the kernel into this file,
            #     we delete the file containing the kernel
            os.remove(os.path.join(path,filename))
            outstr += "} // END FUNCTION\n"
            # Add C code string to dictionary (Python dictionaries are immutable)
            Csrcdict[append_to_make_code_defn_list(filename.replace(".h",".c"))] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.d: Enforcing conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\det{\hat{\gamma}_{ij}}=1$) [Back to top]¶

$$\label{enforcegammahatconstraint}$$

Here we construct the driver function for enforcing the conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint. The BSSN equations are not strongly hyperbolic if this condition is not set.

In [28]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next, the driver for enforcing detgammabar = detgammahat constraint:
    outstr = common_includes + """
void BaikalETK_enforce_detgammahat_constraint(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    #include "enforcedetgammabar_constraint.h"
}
"""
    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("driver_enforcedetgammabar_constraint.c")] = \
        outstr.replace("BaikalETK", ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.e: Diagnostics: Computing Hamiltonian & momentum constraints [Back to top]¶

$$\label{diagnostics}$$

The BSSN Hamiltonian & momentum constraints are useful diagnostics of a numerical-relativity calculation's health, as both should converge to zero with increasing numerical resolution. Here we construct the driver function.

In [29]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    # Next, the driver for computing the BSSN Hamiltonian & momentum constraints
    outstr = """
#include <math.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "SIMD/SIMD_intrinsics.h" // Contains needed definition of REAL_SIMD_ARRAY
#include "finite_difference_functions.h"

void BaikalETK_BSSN_constraints(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    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);
"""
    for filename in BaikalETK_src_filelist:
        if "BSSN_constraints_" in filename:
            array = filename.replace(".","_").split("_")
            outstr += "    if(FD_order == "+str(array[-2])+") {\n"
            outstr += "        #include \""+filename+"\"\n"
            outstr += "    }\n"
    outstr += "}\n"

    # Add C code string to dictionary (Python dictionaries are immutable)
    Csrcdict[append_to_make_code_defn_list("driver_BSSN_constraints.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.f: `driver_BSSN_T4UU()`: Compute $T^{\mu\nu}$ from `TmunuBase`'s $T_{\mu\nu}$ [Back to top]¶

$$\label{t4uu}$$

Here we implement $T^{\mu\nu} = g^{\mu \delta} g^{\nu \gamma} T_{\delta \gamma}.$

In [30]:

%%writefile -a $BaikalETKdir/BaikalETK_C_drivers_codegen.py

    if enable_stress_energy_source_terms == True:
        # Declare T4DD as a set of gridfunctions. These won't
        #    actually appear in interface.ccl, as interface.ccl
        #    was set above. Thus before calling the code output
        #    by FD_outputC(), we'll have to set pointers
        #    to the actual gridfunctions they reference.
        #    (In this case the eTab's.)
        T4DD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","T4DD","sym01",DIM=4)
        import BSSN.ADMBSSN_tofrom_4metric as AB4m
        AB4m.g4UU_ito_BSSN_or_ADM("BSSN")

        T4UUraised = ixp.zerorank2(DIM=4)
        for mu in range(4):
            for nu in range(4):
                for delta in range(4):
                    for gamma in range(4):
                        T4UUraised[mu][nu] += AB4m.g4UU[mu][delta]*AB4m.g4UU[nu][gamma]*T4DD[delta][gamma]

        T4UU_expressions = [
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU00"),rhs=T4UUraised[0][0]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU01"),rhs=T4UUraised[0][1]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU02"),rhs=T4UUraised[0][2]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU03"),rhs=T4UUraised[0][3]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU11"),rhs=T4UUraised[1][1]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU12"),rhs=T4UUraised[1][2]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU13"),rhs=T4UUraised[1][3]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU22"),rhs=T4UUraised[2][2]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU23"),rhs=T4UUraised[2][3]),
            lhrh(lhs=gri.gfaccess("in_gfs","T4UU33"),rhs=T4UUraised[3][3])]

        outCparams = "outCverbose=False,includebraces=False,preindent=2,enable_SIMD=True"
        T4UUstr = fin.FD_outputC("returnstring",T4UU_expressions, outCparams)
        T4UUstr_loop = lp.loop(["i2","i1","i0"],["0","0","0"],["cctk_lsh[2]","cctk_lsh[1]","cctk_lsh[0]"],
                               ["1","1","SIMD_width"],["#pragma omp parallel for","",""],"",T4UUstr)

        outstr = """
#include <math.h>

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

#include "SIMD/SIMD_intrinsics.h"

void BaikalETK_driver_BSSN_T4UU(CCTK_ARGUMENTS) {
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    const CCTK_REAL *restrict T4DD00GF = eTtt;
    const CCTK_REAL *restrict T4DD01GF = eTtx;
    const CCTK_REAL *restrict T4DD02GF = eTty;
    const CCTK_REAL *restrict T4DD03GF = eTtz;
    const CCTK_REAL *restrict T4DD11GF = eTxx;
    const CCTK_REAL *restrict T4DD12GF = eTxy;
    const CCTK_REAL *restrict T4DD13GF = eTxz;
    const CCTK_REAL *restrict T4DD22GF = eTyy;
    const CCTK_REAL *restrict T4DD23GF = eTyz;
    const CCTK_REAL *restrict T4DD33GF = eTzz;
"""+T4UUstr_loop+"""
}\n"""

        # Add C code string to dictionary (Python dictionaries are immutable)
        Csrcdict[append_to_make_code_defn_list("driver_BSSN_T4UU.c")] = outstr.replace("BaikalETK",ThornName)

Appending to BaikalETK_validate/BaikalETK_C_drivers_codegen.py

Step 4.g: Output all C driver functions needed for `Baikal`/`BaikalVacuum` [Back to top]¶

$$\label{outcdrivers}$$

First we call the above functions (output above to the BaikalETK_validate.BaikalETK_C_drivers_codegen Python module) to store all needed driver C files to a Python dictionary, then we simply outputs the dictionary to the appropriate files.

In [31]:

import BaikalETK_validate.BaikalETK_C_drivers_codegen as driver

# The following Python dictionaries consist of a key, which is the filename
#    in the thorn's src/ directory (e.g., "driver_BSSN_constraints.c"),
#    and a value, which is the corresponding source code, stored as a
#    Python string.
Vac_Csrcdict = {}
Reg_Csrcdict = {}

# We'll need lists of gridfunctions for these driver functions
evol_gfs_list    = cclgen.evol_gfs_list
aux_gfs_list     = cclgen.aux_gfs_list
auxevol_gfs_list = cclgen.auxevol_gfs_list

# Generate driver codes for Baikal thorn (i.e., populate the Reg_Csrcdict dictionary)
driver.driver_C_codes(Reg_Csrcdict, "Baikal",
                               cclgen.rhs_list,cclgen.evol_gfs_list,cclgen.aux_gfs_list,cclgen.auxevol_gfs_list,
                               LapseCondition = LapseCondition, enable_stress_energy_source_terms=True)

# Generate driver codes for BaikalVacuum thorn (i.e., populate the Vac_Csrcdict dictionary)
driver.driver_C_codes(Vac_Csrcdict, "BaikalVacuum",
                               cclgen.rhs_list,cclgen.evol_gfs_list,cclgen.aux_gfs_list,cclgen.auxevol_gfs_list,
                               LapseCondition = LapseCondition, enable_stress_energy_source_terms=False)

# Next we output the contents of the Reg_Csrcdict and
#   Vac_Csrcdict dictionaries to files in the respective
#   thorns' directories.
for key,val in Reg_Csrcdict.items():
    with open(os.path.join("Baikal","src",key),"w") as file:
        file.write(val)
for key,val in Vac_Csrcdict.items():
    with open(os.path.join("BaikalVacuum","src",key),"w") as file:
        file.write(val)

Step 4.h: `make.code.defn`: List of all C driver functions needed to compile `BaikalETK` [Back to top]¶

$$\label{makecodedefn}$$

When constructing each C code driver function above, we called the append_to_make_code_defn_list() function, which built a list of each C code driver file. We'll now add each of those files to the make.code.defn file, used by the Einstein Toolkit's build system.

In [32]:

# Finally output the thorns' make.code.defn files, consisting of
#   a list of all C codes in the above dictionaries. This is
#   part of the ETK build system so that these files are output.

def output_make_code_defn(dictionary, ThornName):
    with open(os.path.join(ThornName, "src", "make.code.defn"), "w") as file:
        file.write("""
# Main make.code.defn file for thorn """+ThornName+"""

# Source files in this directory
SRCS =""")
        filestring = ""

        list_of_C_driver_files = list(dictionary.keys())
        for i in range(len(list_of_C_driver_files)):
            filestring += "      "+list_of_C_driver_files[i]
            if i != len(list_of_C_driver_files)-1:
                filestring += " \\\n"
            else:
                filestring += "\n"
        file.write(filestring)

output_make_code_defn(Reg_Csrcdict, "Baikal")
output_make_code_defn(Vac_Csrcdict, "BaikalVacuum")

print("Finished generating Baikal and BaikalVacuum Einstein Toolkit thorns!")

Finished generating Baikal and BaikalVacuum Einstein Toolkit thorns!

Step 5: Code validation [Back to top]¶

$$\label{code_validation}$$

Here we will show plots demonstrating good agreement between BaikalETK and, e.g., McLachlan (another, trusted ETK thorn).

Step 5.a: Validation against BaikalETK/ Python modules [Back to top]¶

$$\label{self_validation}$$

As a self-validation check, we verify that the ETK thorns just generated here perfectly match those generated by the BaikalETK.BaikalETK_main_codegen_driver Python module.

In [33]:

# First move the Baikal and Baikal_vacuum thorns just generated by this notebook to
#  the validate/ subdirectory:
shutil.rmtree("validate", ignore_errors=True)
cmd.mkdir("validate")
shutil.move("Baikal","validate")
shutil.move("BaikalVacuum","validate")

# Next generate Baikal and BaikalVacuum directly from the Python module.
#  These will output to the Baikal and BaikalVacuum directories.
!python3 BaikalETK/BaikalETK_main_codegen_driver.py

***************************************
Starting parallel C kernel codegen...
***************************************
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for Ricci tensor...
Generating symbolic expressions for Ricci tensor...
Generating symbolic expressions for BSSN RHSs...
Generating symbolic expressions for Ricci tensor...
Generating symbolic expressions for Ricci tensor...
Generating symbolic expressions for Ricci tensor...
Generating optimized C code (FD_order=4) for gamma constraint. May take a while, depending on CoordSystem.
Generating optimized C code (FD_order=4) for gamma constraint. May take a while, depending on CoordSystem.
(BENCH) Finished gamma constraint C codegen (FD_order=4) in 0.07429814338684082 seconds.
(BENCH) Finished gamma constraint C codegen (FD_order=4) in 0.07741808891296387 seconds.
(BENCH) Finished BSSN RHS symbolic expressions in 0.7707383632659912 seconds.
Generating C code for BSSN RHSs (FD_order=6,Tmunu=False) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 0.7784299850463867 seconds.
Generating C code for BSSN RHSs (FD_order=4,Tmunu=False) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 0.7759640216827393 seconds.
Generating C code for BSSN RHSs (FD_order=8,Tmunu=False) in Cartesian coordinates.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
(BENCH) Finished Ricci symbolic expressions in 1.3086082935333252 seconds.
Generating C code for Ricci tensor (FD_order=6) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 1.2989954948425293 seconds.
Generating C code for Ricci tensor (FD_order=4) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 1.3120527267456055 seconds.
Generating C code for Ricci tensor (FD_order=8) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 1.3443379402160645 seconds.
Generating C code for Ricci tensor (FD_order=4) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 1.3933684825897217 seconds.
Generating C code for BSSN RHSs (FD_order=4,Tmunu=True) in Cartesian coordinates.
(BENCH) Finished BSSN RHS symbolic expressions in 1.3973097801208496 seconds.
Generating C code for BSSN RHSs (FD_order=2,Tmunu=True) in Cartesian coordinates.
(BENCH) Finished Ricci symbolic expressions in 1.5551927089691162 seconds.
Generating C code for Ricci tensor (FD_order=2) in Cartesian coordinates.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.
(BENCH) Finished BSSN_RHS C codegen (FD_order=4,Tmunu=False) in 9.570033073425293 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=6,Tmunu=False) in 10.808182716369629 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=2,Tmunu=True) in 11.086544752120972 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=4,Tmunu=True) in 11.621188163757324 seconds.
(BENCH) Finished BSSN_RHS C codegen (FD_order=8,Tmunu=False) in 13.593563318252563 seconds.
(BENCH) Finished Ricci C codegen (FD_order=4) in 16.840373754501343 seconds.
(BENCH) Finished Ricci C codegen (FD_order=4) in 16.979881763458252 seconds.
(BENCH) Finished Ricci C codegen (FD_order=6) in 17.595633506774902 seconds.
(BENCH) Finished Ricci C codegen (FD_order=2) in 18.98304510116577 seconds.
(BENCH) Finished Ricci C codegen (FD_order=8) in 19.302734851837158 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=6,Tmunu=False) in 20.034168004989624 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=2,Tmunu=True) in 19.196338176727295 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=4,Tmunu=True) in 18.98688793182373 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=4,Tmunu=False) in 21.114547729492188 seconds.
(BENCH) Finished Hamiltonian & momentum constraint C codegen (FD_order=8,Tmunu=False) in 21.257333755493164 seconds.
Finished C kernel codegen for Baikal and BaikalVacuum in 22.234663009643555 seconds.
Finished generating Baikal and BaikalVacuum Einstein Toolkit thorns!

In [34]:

# Then compare all files generated by this notebook
#      (output moved in previous code cell to validate/)
#      and the separate Python module (output to Baikal
#      and BaikalVacuum).
import difflib

def compare_two_files(filepath1,filepath2):
    with open(filepath1) as file1, open(filepath2) as file2:
        # Read the lines of each file
        file1_lines = file1.readlines()
        file2_lines = file2.readlines()
        num_diffs = 0
        file1_lines_noleadingwhitespace = []
        for line in file1_lines:
            if line.strip() == "": # If the line contains only whitespace, remove all leading whitespace
                file1_lines_noleadingwhitespace.append(line.lstrip())
            else:
                file1_lines_noleadingwhitespace.append(line)
        file2_lines_noleadingwhitespace = []
        for line in file2_lines:
            if line.strip() == "": # If the line contains only whitespace, remove all leading whitespace
                file2_lines_noleadingwhitespace.append(line.lstrip())
            else:
                file2_lines_noleadingwhitespace.append(line)
        for line in difflib.unified_diff(file1_lines_noleadingwhitespace, file2_lines_noleadingwhitespace,
                                         fromfile=filepath1,
                                         tofile  =filepath2):
            sys.stdout.writelines(line)
            num_diffs = num_diffs + 1
        if num_diffs == 0:
            print("PASSED: "+filepath1+" matches trusted version")
        else:
            print("FAILED (see diff above): "+filepath1+" does NOT match trusted version")

import os
for dirpath, dirnames, filenames in os.walk("validate"):
    if len(filenames) > 1:
        filenames.sort()
        for file in [os.path.join(dirpath, file) for file in filenames]:
            compare_two_files(file,file.replace(os.path.join("validate/"),""))

# print("Ignoring lines with only whitespace:")
# for file in ["BaikalETK_C_drivers_codegen.py","BaikalETK_C_kernels_codegen.py","BaikalETK_ETK_ccl_files_codegen.py"]:
#     compare_two_files(BaikalETKdir,"BaikalETK",file)

PASSED: validate/BaikalVacuum/interface.ccl matches trusted version
PASSED: validate/BaikalVacuum/param.ccl matches trusted version
PASSED: validate/BaikalVacuum/schedule.ccl matches trusted version
PASSED: validate/BaikalVacuum/src/ADM_to_BSSN.c matches trusted version
PASSED: validate/BaikalVacuum/src/ADM_to_BSSN__compute_lambdaU_FD_order_4.h matches trusted version
PASSED: validate/BaikalVacuum/src/ADM_to_BSSN__compute_lambdaU_FD_order_6.h matches trusted version
PASSED: validate/BaikalVacuum/src/ADM_to_BSSN__compute_lambdaU_FD_order_8.h matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_RHSs_enable_Tmunu_False_FD_order_4.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_RHSs_enable_Tmunu_False_FD_order_6.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_RHSs_enable_Tmunu_False_FD_order_8.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_Ricci_FD_order_4.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_Ricci_FD_order_6.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_Ricci_FD_order_8.c matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_constraints_enable_Tmunu_False_FD_order_4.h matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_constraints_enable_Tmunu_False_FD_order_6.h matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_constraints_enable_Tmunu_False_FD_order_8.h matches trusted version
PASSED: validate/BaikalVacuum/src/BSSN_to_ADM.c matches trusted version
PASSED: validate/BaikalVacuum/src/Banner.c matches trusted version
PASSED: validate/BaikalVacuum/src/BoundaryCondition_NewRad.c matches trusted version
PASSED: validate/BaikalVacuum/src/BoundaryConditions.c matches trusted version
PASSED: validate/BaikalVacuum/src/MoL_registration.c matches trusted version
PASSED: validate/BaikalVacuum/src/RegisterSlicing.c matches trusted version
PASSED: validate/BaikalVacuum/src/Symmetry_registration_oldCartGrid3D.c matches trusted version
PASSED: validate/BaikalVacuum/src/driver_BSSN_constraints.c matches trusted version
PASSED: validate/BaikalVacuum/src/driver_enforcedetgammabar_constraint.c matches trusted version
PASSED: validate/BaikalVacuum/src/driver_pt1_BSSN_Ricci.c matches trusted version
PASSED: validate/BaikalVacuum/src/driver_pt2_BSSN_RHSs.c matches trusted version
PASSED: validate/BaikalVacuum/src/enforcedetgammabar_constraint.h matches trusted version
PASSED: validate/BaikalVacuum/src/finite_difference_functions.h matches trusted version
PASSED: validate/BaikalVacuum/src/floor_the_lapse.c matches trusted version
PASSED: validate/BaikalVacuum/src/make.code.defn matches trusted version
PASSED: validate/BaikalVacuum/src/zero_rhss.c matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_inner_read0.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_inner_read1.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_inner_read2.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_outer_read0.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_outer_read1.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__SIMD_outer_read2.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__declare.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__define.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__freemem.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__malloc.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__read0.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__read1.h matches trusted version
PASSED: validate/BaikalVacuum/src/rfm_files/rfm_struct__read2.h matches trusted version
PASSED: validate/Baikal/interface.ccl matches trusted version
PASSED: validate/Baikal/param.ccl matches trusted version
PASSED: validate/Baikal/schedule.ccl matches trusted version
PASSED: validate/Baikal/src/ADM_to_BSSN.c matches trusted version
PASSED: validate/Baikal/src/ADM_to_BSSN__compute_lambdaU_FD_order_2.h matches trusted version
PASSED: validate/Baikal/src/ADM_to_BSSN__compute_lambdaU_FD_order_4.h matches trusted version
PASSED: validate/Baikal/src/BSSN_RHSs_enable_Tmunu_True_FD_order_2.c matches trusted version
PASSED: validate/Baikal/src/BSSN_RHSs_enable_Tmunu_True_FD_order_4.c matches trusted version
PASSED: validate/Baikal/src/BSSN_Ricci_FD_order_2.c matches trusted version
PASSED: validate/Baikal/src/BSSN_Ricci_FD_order_4.c matches trusted version
PASSED: validate/Baikal/src/BSSN_constraints_enable_Tmunu_True_FD_order_2.h matches trusted version
PASSED: validate/Baikal/src/BSSN_constraints_enable_Tmunu_True_FD_order_4.h matches trusted version
PASSED: validate/Baikal/src/BSSN_to_ADM.c matches trusted version
PASSED: validate/Baikal/src/Banner.c matches trusted version
PASSED: validate/Baikal/src/BoundaryCondition_NewRad.c matches trusted version
PASSED: validate/Baikal/src/BoundaryConditions.c matches trusted version
PASSED: validate/Baikal/src/MoL_registration.c matches trusted version
PASSED: validate/Baikal/src/RegisterSlicing.c matches trusted version
PASSED: validate/Baikal/src/Symmetry_registration_oldCartGrid3D.c matches trusted version
PASSED: validate/Baikal/src/driver_BSSN_T4UU.c matches trusted version
PASSED: validate/Baikal/src/driver_BSSN_constraints.c matches trusted version
PASSED: validate/Baikal/src/driver_enforcedetgammabar_constraint.c matches trusted version
PASSED: validate/Baikal/src/driver_pt1_BSSN_Ricci.c matches trusted version
PASSED: validate/Baikal/src/driver_pt2_BSSN_RHSs.c matches trusted version
PASSED: validate/Baikal/src/enforcedetgammabar_constraint.h matches trusted version
PASSED: validate/Baikal/src/finite_difference_functions.h matches trusted version
PASSED: validate/Baikal/src/floor_the_lapse.c matches trusted version
PASSED: validate/Baikal/src/make.code.defn matches trusted version
PASSED: validate/Baikal/src/zero_rhss.c matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_inner_read0.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_inner_read1.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_inner_read2.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_outer_read0.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_outer_read1.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__SIMD_outer_read2.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__declare.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__define.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__freemem.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__malloc.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__read0.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__read1.h matches trusted version
PASSED: validate/Baikal/src/rfm_files/rfm_struct__read2.h matches trusted version

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

In [35]:

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

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

BaikalETK: NRPy+-Based BSSN Solvers for the Einstein Toolkit¶