Notebook

`GiRaFFE_NRPy`: Main Driver¶

Author: Patrick Nelson¶

Notebook Status: Validation in progress

Validation Notes: This code assembles the various parts needed for GRFFE evolution in order.

NRPy+ Source Code for this module:¶

GiRaFFE_NRPy/GiRaFFE_NRPy_Main_Driver.py

Other critical files (in alphabetical order):¶

GiRaFFE_NRPy/GiRaFFE_NRPy_Afield_flux_handwritten.py [tutorial] Generates the expressions to find the flux term of the induction equation.
GiRaFFE_NRPy/GiRaFFE_NRPy_A2B.py [tutorial] Generates the driver to compute the magnetic field from the vector potential/
GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py [tutorial] Generates the code to apply boundary conditions to the vector potential, scalar potential, and three-velocity.
GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py [tutorial] Generates the conservative-to-primitive and primitive-to-conservative solvers.
GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py [tutorial] Generates code to interpolate metric gridfunctions to cell faces.
GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py [tutorial] Genearates code to reconstruct primitive variables on cell faces.
GiRaFFE_NRPy/GiRaFFE_NRPy_Source_Terms.py [tutorial] Genearates code to compute the $\tilde{S}_i$ source term.
GiRaFFE_NRPy/Stilde_flux.py [tutorial] Generates the expressions to find the flux term of the Poynting flux evolution equation.
../GRFFE/equations.py [tutorial] Generates code necessary to compute the source terms.
../GRHD/equations.py [tutorial] Generates code necessary to compute the source terms.

Introduction:¶

Having written all the various algorithms that will go into evolving the GRFFE equations forward through time, we are ready to write a start-to-finish module to do so. However, to help keep things more organized, we will first create a dedicated module to assemble the various functions we need to run, in order, to perform the evolution. This will reduce the length of the standalone C code, improving that notebook's readability.

Table of Contents¶

$$\label{prelim}$$

During a given RK substep, we will perform the following steps in this order, based on the order used in the original GiRaFFE: 0. Step 0: Preliminaries

Step 1: Calculate the right-hand sides
1. Step 1.a: Calculate the portion of the gauge terms for $A_k$, $(\alpha \Phi - \beta^j A_j)$ and $\Phi$, $(\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi])$ inside the parentheses to be finite-differenced.
  1. GRFFE/equations.py, GRHD/equations.py
2. Step 1.b: Calculate the source terms of $\partial_t A_i$, $\partial_t \tilde{S}_i$, and $\partial_t [\sqrt{\gamma} \Phi]$ right-hand sides
  1. GRFFE/equations.py, GRHD/equations.py, GiRaFFE_NRPy/GiRaFFE_NRPy_Source_Terms
3. Step 1.c: Calculate the Flux terms
  1. In each direction:
  2. Interpolate the metric gridfunctions to cell faces
    1. GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py
  3. Reconstruct primitives $\bar{v}^i$ and $B^i$ on cell faces with the piecewise-parabolic method
    1. GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py
  4. Compute the fluxes of $\tilde{S}_i$ and $A_i$ and add the appropriate combinations to the evolution equation right-hand sides
    1. GiRaFFE_NRPy/Stilde_flux.py, GiRaFFE_NRPy/GiRaFFE_NRPy_Afield_flux_handwritten.py
Step 2: Recover the primitive variables and apply boundary conditions (post-step)
1. Step 2.a: Apply boundary conditions to $A_i$ and $\sqrt{\gamma} \Phi$
  1. GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py
2. Step 2.b: Compute $B^i$ from $A_i$
  1. GiRaFFE_NRPy/GiRaFFE_NRPy_A2B.py
3. Step 2.c: Run the Conservative-to-Primitive solver
  1. This applies fixes to $\tilde{S}_i$, then computes $\bar{v}^i$. A current sheet prescription is then applied to $\bar{v}^i$, and $\tilde{S}_i$ is recomputed to be consistent.
  2. GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py
4. Step 2.d: Apply outflow boundary conditions to $\bar{v}^i$
  1. GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py
Step 3: Write out the C code function
Step 3: Self-Validation against GiRaFFE_NRPy_Main_Drive.py
Step 5: Output this notebook to $\LaTeX$-formatted PDF file

Step 0: Preliminaries [Back to top]¶

$$\label{prelim}$$

We begin by importing the NRPy+ core functionality. We also import the GRHD module and the GRFFE module.

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)

from outputC import outCfunction, 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 cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface

thismodule = "GiRaFFE_NRPy_Main_Driver"

par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",2)

out_dir = os.path.join("GiRaFFE_standalone_Ccodes")
cmd.mkdir(out_dir)

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.

# Default Kreiss-Oliger dissipation strength
default_KO_strength = 0.1
diss_strength = par.Cparameters("REAL", thismodule, "diss_strength", default_KO_strength)

outCparams = "outCverbose=False,CSE_sorting=none"

Step 1: Calculate the right-hand sides [Back to top]¶

$$\label{rhs}$$

In the method of lines using Runge-Kutta methods, each timestep involves several "RK substeps" during which we will run the same set of function calls. These can be divided into two groups: one in which the RHSs themselves are calculated, and a second in which boundary conditions are applied and auxiliary variables updated (the post-step). Here, we focus on that first group.

Step 1.a: Calculate the portion of the gauge terms for $A_k$, $(\alpha \Phi - \beta^j A_j)$ and $\Phi$, $(\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi])$ inside the parentheses to be finite-differenced. [Back to top]¶

$$\label{operand}$$

The gauge terms of our evolution equations consist of two derivative terms: the Lorentz gauge term of $\partial_t A_k$, which is $\partial_k (\alpha \Phi - \beta^j A_j)$ and the non-damping, flux-like term of $\partial_t [\psi^6 \Phi]$, which is $\partial_j (\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi])$. We can save some effort and execution time (at the cost of memory needed) by computing the derivative operands, storing them, and then finite-differencing that stored variable. For more information, see the notebook for the implementation and the validation, as well as Tutorial-GRFFE_Equations-Cartesian and Tutorial-GRHD_Equations-Cartesian for the terms themselves.

In [2]:

import GRHD.equations as GRHD    # NRPy+: Generate general relativistic hydrodynamics equations
import GRFFE.equations as GRFFE  # NRPy+: Generate general relativisitic force-free electrodynamics equations

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")
AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
psi6Phi = gri.register_gridfunctions("EVOL","psi6Phi")
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD")

ixp.register_gridfunctions_for_single_rank1("AUXEVOL","PhievolParenU",DIM=3)
gri.register_gridfunctions("AUXEVOL","AevolParen")

GRHD.compute_sqrtgammaDET(gammaDD)
GRFFE.compute_AD_source_term_operand_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)
GRFFE.compute_psi6Phi_rhs_flux_term_operand(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)

parens_to_print = [\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","AevolParen"),rhs=GRFFE.AevolParen),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU0"),rhs=GRFFE.PhievolParenU[0]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU1"),rhs=GRFFE.PhievolParenU[1]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU2"),rhs=GRFFE.PhievolParenU[2]),\
                  ]

subdir = "RHSs"
cmd.mkdir(os.path.join(out_dir, subdir))
desc = "Calculate quantities to be finite-differenced for the GRFFE RHSs"
name = "calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs"
outCfunction(
    outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *restrict params,const REAL *restrict in_gfs,REAL *restrict auxevol_gfs",
    body     = fin.FD_outputC("returnstring",parens_to_print,params=outCparams),
    loopopts ="AllPoints",
    rel_path_to_Cparams=os.path.join("../"))

Output C function calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs() to file GiRaFFE_standalone_Ccodes/RHSs/calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs.h

Step 1.b: Calculate the source terms of $\partial_t A_i$, $\partial_t \tilde{S}_i$, and $\partial_t [\sqrt{\gamma} \Phi]$ right-hand sides [Back to top]¶

$$\label{source}$$

With the operands of the gradient of divergence operators stored in memory from the previous step, we can now calculate the terms on the RHS of $A_i$ and $[\sqrt{\gamma} \Phi]$ that involve the derivatives of those terms. We also compute the other term in the RHS of $[\sqrt{\gamma} \Phi]$, which is a straightforward damping term.

In [3]:

xi_damping = par.Cparameters("REAL",thismodule,"xi_damping",0.1)
GRFFE.compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)

AevolParen_dD = ixp.declarerank1("AevolParen_dD",DIM=3)
PhievolParenU_dD = ixp.declarerank2("PhievolParenU_dD","nosym",DIM=3)

A_rhsD = ixp.zerorank1()
psi6Phi_rhs = GRFFE.psi6Phi_damping

for i in range(3):
    A_rhsD[i] += -AevolParen_dD[i]
    psi6Phi_rhs += -PhievolParenU_dD[i][i]

# Add Kreiss-Oliger dissipation to the GRFFE RHSs:
# psi6Phi_dKOD = ixp.declarerank1("psi6Phi_dKOD")
# AD_dKOD    = ixp.declarerank2("AD_dKOD","nosym")
# for i in range(3):
#     psi6Phi_rhs += diss_strength*psi6Phi_dKOD[i]*rfm.ReU[i] # ReU[i] = 1/scalefactor_orthog_funcform[i]
#     for j in range(3):
#         A_rhsD[j] += diss_strength*AD_dKOD[j][i]*rfm.ReU[i] # ReU[i] = 1/scalefactor_orthog_funcform[i]

RHSs_to_print = [\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD0"),rhs=A_rhsD[0]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD1"),rhs=A_rhsD[1]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD2"),rhs=A_rhsD[2]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","psi6Phi"),rhs=psi6Phi_rhs),\
                ]

desc = "Calculate AD gauge term and psi6Phi RHSs"
name = "calculate_AD_gauge_psi6Phi_RHSs"
source_Ccode = outCfunction(
    outfile  = "returnstring", desc=desc, name=name,
    params   ="const paramstruct *params,const REAL *in_gfs,const REAL *auxevol_gfs,REAL *rhs_gfs",
    body     = fin.FD_outputC("returnstring",RHSs_to_print,params=outCparams),
    loopopts ="InteriorPoints",
    rel_path_to_Cparams=os.path.join("../")).replace("= NGHOSTS","= NGHOSTS_A2B").replace("NGHOSTS+Nxx0","Nxx_plus_2NGHOSTS0-NGHOSTS_A2B").replace("NGHOSTS+Nxx1","Nxx_plus_2NGHOSTS1-NGHOSTS_A2B").replace("NGHOSTS+Nxx2","Nxx_plus_2NGHOSTS2-NGHOSTS_A2B")
# Note the above .replace() functions. These serve to expand the loop range into the ghostzones, since
# the second-order FD needs fewer than some other algorithms we use do.
with open(os.path.join(out_dir,subdir,name+".h"),"w") as file:
    file.write(source_Ccode)

We also need to compute the source term of the $\tilde{S}_i$ evolution equation. This term involves derivatives of the four metric, so we can save some effort here by taking advantage of the interpolations done of the metric gridfunctions to the cell faces, which will allow us to take a finite-difference derivative with the accuracy of a higher order and the computational cost of a lower order. However, it will require some more complicated coding, detailed in Tutorial-GiRaFFE_NRPy-Source_Terms

In [4]:

import GiRaFFE_NRPy.GiRaFFE_NRPy_Source_Terms as source
# Declare this symbol:
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")

source.write_out_functions_for_StildeD_source_term(os.path.join(out_dir,subdir),outCparams,gammaDD,betaU,alpha,
                                                   ValenciavU,BU,sqrt4pi)

Output C function calculate_StildeD0_source_term() to file GiRaFFE_standalone_Ccodes/RHSs/calculate_StildeD0_source_term.h
Output C function calculate_StildeD1_source_term() to file GiRaFFE_standalone_Ccodes/RHSs/calculate_StildeD1_source_term.h
Output C function calculate_StildeD2_source_term() to file GiRaFFE_standalone_Ccodes/RHSs/calculate_StildeD2_source_term.h

Step 1.c: Calculate the Flux terms [Back to top]¶

$$\label{flux}$$

Now, we will compute the flux terms of $\partial_t A_i$ and $\partial_t \tilde{S}_i$. To do so, we will first need to interpolate the metric gridfunctions to cell faces and to reconstruct the primitives on the cell faces using the code detailed in Tutorial-GiRaFFE_NRPy-Metric_Face_Values and in Tutorial-GiRaFFE_NRPy-PPM.

In [5]:

subdir = "FCVAL"
cmd.mkdir(os.path.join(out_dir, subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL
FCVAL.GiRaFFE_NRPy_FCVAL(os.path.join(out_dir,subdir))

In [6]:

subdir = "PPM"
cmd.mkdir(os.path.join(out_dir, subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM
PPM.GiRaFFE_NRPy_PPM(os.path.join(out_dir,subdir))

Here, we will write the function to compute the electric field contribution to the induction equation RHS. This is coded with documentation in Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten. The recontructed values in the $i^{\rm th}$ direction will contribute to the $j^{\rm th}$ and $k^{\rm th}$ component of the electric field. That is, in Cartesian coordinates, the component $x$ of the electric field will be the average of the values computed on the cell faces in the $\pm y$- and $\pm z$-directions, and so forth for the other components. However, all of these can be written as only a single function as long as we appropriately pass cyclical permutations of the inputs.

In [7]:

import GiRaFFE_NRPy.GiRaFFE_NRPy_Afield_flux_handwritten as Af

# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")

# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)

subdir = "RHSs"
Af.GiRaFFE_NRPy_Afield_flux(os.path.join(out_dir,subdir))

We must do something similar here, albeit a bit simpler. For instance, the $x$ component of $\partial_t \tilde{S}_i$ will be a finite difference of the flux throught the faces in the $\pm x$ direction; for further detail, see Tutorial-GiRaFFE_NRPy-Stilde_flux.

In [8]:

ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")

import GiRaFFE_NRPy.Stilde_flux as Sf
Sf.generate_C_code_for_Stilde_flux(os.path.join(out_dir,subdir), True, alpha_face,gamma_faceDD,beta_faceU,
                                   Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi)

Output C function calculate_Stilde_rhsD() to file GiRaFFE_standalone_Ccodes/RHSs/calculate_Stilde_rhsD.h

Step 2: Recover the primitive variables and apply boundary conditions [Back to top]¶

$$\label{poststep}$$

With the RHSs computed, we can now recover the primitive variables, which are the Valencia three-velocity $\bar{v}^i$ and the magnetic field $B^i$. We can also apply boundary conditions to the vector potential and velocity. By doing this at each RK substep, we can help ensure the accuracy of the following substeps.

Step 2.a: Apply boundary conditions to $A_i$ and $\sqrt{\gamma} \Phi$ [Back to top]¶

$$\label{potential_bc}$$

First, we will apply boundary conditions to the vector potential, $A_i$, and the scalar potential $\sqrt{\gamma} \Phi$. The file we generate here contains both functions we need for BCs, as documented in Tutorial-GiRaFFE_NRPy-BCs.

In [9]:

subdir = "boundary_conditions"
cmd.mkdir(os.path.join(out_dir,subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_BCs as BC
BC.GiRaFFE_NRPy_BCs(os.path.join(out_dir,subdir))

Step 2.b: Compute $B^i$ from $A_i$ [Back to top]¶

$$\label{a2b}$$

Now, we will calculate the magnetic field as the curl of the vector potential at all points in our domain; this requires care to be taken in the ghost zones, which is detailed in Tutorial-GiRaFFE_NRPy-A2B.

In [10]:

subdir = "A2B"
cmd.mkdir(os.path.join(out_dir,subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_A2B as A2B
A2B.GiRaFFE_NRPy_A2B(os.path.join(out_dir,subdir),gammaDD,AD,BU)

Step 2.c: Run the Conservative-to-Primitive solver [Back to top]¶

$$\label{c2p}$$

With these functions, we apply fixes to the Poynting flux, and use that to update the three-velocity. Then, we apply our current sheet prescription to the velocity, and recompute the Poynting flux to agree with the now-fixed velocity. More detail can be found in Tutorial-GiRaFFE_NRPy-C2P_P2C.

In [11]:

import GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C as C2P_P2C
C2P_P2C.GiRaFFE_NRPy_C2P(StildeD,BU,gammaDD,betaU,alpha)

values_to_print = [
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=C2P_P2C.outStildeD[0]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=C2P_P2C.outStildeD[1]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=C2P_P2C.outStildeD[2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=C2P_P2C.ValenciavU[0]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=C2P_P2C.ValenciavU[1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=C2P_P2C.ValenciavU[2])
                  ]

subdir = "C2P"
cmd.mkdir(os.path.join(out_dir,subdir))
desc = "Apply fixes to \tilde{S}_i and recompute the velocity to match with current sheet prescription."
name = "GiRaFFE_NRPy_cons_to_prims"
outCfunction(
    outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *in_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params=outCparams),
    loopopts ="AllPoints,Read_xxs",
    rel_path_to_Cparams=os.path.join("../"))

Output C function GiRaFFE_NRPy_cons_to_prims() to file GiRaFFE_standalone_Ccodes/C2P/GiRaFFE_NRPy_cons_to_prims.h

In [12]:

# TINYDOUBLE = par.Cparameters("REAL",thismodule,"TINYDOUBLE",1e-100)

C2P_P2C.GiRaFFE_NRPy_P2C(gammaDD,betaU,alpha,  ValenciavU,BU, sqrt4pi)

values_to_print = [
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=C2P_P2C.StildeD[0]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=C2P_P2C.StildeD[1]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=C2P_P2C.StildeD[2]),
                  ]

desc = "Recompute StildeD after current sheet fix to Valencia 3-velocity to ensure consistency between conservative & primitive variables."
name = "GiRaFFE_NRPy_prims_to_cons"
outCfunction(
    outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *auxevol_gfs,REAL *in_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params=outCparams),
    loopopts ="AllPoints",
    rel_path_to_Cparams=os.path.join("../"))

Output C function GiRaFFE_NRPy_prims_to_cons() to file GiRaFFE_standalone_Ccodes/C2P/GiRaFFE_NRPy_prims_to_cons.h

Step 2.d: Apply outflow boundary conditions to $\bar{v}^i$ [Back to top]¶

$$\label{velocity_bc}$$

Now, we can apply outflow boundary conditions to the Valencia three-velocity. This specific type of boundary condition helps avoid numerical error "flowing" into our grid.

This function has already been generated above.

Step 3: Write out the C code function [Back to top]¶

$$\label{write_out}$$

Now, we have generated all the functions we will need for the GiRaFFE evolution. So, we will now assemble our evolution driver. This file will first #include all of the files we just generated for easy access. Then, we will write a function that calls these functions in the correct order, iterating over the flux directions as necessary.

In [13]:

%%writefile $out_dir/GiRaFFE_NRPy_Main_Driver.h
// Structure to track ghostzones for PPM:
typedef struct __gf_and_gz_struct__ {
  REAL *gf;
  int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;
// Some additional constants needed for PPM:
const int VX=0,VY=1,VZ=2,BX=3,BY=4,BZ=5;
const int NUM_RECONSTRUCT_GFS = 6;

// Include ALL functions needed for evolution
#include "RHSs/calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs.h"
#include "RHSs/calculate_AD_gauge_psi6Phi_RHSs.h"
#include "PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
#include "FCVAL/interpolate_metric_gfs_to_cell_faces.h"
#include "RHSs/calculate_StildeD0_source_term.h"
#include "RHSs/calculate_StildeD1_source_term.h"
#include "RHSs/calculate_StildeD2_source_term.h"
#include "RHSs/calculate_E_field_flat_all_in_one.h"
#include "RHSs/calculate_Stilde_flux_D0.h"
#include "RHSs/calculate_Stilde_flux_D1.h"
#include "RHSs/calculate_Stilde_flux_D2.h"
#include "RHSs/calculate_Stilde_rhsD.h"
#include "boundary_conditions/GiRaFFE_boundary_conditions.h"
#include "A2B/driver_AtoB.h"
#include "C2P/GiRaFFE_NRPy_cons_to_prims.h"
#include "C2P/GiRaFFE_NRPy_prims_to_cons.h"

void GiRaFFE_NRPy_RHSs(const paramstruct *restrict params,REAL *restrict auxevol_gfs,const REAL *restrict in_gfs,REAL *restrict rhs_gfs) {
#include "set_Cparameters.h"
    // First thing's first: initialize the RHSs to zero!
#pragma omp parallel for
    for(int ii=0;ii<Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*NUM_EVOL_GFS;ii++) {
        rhs_gfs[ii] = 0.0;
    }
    // Next calculate the easier source terms that don't require flux directions
    // This will also reset the RHSs for each gf at each new timestep.
    calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs(params,in_gfs,auxevol_gfs);
    calculate_AD_gauge_psi6Phi_RHSs(params,in_gfs,auxevol_gfs,rhs_gfs);

    // Now, we set up a bunch of structs of pointers to properly guide the PPM algorithm.
    // They also count the number of ghostzones available.
    gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];
    int which_prims_to_reconstruct[NUM_RECONSTRUCT_GFS],num_prims_to_reconstruct;
    const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;

    REAL *temporary = auxevol_gfs + Nxxp2NG012*AEVOLPARENGF; //We're not using this anymore
    // This sets pointers to the portion of auxevol_gfs containing the relevant gridfunction.
    int ww=0;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU0GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF;
    ww++;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU1GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF;
    ww++;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU2GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF;
    ww++;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU0GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU0GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU0GF;
    ww++;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU1GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU1GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU1GF;
    ww++;
    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU2GF;
      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU2GF;
      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU2GF;
    ww++;

    // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:
    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }
    // Left/right variables are not yet defined, yet we set the # of gz's to zero by default:
    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }

    ww=0;
    which_prims_to_reconstruct[ww]=VX; ww++;
    which_prims_to_reconstruct[ww]=VY; ww++;
    which_prims_to_reconstruct[ww]=VZ; ww++;
    which_prims_to_reconstruct[ww]=BX; ww++;
    which_prims_to_reconstruct[ww]=BY; ww++;
    which_prims_to_reconstruct[ww]=BZ; ww++;
    num_prims_to_reconstruct=ww;

    // In each direction, perform the PPM reconstruction procedure.
    // Then, add the fluxes to the RHS as appropriate.
    for(int flux_dirn=0;flux_dirn<3;flux_dirn++) {
        // In each direction, interpolate the metric gfs (gamma,beta,alpha) to cell faces.
        interpolate_metric_gfs_to_cell_faces(params,auxevol_gfs,flux_dirn+1);
        // Then, reconstruct the primitive variables on the cell faces.
        // This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
        reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
                                                which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
        // For example, if flux_dirn==0, then at gamma_faceDD00(i,j,k) represents gamma_{xx}
        // at (i-1/2,j,k), Valenciav_lU0(i,j,k) is the x-component of the velocity at (i-1/2-epsilon,j,k),
        // and Valenciav_rU0(i,j,k) is the x-component of the velocity at (i-1/2+epsilon,j,k).

        if(flux_dirn==0) {
            // Next, we calculate the source term for StildeD. Again, this also resets the rhs_gfs array at
            // each new timestep.
            calculate_StildeD0_source_term(params,auxevol_gfs,rhs_gfs);
            // Now, compute the electric field on each face of a cell and add it to the RHSs as appropriate
            //calculate_E_field_D0_right(params,auxevol_gfs,rhs_gfs);
            //calculate_E_field_D0_left(params,auxevol_gfs,rhs_gfs);
            // Finally, we calculate the flux of StildeD and add the appropriate finite-differences
            // to the RHSs.
            calculate_Stilde_flux_D0(params,auxevol_gfs,rhs_gfs);
        }
        else if(flux_dirn==1) {
            calculate_StildeD1_source_term(params,auxevol_gfs,rhs_gfs);
            //calculate_E_field_D1_right(params,auxevol_gfs,rhs_gfs);
            //calculate_E_field_D1_left(params,auxevol_gfs,rhs_gfs);
            calculate_Stilde_flux_D1(params,auxevol_gfs,rhs_gfs);
        }
        else {
            calculate_StildeD2_source_term(params,auxevol_gfs,rhs_gfs);
            //calculate_E_field_D2_right(params,auxevol_gfs,rhs_gfs);
            //calculate_E_field_D2_left(params,auxevol_gfs,rhs_gfs);
            calculate_Stilde_flux_D2(params,auxevol_gfs,rhs_gfs);
        }
        calculate_Stilde_rhsD(flux_dirn+1,params,auxevol_gfs,rhs_gfs);
        for(int count=0;count<=1;count++) {
            // This function is written to be general, using notation that matches the forward permutation added to AD2,
            // i.e., [F_HLL^x(B^y)]_z corresponding to flux_dirn=0, count=1.
            // The SIGN parameter is necessary because
            // -E_z(x_i,y_j,z_k) = 0.25 ( [F_HLL^x(B^y)]_z(i+1/2,j,k)+[F_HLL^x(B^y)]_z(i-1/2,j,k)
            //                           -[F_HLL^y(B^x)]_z(i,j+1/2,k)-[F_HLL^y(B^x)]_z(i,j-1/2,k) )
            // Note the negative signs on the reversed permutation terms!

            // By cyclically permuting with flux_dirn, we
            // get contributions to the other components, and by incrementing count, we get the backward permutations:
            // Let's suppose flux_dirn = 0. Then we will need to update Ay (count=0) and Az (count=1):
            //     flux_dirn=count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+0)%3=AD1GF <- Updating Ay!
            //        (flux_dirn)%3 = (0)%3 = 0               Vx
            //        (flux_dirn-count+2)%3 = (0-0+2)%3 = 2   Vz .  Inputs Vx, Vz -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
            //     flux_dirn=0,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+1)%3=AD2GF <- Updating Az!
            //        (flux_dirn)%3 = (0)%3 = 0               Vx
            //        (flux_dirn-count+2)%3 = (0-1+2)%3 = 1   Vy .  Inputs Vx, Vy -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
            // Let's suppose flux_dirn = 1. Then we will need to update Az (count=0) and Ax (count=1):
            //     flux_dirn=1,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+0)%3=AD2GF <- Updating Az!
            //        (flux_dirn)%3 = (1)%3 = 1               Vy
            //        (flux_dirn-count+2)%3 = (1-0+2)%3 = 0   Vx .  Inputs Vy, Vx -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
            //     flux_dirn=count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+1)%3=AD0GF <- Updating Ax!
            //        (flux_dirn)%3 = (1)%3 = 1               Vy
            //        (flux_dirn-count+2)%3 = (1-1+2)%3 = 2   Vz .  Inputs Vy, Vz -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
            // Let's suppose flux_dirn = 2. Then we will need to update Ax (count=0) and Ay (count=1):
            //     flux_dirn=2,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+0)%3=AD0GF <- Updating Ax!
            //        (flux_dirn)%3 = (2)%3 = 2               Vz
            //        (flux_dirn-count+2)%3 = (2-0+2)%3 = 1   Vy .  Inputs Vz, Vy -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
            //     flux_dirn=2,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+1)%3=AD1GF <- Updating Ay!
            //        (flux_dirn)%3 = (2)%3 = 2               Vz
            //        (flux_dirn-count+2)%3 = (2-1+2)%3 = 0   Vx .  Inputs Vz, Vx -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
            calculate_E_field_flat_all_in_one(params,
              &auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn-count+2)%3, 0)],
              &auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn-count+2)%3, 0)],
              &auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn-count+2)%3, 0)],
              &auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn-count+2)%3, 0)],
              &auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn-count+2)%3, 0)],
              &auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn-count+2)%3, 0)],
              &rhs_gfs[IDX4ptS(AD0GF+(flux_dirn+1+count)%3,0)], 2.0*((REAL)count)-1.0, flux_dirn);
        }
    }
}

void GiRaFFE_NRPy_post_step(const paramstruct *restrict params,REAL *xx[3],REAL *restrict auxevol_gfs,REAL *restrict evol_gfs,const int n) {
    // First, apply BCs to AD and psi6Phi. Then calculate BU from AD
    apply_bcs_potential(params,evol_gfs);
    driver_A_to_B(params,evol_gfs,auxevol_gfs);
    //override_BU_with_old_GiRaFFE(params,auxevol_gfs,n);
    // Apply fixes to StildeD, then recompute the velocity at the new timestep.
    // Apply the current sheet prescription to the velocities
    GiRaFFE_NRPy_cons_to_prims(params,xx,auxevol_gfs,evol_gfs);
    // Then, recompute StildeD to be consistent with the new velocities
    //GiRaFFE_NRPy_prims_to_cons(params,auxevol_gfs,evol_gfs);
    // Finally, apply outflow boundary conditions to the velocities.
    apply_bcs_velocity(params,auxevol_gfs);
}

Overwriting GiRaFFE_standalone_Ccodes/GiRaFFE_NRPy_Main_Driver.h

Step 4: Self-Validation against `GiRaFFE_NRPy_Main_Drive.py` [Back to top]¶

$$\label{code_validation}$$

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

that were generated in this tutorial and
those that are generated in the module GiRaFFE_NRPy_Main_Driver.py

In [14]:

gri.glb_gridfcs_list = []
# Define the directory that we wish to validate against:
valdir = os.path.join("GiRaFFE_validation_Ccodes")
cmd.mkdir(valdir)

import GiRaFFE_NRPy.GiRaFFE_NRPy_Main_Driver as md
md.GiRaFFE_NRPy_Main_Driver_generate_all(valdir)

Output C function calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs() to file GiRaFFE_validation_Ccodes/RHSs/calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs.h
Output C function calculate_StildeD0_source_term() to file GiRaFFE_validation_Ccodes/RHSs/calculate_StildeD0_source_term.h
Output C function calculate_StildeD1_source_term() to file GiRaFFE_validation_Ccodes/RHSs/calculate_StildeD1_source_term.h
Output C function calculate_StildeD2_source_term() to file GiRaFFE_validation_Ccodes/RHSs/calculate_StildeD2_source_term.h
Output C function calculate_Stilde_rhsD() to file GiRaFFE_validation_Ccodes/RHSs/calculate_Stilde_rhsD.h
Output C function GiRaFFE_NRPy_cons_to_prims() to file GiRaFFE_validation_Ccodes/C2P/GiRaFFE_NRPy_cons_to_prims.h
Output C function GiRaFFE_NRPy_prims_to_cons() to file GiRaFFE_validation_Ccodes/C2P/GiRaFFE_NRPy_prims_to_cons.h

With both sets of codes generated, we can now compare them against each other.

In [15]:

import difflib
import sys

print("Printing difference between original C code and this code...")
# Open the files to compare
files = ["GiRaFFE_NRPy_Main_Driver.h",
         "RHSs/calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs.h",
         "RHSs/calculate_AD_gauge_psi6Phi_RHSs.h",
         "PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c",
         "PPM/loop_defines_reconstruction_NRPy.h",
         "FCVAL/interpolate_metric_gfs_to_cell_faces.h",
         "RHSs/calculate_StildeD0_source_term.h",
         "RHSs/calculate_StildeD1_source_term.h",
         "RHSs/calculate_StildeD2_source_term.h",
         "RHSs/calculate_E_field_flat_all_in_one.h",
         "RHSs/calculate_Stilde_flux_D0.h",
         "RHSs/calculate_Stilde_flux_D1.h",
         "RHSs/calculate_Stilde_flux_D2.h",
         "boundary_conditions/GiRaFFE_boundary_conditions.h",
         "A2B/driver_AtoB.h",
         "C2P/GiRaFFE_NRPy_cons_to_prims.h",
         "C2P/GiRaFFE_NRPy_prims_to_cons.h"]

for file in files:
    print("Checking file " + file)
    with open(os.path.join(valdir,file)) as file1, open(os.path.join(out_dir,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(out_dir,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)

Printing difference between original C code and this code...
Checking file GiRaFFE_NRPy_Main_Driver.h
No difference. TEST PASSED!
Checking file RHSs/calculate_AD_gauge_term_psi6Phi_flux_term_for_RHSs.h
No difference. TEST PASSED!
Checking file RHSs/calculate_AD_gauge_psi6Phi_RHSs.h
No difference. TEST PASSED!
Checking file PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c
No difference. TEST PASSED!
Checking file PPM/loop_defines_reconstruction_NRPy.h
No difference. TEST PASSED!
Checking file FCVAL/interpolate_metric_gfs_to_cell_faces.h
No difference. TEST PASSED!
Checking file RHSs/calculate_StildeD0_source_term.h
No difference. TEST PASSED!
Checking file RHSs/calculate_StildeD1_source_term.h
No difference. TEST PASSED!
Checking file RHSs/calculate_StildeD2_source_term.h
No difference. TEST PASSED!
Checking file RHSs/calculate_E_field_flat_all_in_one.h
No difference. TEST PASSED!
Checking file RHSs/calculate_Stilde_flux_D0.h
No difference. TEST PASSED!
Checking file RHSs/calculate_Stilde_flux_D1.h
No difference. TEST PASSED!
Checking file RHSs/calculate_Stilde_flux_D2.h
No difference. TEST PASSED!
Checking file boundary_conditions/GiRaFFE_boundary_conditions.h
No difference. TEST PASSED!
Checking file A2B/driver_AtoB.h
No difference. TEST PASSED!
Checking file C2P/GiRaFFE_NRPy_cons_to_prims.h
No difference. TEST PASSED!
Checking file C2P/GiRaFFE_NRPy_prims_to_cons.h
No difference. TEST PASSED!

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

In [16]:

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

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

GiRaFFE_NRPy: Main Driver¶