Notebook

Interpolating Metric Quantities on Cell Faces¶

Author: Patrick Nelson¶

Notebook Status: Validated

Validation Notes: This module will be self-validated against its module and will also be validated against the corresponding algorithm in the old GiRaFFE code in this tutorial.

This module presents the functionality of GiRaFFE_NRPy_Metric_Face_Values.py .¶

This notebook presents the macros from the original GiRaFFE that provide the values of the metric gridfunctions interpolated to the cell faces along with the code needed to implement this in NRPy.

Table of Contents¶

$\label{toc}$

This notebook is organized as follows

Step 0: Preliminaries
Step 1: The Interpolator
1. Step 1.a: Interpolator coefficients and definition
2. Step 1.b: Create an array to easily define the gridfunctions we want to interpolate.
3. Step 1.c: Define the loop parameters and body
Step 2: Code Validation against original C code
Step 3: Output this notebook to $\LaTeX$ -formatted PDF file

Step 0: Preliminaries [Back to top]¶

$\label{prelim}$

This first block of code just sets up a subdirectory within GiRaFFE_standalone_Ccodes/ to which we will write the C code and adds core NRPy+ functionality to sys.path.

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 # NRPy+: Core C code output module
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface
Ccodesdir = "GiRaFFE_standalone_Ccodes/FCVAL"
cmd.mkdir(os.path.join(Ccodesdir))

Step 1: The Interpolator [Back to top]¶

$\label{interpolator}$

Here, we we will write the code necessary to interpolate the metric gridfunction $\alpha, \beta^i, \gamma_{ij}$ onto the cell faces. These values will be necessary to compute fluxes of the Poynting vector and electric field through those faces.

Step 1.a: Interpolator coefficients and definition [Back to top]¶

$\label{macros}$

First, we will define the functional form of our interpolator. At some point on our grid $i$ , we will calculate the value of some gridfunction $Q$ at position $i-\tfrac{1}{2}$ with

$Q_{i-1/2} = a_{i-2} Q_{i-2} +a_{i-1} Q_{i-1} +a_{i} Q_{i} +a_{i+1} Q_{i+1}$ and the coefficients we will use for it,

$\begin{align} a_{i-2} &= -0.0625 \\ a_{i-1} &= 0.5625 \\ a_{i} &= 0.5625 \\ a_{i+1} &= -0.0625 \\ \end{align}$ .

In [2]:

%%writefile $Ccodesdir/interpolate_metric_gfs_to_cell_faces.h
// Side note: the following values could be used for cell averaged gfs:
//     am2=-1.0/12.0, am1=7.0/12.0, a0=7.0/12.0, a1=-1.0/12.0
// However, since the metric gfs store the grid point values instead of the cell average,
//     the following coefficients should be used:
//     am2 = -1/16, am1 = 9/16, a0 = 9/16, a1 = -1/16
// This will yield the third-order-accurate face values at m-1/2,
//      using values specified at {m-2,m-1,m,m+1}
#define AM2 -0.0625
#define AM1  0.5625
#define A0   0.5625
#define A1  -0.0625
#define COMPUTE_FCVAL(METRICm2,METRICm1,METRIC,METRICp1) (AM2*(METRICm2) + AM1*(METRICm1) + A0*(METRIC) + A1*(METRICp1))

Overwriting GiRaFFE_standalone_Ccodes/FCVAL/interpolate_metric_gfs_to_cell_faces.h

Step 1.b: Create an array to easily define the gridfunctions we want to interpolate. [Back to top]¶

$\label{gf_struct}$

We will need to apply this interpolation to each gridpoint for several gridfunctions: the lapse $\alpha$ , the shift $\beta^i$ , and the three-metric $\gamma_{ij}$ . Consider that in NRPy+, each gridfunction is assigned an integer identifier with the C macro #define. So, the simplest (and shortest to write!) way to ensure we hit each of these is to create arrays that list each of these identifiers in order, so we can always hit the write gridfunction no matter where each gridfunction lies in the list. We use two arrays; the first identifies the usual gridfunctions from which we will read, and the second identifies the face-sampled gridfunctions to which we will write.

In [3]:

%%writefile -a $Ccodesdir/interpolate_metric_gfs_to_cell_faces.h

const int metric_gfs_list[10] = {GAMMADD00GF,
                                 GAMMADD01GF,
                                 GAMMADD02GF,
                                 GAMMADD11GF,
                                 GAMMADD12GF,
                                 GAMMADD22GF,
                                 BETAU0GF,
                                 BETAU1GF,
                                 BETAU2GF,
                                 ALPHAGF};

const int metric_gfs_face_list[10] = {GAMMA_FACEDD00GF,
                                      GAMMA_FACEDD01GF,
                                      GAMMA_FACEDD02GF,
                                      GAMMA_FACEDD11GF,
                                      GAMMA_FACEDD12GF,
                                      GAMMA_FACEDD22GF,
                                      BETA_FACEU0GF,
                                      BETA_FACEU1GF,
                                      BETA_FACEU2GF,
                                      ALPHA_FACEGF};

const int num_metric_gfs = 10;

Appending to GiRaFFE_standalone_Ccodes/FCVAL/interpolate_metric_gfs_to_cell_faces.h

Step 1.c: Define the loop parameters and body [Back to top]¶

$\label{loops}$

Next, we will write the function that loops over the entire grid. One additional parameter to consider here is the direction in which we need to do the interpolation. This direction exactly corresponds to the parameter flux_dirn used in the calculation of the flux of the Poynting vector and electric field.

The outermost loop will iterate over the gridfunctions we listed above. Nested inside of that, there will be three loops that go through the grid in the usual way. However, the upper bound will be a little unusual. Instead of covering all points or all interior points, we will write these loops to cover all interior points and one extra past that. This is because we have defined our interpolation on the $i-\tfrac{1}{2}$ face of a cell, but any given calculation will require both that and an interpolation on the $i+\tfrac{1}{2}$ face as well.

In [4]:

# FIXME: The old oop bounds are NGHOSTS to Nxx+NGHOSTS+1, but converting reconstructed velocities to
#        and from Valencia requires extra points, so we crank it to maximum.
desc = "Interpolate metric gridfunctions to cell faces"
name = "interpolate_metric_gfs_to_cell_faces"
interp_Cfunc = outCfunction(
    outfile  = "returnstring", desc=desc, name=name,
    params   ="const paramstruct *params,REAL *auxevol_gfs,const int flux_dirn",
    preloop  ="""    int in_gf,out_gf;
    REAL Qm2,Qm1,Qp0,Qp1;

""" ,
    body     ="""    for(int gf = 0;gf < num_metric_gfs;gf++) {
        in_gf  = metric_gfs_list[gf];
        out_gf = metric_gfs_face_list[gf];
        for (int i2 = 2;i2 < Nxx_plus_2NGHOSTS2-1;i2++) {
            for (int i1 = 2;i1 < Nxx_plus_2NGHOSTS1-1;i1++) {
                for (int i0 = 2;i0 < Nxx_plus_2NGHOSTS0-1;i0++) {
                    Qm2 = auxevol_gfs[IDX4S(in_gf,i0-2*kronecker_delta[flux_dirn][0],i1-2*kronecker_delta[flux_dirn][1],i2-2*kronecker_delta[flux_dirn][2])];
                    Qm1 = auxevol_gfs[IDX4S(in_gf,i0-kronecker_delta[flux_dirn][0],i1-kronecker_delta[flux_dirn][1],i2-kronecker_delta[flux_dirn][2])];
                    Qp0 = auxevol_gfs[IDX4S(in_gf,i0,i1,i2)];
                    Qp1 = auxevol_gfs[IDX4S(in_gf,i0+kronecker_delta[flux_dirn][0],i1+kronecker_delta[flux_dirn][1],i2+kronecker_delta[flux_dirn][2])];
                    auxevol_gfs[IDX4S(out_gf,i0,i1,i2)] = COMPUTE_FCVAL(Qm2,Qm1,Qp0,Qp1);
                }
            }
        }
    }
""",
    rel_path_to_Cparams=os.path.join("../"))

with open(os.path.join(Ccodesdir,"interpolate_metric_gfs_to_cell_faces.h"),"a") as file:
    file.write(interp_Cfunc)

Step 2: Code Validation against `GiRaFFE_NRPy_FCVAL.py` [Back to top]¶

$\label{code_validation}$

Now, we will confirm that the code we have written here is the same as that generated by the module GiRaFFE_NRPy_FCVAL.py.

In [5]:

# Define the directory that we wish to validate against:
valdir = "GiRaFFE_NRPy/GiRaFFE_Ccode_library/FCVAL/"

import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL
FCVAL.GiRaFFE_NRPy_FCVAL(valdir)

import difflib
import sys

print("Printing difference between original C code and this code...")
# Open the files to compare
file = "interpolate_metric_gfs_to_cell_faces.h"

print("Checking file " + file)
with open(os.path.join(valdir,file)) as file1, open(os.path.join(Ccodesdir,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(Ccodesdir+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 interpolate_metric_gfs_to_cell_faces.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-FCVAL.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [6]:

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

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