Notebook

Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: Piecewise-Parabolic Method¶

Author: Patrick Nelson¶

Edits by Terrence Pierre Jacques¶

This module Validates the PPM routine for `GiRaFFE`.¶

Notebook Status: Validated

Validation Notes: This module validates the routines in Tutorial-GiRaFFE_NRPy-PPM.

NRPy+ Source Code for this module:¶

GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py [tutorial] Reconstruct primitive variables at cell interfaces

Introduction:¶

This notebook validates the PPM routine for use in GiRaFFE_NRPy. This code was adapted from the corresonding algorithm used in the original GiRaFFE to work with the NRPy+ infrastructure. So, we will test it against that code by considering an analytic form of the three-velocity and magnetic fields over a small grid. We will also repeat the test with an added sharp feature, which is an important test for this algorithm since it is specifically designed to handle them.

When this notebook is run, the significant digits of agreement between the old GiRaFFE and new GiRaFFE_NRPy versions of the algorithm will be printed to the screen right after the code is run here.

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows:

Step 1: Set up core functions and parameters for unit testing the PPM algorithm
1. Step 1.a Valencia three-velocity
2. Step 1.b Magnetic field
3. Step 1.c Generate C functions to write the test data
4. Step 1.d Set free parameters in the code
5. Step 1.e Download GiRaFFE files from Bitbucket
6. Step 1.f Generate GiRaFFE_NRPy Files
Step 2: PPM_unit_test.C: The Main C Code
1. Step 2.a: Compile and run the code
Step 3: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Set up core functions and parameters for unit testing the PPM algorithm [Back to top]¶

$$\label{setup}$$

We'll start by appending the relevant paths to sys.path so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. We must also set the desired finite differencing order.

In [1]:

import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions
# First, we'll add the parent directory to the list of directories Python will check for modules.
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)
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, nrpyAbs # NRPy+: Core C code output module
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import grid as gri               # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm   # NRPy+: Reference metric support
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface

Ccodesdir = "Start-to-Finish-UnitTests/PPM_UnitTest/"

# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)

outdir = os.path.join(Ccodesdir,"output/")
cmd.mkdir(outdir)

thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-PPM"

Use_Shock_Data = True

Step 1.a: Valencia three-velocity [Back to top]¶

$$\label{velocity}$$

Here, we'll generate some functions for the velocity. Let's choose arctangents, since those have asymptotes that can be easily manipulated to prevent accidentally setting superluminal speeds. \begin{align} \bar{v}^x &= \frac{2}{\pi} \arctan(ax + by + cz) \\ \bar{v}^y &= \frac{2}{\pi} \arctan(bx + cy + az) \\ \bar{v}^z &= \frac{2}{\pi} \arctan(cx + ay + bz) \\ \end{align} If we want to add a jump at the origin, we can simply add $\max(0,x)$ to the argument of the arctangent. This will add a shock in the $x$-direction. The maximum will be described without the use of if statements as $$ \max(a,b) = \tfrac{1}{2} \left( a+b + \lvert a-b \rvert \right). $$

In [2]:

def max_noif(a,b):
    return sp.Rational(1,2)*(a+b+nrpyAbs(a-b))

a,b,c = par.Cparameters("REAL",thismodule,["a","b","c"],1e300) # Note that this default value allows us to set
                                                               # these directly in the C code
M_PI  = par.Cparameters("#define",thismodule,["M_PI"], "")

par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
x = rfm.xx_to_Cart[0]
y = rfm.xx_to_Cart[1]
z = rfm.xx_to_Cart[2]

args = ixp.zerorank1()
args[0] = a*x + b*y + c*z
args[1] = b*x + c*y + a*z
args[2] = c*x + a*y + b*z
if Use_Shock_Data:
    for i in range(3):
        args[i] += max_noif(0,5*x)

ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
for i in range(3):
    ValenciavU[i] = (sp.sympify(2)/M_PI)*sp.atan(args[i])

Step 1.b: Magnetic field [Back to top]¶

$$\label{magnetic}$$

We'll also need some functions for the magnetic field. Exponentials sound fun. \begin{align} B^x &= \exp(ey+fz) \\ B^y &= \exp(fz+dx) \\ B^z &= \exp(dx+ey) \\ \end{align} In this case, we'll add $\max{(0,x)}$ to the field to add the jump.

In [3]:

d,e,f = par.Cparameters("REAL",thismodule,["d","e","f"],1e300) # Note that this default value allows us to set
                                                               # these directly in the C code
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
BU[0] = sp.exp(e*y+f*z)
BU[1] = sp.exp(f*z+d*x)
BU[2] = sp.exp(d*x+e*y)
if Use_Shock_Data:
    for i in range(3):
        BU[i] += max_noif(0,5*x)

Step 1.c: Generate C functions to write the test data [Back to top]¶

$$\label{functions}$$

In [4]:

BU_to_print = [\
                lhrh(lhs=gri.gfaccess("out_gfs","BU0"),rhs=BU[0]),\
                lhrh(lhs=gri.gfaccess("out_gfs","BU1"),rhs=BU[1]),\
                lhrh(lhs=gri.gfaccess("out_gfs","BU2"),rhs=BU[2]),\
               ]

desc = "Calculate sample magnetic field data"
name = "calculate_BU"
outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
    body     = fin.FD_outputC("returnstring",BU_to_print,params="outCverbose=False"),
    loopopts="AllPoints,Read_xxs")

ValenciavU_to_print = [\
                       lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU0"),rhs=ValenciavU[0]),\
                       lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU1"),rhs=ValenciavU[1]),\
                       lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU2"),rhs=ValenciavU[2]),\
                      ]

desc = "Calculate sample velocity data"
name = "calculate_ValenciavU"
outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
    body     = fin.FD_outputC("returnstring",ValenciavU_to_print,params="outCverbose=False"),
    loopopts="AllPoints,Read_xxs")

Output C function calculate_BU() to file Start-to-Finish-UnitTests/PPM_UnitTest/calculate_BU.h
Output C function calculate_ValenciavU() to file Start-to-Finish-UnitTests/PPM_UnitTest/calculate_ValenciavU.h

Step 1.d: Set free parameters in the code [Back to top]¶

$$\label{free_parameters}$$

We also need to create the files that interact with NRPy's C parameter interface.

In [5]:

# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h

# Step 3.d.ii: Set free_parameters.h
with open(os.path.join(Ccodesdir,"free_parameters.h"),"w") as file:
    file.write("""
// Override parameter defaults with values based on command line arguments and NGHOSTS.
// We'll use this grid. It has one point and one ghost zone.
const int NGHOSTS = 3;
params.Nxx0 = 1;
params.Nxx1 = 1;
params.Nxx2 = 1;
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;
// Step 0d: Set up space and time coordinates
// Step 0d.i: Declare \Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:
const REAL xxmin[3] = {-1.0,-1.0,-1.0};
const REAL xxmax[3] = { 1.0, 1.0, 1.0};

params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx_plus_2NGHOSTS0-1.0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx_plus_2NGHOSTS1-1.0);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx_plus_2NGHOSTS2-1.0);
//printf("dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\n",params.dxx0,params.dxx1,params.dxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;
\n""")

# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))

Step 1.e: Download `GiRaFFE` files from Bitbucket [Back to top]¶

$$\label{download}$$

Here, we download the relevant portion of the original GiRaFFE code from Bitbucket, which is all contained in the files loop_defines_reconstruction.h and reconstruct_set_of_prims_PPM_GRFFE.C. Normally, these files depend on things declared in GiRaFFE_headers.h, but we've already defined what we need from that file as part of getting our new version working. Some of the numbers like VX=0, etc.,are slightly different, only because we are testing the limited subset of the original's functionality that we need for the new version.

In [6]:

# First download the original IllinoisGRMHD source code
import urllib
original_file_url  = ["https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/loop_defines_reconstruction.h",\
                      "https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/reconstruct_set_of_prims_PPM_GRFFE.C",\
                     ]
original_file_name = ["loop_defines_reconstruction.h",\
                      "reconstruct_set_of_prims_PPM_GRFFE.C",\
                     ]

for i in range(len(original_file_url)):
    original_file_path = os.path.join(Ccodesdir,original_file_name[i])

    # Then download the original IllinoisGRMHD source code
    # We try it here in a couple of ways in an attempt to keep
    # the code more portable
    try:
        original_file_code = urllib.request.urlopen(original_file_url[i]).read().decode('utf-8')
    except:
        original_file_code = urllib.urlopen(original_file_url[i]).read().decode('utf-8')

    # Write down the file the original IllinoisGRMHD source code
    with open(original_file_path,"w") as file:
        file.write(original_file_code)

Step 1.f: Generate `GiRaFFE_NRPy` Files [Back to top]¶

$$\label{module}$$

Here, we generate the functions we want to test by calling the function found here and documented in this tutorial.

In [7]:

import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM
PPM.GiRaFFE_NRPy_PPM(os.path.join(Ccodesdir,"PPM"))

Step 2: `PPM_unit_test.C`: The Main C Code [Back to top]¶

$$\label{mainc}$$

Now that we have our three-velocity and magnetic field to set up data, we will start writing our unit test. We'll also import common C functionality, define REAL, the number of ghost zones, and the faces, and set the standard macros for NRPy+ style memory access. We then compare the output between NRPy+ and ETK using the function reconstruct_set_of_prims_PPM_GRFFE_NRPy() and reconstruct_set_of_prims_PPM_GRFFE().

In [8]:

%%writefile $Ccodesdir/PPM_unit_test.C
// These are common packages that we are likely to need.
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h" // Needed for strncmp, etc.
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#include <time.h>   // Needed to set a random seed.

#define REAL double
#include "declare_Cparameters_struct.h"

REAL a,b,c,d,e,f;

// Standard NRPy+ memory access:
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * (k) ) )
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )

#define VALENCIAVU0GF 0
#define VALENCIAVU1GF 1
#define VALENCIAVU2GF 2
#define BU0GF 3
#define BU1GF 4
#define BU2GF 5
#define VALENCIAV_RU0GF 6
#define VALENCIAV_RU1GF 7
#define VALENCIAV_RU2GF 8
#define B_RU0GF 9
#define B_RU1GF 10
#define B_RU2GF 11
#define VALENCIAV_LU0GF 12
#define VALENCIAV_LU1GF 13
#define VALENCIAV_LU2GF 14
#define B_LU0GF 15
#define B_LU1GF 16
#define B_LU2GF 17
#define NUM_AUXEVOL_GFS 18

// Some specific definitions needed for this file
typedef struct __gf_and_gz_struct__ {
  REAL *gf;
  int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;

const int VX=0,VY=1,VZ=2,BX=3,BY=4,BZ=5;
const int NUM_RECONSTRUCT_GFS = 6;
const int MAXNUMVARS = NUM_RECONSTRUCT_GFS; // For the CCTK version

#include "PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
#include "PPM/loop_defines_reconstruction_NRPy.h"

#include "calculate_BU.h"
#include "calculate_ValenciavU.h"

// Some needed workarounds to get the ETK version of the code to work
#define CCTK_REAL double
#define DECLARE_CCTK_PARAMETERS //'
#define CCTK_GFINDEX3D(thing,i,j,k) ( (i) + cctk_lsh[0] * ( (j) + cctk_lsh[1] * (k) ) )
#define VERR_DEF_PARAMS 0
#define CCTK_VError(dummy,string,a,b,c,d,e,f,g,h) printf(string,a,b,c,d,e,f,g,h)

struct cGH{};
const cGH *cctkGH;

#include "reconstruct_set_of_prims_PPM_GRFFE.C"
#include "loop_defines_reconstruction.h"

void compare_ppm_calc(int flux_dirn){
      paramstruct params;
      #include "set_Cparameters_default.h"

          // Step 0c: Set free parameters, overwriting Cparameters defaults
          //          by hand or with command-line input, as desired.
      #include "free_parameters.h"
      #include "set_Cparameters-nopointer.h"

    const int cctk_lsh[3] = {Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2};

    // We'll define our grid slightly different from how we normally would. We let our outermost
    // ghostzones coincide with xxmin and xxmax instead of the interior of the grid. This means
    // that the ghostzone points will have identical positions so we can do convergence tests of them.    // Step 0d.ii: Set up uniform coordinate grids
    REAL *xx[3];
    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);
    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);
    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);
    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + ((REAL)(j))*dxx0;
    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + ((REAL)(j))*dxx1;
    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + ((REAL)(j))*dxx2;

    //for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) printf("xx[0][%d] = %.15e\n",i,xx[0][i]);

    // This is the array to which we'll write the NRPy+ variables.
    REAL *auxevol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);

    a = (double)(rand()%10-5);
    b = (double)(rand()%10-5);
    c = (double)(rand()%10-5);
    d = (double)(rand()%10-5);
    e = (double)(rand()%10-5);
    f = (double)(rand()%10-5);

    calculate_BU(&params,xx,auxevol_gfs);
    calculate_ValenciavU(&params,xx,auxevol_gfs);

    gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];

    const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
    REAL temporary[Nxxp2NG012];

    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; }

    int which_prims_to_reconstruct[NUM_RECONSTRUCT_GFS] = {VX, VY, VZ, BX, BY, BZ};
    int num_prims_to_reconstruct = NUM_RECONSTRUCT_GFS;

    // Outputs to store data temporarily for comparison:
    REAL vxrC,vxrN,vyrC,vyrN,vzrC,vzrN;
    REAL vxlC,vxlN,vylC,vylN,vzlC,vzlN;
    REAL BxrC,BxrN,ByrC,ByrN,BzrC,BzrN;
    REAL BxlC,BxlN,BylC,BylN,BzlC,BzlN;
    // There is one relevant gridpoint at which to compare: the center, or indices 3,3,3.
    int all_agree = 1;

    // 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, num_prims_to_reconstruct,
                                            which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
    vxrN = auxevol_gfs[IDX4S(VALENCIAV_RU0GF,3,3,3)];
    vyrN = auxevol_gfs[IDX4S(VALENCIAV_RU1GF,3,3,3)];
    vzrN = auxevol_gfs[IDX4S(VALENCIAV_RU2GF,3,3,3)];
    vxlN = auxevol_gfs[IDX4S(VALENCIAV_LU0GF,3,3,3)];
    vylN = auxevol_gfs[IDX4S(VALENCIAV_LU1GF,3,3,3)];
    vzlN = auxevol_gfs[IDX4S(VALENCIAV_LU2GF,3,3,3)];
    BxrN = auxevol_gfs[IDX4S(B_RU0GF,3,3,3)];
    ByrN = auxevol_gfs[IDX4S(B_RU1GF,3,3,3)];
    BzrN = auxevol_gfs[IDX4S(B_RU2GF,3,3,3)];
    BxlN = auxevol_gfs[IDX4S(B_LU0GF,3,3,3)];
    BylN = auxevol_gfs[IDX4S(B_LU1GF,3,3,3)];
    BzlN = auxevol_gfs[IDX4S(B_LU2GF,3,3,3)];
    reconstruct_set_of_prims_PPM_GRFFE(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                                       in_prims,out_prims_r,out_prims_l,temporary);
    vxrC = auxevol_gfs[IDX4S(VALENCIAV_RU0GF,3,3,3)];
    vyrC = auxevol_gfs[IDX4S(VALENCIAV_RU1GF,3,3,3)];
    vzrC = auxevol_gfs[IDX4S(VALENCIAV_RU2GF,3,3,3)];
    vxlC = auxevol_gfs[IDX4S(VALENCIAV_LU0GF,3,3,3)];
    vylC = auxevol_gfs[IDX4S(VALENCIAV_LU1GF,3,3,3)];
    vzlC = auxevol_gfs[IDX4S(VALENCIAV_LU2GF,3,3,3)];
    BxrC = auxevol_gfs[IDX4S(B_RU0GF,3,3,3)];
    ByrC = auxevol_gfs[IDX4S(B_RU1GF,3,3,3)];
    BzrC = auxevol_gfs[IDX4S(B_RU2GF,3,3,3)];
    BxlC = auxevol_gfs[IDX4S(B_LU0GF,3,3,3)];
    BylC = auxevol_gfs[IDX4S(B_LU1GF,3,3,3)];
    BzlC = auxevol_gfs[IDX4S(B_LU2GF,3,3,3)];

    printf("Checking results in the x-direction...\n");
    if(vxrN!=vxrC) {
        printf("vxr agrees to %.3e significant digits.\n",1-log10(2*fabs(vxrN-vxrC)/(fabs(vxrN)+fabs(vxrC))));
        all_agree=0;
    }
    if(vyrN!=vyrC) {
        printf("vyr agrees to %.3e significant digits.\n",1-log10(2*fabs(vyrN-vyrC)/(fabs(vyrN)+fabs(vyrC))));
        all_agree=0;
    }
    if(vzrN!=vzrC) {
        printf("vzr agrees to %.3e significant digits.\n",1-log10(2*fabs(vzrN-vzrC)/(fabs(vzrN)+fabs(vzrC))));
        all_agree=0;
    }
    if(vxlN!=vxlC) {
        printf("vxl agrees to %.3e significant digits.\n",1-log10(2*fabs(vxlN-vxlC)/(fabs(vxlN)+fabs(vxlC))));
        all_agree=0;
    }
    if(vylN!=vylC) {
        printf("vyl agrees to %.3e significant digits.\n",1-log10(2*fabs(vylN-vylC)/(fabs(vylN)+fabs(vylC))));
        all_agree=0;
    }
    if(vzlN!=vzlC) {
        printf("vzl agrees to %.3e significant digits.\n",1-log10(2*fabs(vzlN-vzlC)/(fabs(vzlN)+fabs(vzlC))));
        all_agree=0;
    }
    if(BxrN!=BxrC) {
        printf("Bxr agrees to %.3e significant digits.\n",1-log10(2*fabs(BxrN-BxrC)/(fabs(BxrN)+fabs(BxrC))));
        all_agree=0;
    }
    if(ByrN!=ByrC) {
        printf("Byr agrees to %.3e significant digits.\n",1-log10(2*fabs(ByrN-ByrC)/(fabs(ByrN)+fabs(ByrC))));
        all_agree=0;
    }
    if(BzrN!=BzrC) {
        printf("Bzr agrees to %.3e significant digits.\n",1-log10(2*fabs(BzrN-BzrC)/(fabs(BzrN)+fabs(BzrC))));
        all_agree=0;
    }
    if(BxlN!=BxlC) {
        printf("Bxl agrees to %.3e significant digits.\n",1-log10(2*fabs(BxlN-BxlC)/(fabs(BxlN)+fabs(BxlC))));
        all_agree=0;
    }
    if(BylN!=BylC) {
        printf("Byl agrees to %.3e significant digits.\n",1-log10(2*fabs(BylN-BylC)/(fabs(BylN)+fabs(BylC))));
        all_agree=0;
    }
    if(BzlN!=BzlC) {
        printf("Bzl agrees to %.3e significant digits.\n",1-log10(2*fabs(BzlN-BzlC)/(fabs(BzlN)+fabs(BzlC))));
        all_agree=0;
    }
    if (all_agree) {printf("All quantities are bit-wise identical in this flux direction!\n");}
    else {printf("All other quantities are bit-wise identical in this flux direction!\n");}
}    

int main() {
    for(int flux_dirn=1; flux_dirn<4; flux_dirn++) compare_ppm_calc(flux_dirn);
}

Writing Start-to-Finish-UnitTests/PPM_UnitTest//PPM_unit_test.C

Step 2.a: Compile and run the code [Back to top]¶

$$\label{compile_run}$$

Now that we have our file, we can compile it and run the executable.

In [9]:

import time
results_file = "out.txt"

print("Now compiling, should take ~2 seconds...\n")
start = time.time()
cmd.C_compile(os.path.join(Ccodesdir,"PPM_unit_test.C"), os.path.join(outdir,"PPM_unit_test"))
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")

# Change to output directory
os.chdir(outdir)

print("Now running...\n")
start = time.time()
cmd.Execute(os.path.join("PPM_unit_test"),file_to_redirect_stdout=results_file)
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")

Now compiling, should take ~2 seconds...

Compiling executable...
(EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Start-to-Finish-UnitTests/PPM_UnitTest/PPM_unit_test.C -o Start-to-Finish-UnitTests/PPM_UnitTest/output/PPM_unit_test -lm`...
(BENCH): Finished executing in 0.6103384494781494 seconds.
Finished compilation.
Finished in 0.6203911304473877 seconds.


Now running...

(EXEC): Executing `taskset -c 0,1,2,3,4,5 ./PPM_unit_test `...
(BENCH): Finished executing in 0.21147823333740234 seconds.
Finished in 0.22610020637512207 seconds.

Here, we will check the output and make sure all the quantities agree.

In [10]:

with open(results_file,"r") as file:
    for i in range(6):
        output = file.readline()
        print(output)
        if i%2==1 and output!="All quantities are bit-wise identical in this flux direction!\n":
            # If this isn't (0 indexed) line 1,3,5 of this file, something went wrong!
            sys.exit(1)

Checking results in the x-direction...

All quantities are bit-wise identical in this flux direction!

Checking results in the x-direction...

All quantities are bit-wise identical in this flux direction!

Checking results in the x-direction...

All quantities are bit-wise identical in this flux direction!

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

In [11]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface

# Change to NRPy directory
os.chdir("../../../")

cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-PPM",location_of_template_file=os.path.join(".."))

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

Start-to-Finish Example: Unit Testing GiRaFFE_NRPy: Piecewise-Parabolic Method¶