GiRaFFE_NRPy
1D tests¶Here we use NRPy+ to generate the C source code necessary to set up initial data for an Alfvén wave (see the original GiRaFFE paper). Then we use it to generate the RHS expressions for Method of Lines time integration based on the explicit Runge-Kutta fourth-order scheme (RK4).
This notebook is organized as follows
GiRaFFEfood_NRPy
initial data modulesfree_parameters.h
GiRaFFE_NRPy_standalone.c
: The Main C Codeimport 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)
# Step P1: Import needed NRPy+ core modules:
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 cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
# Step P2: Create C code output directory:
Ccodesdir = os.path.join("GiRaFFE_staggered_1D_Tests_standalone_Ccodes/")
# 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
# !rm -r ScalarWaveCurvilinear_Playground_Ccodes
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)
# Step P3: Create executable output directory:
outdir = os.path.join(Ccodesdir,"output/")
cmd.mkdir(outdir)
# Step P5: Set timestepping algorithm (we adopt the Method of Lines)
REAL = "double" # Best to use double here.
default_CFL_FACTOR= 0.5 # (GETS OVERWRITTEN WHEN EXECUTED.) In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower.
# Step P6: Set the finite differencing order to 2.
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",4)
thismodule = "Start_to_Finish-GiRaFFE_NRPy-1D_tests"
TINYDOUBLE = par.Cparameters("REAL", thismodule, "TINYDOUBLE", 1e-100)
import GiRaFFE_NRPy.GiRaFFE_NRPy_Main_Driver_staggered as md
# par.set_paramsvals_value("GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C::enforce_speed_limit_StildeD = False")
par.set_paramsvals_value("GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C::enforce_current_sheet_prescription = False")
We will first write the C codes needed for GRFFE evolution. We have already written a module to generate all these codes and call the functions in the appropriate order, so we will import that here. We will take the slightly unusual step of doing this before we generate the initial data functions because the main driver module will register all the gridfunctions we need. It will also generate functions that, in addition to their normal spot in the MoL timestepping, will need to be called during the initial data step to make sure all the variables are appropriately filled in.
All of this is handled with a single call to GiRaFFE_NRPy_Main_Driver_generate_all()
, which will register gridfunctions, write all the C code kernels, and write the C code functions to call those.
md.GiRaFFE_NRPy_Main_Driver_generate_all(Ccodesdir)
Output C function calculate_StildeD0_source_term() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/RHSs/calculate_StildeD0_source_term.h Output C function calculate_StildeD1_source_term() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/RHSs/calculate_StildeD1_source_term.h Output C function calculate_StildeD2_source_term() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/RHSs/calculate_StildeD2_source_term.h Output C function calculate_Stilde_rhsD() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/RHSs/calculate_Stilde_rhsD.h Output C function GiRaFFE_NRPy_cons_to_prims() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/C2P/GiRaFFE_NRPy_cons_to_prims.h Output C function GiRaFFE_NRPy_prims_to_cons() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/C2P/GiRaFFE_NRPy_prims_to_cons.h
RK_method = "RK4"
# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.
# As described above the Table of Contents, this is a 3-step process:
# 3.A: Evaluate RHSs (RHS_string)
# 3.B: Apply boundary conditions (post_RHS_string, pt 1)
import MoLtimestepping.C_Code_Generation as MoL
from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict
RK_order = Butcher_dict[RK_method][1]
cmd.mkdir(os.path.join(Ccodesdir,"MoLtimestepping/"))
MoL.MoL_C_Code_Generation(RK_method,
RHS_string = """
GiRaFFE_NRPy_RHSs(¶ms,auxevol_gfs,RK_INPUT_GFS,RK_OUTPUT_GFS);""",
post_RHS_string = """
GiRaFFE_NRPy_post_step(¶ms,xx,auxevol_gfs,RK_OUTPUT_GFS,n+1);\n""",
outdir = os.path.join(Ccodesdir,"MoLtimestepping/"))
GiRaFFEfood_NRPy
initial data modules [Back to top]¶With the preliminaries out of the way, we will write the C functions to set up initial data. There are two categories of initial data that must be set: the spacetime metric variables, and the GRFFE plasma variables. We will set up the spacetime first, namely the Minkowski spacetime.
gammaDD = ixp.zerorank2(DIM=3)
for i in range(3):
for j in range(3):
if i==j:
gammaDD[i][j] = sp.sympify(1) # else: leave as zero
betaU = ixp.zerorank1() # All should be 0
alpha = sp.sympify(1)
# Description and options for this initial data
desc = "Generate a flat spacetime metric."
loopopts_id ="AllPoints" # we don't need to read coordinates for flat spacetime.
# For testing: Also set inverse metric:
gammaUU, unused_gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
name = "set_initial_spacetime_metric_data"
values_to_print = [
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD00"),rhs=gammaDD[0][0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD01"),rhs=gammaDD[0][1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD02"),rhs=gammaDD[0][2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD11"),rhs=gammaDD[1][1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD12"),rhs=gammaDD[1][2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD22"),rhs=gammaDD[2][2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU0"),rhs=betaU[0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU1"),rhs=betaU[1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU2"),rhs=betaU[2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","alpha"),rhs=alpha)
]
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",values_to_print,params="outCverbose=False"),
loopopts = loopopts_id)
Output C function set_initial_spacetime_metric_data() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/set_initial_spacetime_metric_data.h
Now, we will write out the initials data function for the GRFFE variables.
initial_data_dir = os.path.join(Ccodesdir,"InitialData/")
cmd.mkdir(initial_data_dir)
ID_opts = ["AlfvenWave", "ThreeAlfvenWaves", "DegenAlfvenWave", "FastWave", "FFEBD"]
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gid
for initial_data in ID_opts:
if initial_data=="AlfvenWave":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True)
desc = "Generate Alfven wave 1D initial data for GiRaFFEfood_NRPy."
elif initial_data=="ThreeAlfvenWaves":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = "ThreeWaves", stagger_enable = True)
desc = "Generate three Alfven wave 1D initial data for GiRaFFEfood_NRPy."
elif initial_data=="DegenAlfvenWave":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True)
desc = "Generate degenerate Alfven wave 1D initial data for GiRaFFEfood_NRPy."
elif initial_data=="FastWave":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True)
desc = "Generate fast wave 1D initial data for GiRaFFEfood_NRPy."
elif initial_data=="FFEBD":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = "FFE_Breakdown", stagger_enable = True)
desc = "Generate FFE breakdown 1D initial data for GiRaFFEfood_NRPy."
name = initial_data
values_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","AD0"),rhs=gid.AD[0]),
lhrh(lhs=gri.gfaccess("out_gfs","AD1"),rhs=gid.AD[1]),
lhrh(lhs=gri.gfaccess("out_gfs","AD2"),rhs=gid.AD[2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=gid.ValenciavU[0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=gid.ValenciavU[1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=gid.ValenciavU[2]),
# lhrh(lhs=gri.gfaccess("auxevol_gfs","BU0"),rhs=gid.BU[0]),
# lhrh(lhs=gri.gfaccess("auxevol_gfs","BU1"),rhs=gid.BU[1]),
# lhrh(lhs=gri.gfaccess("auxevol_gfs","BU2"),rhs=gid.BU[2]),
lhrh(lhs=gri.gfaccess("out_gfs","psi6Phi"),rhs=sp.sympify(0))
]
outCfunction(
outfile = os.path.join(initial_data_dir,name+".c"), 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"),
rel_path_to_Cparams='../',
loopopts ="AllPoints,Read_xxs")
inital_data_body = """
const char *option1 = "AlfvenWave";
const char *option2 = "ThreeAlfvenWaves";
const char *option3 = "DegenAlfvenWave";
const char *option4 = "FastWave";
const char *option5 = "FFEBD";
if (strcmp(initial_data_option, option1) == 0) {
AlfvenWave(params, xx, auxevol_gfs, out_gfs);
}
else if (strcmp(initial_data_option, option2) == 0) {
ThreeAlfvenWaves(params, xx, auxevol_gfs, out_gfs);
}
else if (strcmp(initial_data_option, option3) == 0) {
DegenAlfvenWave(params, xx, auxevol_gfs, out_gfs);
}
else if (strcmp(initial_data_option, option4) == 0) {
FastWave(params, xx, auxevol_gfs, out_gfs);
}
else if (strcmp(initial_data_option, option5) == 0) {
FFEBD(params, xx, auxevol_gfs, out_gfs);
}
else {
printf("ERROR: Invalid choice of initial data.");
exit(1);
}
"""
name = "initial_data"
desc = "Main initial data function."
includes = ["AlfvenWave.c", "ThreeAlfvenWaves.c", "DegenAlfvenWave.c", "FastWave.c", "FFEBD.c"]
outCfunction(
outfile = os.path.join(initial_data_dir,name+".h"), desc=desc, name=name,
params ="const char *initial_data_option, const paramstruct *restrict params,REAL *xx[3],REAL *restrict auxevol_gfs,REAL *restrict out_gfs",
body = inital_data_body,
includes = includes,
prefunc="#include <string.h>",
rel_path_to_Cparams='../',
loopopts ="")
Output C function AlfvenWave() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/AlfvenWave.c Output C function ThreeAlfvenWaves() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/ThreeAlfvenWaves.c Output C function DegenAlfvenWave() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/DegenAlfvenWave.c Output C function FastWave() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/FastWave.c Output C function FFEBD() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/FFEBD.c Output C function initial_data() to file GiRaFFE_staggered_1D_Tests_standalone_Ccodes/InitialData/initial_data.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 initial data parameters, as well as grid domain & reference metric parameters, applying domain_size
and sinh_width
/SymTP_bScale
(if applicable) as set above
# Step 3.e: Output C codes needed for declaring and setting Cparameters; also set free_parameters.h
# Step 3.e.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))
# Step 3.e.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.
params.Nxx0 = atoi(argv[1]);
params.Nxx1 = atoi(argv[2]);
params.Nxx2 = atoi(argv[3]);
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.3255,-0.085,-0.085};
const REAL xxmax[3] = { 1.6745, 0.115, 0.115};
params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0+1);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1+1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2+1);
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;
const int poison_grids = 0;
// Standard GRFFE parameters:
params.GAMMA_SPEED_LIMIT = 2000.0;
params.diss_strength = 0.1;
""")
Next apply singular, curvilinear coordinate boundary conditions as documented in the corresponding NRPy+ tutorial notebook
...But, for the moment, we're actually just using this because it writes the file gridfunction_defines.h
.
import CurviBoundaryConditions.CurviBoundaryConditions as cbcs
cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,"boundary_conditions/"),Cparamspath=os.path.join("../"),enable_copy_of_static_Ccodes=False)
Wrote to file "GiRaFFE_staggered_1D_Tests_standalone_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h" Evolved parity: ( AD0:1, AD1:2, AD2:3, StildeD0:1, StildeD1:2, StildeD2:3, psi6Phi:0 ) AuxEvol parity: ( BU0:1, BU1:2, BU2:3, B_lU0:1, B_lU1:2, B_lU2:3, B_rU0:1, B_rU1:2, B_rU2:3, BstaggerU0:1, BstaggerU1:2, BstaggerU2:3, Bstagger_lU0:1, Bstagger_lU1:2, Bstagger_lU2:3, Bstagger_rU0:1, Bstagger_rU1:2, Bstagger_rU2:3, Stilde_flux_HLLED0:1, Stilde_flux_HLLED1:2, Stilde_flux_HLLED2:3, ValenciavU0:1, ValenciavU1:2, ValenciavU2:3, Valenciav_lU0:1, Valenciav_lU1:2, Valenciav_lU2:3, Valenciav_llU0:1, Valenciav_llU1:2, Valenciav_llU2:3, Valenciav_lrU0:1, Valenciav_lrU1:2, Valenciav_lrU2:3, Valenciav_rU0:1, Valenciav_rU1:2, Valenciav_rU2:3, Valenciav_rlU0:1, Valenciav_rlU1:2, Valenciav_rlU2:3, Valenciav_rrU0:1, Valenciav_rrU1:2, Valenciav_rrU2:3, alpha:0, alpha_face:0, betaU0:1, betaU1:2, betaU2:3, beta_faceU0:1, beta_faceU1:2, beta_faceU2:3, cmax_x:0, cmax_y:0, cmax_z:0, cmin_x:0, cmin_y:0, cmin_z:0, gammaDD00:4, gammaDD01:5, gammaDD02:6, gammaDD11:7, gammaDD12:8, gammaDD22:9, gamma_faceDD00:4, gamma_faceDD01:5, gamma_faceDD02:6, gamma_faceDD11:7, gamma_faceDD12:8, gamma_faceDD22:9, psi6_temp:0, psi6center:0 ) Wrote to file "GiRaFFE_staggered_1D_Tests_standalone_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h"
# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),
# and set the CFL_FACTOR (which can be overwritten at the command line)
with open(os.path.join(Ccodesdir,"GiRaFFE_NRPy_REAL__NGHOSTS__CFL_FACTOR.h"), "w") as file:
file.write("""
// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER
#define NGHOSTS """+str(3)+"""
#define NGHOSTS_A2B """+str(2)+"""
// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point
// numbers are stored to at least ~16 significant digits
#define REAL """+REAL+"""
// Part P0.c: Set the CFL Factor. Can be overwritten at command line.
REAL CFL_FACTOR = """+str(default_CFL_FACTOR)+";")
%%writefile $Ccodesdir/GiRaFFE_NRPy_standalone.c
// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.
#include "GiRaFFE_NRPy_REAL__NGHOSTS__CFL_FACTOR.h"
#include "declare_Cparameters_struct.h"
const int NSKIP_1D_OUTPUT = 1;
// Step P1: Import needed header files
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884L
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.707106781186547524400844362104849039L
#endif
// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of
// data in a 1D array. In this case, consecutive values of "i"
// (all other indices held to a fixed value) are consecutive in memory, where
// consecutive values of "j" (fixing all other indices) are separated by
// Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of
// "k" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.
#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) ) ) )
#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)
#define LOOP_ALL_GFS_GPS(ii) _Pragma("omp parallel for") \
for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)
// Step P3: Set gridfunction macros
#include "boundary_conditions/gridfunction_defines.h"
// Step P4: Include the RHS, BC, and primitive recovery functions
#include "GiRaFFE_NRPy_Main_Driver.h"
// Step P5: Include the initial data functions
#include "set_initial_spacetime_metric_data.h"
#include "InitialData/initial_data.h"
// main() function:
// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates
// Step 1: Set up scalar wave initial data
// Step 2: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,
// applying quadratic extrapolation outer boundary conditions.
// Step 3: Output relative error between numerical and exact solution.
// Step 4: Free all allocated memory
int main(int argc, const char *argv[]) {
paramstruct params;
#include "set_Cparameters_default.h"
// Step 0a: Read command-line input, error out if nonconformant
if(argc != 5 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {
printf("Error: Expected three command-line arguments: ./GiRaFFE_NRPy_standalone [Nx] [Ny] [Nz],\n");
printf("where Nx is the number of grid points in the x direction, and so forth.\n");
printf("Nx,Ny,Nz MUST BE larger than NGHOSTS (= %d)\n",NGHOSTS);
exit(1);
}
// 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"
// ... and then set up the numerical grid structure in time:
const REAL t_final = 2.0;
const REAL CFL_FACTOR = 0.5; // Set the CFL Factor
// Step 0c: Allocate memory for gridfunctions
const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
// Step 0k: Allocate memory for gridfunctions
#include "MoLtimestepping/RK_Allocate_Memory.h"
REAL *restrict auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);
REAL *evol_gfs_exact = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot);
REAL *auxevol_gfs_exact = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);
// For debugging, it can be useful to set everything to NaN initially.
if(poison_grids) {
for(int ii=0;ii<NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
y_n_gfs[ii] = 1.0/0.0;
y_nplus1_running_total_gfs[ii] = 1.0/0.0;
//k_odd_gfs[ii] = 1.0/0.0;
//k_even_gfs[ii] = 1.0/0.0;
diagnostic_output_gfs[ii] = 1.0/0.0;
evol_gfs_exact[ii] = 1.0/0.0;
}
for(int ii=0;ii<NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
auxevol_gfs[ii] = 1.0/0.0;
auxevol_gfs_exact[ii] = 1.0/0.0;
}
}
// Step 0d: Set up coordinates: Set dx, and then dt based on dx_min and CFL condition
// This is probably already defined above, but just in case...
#ifndef MIN
#define MIN(A, B) ( ((A) < (B)) ? (A) : (B) )
#endif
REAL dt = CFL_FACTOR * MIN(dxx0,MIN(dxx1,dxx2)); // CFL condition
int Nt = (int)(t_final / dt + 0.5); // The number of points in time.
//Add 0.5 to account for C rounding down integers.
// 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+1)*dxx0;
for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS+1)*dxx1;
for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS+1)*dxx2;
// Step 1: Set up initial data to be exact solution at time=0:
REAL time = 0.0;
set_initial_spacetime_metric_data(¶ms, xx, auxevol_gfs);
const char *initial_data_option = argv[4];
initial_data(initial_data_option, ¶ms, xx, auxevol_gfs, y_n_gfs);
// Fill in the remaining quantities
GiRaFFE_compute_B_and_Bstagger_from_A(¶ms,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD00GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD01GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD02GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD11GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD12GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD22GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*PSI6_TEMPGF, /* Temporary storage,overwritten */
y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD0GF,
y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD1GF,
y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD2GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU0GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU1GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU2GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU0GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU1GF,
auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU2GF);
//override_BU_with_old_GiRaFFE(¶ms,auxevol_gfs,0);
GiRaFFE_NRPy_prims_to_cons(¶ms,auxevol_gfs,y_n_gfs);
// Extra stack, useful for debugging:
GiRaFFE_NRPy_cons_to_prims(¶ms,xx,auxevol_gfs,y_n_gfs);
for(int n=0;n<=Nt;n++) { // Main loop to progress forward in time.
//for(int n=0;n<=1;n++) { // Main loop to progress forward in time.
// Step 1a: Set current time to correct value & compute exact solution
time = ((REAL)n)*dt;
/* Step 2: Validation: Output relative error between numerical and exact solution, */
if(time == 0.0 || time == 0.5 || time == 1.0 || time == 2.0 || time == 0.02 || time == 0.56) {
// Step 2c: Output relative error between exact & numerical at center of grid.
const int i0mid=Nxx_plus_2NGHOSTS0/2;
const int i1mid=Nxx_plus_2NGHOSTS1/2;
const int i2mid=Nxx_plus_2NGHOSTS2/2;
char filename[100];
sprintf(filename,"out%d__%s-%08d.txt", Nxx0, initial_data_option, n);
FILE *out2D = fopen(filename, "w");
for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
const int idx = IDX3S(i0,i1mid,i2mid);
fprintf(out2D,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",
xx[0][i0],
auxevol_gfs[IDX4ptS(BU0GF,idx)],auxevol_gfs[IDX4ptS(BU1GF,idx)],auxevol_gfs[IDX4ptS(BU2GF,idx)],
y_n_gfs[IDX4ptS(AD0GF,idx)],y_n_gfs[IDX4ptS(AD1GF,idx)],y_n_gfs[IDX4ptS(AD2GF,idx)],
y_n_gfs[IDX4ptS(STILDED0GF,idx)],y_n_gfs[IDX4ptS(STILDED1GF,idx)],y_n_gfs[IDX4ptS(STILDED2GF,idx)],
auxevol_gfs[IDX4ptS(VALENCIAVU0GF,idx)],auxevol_gfs[IDX4ptS(VALENCIAVU1GF,idx)],auxevol_gfs[IDX4ptS(VALENCIAVU2GF,idx)],
y_n_gfs[IDX4ptS(PSI6PHIGF,idx)], time);
}
fclose(out2D);
// For convergence testing, we'll shift the grid x -> x-1 and output initial data again, giving the exact solution.
LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,1,0,1) {
xx[0][i0] += -mu_AW*time;
//xx[0][i0] += -time;
}
set_initial_spacetime_metric_data(¶ms,xx,auxevol_gfs_exact);
initial_data(initial_data_option, ¶ms,xx,auxevol_gfs_exact,evol_gfs_exact);
// Fill in the remaining quantities
//driver_A_to_B(¶ms,evol_gfs_exact,auxevol_gfs_exact);
GiRaFFE_NRPy_prims_to_cons(¶ms,auxevol_gfs_exact,evol_gfs_exact);
// And now, we'll set the grid back to rights.
LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,1,0,1) {
xx[0][i0] -= -mu_AW*time;
//xx[0][i0] -= -time;
}
sprintf(filename,"out%d-%08d_exact.txt",Nxx0,n);
FILE *out2D_exact = fopen(filename, "w");
for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
const int idx = IDX3S(i0,i1mid,i2mid);
fprintf(out2D_exact,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",
xx[0][i0],
auxevol_gfs_exact[IDX4ptS(BU0GF,idx)],auxevol_gfs_exact[IDX4ptS(BU1GF,idx)],auxevol_gfs_exact[IDX4ptS(BU2GF,idx)],
evol_gfs_exact[IDX4ptS(AD0GF,idx)],evol_gfs_exact[IDX4ptS(AD1GF,idx)],evol_gfs_exact[IDX4ptS(AD2GF,idx)],
evol_gfs_exact[IDX4ptS(STILDED0GF,idx)],evol_gfs_exact[IDX4ptS(STILDED1GF,idx)],evol_gfs_exact[IDX4ptS(STILDED2GF,idx)],
auxevol_gfs_exact[IDX4ptS(VALENCIAVU0GF,idx)],auxevol_gfs_exact[IDX4ptS(VALENCIAVU1GF,idx)],auxevol_gfs_exact[IDX4ptS(VALENCIAVU2GF,idx)],
evol_gfs_exact[IDX4ptS(PSI6PHIGF,idx)]);
}
fclose(out2D_exact);
}
// Step 3: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,
// applying quadratic extrapolation outer boundary conditions.
// Step 3.b: Step forward one timestep (t -> t+dt) in time using
// chosen RK-like MoL timestepping algorithm
#include "MoLtimestepping/RK_MoL.h"
} // End main loop to progress forward in time.
// Step 4: Free all allocated memory
#include "MoLtimestepping/RK_Free_Memory.h"
free(auxevol_gfs);
free(auxevol_gfs_exact);
free(evol_gfs_exact);
for(int i=0;i<3;i++) free(xx[i]);
return 0;
}
Writing GiRaFFE_staggered_1D_Tests_standalone_Ccodes//GiRaFFE_NRPy_standalone.c
To aid in the cross-platform-compatible (with Windows, MacOS, & Linux) compilation and execution, we make use of cmdline_helper
(Tutorial).
cmd.C_compile(os.path.join(Ccodesdir,"GiRaFFE_NRPy_standalone.c"),
os.path.join(Ccodesdir,"output","GiRaFFE_NRPy_standalone"),compile_mode="optimized")
# Change to output directory
os.chdir(outdir)
# Clean up existing output files
cmd.delete_existing_files("out*.txt")
cmd.delete_existing_files("out*.png")
# ID options are: "AlfvenWave", "ThreeAlfvenWaves", "DegenAlfvenWave", "FastWave", "FFEBD"
for opt in ID_opts:
cmd.Execute("GiRaFFE_NRPy_standalone", "299 4 4 "+opt, "out_298"+opt+".txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "1280 9 9 "+opt, "out_1280"+opt+".txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "1280 32 32 "+opt, "out_"+opt+".txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "149 9 9 AlfvenWave","out149.txt")
# Return to root directory
os.chdir(os.path.join("../../"))
Compiling executable... (EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops GiRaFFE_staggered_1D_Tests_standalone_Ccodes/GiRaFFE_NRPy_standalone.c -o GiRaFFE_staggered_1D_Tests_standalone_Ccodes/output/GiRaFFE_NRPy_standalone -lm`... (BENCH): Finished executing in 3.215022563934326 seconds. Finished compilation. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 299 4 4 AlfvenWave`... (BENCH): Finished executing in 17.259787559509277 seconds. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 299 4 4 ThreeAlfvenWaves`... (BENCH): Finished executing in 17.4488844871521 seconds. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 299 4 4 DegenAlfvenWave`... (BENCH): Finished executing in 17.24630904197693 seconds. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 299 4 4 FastWave`... (BENCH): Finished executing in 17.657445907592773 seconds. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 299 4 4 FFEBD`... (BENCH): Finished executing in 18.245166063308716 seconds.
Now we plot the data and recreate figure 1 from the GiRaFFE paper. We reconstruct the electric field via
$$ E_i = -\epsilon_{ijk}v^j B^k $$the calc_E
function below. We also calculate the FFE condition $B^2 - E^2$ below using the calc_Bsquared_minus_Esquared
function.
eDDD = ixp.LeviCivitaSymbol_dim3_rank3()
def calc_E(data):
VU0 = data[:, 10]
VU1 = data[:, 11]
VU2 = data[:, 12]
BU0 = data[:, 1]
BU1 = data[:, 2]
BU2 = data[:, 3]
VU = [VU0, VU1, VU2]
BU = [BU0, BU1, BU2]
ED = np.zeros((VU0.size, 3))
for i in range(3):
for j in range(3):
for k in range(3):
ED[:,i] = ED[:,i] - eDDD[i][j][k]*VU[j]*BU[k]
return ED
def calc_Bsquared_minus_Esquared(data):
EU = calc_E(data)
BU0 = data[:, 1]
BU1 = data[:, 2]
BU2 = data[:, 3]
return (BU0**2 + BU1**2 + BU2**2) - (EU[:,0]**2 + EU[:,1]**2 + EU[:,2]**2)
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
# 0 1 2 3 4 5 6 7 8 9 10 11 12 13
labels = ["x","BU0","BU1","BU2","AD0","AD1","AD2","StildeD0","StildeD1","StildeD2","ValenciavU0","ValenciavU1","ValenciavU2", "psi6Phi"]
fig = plt.figure(figsize=(6, 15))
# spec = mpl.gridspec.GridSpec(ncols=6, nrows=2,wspace=0.65, hspace=0.4) # 6 columns evenly divides both 2 & 3
# ax1 = fig.add_subplot(spec[0,0:2]) # row 0 with axes spanning 2 cols on evens
# ax2 = fig.add_subplot(spec[0,2:4])
# ax3 = fig.add_subplot(spec[0,4:])
# ax4 = fig.add_subplot(spec[1,1:3]) # row 0 with axes spanning 2 cols on odds
# ax5 = fig.add_subplot(spec[1,3:5])
gs = gridspec.GridSpec(nrows=5, ncols=1, hspace=0.5)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[2, 0])
ax4 = fig.add_subplot(gs[3, 0])
ax5 = fig.add_subplot(gs[4, 0])
Data_num_Fast_A = np.loadtxt(os.path.join(Ccodesdir,"output","out299__FastWave-00000000.txt"))
Data_num_Fast_B = np.loadtxt(os.path.join(Ccodesdir,"output","out299__FastWave-00000100.txt"))
E_Fast_A = calc_E(Data_num_Fast_A)
E_Fast_B = calc_E(Data_num_Fast_B)
ax1.scatter(Data_num_Fast_A[:,0], np.abs(E_Fast_A[:,2]), s=1,label = 't = 0')
ax1.plot(Data_num_Fast_B[:,0], np.abs(E_Fast_B[:,2]), 'k-', label = 't = 0.5')
ax1.set_xlim(-0.5, 1.5)
ax1.set_ylim(0.6)
ax1.text(0.95, 0.01, 'Fast Wave',
verticalalignment='bottom', horizontalalignment='right',
transform=ax1.transAxes,
color='black', fontsize=14)
ax1.set_xlabel('x')
ax1.set_ylabel(r'$|E^z|$')
ax1.legend()
Data_num_Alf_A = np.loadtxt(os.path.join(Ccodesdir,"output","out299__AlfvenWave-00000000.txt"))
Data_num_Alf_B = np.loadtxt(os.path.join(Ccodesdir,"output","out299__AlfvenWave-00000400.txt"))
ax2.scatter(Data_num_Alf_A[:,0], Data_num_Alf_A[:,3], s=1, label = 't = 0')
ax2.plot(Data_num_Alf_B[:,0], Data_num_Alf_B[:,3], 'k-', label = 't = 2.0')
ax2.set_xlim(-1.5, 1.5)
ax2.set_ylim(1.1)
ax2.text(0.95, 0.01, 'Alfven Wave',
verticalalignment='bottom', horizontalalignment='right',
transform=ax2.transAxes,
color='black', fontsize=14)
ax2.set_xlabel('x')
ax2.set_ylabel(r'$B^z$')
ax2.legend(loc='center right')
Data_num_DegenAlf_A = np.loadtxt(os.path.join(Ccodesdir,"output","out299__DegenAlfvenWave-00000000.txt"))
Data_num_DegenAlf_B = np.loadtxt(os.path.join(Ccodesdir,"output","out299__DegenAlfvenWave-00000200.txt"))
E_DegenAlf_A = calc_E(Data_num_DegenAlf_A)
E_DegenAlf_B = calc_E(Data_num_DegenAlf_B)
ax3.scatter(Data_num_DegenAlf_A[:,0], E_DegenAlf_A[:,1], s=1, label = 't = 0')
ax3.plot(Data_num_DegenAlf_B[:,0], E_DegenAlf_B[:,1], 'k-', label = 't = 1.0')
ax3.set_xlim(-1.5, 1.5)
ax3.set_ylim(-1.35)
ax3.text(0.95, 0.01, 'Deg. Alfven Wave',
verticalalignment='bottom', horizontalalignment='right',
transform=ax3.transAxes,
color='black', fontsize=14)
ax3.set_xlabel('x')
ax3.set_ylabel(r'$E^y$')
ax3.legend()
# Data_num_ThreeAlf_A = np.loadtxt(os.path.join(Ccodesdir,"output","out149__ThreeAlfvenWaves-00000000.txt"))
Data_num_ThreeAlf_B = np.loadtxt(os.path.join(Ccodesdir,"output","out299__ThreeAlfvenWaves-00000112.txt"))
# ax2.plot(Data_num_ThreeAlf_A[:,0], Data_num_ThreeAlf_A[:,2], 'k-')
ax4.scatter(Data_num_ThreeAlf_B[:,0], Data_num_ThreeAlf_B[:,2], s=1, label = 't = 0.56')
ax4.set_xlim(-1.0, 1.0)
# ax4.set_ylim()
ax4.text(0.95, 0.01, 'Three Waves',
verticalalignment='bottom', horizontalalignment='right',
transform=ax4.transAxes,
color='black', fontsize=14)
ax4.set_xlabel('x')
ax4.set_ylabel(r'$B^y$')
ax4.legend(loc='center')
Data_num_FFEBD_A = np.loadtxt(os.path.join(Ccodesdir,"output","out299__FFEBD-00000000.txt"))
Data_num_FFEBD_B = np.loadtxt(os.path.join(Ccodesdir,"output","out299__FFEBD-00000004.txt"))
B2mE2_A = calc_Bsquared_minus_Esquared(Data_num_FFEBD_A)
B2mE2_B = calc_Bsquared_minus_Esquared(Data_num_FFEBD_B)
ax5.scatter(Data_num_FFEBD_A[:,0], B2mE2_A, s=1, label = 't = 0')
ax5.plot(Data_num_FFEBD_B[:,0], B2mE2_B, 'k-', label = 't = 0.02')
ax5.set_xlim(-0.4, 0.6)
ax5.text(0.95, 0.01, 'FFE Breakdown',
verticalalignment='bottom', horizontalalignment='right',
transform=ax5.transAxes,
color='black', fontsize=14)
ax5.set_xlabel('x')
ax5.set_ylabel(r'$B^2 - E^2$')
ax5.legend()
plt.savefig(os.path.join(Ccodesdir,"output","NRPy-GiRaFFE"), dpi=800, bbox_inches="tight")
plt.close(fig)
img1 = plt.imread(os.path.join(Ccodesdir,"output","NRPy-GiRaFFE.png"))
img2 = plt.imread(os.path.join("GiRaFFE_NRPy/example_par_files/figure1_GiRaFFE_paper.png"))
NUM_ROWS = 1
IMGs_IN_ROW = 2
f, ax = plt.subplots(NUM_ROWS, IMGs_IN_ROW, figsize=(28,18))
plt.subplots_adjust(wspace=0.05)
plt.axis('off')
ax[0].imshow(img1)
ax[1].imshow(img2)
ax[0].set_title('image 1')
ax[1].set_title('image 2')
# title = 'side by side view of images'
# f.suptitle(title, fontsize=16)
plt.tight_layout()
# plt.xticks([])
# plt.yticks([])
plt.show()
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-GiRaFFE_NRPy-1D_tests-staggered.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-GiRaFFE_NRPy-1D_tests-staggered",location_of_template_file=os.path.join(".."))
Created Tutorial-Start_to_Finish-GiRaFFE_NRPy-1D_tests-staggered.tex, and compiled LaTeX file to PDF file Tutorial-Start_to_Finish- GiRaFFE_NRPy-1D_tests-staggered.pdf