GiRaFFEfood_NRPy
against original, trusted GiRaFFEfood
:¶Notebook Status: Validated
Validation Notes: This module validates all expressions used to set up initial data in
against the C-code implementation of these expressions found in the original (trusted) GiRaFFEfood
Einstein Toolkit thorn, and confirms roundoff-level agreement.
This notebook validates the initial data routines that we will use for GiRaFFE_NRPy
, collectively referred to as GiRaFFEfood_NRPy
. To do so, we will generate the initial data with both our code and the original GiRaFFEfood
code. Then, we will directly compare the velocities and show round-off level agreement between the two.
When this notebook is run, the significant digits of agreement between the old GiRaFFE
and new GiRaFFE_NRPy
versions of the algorithm will be evaluated. If the agreement falls below a thresold, the point, quantity, and level of agreement are reported here.
This notebook is organized as follows
GiRaFFEfood_NRPy_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.
# There are several initial data routines we need to test. We'll control which one we use with a string option
initial_data = "AlignedRotator" # Valid options: "AllTests", "AlignedRotator", "AlfvenWave", "FastWave",
# "DegenAlfvenWave", "ThreeWaves", "FFE_Breakdown"
import os, sys, shutil # 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("..", "Deprecated")
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 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
out_dir = "Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests/"
# 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(out_dir, ignore_errors=True)
cmd.mkdir(out_dir)
thismodule = "Start_to_Finish_UnitTest-GiRaFFEfood_NRPy"
# Register the gridfunctions we need for this function
AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
# gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU")
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gf
gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True)
if initial_data=="AlfvenWave":
desc = "Generate Alfven wave 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="FastWave":
desc = "Generate fast wave 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="DegenAlfvenWave":
desc = "Generate degenerate Alfven wave 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="ThreeWaves":
desc = "Generate three waves 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="FFE_Breakdown":
desc = "Generate FFE breakdown 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="AlignedRotator":
desc = "Generate aligned rotator initial test data for GiRaFFEfood_NRPy."
else:
print("Unsupported Initial Data string "+initial_data+"! Supported ID: AllTests, AlfvenWave, FastWave, DegenAlfvenWave, ThreeWaves, FFE_Breakdown, AlignedRotator, or ExactWald")
name = "GiRaFFE_NRPy_initial_data"
values_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","AD0"),rhs=gf.AD[0]),
lhrh(lhs=gri.gfaccess("out_gfs","AD1"),rhs=gf.AD[1]),
lhrh(lhs=gri.gfaccess("out_gfs","AD2"),rhs=gf.AD[2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=gf.ValenciavU[0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=gf.ValenciavU[1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=gf.ValenciavU[2])
]
outCfunction(
outfile = os.path.join(out_dir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *out_gfs",
body = fin.FD_outputC("returnstring",values_to_print,params="outCverbose=False"),
loopopts ="AllPoints,Read_xxs")
Output C function GiRaFFE_NRPy_initial_data() to file Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests/GiRaFFE_NRPy_initial_data.h
# First download the original GiRaFFE source code
import urllib
original_file_url = [
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFEfood/src/AlfvenWave.cc",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFEfood/src/FastWave.cc",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/b826f6578b0c2b5a43ad9171e65a1b0af88d8b77/GiRaFFEfood/src/DegenAlfvenWave.cc",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/b826f6578b0c2b5a43ad9171e65a1b0af88d8b77/GiRaFFEfood/src/ThreeAlfvenWave.cc",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/b826f6578b0c2b5a43ad9171e65a1b0af88d8b77/GiRaFFEfood/src/FFEBreakdown.cc",
"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFEfood/src/AlignedRotator.cc",
]
original_file_name = [
"AlfvenWave.cc",
"FastWave.cc",
"DegenAlfvenWave.cc",
"ThreeAlfvenWave.cc",
"FFEBreakdown.cc",
"AlignedRotator.cc",
]
for i in range(len(original_file_url)):
original_file_path = os.path.join(out_dir,original_file_name[i])
# Then download the original GiRaFFE 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 GiRaFFE 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
# 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 = 3;
// Set free-parameter values for the initial data.
// Override parameter defaults with values based on command line arguments and NGHOSTS.
const int Nx0x1x2 = 5;
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] = {-1.0,-1.0,-1.0};
const REAL xxmax[3] = { 1.0, 1.0, 1.0};
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;
\n""")
# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
import deprecated_NRPy_param_funcs as evil_par
evil_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("""//""")
GiRaFFEfood_NRPy_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/GiRaFFEfood_NRPy_unit_test.C
// These are common packages that we are likely to need.
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <string> // 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) ) ) )
// Standard formula to calculate significant digits of agreement:
//#define SDA(a,b) 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b)))
REAL SDA(REAL a, REAL b) {
return (a<1.0e-15 && b<1.0e-15) ? 16.0 : 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b)));
}
// Memory access definitions for NRPy+
#define BU0GF 0
#define BU1GF 1
#define BU2GF 2
#define VALENCIAVU0GF 3
#define VALENCIAVU1GF 4
#define VALENCIAVU2GF 5
#define BETAU0GF 6
#define BETAU1GF 7
#define BETAU2GF 8
#define ALPHAGF 9
#define NUM_AUXEVOL_GFS 10
#define AD0GF 0
#define AD1GF 1
#define AD2GF 2
#define NUM_EVOL_GFS 3
// Include the functions that we want to test:
#include "GiRaFFE_NRPy_initial_data.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] = {11,11,11};
const int grid_size = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];
const int grid_size_3_comp = grid_size*3;
CCTK_REAL *Avec = (REAL *)malloc(sizeof(REAL) * grid_size_3_comp);
CCTK_REAL *vel = (REAL *)malloc(sizeof(REAL) * grid_size_3_comp);
CCTK_REAL *Ax = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *Ay = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *Az = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *vx = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *vy = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *vz = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *Bx = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *By = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *Bz = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *x = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *y = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *z = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *r = (REAL *)malloc(sizeof(REAL) * grid_size);;
CCTK_REAL *alp;
CCTK_REAL *betax;
CCTK_REAL *betay;
CCTK_REAL *betaz;
// We need to declare these to compile functions we won't call:
int Compute_Exact_Every;
int cctk_iteration;
CCTK_REAL *delpsi6phi;
CCTK_REAL *psi6phi;
CCTK_REAL *delAx;
CCTK_REAL *delAy;
CCTK_REAL *delAz;
CCTK_REAL *exactBx;
CCTK_REAL *exactBy;
CCTK_REAL *exactBz;
CCTK_REAL *delBx;
CCTK_REAL *delBy;
CCTK_REAL *delBz;
CCTK_REAL *exactVx;
CCTK_REAL *exactVy;
CCTK_REAL *exactVz;
CCTK_REAL *delvx;
CCTK_REAL *delvy;
CCTK_REAL *delvz;
CCTK_REAL cctk_time = 0.0;
CCTK_REAL *exactBx_ThreeWaves;
CCTK_REAL *exactBy_ThreeWaves;
CCTK_REAL *exactBz_ThreeWaves;
CCTK_REAL *delBx_ThreeWaves;
CCTK_REAL *delBy_ThreeWaves;
CCTK_REAL *delBz_ThreeWaves;
CCTK_REAL *mhd_st_x;
CCTK_REAL *mhd_st_y;
CCTK_REAL *mhd_st_z;
CCTK_REAL *Ex;
CCTK_REAL *Ey;
CCTK_REAL *Ez;
CCTK_REAL *B2mE2;
// Set constants to default for comparison
CCTK_REAL wave_speed = -0.5;
CCTK_REAL Omega_aligned_rotator = 1e3;
CCTK_REAL R_NS_aligned_rotator = 1.0;
CCTK_REAL B_p_aligned_rotator = 1e-5;
CCTK_REAL Wald_B0 = 1.0;
// Define dz in CCTK
CCTK_REAL cactus_dxx[3];
#define CCTK_DELTA_SPACE(i) cactus_dxx[i]
// Dummy ETK function:
#define CCTK_GFINDEX3D(cctkGH,i,j,k) (i) + cctk_lsh[0] * ( (j) + cctk_lsh[1] * (k) )
#define CCTK_GFINDEX4D(cctkGH,i,j,k,g) \
( (i) + cctk_lsh[0] * ( (j) + cctk_lsh[1] * ( (k) + cctk_lsh[2] * (g) ) ) )
#define CCTK_VInfo(...) //
//#define CCTK_VWarn(...) //
#include "AlfvenWave.cc"
#include "FastWave.cc"
#include "AlignedRotator.cc"
#include "DegenAlfvenWave.cc"
#include "ThreeAlfvenWave.cc"
#include "FFEBreakdown.cc"
int main(int argc, char** argv) {
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_dxx[0] = dxx0;
cactus_dxx[1] = dxx1;
cactus_dxx[2] = dxx2;
// 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] + (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;
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++)
for(int j=0;j<Nxx_plus_2NGHOSTS1;j++)
for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
int index = CCTK_GFINDEX3D(cctkGH,i,j,k);
x[index] = xx[0][i];
y[index] = xx[1][j];
z[index] = xx[2][k];
r[index] = sqrt(x[index]*x[index] + y[index]*y[index] + z[index]*z[index]);
}
//for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) printf("x[%d] = %.5e\n",j,xx[0][j]);
// 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);
REAL *evol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
GiRaFFE_NRPy_initial_data(¶ms,xx,auxevol_gfs,evol_gfs);
if(atoi(argv[1])==0) GiRaFFEfood_AlfvenWave();
else if(atoi(argv[1])==1) GiRaFFEfood_AlignedRotator();
else if(atoi(argv[1])==2) GiRaFFEfood_FastWave();
else if(atoi(argv[1])==3) GiRaFFEfood_DegenAlfvenWave();
else if(atoi(argv[1])==4) GiRaFFEfood_ThreeAlfvenWave();
else if(atoi(argv[1])==5) GiRaFFEfood_FFEBreakdown();
int all_agree = 1;
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++){
if(SDA(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)])<10.0){
printf("Quantity ValenciavU0 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)]),i0,i1,i2);
all_agree=0;
}
if(SDA(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)])<10.0){
printf("Quantity ValenciavU1 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]),i0,i1,i2);
all_agree=0;
}
if(SDA(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)])<10.0){
printf("Quantity ValenciavU2 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]),i0,i1,i2);
all_agree=0;
}
//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",vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]);
}
}
}
// Shift the grid to compare A_x
//for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] += 0.5*dxx1;
//for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] += 0.5*dxx2;
GiRaFFE_NRPy_initial_data(¶ms,xx,auxevol_gfs,evol_gfs);
if(atoi(argv[1])==0) GiRaFFEfood_AlfvenWave();
else if(atoi(argv[1])==1) GiRaFFEfood_AlignedRotator();
else if(atoi(argv[1])==2) GiRaFFEfood_FastWave();
else if(atoi(argv[1])==3) GiRaFFEfood_DegenAlfvenWave();
else if(atoi(argv[1])==4) GiRaFFEfood_ThreeAlfvenWave();
else if(atoi(argv[1])==5) GiRaFFEfood_FFEBreakdown();
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++){
if(SDA(evol_gfs[IDX4S(AD0GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)])<10.0){
printf("Quantity AD0 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(evol_gfs[IDX4S(AD0GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)]),i0,i1,i2);
all_agree=0;
}
}
}
}
// Shift the grid to compare A_y
//for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] += 0.5*dxx0;
//for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] -= 0.5*dxx1;
GiRaFFE_NRPy_initial_data(¶ms,xx,auxevol_gfs,evol_gfs);
if(atoi(argv[1])==0) GiRaFFEfood_AlfvenWave();
else if(atoi(argv[1])==1) GiRaFFEfood_AlignedRotator();
else if(atoi(argv[1])==2) GiRaFFEfood_FastWave();
else if(atoi(argv[1])==3) GiRaFFEfood_DegenAlfvenWave();
else if(atoi(argv[1])==4) GiRaFFEfood_ThreeAlfvenWave();
else if(atoi(argv[1])==5) GiRaFFEfood_FFEBreakdown();
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++){
if(SDA(evol_gfs[IDX4S(AD1GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)])<10.0){
printf("Quantity AD1 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(evol_gfs[IDX4S(AD1GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]),i0,i1,i2);
all_agree=0;
printf("NRPy: %.15e\n",evol_gfs[IDX4S(AD1GF, i0,i1,i2)]);
printf("CCTK: %.15e\n\n",Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]);
}
}
}
}
// Shift the grid to compare A_z
//for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] += 0.5*dxx1;
//for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] -= 0.5*dxx2;
GiRaFFE_NRPy_initial_data(¶ms,xx,auxevol_gfs,evol_gfs);
if(atoi(argv[1])==0) GiRaFFEfood_AlfvenWave();
else if(atoi(argv[1])==1) GiRaFFEfood_AlignedRotator();
else if(atoi(argv[1])==2) GiRaFFEfood_FastWave();
else if(atoi(argv[1])==3) GiRaFFEfood_DegenAlfvenWave();
else if(atoi(argv[1])==4) GiRaFFEfood_ThreeAlfvenWave();
else if(atoi(argv[1])==5) GiRaFFEfood_FFEBreakdown();
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++){
if(SDA(evol_gfs[IDX4S(AD2GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)])<10.0){
printf("Quantity AD2 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\n",
SDA(evol_gfs[IDX4S(AD2GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]),i0,i1,i2);
all_agree=0;
}
}
}
}
//printf("NRPy: %.15e,%.15e,%.15e\n",evol_gfs[IDX4S(AD0GF, i0,i1,i2)],evol_gfs[IDX4S(AD1GF, i0,i1,i2)],evol_gfs[IDX4S(AD2GF, i0,i1,i2)]);
//printf("CCTK: %.15e,%.15e,%.15e\n",Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]);
if(all_agree) printf("All quantities agree at all points!\n");
}
Writing Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests//GiRaFFEfood_NRPy_unit_test.C
import time
print("Now compiling, should take ~2 seconds...\n")
start = time.time()
# cmd.C_compile(os.path.join(out_dir,"GiRaFFEfood_NRPy_unit_test.C"), os.path.join(out_dir,"GiRaFFEfood_NRPy_unit_test"))
!g++ -Ofast -fopenmp -march=native -funroll-loops Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests/GiRaFFEfood_NRPy_unit_test.C -o Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests/GiRaFFEfood_NRPy_unit_test -lstdc++
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")
results_file = "out_GiRaFFEfood_NRPy_test.txt"
# os.chdir(out_dir)
os.chdir(out_dir)
# cmd.Execute(os.path.join("GiRaFFEfood_NRPy_unit_test"))
if initial_data=="AlfvenWave":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","0",results_file)
elif initial_data=="AlignedRotator":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","1",results_file)
elif initial_data=="FastWave":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","2",results_file)
elif initial_data=="DegenAlfvenWave":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","3",results_file)
elif initial_data=="ThreeWaves":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","4",results_file)
elif initial_data=="FFE_Breakdown":
cmd.Execute("GiRaFFEfood_NRPy_unit_test","5",results_file)
os.chdir(os.path.join("../"))
Now compiling, should take ~2 seconds... Finished in 1.018186092376709 seconds. (EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./GiRaFFEfood_NRPy_unit_test 1`... (BENCH): Finished executing in 0.20 seconds.
Here, we add some emergency brakes so that if the output from the test isn't good, we throw an error to stop the notebook dead in its tracks. This way, our automatic testing infrastructure can let us know if something goes wrong. We will also print the output from the test for convenience's sake.
with open(os.path.join(out_dir,results_file),"r") as file:
output = file.readline()
print(output)
if output!="All quantities agree at all points!\n": # If this isn't the first line of this file, something went wrong!
sys.exit(1)
All quantities agree at all points!
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-GiRaFFEfood_NRPy.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-GiRaFFEfood_NRPy-Flat_Spacetime_Tests")
Created Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy- Flat_Spacetime_Tests.tex, and compiled LaTeX file to PDF file Tutorial- Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Flat_Spacetime_Tests.pdf