Notebook

`GiRaFFE_NRPy`: Source Terms¶

Author: Patrick Nelson¶

Notebook Status: Validated

Validation Notes: This code produces the expected results for generated functions.

This module presents the functionality of GiRaFFE_NRPy_Source_Terms.py.¶

Introduction:¶

This writes and documents the C code that GiRaFFE_NRPy uses to compute the source terms for the right-hand sides of the evolution equations for the unstaggered prescription.

The equations themselves are already coded up in other functions; however, for the $\tilde{S}_i$ source term, we will need derivatives of the metric. It will be most efficient and accurate to take them using the interpolated metric values that we will have calculated anyway; however, we will need to write our derivatives in a nonstandard way within NRPy+ in order to take advantage of this, writing our own code for memory access.

Table of Contents¶

$\label{toc}$

This notebook is organized as follows

Step 1: The $\tilde{S}_i$ source term
Step 2: Code Validation against original C code
Step 3: Output this notebook to $\LaTeX$ -formatted PDF file

In [1]:

# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

import cmdline_helper as cmd
outdir = os.path.join("GiRaFFE_NRPy","GiRaFFE_Ccode_validation","RHSs")
cmd.mkdir(outdir)

Step 1: The $\tilde{S}_i$ source term [Back to top]¶

$\label{stilde_source}$

We start in the usual way - import the modules we need. We will also import the Levi-Civita symbol from indexedexp.py and use it to set the Levi-Civita tensor $\epsilon^{ijk} = [ijk]/\sqrt{\gamma}$ .

In [2]:

# Step 1: The StildeD RHS *source* term
from outputC import outputC, outCfunction, add_to_Cfunction_dict # NRPy+: Core C code output module
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import GRHD.equations as GRHD    # NRPy+: Generate general relativistic hydrodynamics equations
import GRFFE.equations as GRFFE  # NRPy+: Generate general relativistic force-free electrodynamics equations

thismodule = "GiRaFFE_NRPy_Source_Terms"

def generate_memory_access_code():
    # There are several pieces of C code that we will write ourselves because we need to do things
    # a little bit outside of what NRPy+ is built for.
    # First, we will write general memory access. We will read in values from memory at a given point
    # for each quantity we care about.
    global general_access
    general_access = ""
    for var in ["GAMMADD00", "GAMMADD01", "GAMMADD02",
                "GAMMADD11", "GAMMADD12", "GAMMADD22",
                "BETAU0", "BETAU1", "BETAU2","ALPHA",
                "BU0","BU1","BU2",
                "VALENCIAVU0","VALENCIAVU1","VALENCIAVU2"]:
        lhsvar = var.lower().replace("dd","DD").replace("u","U").replace("bU","BU").replace("valencia","Valencia")
        # e.g.,
        # const REAL gammaDD00dD0 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)];
        general_access += "const REAL "+lhsvar+" = auxevol_gfs[IDX4S("+var+"GF,i0,i1,i2)];\n"

    # This quick function returns a nearby point for memory access. We need this because derivatives are not local operations.
    def idxp1(dirn):
        if dirn==0:
            return "i0+1,i1,i2"
        if dirn==1:
            return "i0,i1+1,i2"
        if dirn==2:
            return "i0,i1,i2+1"

    # Next we evaluate needed derivatives of the metric, based on their values at cell faces
    global metric_deriv_access
    metric_deriv_access = []
    for dirn in range(3):
        metric_deriv_access.append("")
        for var in ["GAMMA_FACEDDdD00", "GAMMA_FACEDDdD01", "GAMMA_FACEDDdD02",
                    "GAMMA_FACEDDdD11", "GAMMA_FACEDDdD12", "GAMMA_FACEDDdD22",
                    "BETA_FACEUdD0", "BETA_FACEUdD1", "BETA_FACEUdD2","ALPHA_FACEdD"]:
            lhsvar = var.lower().replace("dddd","DDdD").replace("udd","UdD").replace("dd","dD").replace("u","U").replace("_face","")
            rhsvar = var.replace("dD","")
            # e.g.,
            # const REAL gammaDDdD000 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx0;
            metric_deriv_access[dirn] += "const REAL "+lhsvar+str(dirn)+" = (auxevol_gfs[IDX4S("+rhsvar+"GF,"+idxp1(dirn)+")]-auxevol_gfs[IDX4S("+rhsvar+"GF,i0,i1,i2)])/dxx"+str(dirn)+";\n"
        metric_deriv_access[dirn] += "REAL Stilde_rhsD"+str(dirn)+";\n"

    # This creates the C code that writes to the Stilde_rhs direction specified.
    global write_final_quantity
    write_final_quantity = []
    for dirn in range(3):
        write_final_quantity.append("")
        write_final_quantity[dirn] += "rhs_gfs[IDX4S(STILDED"+str(dirn)+"GF,i0,i1,i2)] += Stilde_rhsD"+str(dirn)+";"

def write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi):
    generate_memory_access_code()
    # First, we declare some dummy tensors that we will use for the codegen.
    gammaDDdD  = ixp.declarerank3("gammaDDdD","sym01",DIM=3)
    betaUdD = ixp.declarerank2("betaUdD","nosym",DIM=3)
    alphadD = ixp.declarerank1("alphadD",DIM=3)

    # We need to rerun a few of these functions with the reset lists to make sure these functions
    # don't cheat by using analytic expressions
    GRHD.compute_sqrtgammaDET(gammaDD)
    GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)
    GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)
    GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
    GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)
    GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)
    GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)
    for i in range(3):
        desc = "Adds the source term to StildeD"+str(i)+"."
        name = "calculate_StildeD"+str(i)+"_source_term"
        outCfunction(
            outfile  = os.path.join(outdir,name+".h"), desc=desc, name=name,
            params   ="const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs",
            body     = general_access \
                      +metric_deriv_access[i]\
                      +outputC(GRHD.S_tilde_source_termD[i],"Stilde_rhsD"+str(i),"returnstring",params=outCparams)\
                      +write_final_quantity[i],
            loopopts ="InteriorPoints",
            rel_path_to_Cparams=os.path.join("../"))

def add_to_Cfunction_dict__functions_for_StildeD_source_term(outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi,
                                                             includes=None, rel_path_to_Cparams=os.path.join("../"),
                                                             path_from_rootsrcdir_to_this_Cfunc=os.path.join("RHSs/")):
    generate_memory_access_code()
    # First, we declare some dummy tensors that we will use for the codegen.
    gammaDDdD  = ixp.declarerank3("gammaDDdD","sym01",DIM=3)
    betaUdD = ixp.declarerank2("betaUdD","nosym",DIM=3)
    alphadD = ixp.declarerank1("alphadD",DIM=3)

    # We need to rerun a few of these functions with the reset lists to make sure these functions
    # don't cheat by using analytic expressions
    GRHD.compute_sqrtgammaDET(gammaDD)
    GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)
    GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)
    GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
    GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)
    GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)
    GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)
    for i in range(3):
        desc = "Adds the source term to StildeD"+str(i)+"."
        name = "calculate_StildeD"+str(i)+"_source_term"
        params   ="const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs"
        body     = general_access \
                  +metric_deriv_access[i]\
                  +outputC(GRHD.S_tilde_source_termD[i],"Stilde_rhsD"+str(i),"returnstring",params=outCparams)\
                  +write_final_quantity[i]
        loopopts ="InteriorPoints"
        add_to_Cfunction_dict(
            includes=includes,
            desc=desc,
            name=name, params=params,
            body=body, loopopts=loopopts,
            path_from_rootsrcdir_to_this_Cfunc=path_from_rootsrcdir_to_this_Cfunc,
            rel_path_to_Cparams=rel_path_to_Cparams)

Step 2: Code Validation against original C code [Back to top]¶

$\label{code_validation}$

To validate the code in this tutorial we check for agreement between the files

that were written in this tutorial and
those that are stored in GiRaFFE_NRPy/GiRaFFE_Ccode_library or generated by GiRaFFE_NRPy_A2B.py

In [3]:

# Declare gridfunctions necessary to generate the C code:
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01",DIM=3)
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU",DIM=3)
alpha = gri.register_gridfunctions("AUXEVOL","alpha",DIM=3)
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU",DIM=3)
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU",DIM=3)
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD",DIM=3)
# Declare this symbol:
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")

# First, we generate the file using the functions written in this notebook:
outCparams = "outCverbose=False"
write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)

# Define the directory that we wish to validate against:
valdir = os.path.join("GiRaFFE_NRPy","GiRaFFE_Ccode_library","RHSs")
cmd.mkdir(valdir)

import GiRaFFE_NRPy.GiRaFFE_NRPy_Source_Terms as source
source.write_out_functions_for_StildeD_source_term(valdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)

import difflib
import sys

print("Printing difference between original C code and this code...")
# Open the files to compare
files = ["calculate_StildeD0_source_term.h","calculate_StildeD1_source_term.h","calculate_StildeD2_source_term.h"]

for file in files:
    print("Checking file " + file)
    with open(os.path.join(valdir,file)) as file1, open(os.path.join(outdir,file)) as file2:
        # Read the lines of each file
        file1_lines = file1.readlines()
        file2_lines = file2.readlines()
        num_diffs = 0
        for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir+file), tofile=os.path.join(outdir+file)):
            sys.stdout.writelines(line)
            num_diffs = num_diffs + 1
        if num_diffs == 0:
            print("No difference. TEST PASSED!")
        else:
            print("ERROR: Disagreement found with .py file. See differences above.")
            sys.exit(1)

Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_validation\RHSs\calculate_StildeD0_source_term.h
Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_validation\RHSs\calculate_StildeD1_source_term.h
Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_validation\RHSs\calculate_StildeD2_source_term.h
Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_library\RHSs\calculate_StildeD0_source_term.h
Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_library\RHSs\calculate_StildeD1_source_term.h
Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\GiRaFFE_Ccode_library\RHSs\calculate_StildeD2_source_term.h
Printing difference between original C code and this code...
Checking file calculate_StildeD0_source_term.h
No difference. TEST PASSED!
Checking file calculate_StildeD1_source_term.h
No difference. TEST PASSED!
Checking file calculate_StildeD2_source_term.h
No difference. TEST PASSED!

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

$\label{latex_pdf_output}$

The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$ -formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-GiRaFFE_NRPy_C_code_library-Source_Terms (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [4]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFE_NRPy-Source_Terms",location_of_template_file=os.path.join(".."))

Notebook output to PDF is only supported on Linux systems, with pdflatex installed.

GiRaFFE_NRPy: Source Terms¶