GiRaFFE_NRPy
: Fluxes of $\tilde{S}_i$¶This notebook validates our new, NRPyfied HLLE solver against the function from the original GiRaFFE
that calculates the flux for $\tilde{S}_i$ according to the method of Harten, Lax, von Leer, and Einfeldt (HLLE), assuming that we have calculated the values of the flux on the cell faces according to the piecewise-parabolic method (PPM) of Colella and Woodward (1984), modified for the case of GRFFE.
Module Status: Validated: This code has passed unit tests against the original GiRaFFE
version.
Validation Notes: This demonstrates the validation of Tutorial-GiRaFFE_NRPy-Stilde-flux.
It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. Here, we expect our functions GRFFE__S_i__flux_in_dir_*.h
to produce identical output to the function GRFFE__S_i__flux.C
in the original GiRaFFE
. It should be noted that the two codes handle the parameter flux_dirn
(the direction in which the code is presently calculating the flux through the cell) differently; in the original GiRaFFE
, the function GRFFE__S_i__flux()
expects a parameter flux_dirn
with value 1, 2, or 3, corresponding to the functions GRFFE__S_i__flux_in_dir_0()
, GRFFE__S_i__flux_in_dir_1()
, and GRFFE__S_i__flux_in_dir_2()
, respectively, in GiRaFFE_NRPy
.
This notebook validates the NRPyfied Stilde_flux solver against the original GiRaFFE
code. This will be done at a point with a random but realistic spacetime and a variety of magnetic fields and Valencia velocities to test edge cases.
We'll write this in C because the codes we want to test are already written that way, and we would like to avoid modifying the files as much as possible. To do so, we will print the C code to a file. We will begin by including core functionality. We will also define standard parameters needed for GRFFE and NRPy+.
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.
This notebook is organized as follows
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 will also declare the gridfunctions that are needed for this portion of the code.
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
The last things NRPy+ will require are the definition of type REAL
and, of course, the functions we are testing. These files are generated on the fly.
import shutil, os, sys # Standard Python modules for multiplatform OS-level functions
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
out_dir = "Validation"
cmd.mkdir(out_dir)
thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-Stilde_flux"
import GiRaFFE_NRPy.Stilde_flux as Sf
# 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,gammadet_face = gri.register_gridfunctions("AUXEVOL",["alpha_face","gammadet_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)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")
gri.register_gridfunctions("AUXEVOL","cmax_x")
gri.register_gridfunctions("AUXEVOL","cmin_x")
gri.register_gridfunctions("AUXEVOL","cmax_y")
gri.register_gridfunctions("AUXEVOL","cmin_y")
gri.register_gridfunctions("AUXEVOL","cmax_z")
gri.register_gridfunctions("AUXEVOL","cmin_z")
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
# And the function to which we'll write the output data:
Stilde_fluxD = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_fluxD",DIM=3)
# And now, we'll write the files
# In practice, the C functions should only loop over the interior; here, however, we can loop over
# all points and set fewer parameters.
subdir = "RHSs"
cmd.mkdir(os.path.join(out_dir,subdir))
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,
outCparams = "outCverbose=False,CSE_sorting=none",write_cmax_cmin=True)
Output C function calculate_Stilde_rhsD() to file Validation/RHSs/calculate_Stilde_rhsD.h
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.
# 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.
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;
// Standard GRFFE parameters:
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))
GiRaFFE
for comparison [Back to top]¶We'll also need to download the files in question from the GiRaFFE
bitbucket repository. This code was originally written by Leo Werneck in the IllinoisGRMHD documentation; we have modified it to download the files we want. Of note is the addition of the for
loop since we need three files (The function GRFFE__S_i__flux()
depends on two other files for headers and functions).
# 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/5611b2f0b17135538c9d9d17c7da062abe0401b6/GiRaFFE/src/GRFFE__S_i__flux.C"
]
original_file_name = ["GiRaFFE_headers.h",
"inlined_functions.C",
"GRFFE__S_i__flux.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)
Stilde_flux_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.Now we can write our C code. First, we will import our usual libraries and define the various constants and macros we need, taking care to imitate CCTK functionality wherever necessary. In the main function, we will fill all relevant arrays with (appropriate) random values. That is, if a certain gridfunction should never be negative, we will make sure to only generate positive numbers for it. We must also contend with the fact that in NRPy+, we chose to use the Valencia 3-velocity $v^i_{(n)}$, while in ETK, we used the drift velocity $v^i$; the two are related by $$v^i = \alpha v^i_{(n)} - \beta^i.$$
%%writefile $out_dir/Stilde_flux_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.
// Standard GRFFE parameters:
const double GAMMA_SPEED_LIMIT = 2000.0;
#define NGHOSTS 1
// Standard NRPy+ memory access:
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )
// Standard CCTK memory access:
#define CCTK_GFINDEX3D(thing,i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * (k) ) )
// Let's also #define the NRPy+ gridfunctions
#define ALPHA_FACEGF 0
#define GAMMA_FACEDD00GF 1
#define GAMMA_FACEDD01GF 2
#define GAMMA_FACEDD02GF 3
#define GAMMA_FACEDD11GF 4
#define GAMMA_FACEDD12GF 5
#define GAMMA_FACEDD22GF 6
#define BETA_FACEU0GF 7
#define BETA_FACEU1GF 8
#define BETA_FACEU2GF 9
#define VALENCIAV_RU0GF 10
#define VALENCIAV_RU1GF 11
#define VALENCIAV_RU2GF 12
#define B_RU0GF 13
#define B_RU1GF 14
#define B_RU2GF 15
#define VALENCIAV_LU0GF 16
#define VALENCIAV_LU1GF 17
#define VALENCIAV_LU2GF 18
#define B_LU0GF 19
#define B_LU1GF 20
#define B_LU2GF 21
#define CMAX_XGF 22
#define CMIN_XGF 23
#define CMAX_YGF 24
#define CMIN_YGF 25
#define CMAX_ZGF 26
#define CMIN_ZGF 27
#define STILDE_FLUX_HLLED0GF 28
#define STILDE_FLUX_HLLED1GF 29
#define STILDE_FLUX_HLLED2GF 30
#define NUM_AUXEVOL_GFS 31
#define STILDED0GF 0
#define STILDED1GF 1
#define STILDED2GF 2
#define NUM_EVOL_GFS 3
// The NRPy+ versions of the function. These should require relatively little modification.
// We will need this define, though:
#define REAL double
#include "declare_Cparameters_struct.h"
#include "RHSs/calculate_Stilde_flux_D0.h"
#include "RHSs/calculate_Stilde_flux_D1.h"
#include "RHSs/calculate_Stilde_flux_D2.h"
// Some needed workarounds to get the ETK version of the code to work
#define CCTK_REAL double
#define DECLARE_CCTK_PARAMETERS //
struct cGH{};
const cGH cctkGH;
// And include the code we want to test against
#include "GiRaFFE_headers.h"
#include "inlined_functions.C"
#include "GRFFE__S_i__flux.C"
// Include this one later because it requries kronecker_delta from GiRaFFE_headers.h
#include "RHSs/calculate_Stilde_rhsD.h"
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"
// This is the array to which we'll write the NRPy+ variables.
int total_num_gridpoints = NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0;
REAL *auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * total_num_gridpoints);
REAL *rhs_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * total_num_gridpoints);
// These are the arrays to which we will write the ETK variables.
CCTK_REAL dxi[4] = {0,invdx0,invdx1,invdx2}; // Note that we have to match flux_dirn=1:3
CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX];
CCTK_REAL Ur[MAXNUMVARS];
CCTK_REAL Ul[MAXNUMVARS];
CCTK_REAL FACEVAL[NUMVARS_FOR_METRIC_FACEVALS];
CCTK_REAL cmax[total_num_gridpoints]; CCTK_REAL cmin[total_num_gridpoints];
CCTK_REAL st_x_flux[total_num_gridpoints]; CCTK_REAL st_y_flux[total_num_gridpoints]; CCTK_REAL st_z_flux[total_num_gridpoints];
CCTK_REAL st_x_rhs[total_num_gridpoints]; CCTK_REAL st_y_rhs[total_num_gridpoints]; CCTK_REAL st_z_rhs[total_num_gridpoints];
for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++) {
for(int gf = 0; gf < NUM_EVOL_GFS; gf++) {
// Initialize the RHSs to 0:
rhs_gfs[IDX4S(gf,i0,i1,i2)] = 0;
}
st_x_rhs[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;
st_y_rhs[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;
st_z_rhs[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;
}
// Note that we use the CCTK-style flux_dirn indexing, and run from 1 to 3. This is for compatibility with
// the old code, since GiRaFFE_NRPy uses this integer for fewer things.
for(int flux_dirn = 1; flux_dirn <= 3; flux_dirn++) {
for (int i2 = NGHOSTS; i2 < NGHOSTS+Nxx2+1; i2++) for (int i1 = NGHOSTS; i1 < NGHOSTS+Nxx1+1; i1++) for (int i0 = NGHOSTS; i0 < NGHOSTS+Nxx0+1; i0++) {
// for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++) {
// 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
//const long int seed = 1574393335;
//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(ALPHA_FACEGF, i0,i1,i2)] =1.0+(double)rand()/RAND_MAX-0.2-0.1;
const double alpha = auxevol_gfs[IDX4S(ALPHA_FACEGF, i0,i1,i2)];
METRIC_LAP_PSI4[LAPSE] = alpha;
//METRIC_LAP_PSI4[LAPM1] = METRIC_LAP_PSI4[LAPSE]-1;
auxevol_gfs[IDX4S(GAMMA_FACEDD00GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMA_FACEDD01GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMA_FACEDD02GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMA_FACEDD11GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMA_FACEDD12GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
auxevol_gfs[IDX4S(GAMMA_FACEDD22GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;
// Generated by NRPy+:
const double gammaDD00 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF, i0,i1,i2)];
const double gammaDD01 = auxevol_gfs[IDX4S(GAMMA_FACEDD01GF, i0,i1,i2)];
const double gammaDD02 = auxevol_gfs[IDX4S(GAMMA_FACEDD02GF, i0,i1,i2)];
const double gammaDD11 = auxevol_gfs[IDX4S(GAMMA_FACEDD11GF, i0,i1,i2)];
const double gammaDD12 = auxevol_gfs[IDX4S(GAMMA_FACEDD12GF, i0,i1,i2)];
const double gammaDD22 = auxevol_gfs[IDX4S(GAMMA_FACEDD22GF, 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;
METRIC_LAP_PSI4[PSI6] = sqrt(tmp5);
METRIC_LAP_PSI4[PSI2] = pow(METRIC_LAP_PSI4[PSI6],1.0/3.0);
METRIC_LAP_PSI4[PSI4] = METRIC_LAP_PSI4[PSI2]*METRIC_LAP_PSI4[PSI2];
const double Psim4 = 1.0/METRIC_LAP_PSI4[PSI4];
METRIC_LAP_PSI4[PSIM4] = Psim4;
// Copied from the ETK implementation
CCTK_REAL gtxxL = gammaDD00*Psim4;
CCTK_REAL gtxyL = gammaDD01*Psim4;
CCTK_REAL gtxzL = gammaDD02*Psim4;
CCTK_REAL gtyyL = gammaDD11*Psim4;
CCTK_REAL gtyzL = gammaDD12*Psim4;
CCTK_REAL gtzzL = gammaDD22*Psim4;
/*********************************
* Apply det gtij = 1 constraint *
*********************************/
const CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL -
gtxzL * gtyyL * gtxzL - gtxyL * gtxyL * gtzzL - gtxxL * gtyzL * gtyzL;
/*const CCTK_REAL gtijdet_Fm1o3 = fabs(1.0/cbrt(gtijdet));
gtxxL = gtxxL * gtijdet_Fm1o3;
gtxyL = gtxyL * gtijdet_Fm1o3;
gtxzL = gtxzL * gtijdet_Fm1o3;
gtyyL = gtyyL * gtijdet_Fm1o3;
gtyzL = gtyzL * gtijdet_Fm1o3;
gtzzL = gtzzL * gtijdet_Fm1o3;*/
FACEVAL[GXX] = gtxxL;
FACEVAL[GXY] = gtxyL;
FACEVAL[GXZ] = gtxzL;
FACEVAL[GYY] = gtyyL;
FACEVAL[GYZ] = gtyzL;
FACEVAL[GZZ] = gtzzL;
FACEVAL[GUPXX] = ( gtyyL * gtzzL - gtyzL * gtyzL )/gtijdet;
FACEVAL[GUPYY] = ( gtxxL * gtzzL - gtxzL * gtxzL )/gtijdet;
FACEVAL[GUPZZ] = ( gtxxL * gtyyL - gtxyL * gtxyL )/gtijdet;
/*auxevol_gfs[IDX4S(GAMMA_FACEUU00GF, i0,i1,i2)] = FACEVAL[GUPXX];
auxevol_gfs[IDX4S(GAMMA_FACEUU11GF, i0,i1,i2)] = FACEVAL[GUPYY];
auxevol_gfs[IDX4S(GAMMA_FACEUU22GF, i0,i1,i2)] = FACEVAL[GUPZZ];*/
auxevol_gfs[IDX4S(BETA_FACEU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
const double betax = auxevol_gfs[IDX4S(BETA_FACEU0GF, i0,i1,i2)];
FACEVAL[SHIFTX] = betax;
auxevol_gfs[IDX4S(BETA_FACEU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
const double betay = auxevol_gfs[IDX4S(BETA_FACEU1GF, i0,i1,i2)];
FACEVAL[SHIFTY] = betay;
auxevol_gfs[IDX4S(BETA_FACEU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;
const double betaz = auxevol_gfs[IDX4S(BETA_FACEU2GF, i0,i1,i2)];
FACEVAL[SHIFTZ] = betaz;
/* Generate physically meaningful speeds */
auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
/* Superluminal speeds for testing */
/*auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;*/
Ur[VX] = alpha*auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)]-betax;
Ur[VY] = alpha*auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)]-betay;
Ur[VZ] = alpha*auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)]-betaz;
auxevol_gfs[IDX4S(B_RU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Ur[BX_CENTER] = auxevol_gfs[IDX4S(B_RU0GF, i0,i1,i2)];
auxevol_gfs[IDX4S(B_RU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Ur[BY_CENTER] = auxevol_gfs[IDX4S(B_RU1GF, i0,i1,i2)];
// Set Bz to enforce orthogonality of u and B; lower index on drift v^i, then set
// Bz = -(ux*Bx+uy*By)/uz
REAL ux,uy,uz;
ux = gammaDD00*(Ur[VX]+betax) + gammaDD01*(Ur[VY]+betay) + gammaDD02*(Ur[VZ]+betaz);
uy = gammaDD01*(Ur[VX]+betax) + gammaDD11*(Ur[VY]+betay) + gammaDD12*(Ur[VZ]+betaz);
uz = gammaDD02*(Ur[VX]+betax) + gammaDD12*(Ur[VY]+betay) + gammaDD22*(Ur[VZ]+betaz);
auxevol_gfs[IDX4S(B_RU2GF, i0,i1,i2)] = -(ux*Ur[BX_CENTER] + uy*Ur[BY_CENTER])/uz;
Ur[BZ_CENTER] = auxevol_gfs[IDX4S(B_RU2GF, i0,i1,i2)];
/* Generate physically meaningful speeds */
auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
/* Superluminal speeds for testing */
/*auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;
auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;*/
Ul[VX] = alpha*auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)]-betax;
Ul[VY] = alpha*auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)]-betay;
Ul[VZ] = alpha*auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)]-betaz;
auxevol_gfs[IDX4S(B_LU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Ul[BX_CENTER] = auxevol_gfs[IDX4S(B_LU0GF, i0,i1,i2)];
auxevol_gfs[IDX4S(B_LU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;
Ul[BY_CENTER] = auxevol_gfs[IDX4S(B_LU1GF, i0,i1,i2)];
// Set Bz to enforce orthogonality of u and B; lower index on drift v^i, then set
// Bz = -(ux*Bx+uy*By)/uz
ux = gammaDD00*(Ul[VX]+betax) + gammaDD01*(Ul[VY]+betay) + gammaDD02*(Ul[VZ]+betaz);
uy = gammaDD01*(Ul[VX]+betax) + gammaDD11*(Ul[VY]+betay) + gammaDD12*(Ul[VZ]+betaz);
uz = gammaDD02*(Ul[VX]+betax) + gammaDD12*(Ul[VY]+betay) + gammaDD22*(Ul[VZ]+betaz);
auxevol_gfs[IDX4S(B_LU2GF, i0,i1,i2)] = -(ux*Ul[BX_CENTER] + uy*Ul[BY_CENTER])/uz;
Ul[BZ_CENTER] = auxevol_gfs[IDX4S(B_LU2GF, i0,i1,i2)];
// Compute the fluxes for the ETK version
const int index = CCTK_GFINDEX3D(cctkGH,i0,i1,i2);
GRFFE__S_i__flux(i0,i1,i2,flux_dirn,Ul,Ur,FACEVAL,METRIC_LAP_PSI4,cmax[index],cmin[index],st_x_flux[index],st_y_flux[index],st_z_flux[index]);
}
// Compute the fluxes for the NRPy+ version. Note the one-offset between flux_dirn and the function names,
// since we set flux_dirn for backwards compatibility.
if(flux_dirn==1) calculate_Stilde_flux_D0(¶ms,auxevol_gfs,rhs_gfs);
if(flux_dirn==2) calculate_Stilde_flux_D1(¶ms,auxevol_gfs,rhs_gfs);
if(flux_dirn==3) calculate_Stilde_flux_D2(¶ms,auxevol_gfs,rhs_gfs);
calculate_Stilde_rhsD(flux_dirn,¶ms,auxevol_gfs,rhs_gfs);
for(int i0=0;i0<Nxx_plus_2NGHOSTS0-1;i0++) for(int i1=0;i1<Nxx_plus_2NGHOSTS1-1;i1++) for(int i2=0;i2<Nxx_plus_2NGHOSTS2-1;i2++) {
const int index = CCTK_GFINDEX3D(cctkGH,i0,i1,i2);
const int indexp1 = CCTK_GFINDEX3D(cctkGH,i0+kronecker_delta[flux_dirn][0],i1+kronecker_delta[flux_dirn][1],i2+kronecker_delta[flux_dirn][2]);
st_x_rhs[index] += (st_x_flux[index] - st_x_flux[indexp1] ) * dxi[flux_dirn];
st_y_rhs[index] += (st_y_flux[index] - st_y_flux[indexp1] ) * dxi[flux_dirn];
st_z_rhs[index] += (st_z_flux[index] - st_z_flux[indexp1] ) * dxi[flux_dirn];
}
}
#define SDA(a,b) 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b))+1.0e-15)
//printf("Valencia 3-velocity (right): %.4e, %.4e, %.4e\n",auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)]);
//printf("Valencia 3-velocity (left): %.4e, %.4e, %.4e\n\n",auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)]);
// We'll compare the output in-between each
int i0 = 1, i1 = 1, i2 = 1;
const int index = CCTK_GFINDEX3D(cctkGH,i0,i1,i2);
printf("%.1f %.1f %.1f\n",SDA(rhs_gfs[IDX4S(STILDED0GF,i0,i1,i2)],st_x_rhs[index]),
SDA(rhs_gfs[IDX4S(STILDED1GF,i0,i1,i2)],st_y_rhs[index]),
SDA(rhs_gfs[IDX4S(STILDED2GF,i0,i1,i2)],st_z_rhs[index]));
//printf("NRPy+ Results: %.16e, %.16e, %.16e\n",rhs_gfs[IDX4S(STILDED0GF, i0, i1, i2)],rhs_gfs[IDX4S(STILDED1GF, i0,i1,i2)],rhs_gfs[IDX4S(STILDED2GF, i0,i1,i2)]);
//printf("ETK Results: %.16e, %.16e, %.16e\n",st_x_rhs[index],st_y_rhs[index],st_z_rhs[index]);
}
Overwriting Validation/Stilde_flux_unit_test.C
import cmdline_helper as cmd
results_file = "out.txt"
cmd.C_compile(os.path.join(out_dir,"Stilde_flux_unit_test.C"), os.path.join(out_dir,"Stilde_flux_unit_test"))
cmd.Execute(os.path.join(out_dir,"Stilde_flux_unit_test"),file_to_redirect_stdout=os.path.join(out_dir,results_file))
# with open(os.path.join(out_dir,results_file),"r") as file:
# for i in range(18):
# output = file.readline()
# if "NRPy+" in output:
# substrings = output.split(": ")
# substrings = substrings[1].split(", ")
# for j in range(3):
# if float(substrings[j])!=0.0:
# sys.exit(1)
# print(output)
Compiling executable... (EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops Validation/Stilde_flux_unit_test.C -o Validation/Stilde_flux_unit_test -lm`... cc1plus: warning: command line option ‘-std=gnu99’ is valid for C/ObjC but not for C++ (BENCH): Finished executing in 1.0110466480255127 seconds. Finished compilation. (EXEC): Executing `taskset -c 0,1,2,3 ./Validation/Stilde_flux_unit_test `... (BENCH): Finished executing in 0.20968246459960938 seconds.
Below are the numbers we care about. These are the Significant Digits of Agreement between the HLLE fluxes 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 14 significant digits of agreement.
Note that in the case of very high velocities, numerical error will accumulate and reduce agreement significantly due to a catastrophic cancellation.
print("Now comparing the significant digits of agreements in each direction...")
with open(os.path.join(out_dir,results_file),"r") as file:
output = file.readline()
output = output.split()
for i in range(3):
print("SDA in Stilde_rhsD"+str(i)+" = "+output[i])
if float(output[i])<14:
sys.exit(1)
Now comparing the significant digits of agreements in each direction... SDA in Stilde_rhsD0 = 15.6 SDA in Stilde_rhsD1 = 15.6 SDA in Stilde_rhsD2 = 15.8
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-Stilde_flux.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-Stilde_flux",location_of_template_file=os.path.join(".."))
Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Stilde_flux.tex, and compiled LaTeX file to PDF file Tutorial-Start_to_Finish_UnitTest- GiRaFFE_NRPy-Stilde_flux.pdf