GiRaFFE_NRPy
: Conservative-to-Primitive and Primitive-to-Conservative Solvers¶Notebook Status: Validated
Validation Notes: This module will validate the routines in Tutorial-GiRaFFE_NRPy-C2P_P2C.
This notebook validates the NRPyfied C2P and P2C solvers against the original GiRaFFE
code. This will be done at a point with a random but realistic spacetime and a variety of Poynting fluxes and Valencia velocities to test edge cases.
When this notebook is run, the significant digits of agreement between the old GiRaFFE
and new GiRaFFE_NRPy
versions of the algorithm outputs of Poynting flux and Valencia three-velocity will be printed to the screen right after the code is run here.
This notebook is organized as follows
C2P_P2C_unit_test.c
: The Main C CodeWe'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 will also declare the gridfunctions that are needed for this portion of the code.
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 * # 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 cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
out_dir = "Validation/"
cmd.mkdir(out_dir)
thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C"
# Register the gridfunctions we need for this function
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU")
First, we'll use NRPy+ to build the C function that will apply fixes to $\tilde{S}_i$ and recompute the velocity to match, along with the current sheet prescription. Note that the NRPy+ version of this code also speed-limits the velocities.
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])
]
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,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="outCverbose=False").replace("IDX4","IDX4S"),
loopopts ="AllPoints,Read_xxs")
Output C function GiRaFFE_NRPy_cons_to_prims() to file Validation/GiRaFFE_NRPy_cons_to_prims.h
# Declare this symbol:
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
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=sp.simplify(C2P_P2C.StildeD[0])),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=sp.simplify(C2P_P2C.StildeD[1])),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=sp.simplify(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,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="outCverbose=False").replace("IDX4","IDX4S"),
loopopts ="AllPoints")
Output C function GiRaFFE_NRPy_prims_to_cons() to file Validation/GiRaFFE_NRPy_prims_to_cons.h
# First download the original IllinoisGRMHD source code
import urllib
original_file_url = ["https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/GiRaFFE/src/GiRaFFE_headers.h",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/GiRaFFE/src/inlined_functions.C",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/driver_conserv_to_prims_FFE.C",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/compute_conservatives_FFE.C",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
]
original_file_name = ["GiRaFFE_headers.h",
"inlined_functions.C",
"driver_conserv_to_prims_FFE.C",
"compute_conservatives_FFE.C",
"convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
]
for i in range(len(original_file_url)):
original_file_path = os.path.join(out_dir,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)
free_parameters.h
[Back to top]¶Based on declared NRPy+ Cparameters, first we generate declare_Cparameters_struct.h
, set_Cparameters_default.h
, and set_Cparameters[-SIMD].h
.
Then we output free_parameters.h
, which sets some basic grid parameters as well as the speed limit parameter we need for this function.
# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
# par.generate_Cparameters_Ccodes(os.path.join(out_dir))
# Step 3.d.ii: Set free_parameters.h
with open(os.path.join(out_dir,"free_parameters.h"),"w") as file:
file.write("""
// Set free-parameter values.
const int NGHOSTS = 0;
// Set free-parameter values for the initial data.
// Override parameter defaults with values based on command line arguments and NGHOSTS.
const int Nx0x1x2 = 1;
params.Nxx0 = Nx0x1x2;
params.Nxx1 = Nx0x1x2;
params.Nxx2 = Nx0x1x2;
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] = {0.0,0.0,0.0};
const REAL xxmax[3] = {0.1,0.1,0.1};
params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;
params.GAMMA_SPEED_LIMIT = 2000.0;
\n""")
# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(out_dir))
The original GiRaFFE
code depends on some functionalities of the CCTK. Since we only care about this one small function, we can get around this by creating some nearly-empty, non-functional files that can be included to satisfy the pre-processor without changing functionality. We will later replace what little functionality we need with some basic global variables and macros.
#incldue "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"
with open(os.path.join(out_dir,"cctk.h"),"w") as file:
file.write("""//""")
with open(os.path.join(out_dir,"cctk_Arguments.h"),"w") as file:
file.write("""#define DECLARE_CCTK_ARGUMENTS //
#define CCTK_ARGUMENTS void
""")
with open(os.path.join(out_dir,"cctk_Parameters.h"),"w") as file:
file.write("""#define DECLARE_CCTK_PARAMETERS //
""")
with open(os.path.join(out_dir,"Symmetry.h"),"w") as file:
file.write("""//""")
C2P_P2C_unit_test.c
: The Main C Code [Back to top]¶Now that we have our vector potential and analytic magnetic field to compare against, 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.
%%writefile $out_dir/C2P_P2C_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"
// Standard NRPy+ memory access:
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
// Memory access definitions for NRPy+
#define GAMMADD00GF 0
#define GAMMADD01GF 1
#define GAMMADD02GF 2
#define GAMMADD11GF 3
#define GAMMADD12GF 4
#define GAMMADD22GF 5
#define BETAU0GF 6
#define BETAU1GF 7
#define BETAU2GF 8
#define ALPHAGF 9
#define BU0GF 10
#define BU1GF 11
#define BU2GF 12
#define VALENCIAVU0GF 13
#define VALENCIAVU1GF 14
#define VALENCIAVU2GF 15
#define NUM_AUXEVOL_GFS 16
#define STILDED0GF 0
#define STILDED1GF 1
#define STILDED2GF 2
#define NUM_EVOL_GFS 3
// Include the functions that we want to test:
#include "GiRaFFE_NRPy_cons_to_prims.h"
#include "GiRaFFE_NRPy_prims_to_cons.h"
// Define CCTK macros
#define CCTK_REAL double
#define CCTK_INT int
struct cGH{};
const cGH* cctkGH;
// GiRaFFE parameters in ETK
const CCTK_REAL min_radius_inside_of_which_conserv_to_prims_FFE_and_FFE_evolution_is_DISABLED = -1;
const int current_sheet_null_v = 1;
// More definitions to interface with ETK code:
const int cctk_lsh[3] = {1,1,1};
CCTK_REAL gxx[1];
CCTK_REAL gxy[1];
CCTK_REAL gxz[1];
CCTK_REAL gyy[1];
CCTK_REAL gyz[1];
CCTK_REAL gzz[1];
CCTK_REAL alp[1];
CCTK_REAL gtxx[1];
CCTK_REAL gtxy[1];
CCTK_REAL gtxz[1];
CCTK_REAL gtyy[1];
CCTK_REAL gtyz[1];
CCTK_REAL gtzz[1];
CCTK_REAL gtupxx[1];
CCTK_REAL gtupxy[1];
CCTK_REAL gtupxz[1];
CCTK_REAL gtupyy[1];
CCTK_REAL gtupyz[1];
CCTK_REAL gtupzz[1];
CCTK_REAL phi_bssn[1];
CCTK_REAL psi_bssn[1];
CCTK_REAL lapm1[1];
CCTK_REAL betax[1];
CCTK_REAL betay[1];
CCTK_REAL betaz[1];
CCTK_REAL mhd_st_x[1];
CCTK_REAL mhd_st_y[1];
CCTK_REAL mhd_st_z[1];
CCTK_REAL vx[1];
CCTK_REAL vy[1];
CCTK_REAL vz[1];
CCTK_REAL Bx[1];
CCTK_REAL By[1];
CCTK_REAL Bz[1];
CCTK_REAL x[1];
CCTK_REAL y[1];
CCTK_REAL z[1];
CCTK_REAL r[1];
// Define dz in CCTK
CCTK_REAL cactus_dz;
#define CCTK_DELTA_SPACE(i) cactus_dz
CCTK_REAL GAMMA_SPEED_LIMIT = 2000.0;
// Dummy ETK function:
#define CCTK_GFINDEX3D(cctkGH,i,j,k) 0
#define CCTK_VInfo(a01,a02,a03,a04,a05,a06,a07,a08,a09,a10,a11,a12) //
#define CCTK_VWarn(b01,b02,b03,b04,b05,b06,b07,b08,b09,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25) //
#include "driver_conserv_to_prims_FFE.C"
#include "compute_conservatives_FFE.C"
#include "convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
int main() {
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"
// Now that we've calculated dxx2, we can define a cactus equivalent
cactus_dz = dxx2;
// We'll define all indices to be 0. No need to complicate memory access
const int i0 = 0;
const int i1 = 0;
const int i2 = 0;
// This is the array to which we'll write the NRPy+ variables.
REAL *auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS);
REAL *evol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS);
// Step 0e: Set up cell-centered Cartesian 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] + (j-NGHOSTS)*dxx0;
for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;
for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;
x[0] = xx[0][0];
y[0] = xx[1][0];
z[0] = xx[2][0];
r[0] = sqrt(xx[0][0]*xx[0][0] + xx[1][0]*xx[1][0] + xx[2][0]*xx[2][0]);
// Now, it's time to make the random numbers.
//const long int seed = time(NULL); // seed = 1570632212; is an example of a seed that produces
// bad agreement for high speeds
//srand(seed); // Set the seed
//printf("seed for random number generator = %ld; RECORD IF AGREEMENT IS BAD\\n\\n",seed);
// We take care to make sure the corresponding quantities have the SAME value.
auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
const double alphaL = auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)];
alp[0] = alphaL;
auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
// Generated by NRPy+:
const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];
const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];
const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];
const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];
const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];
const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];
/*
* NRPy+ Finite Difference Code Generation, Step 2 of 1: Evaluate SymPy expressions and write to main memory:
*/
const double tmp0 = gammaDD11*gammaDD22;
const double tmp1 = pow(gammaDD12, 2);
const double tmp2 = gammaDD02*gammaDD12;
const double tmp3 = pow(gammaDD01, 2);
const double tmp4 = pow(gammaDD02, 2);
const double tmp5 = gammaDD00*tmp0 - gammaDD00*tmp1 + 2*gammaDD01*tmp2 - gammaDD11*tmp4 - gammaDD22*tmp3;
const double tmp6 = 1.0/tmp5;
// Set the ETK metric:
gxx[0] = gammaDD00;
gxy[0] = gammaDD01;
gxz[0] = gammaDD02;
gyy[0] = gammaDD11;
gyz[0] = gammaDD12;
gzz[0] = gammaDD22;
auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
const double betaxL = auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)];
betax[0] = betaxL;
auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
const double betayL = auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)];
betay[0] = betayL;
auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.02-0.01;
const double betazL = auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)];
betaz[0] = betazL;
/* Generate physically meaningful speeds */
auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
vx[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-betaxL;
auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
vy[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-betayL;
auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
vz[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-betazL;
/* Superluminal speeds for testing */
/*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
vx[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-betaxL;
auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
vy[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-betayL;
auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
vz[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-betazL;*/
auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Bx[0] = auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)];
auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
By[0] = auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)];
auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Bz[0] = auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)];
evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
mhd_st_x[0] = evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)];
evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
mhd_st_y[0] = evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)];
evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
mhd_st_z[0] = evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)];
GiRaFFE_conserv_to_prims_FFE();
GiRaFFE_NRPy_cons_to_prims(¶ms,xx,auxevol_gfs,evol_gfs);
GiRaFFE_NRPy_prims_to_cons(¶ms,xx,auxevol_gfs,evol_gfs);
//printf("Checking the Poynting Fluxes:\n");
printf("%.1f, %.1f, %.1f\n",1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)]-mhd_st_x[0])/(fabs(evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)])+fabs(mhd_st_x[0]))),
1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)]-mhd_st_y[0])/(fabs(evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)])+fabs(mhd_st_y[0]))),
1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)]-mhd_st_z[0])/(fabs(evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)])+fabs(mhd_st_z[0]))));
//printf("NRPy: %.15e,%.15e,%.15e\n",evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)],evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)],evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)]);
//printf("CCTK: %.15e,%.15e,%.15e\n",mhd_st_x[0],mhd_st_y[0],mhd_st_z[0]);
//printf("\n\n");
//printf("Checking the Valencia Velocities:\n");
printf("%.1f, %.1f, %.1f\n",1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-(vx[0]+betaxL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)])+fabs((vx[0]+betaxL)/alphaL))),
1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-(vy[0]+betayL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)])+fabs((vy[0]+betayL)/alphaL))),
1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-(vz[0]+betazL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)])+fabs((vz[0]+betazL)/alphaL))));
//printf("NRPy: %.15e,%.15e,%.15e\n",auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]);
//printf("CCTK: %.15e,%.15e,%.15e\n",(vx[0]+betaxL)/alphaL,(vy[0]+betayL)/alphaL,(vz[0]+betazL)/alphaL);
}
Writing Validation//C2P_P2C_unit_test.C
import time
results_file = "out.txt"
print("Now compiling, should take ~2 seconds...\n")
start = time.time()
# Can't use C here; must use C++
# This command should work on all systems:
!g++ -Ofast -fopenmp -march=native -funroll-loops Validation/C2P_P2C_unit_test.C -o Validation/C2P_P2C_unit_test -lm
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")
os.chdir(os.path.join("Validation"))
cmd.Execute(os.path.join("C2P_P2C_unit_test"),file_to_redirect_stdout=os.path.join(results_file))
os.chdir(os.path.join(".."))
Now compiling, should take ~2 seconds... Finished in 0.48241138458251953 seconds. (EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./C2P_P2C_unit_test `... (BENCH): Finished executing in 0.20868754386901855 seconds.
Below are the numbers we care about. These are the Significant Digits of Agreement between the corrected StildeD and ValenciavU as computed by NRPy+ and ETK. Each row represents a flux direction; each entry therein corresponds to a component of StildeD. Each pair of outputs should show at least 10 significant digits of agreement.
import numpy as np
with open(os.path.join(out_dir,results_file),"r") as file:
for i in range(2):
output = file.readline()
substrings = output.split(", ")
sda = [float(i) for i in substrings]
enough_digits = np.array(sda)>14
if not enough_digits.all():
sys.exit(1)
While we have chosen to use the Valencia three-velocity in this version of GiRaFFE
, we have also frequently made use of a relationship expressing this in terms of the drift velocity, which was used in the original GiRaFFE
. The usefulness of this relationship to drift velocity extends beyond merely translating the original code. As discussed in Paschalidis, et al., Sec. III.A (just above Eq. 45, with a proof in Appendix A), there is a one-parameter family of velocity definitions that fulfill the GRFFE conditions. The drift velocity sets this parameter to 0, which minimizes the Lorentz factor and guarantees that the four-velocity and magnetic fields are orthogonal to each other. This simplifies the form of $b^\mu$ and quantities that depend on it.
This must be taken into account in developing this unit test, because NRPy+'s GRFFE module defaults to using a definition of $b^\mu$ that does not assume that this criterion is met, while the original GiRaFFE
code assumes this in its C2P and P2C solvers. So, if we do not guarantee that our test data fulfills this criterion, these two different routines will produce different results. We will now go through the derivation of the equation used by GiRaFFE
from first principles to show where this extra term appears.
This is the equation used by GiRaFFE
, pulled from Eqs. 47 and 85 of this paper:
$$\tilde{S}_i = \gamma_{ij} \frac{\bar{v}^j \sqrt{\gamma}B^2}{4 \pi},$$
or $$\tilde{S}_i = \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha},$$
where $\bar{v}^j$ is the Valencia 3-velocity and $v^j$ is the drift velocity.
In IllinoisGRMHD (IGM), the expression used is $\tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i},$ where \begin{align} T^{\mu\nu}_{\rm EM} &= b^2 u^{\mu} u^{\nu} + \frac{1}{2} b^2 g^{\mu\nu} - b^\mu b^\nu \\ b^0 &= \frac{u_j B^j}{\sqrt{4\pi} \alpha} \\ b^i &= \frac{B^i + (u_j B^j) u^i}{\sqrt{4\pi} \alpha u^0} \\ u^i &= u^0 v^i. \end{align}
Now, we'll work to bridge the gap between the expression in terms of the stress-energy tensor and that used in GiRaFFE
. Consider that $T^0_{{\rm EM}i} = g_{i \mu} T^{0 \mu}_{{\rm EM}}.$ We'll thus start with the index lowering operation.
\begin{align}
g_{\mu \xi} T^{\xi \nu}_{\rm EM} &= g_{\mu \xi} \left( b^2 u^{\xi} u^{\nu} + \frac{1}{2} b^2 g^{\xi\nu} - b^\xi b^\nu \right) \\
&= b^2 u_{\mu} u^{\nu} + \frac{1}{2} b^2 \delta^\nu_\mu - b_\mu b^\nu.
\end{align}
After doing so, it is still apparent that we will need $b^2 = g_{\mu\nu} b^\mu b^\nu$, where $$ g_{\mu\nu} = \begin{pmatrix} -\alpha^2 + \beta^k \beta_k & \beta_i \\ \beta_j & \gamma_{ij} \end{pmatrix}. $$ Expanding out the implied sum, we see that \begin{align} b^2 &= g_{00} b^0 b^0 + g_{i0} b^i b^0 + g_{0j} b^0 b^j + g_{ij} b^i b^j \\ &= \left(-\alpha^2 + \beta^k \beta_k\right) b^0 b^0 + \beta_i b^i b^0 + \beta_j b^0 b^j + \gamma_{ij} b^i b^j \\ &= \left(-\alpha^2 + \beta^k \beta_k\right) b^0 b^0 + 2 \beta_i b^i b^0 + \gamma_{ij} b^i b^j \\ \end{align}
Now, it will be useful to plug in the definition of $b^\mu$. Recall that \begin{align} b^0 &= \frac{u_j B^j}{\sqrt{4\pi} \alpha} \\ b^i &= \frac{B^i + (u_j B^j) u^i}{\sqrt{4\pi} \alpha u^0} \\ \end{align}
Then \begin{align} b^2 &= \frac{1}{4 \pi \alpha^2} \left( -\alpha^2 (u_k B^k)^2 + \beta^k \beta_k (u_k B^k)^2 + 2 \beta_i (u_k B^k) (B^i + (u_k B^k) u^i)/u^0 + \gamma_{ij}(B^i + (u_k B^k) u^i)(B^j + (u_k B^k) u^j)/(u^0)^2 \right) \\ \end{align}
This ultimately reduces to (from Paschalidis, et al.): \begin{align} b^2 &= \frac{B^2 + (u_\mu B^\mu)^2}{4 \pi \alpha^2 (u^0)^2} \\ \end{align}
So, we insert this into the stress energy tensor: $$ T^0_{{\rm EM}i} = \frac{B^2 + (u_\mu B^\mu)^2}{4 \pi \alpha^2 (u^0)^2} u_i u^0 - \gamma_{ij} b^j b^0 $$ Note that since $B^\mu$ is purely spatial, $B^0=0$, thus $u_\mu B^\mu = u_k B^k$
\begin{align} T^0_{{\rm EM}i} &= \frac{B^2 + (u_k B^k)^2}{4 \pi \alpha^2 u^0} u_i - \gamma_{ij} \frac{(B^j + (u_k B^k) u^j)(u_k B^k)}{4 \pi \alpha^2 u^0} \\ &= \frac{1}{4 \pi \alpha^2 u^0} \left( B^2 u_i + (u_k B^k)^2 u_i - \gamma_{ij}(B^j + (u_k B^k) u^j)(u_k B^k)\right) \\ &= \frac{1}{4 \pi \alpha^2 u^0} \left( B^2 u_i + (u_k B^k)^2 u_i - (B_i + (u_k B^k) u_i)(u_k B^k)\right) \\ \end{align}Now, we will substitute in $u_j/u_0 = \beta_j + \gamma_{lj} v^l$: \begin{align} &= \frac{1}{4 \pi \alpha^2} \left( B^2 (\beta_i + \gamma_{li} v^l) - B_i (u_k B^k)\right) \\ &= \frac{1}{4 \pi \alpha^2} \left( B^2 (\gamma_{ij} \beta^j + \gamma_{ij} v^j) - B_i (u_k B^k)\right) \\ \end{align}
Now, we multiply by $\alpha \sqrt{\gamma}$: \begin{align} \tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i} &= \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha} - \frac{\sqrt{\gamma}B_i (u_k B^k)}{4 \pi \alpha} \end{align}
We note that there is an extra term here beyond the one that GiRaFFE
uses for its P2C solver. This term exists because we did not make the assumption that $(u_k B^k) = 0$, which follows from the definition of drift velocity (see here, Sec. III.A, just above Eq. 45). Making that assumption, we get the expression used in GiRaFFE
:
\begin{align}
\tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i} &= \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha}
\end{align}
Normally, this relation does not come into play because $u_i \propto v_i \propto \tilde{S}_i$ and $\tilde{S}_i$ is calculated initially as the cross product of the electric and magnetic fields, so is guaranteed to be orthogonal to $B^i$. Here, however, we did nothing to guarantee this beforehand, so the difference shows up. Note also that the entire derivation becomes far easier if you assume this sooner since $b^0 = 0$ and $b^i = B^i / \left(\sqrt{4 \pi} \alpha u^0\right)$
To handle this more smoothly moving forward, we have added an extra function to the GRFFE module that computes $b^\mu$ under the assumption that $u_j B^j = 0$ for use only in GRFFE codes (i.e., without GRHD). This is desirable for more than just this small case since the extra term might cause small errors to propagate.
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-C2P_P2C.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
# cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C",location_of_template_file=os.path.join(".."))