GiRaFFE_NRPy
3D tests¶Here we use NRPy+ to generate the C source code necessary to set up initial data for an Exact Wald (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, add_to_Cfunction_dict # 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_new_way_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",2)
# Step P7: Enable SIMD-optimized code?
# I.e., generate BSSN and Ricci C code kernels using SIMD-vectorized
# compiler intrinsics, which *greatly improve the code's performance*,
# though at the expense of making the C-code kernels less
# human-readable.
# * Important note in case you wish to modify the BSSN/Ricci kernels
# here by adding expressions containing transcendental functions
# (e.g., certain scalar fields):
# Note that SIMD-based transcendental function intrinsics are not
# supported by the default installation of gcc or clang (you will
# need to use e.g., the SLEEF library from sleef.org, for this
# purpose). The Intel compiler suite does support these intrinsics
# however without the need for external libraries.
enable_SIMD = False
# Step 1.b: Enable reference metric precomputation.
enable_rfm_precompute = False
if enable_SIMD and not enable_rfm_precompute:
print("ERROR: SIMD does not currently handle transcendental functions,\n")
print(" like those found in rfmstruct (rfm_precompute).\n")
print(" Therefore, enable_SIMD==True and enable_rfm_precompute==False\n")
print(" is not supported.\n")
sys.exit(1)
# Step 1.c: Enable "FD functions". In other words, all finite-difference stencils
# will be output as inlined static functions. This is essential for
# compiling highly complex FD kernels with using certain versions of GCC;
# GCC 10-ish will choke on BSSN FD kernels at high FD order, sometimes
# taking *hours* to compile. Unaffected GCC versions compile these kernels
# in seconds. FD functions do not slow the code performance, but do add
# another header file to the C source tree.
# With gcc 7.5.0, enable_FD_functions=True decreases performance by 10%
enable_FD_functions = False
thismodule = "Start_to_Finish-GiRaFFE_NRPy-1D_tests"
TINYDOUBLE = par.Cparameters("REAL", thismodule, "TINYDOUBLE", 1e-100)
import GiRaFFE_NRPy.GiRaFFE_NRPy_Main_Driver_staggered_new_way 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.
RK_method = "Euler"
# 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.
# There are several initial data routines we need to test. We'll control which one we use with a string option
initial_data = "ExactWald" # Valid options: "ExactWald", "AlignedRotator"
spacetime = "ShiftedKerrSchild" # Valid options: "ShiftedKerrSchild", "flat"
if spacetime == "ShiftedKerrSchild":
import BSSN.ShiftedKerrSchild as sks
sks.ShiftedKerrSchild(True)
import reference_metric as rfm
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
# Use the Jacobian matrix to transform the vectors to Cartesian coordinates.
par.set_parval_from_str("reference_metric::CoordSystem","Spherical")
rfm.reference_metric()
Jac_dUCart_dDrfmUD,Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()
# Transform the coordinates of the Jacobian matrix from spherical to Cartesian:
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
tmpa,tmpb,tmpc = sp.symbols("tmpa,tmpb,tmpc")
for i in range(3):
for j in range(3):
Jac_dUCart_dDrfmUD[i][j] = Jac_dUCart_dDrfmUD[i][j].subs([(rfm.xx[0],tmpa),(rfm.xx[1],tmpb),(rfm.xx[2],tmpc)])
Jac_dUCart_dDrfmUD[i][j] = Jac_dUCart_dDrfmUD[i][j].subs([(tmpa,rfm.xxSph[0]),(tmpb,rfm.xxSph[1]),(tmpc,rfm.xxSph[2])])
Jac_dUrfm_dDCartUD[i][j] = Jac_dUrfm_dDCartUD[i][j].subs([(rfm.xx[0],tmpa),(rfm.xx[1],tmpb),(rfm.xx[2],tmpc)])
Jac_dUrfm_dDCartUD[i][j] = Jac_dUrfm_dDCartUD[i][j].subs([(tmpa,rfm.xxSph[0]),(tmpb,rfm.xxSph[1]),(tmpc,rfm.xxSph[2])])
gammaSphDD = ixp.zerorank2()
for i in range(3):
for j in range(3):
gammaSphDD[i][j] += sks.gammaSphDD[i][j].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
betaSphU = ixp.zerorank1()
for i in range(3):
betaSphU[i] += sks.betaSphU[i].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
alpha = sks.alphaSph.subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
gammaDD = rfm.basis_transform_tensorDD_from_rfmbasis_to_Cartesian(Jac_dUrfm_dDCartUD, gammaSphDD)
unused_gammaUU,gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)
sqrtgammaDET = sp.sqrt(gammaDET)
betaU = rfm.basis_transform_vectorD_from_rfmbasis_to_Cartesian(Jac_dUrfm_dDCartUD, betaSphU)
loopopts_id ="AllPoints,Read_xxs"
elif spacetime == "flat":
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.
desc = "Generate a spinning black hole with Shifted Kerr Schild metric in a Cartesian basis."
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_new_way_standalone_Ccodes/set_initial_spacetime_metric_data.h
Now, we will write out the initial data function for the GRFFE variables.
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gid
if initial_data=="ExactWald":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True,M=sks.M,KerrSchild_radial_shift=sks.r0,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaDET)
desc = "Generate exact Wald initial test data for GiRaFFEfood_NRPy."
elif initial_data=="SplitMonopole":
gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True,M=sks.M,a=sks.a,KerrSchild_radial_shift=sks.r0,alpha=alpha,betaU=betaSphU,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaSphDET)
desc = "Generate Split Monopole initial test data for GiRaFFEfood_NRPy."
else:
print("Unsupported Initial Data string "+initial_data+"! Supported ID: ExactWald, or SplitMonopole")
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(Ccodesdir,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 initial_data() to file GiRaFFE_staggered_new_way_standalone_Ccodes/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
import deprecated_NRPy_param_funcs as evil_par
evil_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.5,-1.5,-1.5};
const REAL xxmax[3] = { 1.5, 1.5, 1.5};
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 = %.15e,%.15e,%.15e\\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 = 1;
// Standard GRFFE parameters:
params.GAMMA_SPEED_LIMIT = 2000.0;
params.diss_strength = 4.0;
""")
if initial_data=="ExactWald":
with open(os.path.join(Ccodesdir,"free_parameters.h"),"a") as file:
file.write("""params.r0 = 0.4;
params.a = 0.0;
""")
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(Ccodesdir,enable_copy_of_static_Ccodes=False)
Wrote to file "GiRaFFE_staggered_new_way_standalone_Ccodes/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_new_way_standalone_Ccodes/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)+";")
#include "GiRaFFE_NRPy_REAL__NGHOSTS__CFL_FACTOR.h"
#include "declare_Cparameters_struct.h"
def add_to_Cfunction_dict_main__GiRaFFE_NRPy_3D_tests_staggered():
includes = ["NRPy_basic_defines.h", "GiRaFFE_main_defines.h", "NRPy_function_prototypes.h", "time.h", "set_initial_spacetime_metric_data.h", "initial_data.h"]
desc = """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
"""
prefunc = """const int NSKIP_1D_OUTPUT = 1;
const REAL outer_integration_radius = 1.5; const REAL inner_integration_radius = 0.6;
#define LOOP_ALL_GFS_GPS(ii) _Pragma("omp parallel for") \
for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)
"""
c_type = "int"
name = "main"
params = "int argc, const char *argv[]"
body = """
paramstruct params;
#include "set_Cparameters_default.h"
// Step 0a: Read command-line input, error out if nonconformant
if(argc != 4 || 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);
initial_data(¶ms,xx,auxevol_gfs,y_n_gfs);
// Code to perturb the initial data:
/*
for(int ii=0;ii<NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
y_n_gfs[ii] += 1.0e-15;
y_nplus1_running_total_gfs[ii] += 1.0e-15;
//k_odd_gfs[ii] = 1.0/0.0;
//k_even_gfs[ii] = 1.0/0.0;
diagnostic_output_gfs[ii] += 1.0e-15;
evol_gfs_exact[ii] += 1.0e-15;
}
for(int ii=0;ii<NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
auxevol_gfs[ii] += 1.0e-15;
auxevol_gfs_exact[ii] += 1.0e-15;
}
*/
// 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);
GiRaFFE_NRPy_prims_to_cons(¶ms,auxevol_gfs,y_n_gfs);
GiRaFFE_NRPy_cons_to_prims(¶ms,xx,auxevol_gfs,y_n_gfs);
GiRaFFE_NRPy_prims_to_cons(¶ms,auxevol_gfs,y_n_gfs);
GiRaFFE_NRPy_cons_to_prims(¶ms,xx,auxevol_gfs,y_n_gfs);
GiRaFFE_NRPy_prims_to_cons(¶ms,auxevol_gfs,y_n_gfs);
//GiRaFFE_NRPy_cons_to_prims(¶ms,xx,auxevol_gfs,y_n_gfs);
// Fill in the exact gridfunctions the same way as above: In Exact Wald, initial data IS exact solution!
for(int ii=0;ii<NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
evol_gfs_exact[ii] = y_n_gfs[ii];
}
for(int ii=0;ii<NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
auxevol_gfs_exact[ii] = auxevol_gfs[ii];
}
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((n)%NSKIP_1D_OUTPUT ==0) {
// Step 2c: Output quantities along x axis.
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-%08d.txt",Nxx0,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\\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(ALPHAGF,idx)]*auxevol_gfs[IDX4ptS(VALENCIAVU0GF,idx)]-auxevol_gfs[IDX4ptS(BETAU0GF,idx)],
auxevol_gfs[IDX4ptS(ALPHAGF,idx)]*auxevol_gfs[IDX4ptS(VALENCIAVU1GF,idx)]-auxevol_gfs[IDX4ptS(BETAU1GF,idx)],
auxevol_gfs[IDX4ptS(ALPHAGF,idx)]*auxevol_gfs[IDX4ptS(VALENCIAVU2GF,idx)]-auxevol_gfs[IDX4ptS(BETAU2GF,idx)],
y_n_gfs[IDX4ptS(PSI6PHIGF,idx)]);
}
fclose(out2D);
// Compute & output L2 norm. For each quantity Q we care about:
// \\sqrt{\\frac{1}{N} \\sum_{ijk} \\left( Q_{\\rm approx} - Q_{\\rm exact} \\right)^2}
int num_quantities_integrated = 13;
REAL integrated_quantities[num_quantities_integrated];
for(int ii=0; ii<num_quantities_integrated; ii++) {
integrated_quantities[ii] = 0.0;
}
sprintf(filename,"out%d-L2.txt",Nxx0);
FILE *outL2 = fopen(filename, "a");
int num_points = 0;
#ifndef SQR
#define SQR(x) ((x) * (x))
#endif /*SQR*/
for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++) {
REAL xx2 = xx[2][i2];
for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) {
REAL xx1 = xx[1][i1];
for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
REAL xx0 = xx[0][i0];
// Make sure we're inside a thick spherical shell:
if(SQR(xx0)+SQR(xx1)+SQR(xx2) < SQR(outer_integration_radius) &&
SQR(xx0)+SQR(xx1)+SQR(xx2) > SQR(inner_integration_radius)) {
const int idx = IDX3S(i0,i1,i2);
// Compute the sum of the squares of the difference between numerical and exact:
integrated_quantities[0] += SQR(auxevol_gfs[IDX4ptS(BU0GF,idx)] - auxevol_gfs_exact[IDX4ptS(BU0GF,idx)]);
integrated_quantities[1] += SQR(auxevol_gfs[IDX4ptS(BU1GF,idx)] - auxevol_gfs_exact[IDX4ptS(BU1GF,idx)]);
integrated_quantities[2] += SQR(auxevol_gfs[IDX4ptS(BU2GF,idx)] - auxevol_gfs_exact[IDX4ptS(BU2GF,idx)]);
integrated_quantities[3] += SQR(y_n_gfs[IDX4ptS(AD0GF,idx)] - evol_gfs_exact[IDX4ptS(AD0GF,idx)]);
integrated_quantities[4] += SQR(y_n_gfs[IDX4ptS(AD1GF,idx)] - evol_gfs_exact[IDX4ptS(AD1GF,idx)]);
integrated_quantities[5] += SQR(y_n_gfs[IDX4ptS(AD2GF,idx)] - evol_gfs_exact[IDX4ptS(AD2GF,idx)]);
integrated_quantities[6] += SQR(y_n_gfs[IDX4ptS(STILDED0GF,idx)] - evol_gfs_exact[IDX4ptS(STILDED0GF,idx)]);
integrated_quantities[7] += SQR(y_n_gfs[IDX4ptS(STILDED1GF,idx)] - evol_gfs_exact[IDX4ptS(STILDED1GF,idx)]);
integrated_quantities[8] += SQR(y_n_gfs[IDX4ptS(STILDED2GF,idx)] - evol_gfs_exact[IDX4ptS(STILDED2GF,idx)]);
integrated_quantities[9] += SQR(auxevol_gfs[IDX4ptS(VALENCIAVU0GF,idx)] - auxevol_gfs_exact[IDX4ptS(VALENCIAVU0GF,idx)]);
integrated_quantities[10] += SQR(auxevol_gfs[IDX4ptS(VALENCIAVU1GF,idx)] - auxevol_gfs_exact[IDX4ptS(VALENCIAVU1GF,idx)]);
integrated_quantities[11] += SQR(auxevol_gfs[IDX4ptS(VALENCIAVU2GF,idx)] - auxevol_gfs_exact[IDX4ptS(VALENCIAVU2GF,idx)]);
integrated_quantities[12] += SQR(y_n_gfs[IDX4ptS(PSI6PHIGF,idx)] - evol_gfs_exact[IDX4ptS(PSI6PHIGF,idx)]);
num_points++;
}
}
}
}
// Divide by N, take the square root:
for(int ii=0; ii<num_quantities_integrated; ii++) {
integrated_quantities[ii] = sqrt(integrated_quantities[ii]/num_points);
}
fprintf(outL2,"%d %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\\n",
n,
integrated_quantities[0],integrated_quantities[1],integrated_quantities[2],
integrated_quantities[3],integrated_quantities[4],integrated_quantities[5],
integrated_quantities[6],integrated_quantities[7],integrated_quantities[8],
integrated_quantities[9],integrated_quantities[10],integrated_quantities[11],
integrated_quantities[12]);
fclose(outL2);
}
// 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;
"""
add_to_Cfunction_dict(
includes=includes,
desc=desc,
c_type=c_type, name=name, params=params,
prefunc = prefunc, body=body,
rel_path_to_Cparams=os.path.join("."), enableCparameters=False)
import GiRaFFE_NRPy.GiRaFFE_NRPy_staggered_Source_Terms as stgsrc
stgsrc.add_to_Cfunction_dict__GiRaFFE_NRPy_staggered_Source_Terms()
md.add_to_Cfunction_dict__cons_to_prims(md.StildeD,md.BU,md.gammaDD,md.betaU,md.alpha,includes=["NRPy_basic_defines.h","GiRaFFE_basic_defines.h"])
md.add_to_Cfunction_dict__prims_to_cons(md.gammaDD,md.betaU,md.alpha,md.ValenciavU,md.BU,md.sqrt4pi,includes=["NRPy_basic_defines.h","GiRaFFE_basic_defines.h"])
import GiRaFFE_NRPy.GiRaFFE_NRPy_Source_Terms as source
source.add_to_Cfunction_dict__functions_for_StildeD_source_term(md.outCparams,md.gammaDD,md.betaU,md.alpha,
md.ValenciavU,md.BU,md.sqrt4pi,includes=["NRPy_basic_defines.h","GiRaFFE_basic_defines.h"])
import GiRaFFE_NRPy.Stilde_flux as Sf
Sf.add_to_Cfunction_dict__Stilde_flux(includes=["NRPy_basic_defines.h","GiRaFFE_basic_defines.h"], inputs_provided = True, alpha_face=md.alpha_face, gamma_faceDD=md.gamma_faceDD,
beta_faceU=md.beta_faceU, Valenciav_rU=md.Valenciav_rU, B_rU=md.B_rU,
Valenciav_lU=md.Valenciav_lU, B_lU=md.B_lU, sqrt4pi=md.sqrt4pi,write_cmax_cmin=True)
import GiRaFFE_NRPy.GiRaFFE_NRPy_staggered_Afield_flux as Af
Af.add_to_Cfunction_dict__GiRaFFE_NRPy_staggered_Afield_flux()
import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL
FCVAL.add_to_Cfunction_dict__GiRaFFE_NRPy_FCVAL(includes=["NRPy_basic_defines.h","GiRaFFE_basic_defines.h"])
import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM
PPM.add_to_Cfunction_dict__GiRaFFE_NRPy_PPM(Ccodesdir)
import GiRaFFE_NRPy.GiRaFFE_NRPy_staggered_A2B as A2B
A2B.add_to_Cfunction_dict__GiRaFFE_NRPy_staggered_A2B()
import GiRaFFE_NRPy.GiRaFFE_NRPy_BCs as BC
BC.add_to_Cfunction_dict__GiRaFFE_NRPy_BCs()
md.add_to_Cfunction_dict__driver_function()
add_to_Cfunction_dict_main__GiRaFFE_NRPy_3D_tests_staggered()
Now, we will register the remaining C functions and contributions to NRPy_basic_defines.h
, then we output NRPy_basic_defines.h
and NRPy_function_prototypes.h
.
import outputC as outC
outC.outputC_register_C_functions_and_NRPy_basic_defines() # #define M_PI, etc.
# Declare paramstruct, register set_Cparameters_to_default(),
# and output declare_Cparameters_struct.h and set_Cparameters[].h:
outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesdir))
gri.register_C_functions_and_NRPy_basic_defines(enable_griddata_struct=False, enable_bcstruct_in_griddata_struct=False,
enable_rfmstruct=False,
enable_MoL_gfs_struct=False,
extras_in_griddata_struct=None) # #define IDX3S(), etc.
fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=True,
enable_SIMD=enable_SIMD) # #define NGHOSTS, and UPWIND() macro if SIMD disabled
# Output functions for computing all finite-difference stencils.
# Must be called after defining all functions depending on FD stencils.
if enable_FD_functions:
fin.output_finite_difference_functions_h(path=Ccodesdir)
# Call this last: Set up NRPy_basic_defines.h and NRPy_function_prototypes.h.
outC.construct_NRPy_basic_defines_h(Ccodesdir, enable_SIMD=enable_SIMD)
with open(os.path.join(Ccodesdir,"GiRaFFE_basic_defines.h"),"w") as file:
file.write("""#define NGHOSTS_A2B """+str(2)+"\n"+"""extern int kronecker_delta[4][3];
extern int MAXFACE;
extern int NUL;
extern int MINFACE;
// Some additional constants needed for PPM:
extern int VX,VY,VZ,
BX_CENTER,BY_CENTER,BZ_CENTER,BX_STAGGER,BY_STAGGER,BZ_STAGGER,
VXR,VYR,VZR,VXL,VYL,VZL; //<-- Be _sure_ to define MAXNUMVARS appropriately!
extern int NUM_RECONSTRUCT_GFS;
// Structure to track ghostzones for PPM:
typedef struct __gf_and_gz_struct__ {
REAL *gf;
int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;
""")
with open(os.path.join(Ccodesdir,"GiRaFFE_main_defines.h"),"w") as file:
file.write("""#define NGHOSTS_A2B """+str(2)+"\n"+PPM.kronecker_code+"""const int MAXFACE = -1;
const int NUL = +0;
const int MINFACE = +1;
// Some additional constants needed for PPM:
const int VX=0,VY=1,VZ=2,
BX_CENTER=3,BY_CENTER=4,BZ_CENTER=5,BX_STAGGER=6,BY_STAGGER=7,BZ_STAGGER=8,
VXR=9,VYR=10,VZR=11,VXL=12,VYL=13,VZL=14; //<-- Be _sure_ to define MAXNUMVARS appropriately!
const int NUM_RECONSTRUCT_GFS = 15;
// Structure to track ghostzones for PPM:
typedef struct __gf_and_gz_struct__ {
REAL *gf;
int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;
""")
outC.construct_NRPy_function_prototypes_h(Ccodesdir)
cmd.new_C_compile(Ccodesdir, os.path.join("output", "GiRaFFE_NRPy_standalone"),
uses_free_parameters_h=True, compiler_opt_option="fast") # fastdebug or debug also supported
# !gcc -g -O2 -fopenmp GiRaFFE_standalone_Ccodes/GiRaFFE_NRPy_standalone.c -o GiRaFFE_NRPy_standalone -lm
# Change to output directory
os.chdir(outdir)
# Clean up existing output files
cmd.delete_existing_files("out*.txt")
cmd.delete_existing_files("out*.png")
# cmd.Execute(os.path.join(Ccodesdir,"output","GiRaFFE_NRPy_standalone"), "640 16 16", os.path.join(outdir,"out640.txt"))
cmd.Execute("GiRaFFE_NRPy_standalone", "64 64 64","out64.txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "239 15 15","out239.txt")
# !OMP_NUM_THREADS=1 valgrind --track-origins=yes -v ./GiRaFFE_NRPy_standalone 1280 32 32
# Return to root directory
os.chdir(os.path.join("../../"))
(EXEC): Executing `make -j10`... (BENCH): Finished executing in 0.8128435611724854 seconds. Finished compilation. (EXEC): Executing `taskset -c 0,1,2,3 ./GiRaFFE_NRPy_standalone 64 64 64`... (BENCH): Finished executing in 18.257782697677612 seconds.
Now, we will load the data generated by the simulation and plot it in order to test for convergence.
import numpy as np
import matplotlib.pyplot as plt
# 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"]
Data_numer = np.loadtxt(os.path.join(Ccodesdir,"output","out64-00000010.txt"))
# Data_num_2 = np.loadtxt(os.path.join(Ccodesdir,"output","out239-00000080.txt"))
# Data_old = np.loadtxt("/home/penelson/OldCactus/Cactus/exe/ABE-GiRaFFEfood_1D_AlfvenWave/giraffe-grmhd_primitives_bi.x.asc")
# Data_o_2 = np.loadtxt("/home/penelson/OldCactus/Cactus/exe/ABE-GiRaFFEfood_1D_AlfvenWave_2/giraffe-grmhd_primitives_bi.x.asc")
# Data_numer = Data_old[5000:5125,11:15] # The column range is chosen for compatibility with the plotting script.
# Data_num_2 = Data_o_2[19600:19845,11:15] # The column range is chosen for compatibility with the plotting script.
Data_exact = np.loadtxt(os.path.join(Ccodesdir,"output","out64-00000010.txt"))
# Data_exa_2 = np.loadtxt(os.path.join(Ccodesdir,"output","out239-00000080_exact.txt"))
predicted_order = 2.0
column = 3
plt.figure()
# # plt.plot(Data_exact[2:-2,0],np.log2(np.absolute((Data_numer[2:-2,column]-Data_exact[2:-2,column])/\
# # (Data_num_2[2:-2:2,column]-Data_exa_2[2:-2:2,column]))),'.')
plt.plot(Data_exact[:,0],Data_exact[:,column])
plt.plot(Data_exact[:,0],Data_numer[:,column],'.')
# plt.xlim(-0.0,1.0)
# # plt.ylim(-1.0,5.0)
# # plt.ylim(-0.0005,0.0005)
# plt.xlabel("x")
# plt.ylabel("BU2")
plt.show()
# Plotting scripts for comparison with original GiRaFFE:
# old_files = ["",
# "giraffe-grmhd_primitives_bi.x.asc","giraffe-grmhd_primitives_bi.x.asc","giraffe-grmhd_primitives_bi.x.asc",
# "giraffe-em_ax.x.asc","giraffe-em_ay.x.asc","giraffe-em_az.x.asc",
# "giraffe-grmhd_conservatives.x.asc","giraffe-grmhd_conservatives.x.asc","giraffe-grmhd_conservatives.x.asc",
# "giraffe-grmhd_primitives_allbutbi.x.asc","giraffe-grmhd_primitives_allbutbi.x.asc","giraffe-grmhd_primitives_allbutbi.x.asc",
# "giraffe-em_psi6phi.x.asc"]
# column = 6
# column_old = [0,12,13,14,12,12,12,12,13,14,12,13,14,12]
# old_path = "/home/penelson/Cactus/exe/ABE-GiRaFFEfood_ExactWald_newcode_agreement"
# perturb_path = "/home/penelson/Cactus/exe/ABE-GiRaFFEfood_ExactWald_newcode_agreement_pert"
# new_path = os.path.join(Ccodesdir,"output")
# data_old = np.loadtxt(os.path.join(old_path,old_files[column]))
# data_per = np.loadtxt(os.path.join(perturb_path,old_files[column]))
# n=80
# data_old = data_old[n*70:n*70+70,:]# Select only the nth timestep
# data_per = data_per[n*70:n*70+70,:]# Select only the nth timestep
# data_new = np.loadtxt(os.path.join(new_path,"out64-00000080.txt"))
# deltaA_old = data_old[125:250,:] - data_old[0:125,:]
# data_new_t0 = np.loadtxt(os.path.join(new_path,"out119-00000000.txt"))
# deltaA_new = data_new[:,:] - data_new_t0[:,:]
# plt.figure()
# plt.semilogy(data_new[:,0],(data_new[:,column]-data_old[:,column_old[column]])/(data_new[:,column]+data_old[:,column_old[column]]),label="Old vs. New")
# plt.plot(data_new[4:22,0],data_new[4:22,column]-data_old[4:22,column_old[column]])
# plt.plot(data_new[49:-4,0],data_new[49:-4,column]-data_old[49:-4,column_old[column]])
# plt.plot(data_new[:,0],data_old[:,column_old[column]])
# plt.plot(data_new[:,0],data_new[:,column],'.')
# plt.plot(data_new[:22,0],data_old[:22,column_old[column]])
# plt.plot(data_new[49:,0],data_old[49:,column_old[column]])
# plt.plot(data_new[:22,0],data_new[:22,column],'.')
# plt.plot(data_new[49:,0],data_new[49:,column],'.')
# plt.plot(data_new[3:-3,0],deltaA_new[3:-3,column],'.')
# plt.plot(data_new[3:-3,0],deltaA_old[3:-3,column_old[column]]-deltaA_new[3:-3,column])
# plt.xlabel(labels[0])
# plt.ylabel(labels[column])
# plt.show()
# plt.figure()
# plt.semilogy(data_new[:,0],(data_old[:,column_old[column]]-data_per[:,column_old[column]])/(data_old[:,column_old[column]]+data_per[:,column_old[column]]),label="Old vs. Perturbed")
# plt.xlabel(labels[0])
# plt.ylabel("E_rel in "+labels[column])
# plt.legend()
# plt.show()
Here, we will use the output L2 norm data and find its convergence to determine how well our code converges.
By doing so, we show that the truncation error dominates and converges to zero at the expected rate. In the code above, we output the L2 norm within a thick spherical shell; for some quantity $Q$, this is given as $$ | Q_{\rm approx} - Q_{\rm exact}| = \sqrt{\frac{1}{N} \sum_{ijk} \left( Q_{\rm approx} - Q_{\rm exact} \right)^2} $$
Here, we import that data and calculate the convergence in the usual way, $$ k = \log_2 \left( \frac{F - F_1}{F - F_2} \right), $$ where $k$ is the convergence order, $F$ is the exact solution, $F_1$ is the approximate solution on the coarser grid with resolution $\Delta x$, and $F_2$ is the approximate solution on the finer grid with resolution $\Delta x/2$. Note that for the exact solution, the L2 norm is $0$.
# l2_lo = np.loadtxt(os.path.join(new_path,"out64-L2.txt"))
# l2_hi = np.loadtxt(os.path.join(new_path,"out128-L2.txt"))
# # 0 1 2 3 4 5 6 7 8 9 10 11 12 13
# labels = ["n","BU0","BU1","BU2","AD0","AD1","AD2","StildeD0","StildeD1","StildeD2","ValenciavU0","ValenciavU1","ValenciavU2", "psi6Phi"]
# def conv_order(low,high):
# order = np.log2(low/high)
# return order
# for column in range(1,len(labels)):
# print(labels[column]+": "+str(conv_order(l2_lo[1,column],l2_hi[2,column])))
This code will create an animation of the wave over time.
# import matplotlib.pyplot as plt
from matplotlib.pyplot import savefig
from IPython.display import HTML
import matplotlib.image as mgimg
import glob
import sys
from matplotlib import animation
cmd.delete_existing_files("out64-00*.png")
globby = glob.glob(os.path.join(Ccodesdir,'output','out64-00*.txt'))
file_list = []
for x in sorted(globby):
file_list.append(x)
number_of_files = int(len(file_list))
for timestep in range(number_of_files):
fig = plt.figure()
numer_filename = file_list[timestep]
exact_filename = os.path.join(Ccodesdir,'output','out64-00000000.txt')
Numer = np.loadtxt(numer_filename)
Exact = np.loadtxt(exact_filename)
plt.title("Exact Wald")
plt.xlabel("x")
plt.ylabel("BU2")
# plt.xlim(-0.5,0.5)
# plt.ylim(1.0,1.7)
plt.plot(Numer[3:22,0],Numer[3:22,3],'.',label="Numerical")
plt.plot(Numer[49:-3,0],Numer[49:-3,3],'.',label="Numerical")
plt.plot(Exact[3:22,0],Exact[3:22,3],label="Exact")
plt.plot(Exact[49:-3,0],Exact[49:-3,3],label="Exact")
plt.legend()
savefig(numer_filename+".png",dpi=150)
plt.close(fig)
sys.stdout.write("%c[2K" % 27)
sys.stdout.write("Processing file "+numer_filename+"\r")
sys.stdout.flush()
Processing file GiRaFFE_staggered_new_way_standalone_Ccodes/output/out64-00000087.txt
## VISUALIZATION ANIMATION, PART 2: Combine PNGs to generate movie ##
# https://stackoverflow.com/questions/14908576/how-to-remove-frame-from-matplotlib-pyplot-figure-vs-matplotlib-figure-frame
# https://stackoverflow.com/questions/23176161/animating-pngs-in-matplotlib-using-artistanimation
# !rm -f GiRaFFE_NRPy-1D_tests.mp4
cmd.delete_existing_files("GiRaFFE_NRPy-3D_tests.mp4")
fig = plt.figure(frameon=False)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
myimages = []
for i in range(number_of_files):
img = mgimg.imread(file_list[i]+".png")
imgplot = plt.imshow(img)
myimages.append([imgplot])
ani = animation.ArtistAnimation(fig, myimages, interval=100, repeat_delay=1000)
plt.close()
ani.save('GiRaFFE_NRPy-1D_tests.mp4', fps=5,dpi=150)
%%HTML
<video width="480" height="360" controls>
<source src="GiRaFFE_NRPy-1D_tests.mp4" type="video/mp4">
</video>
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFE_NRPy_Main_Driver")
Created Tutorial-GiRaFFE_NRPy_Main_Driver.tex, and compiled LaTeX file to PDF file Tutorial-GiRaFFE_NRPy_Main_Driver.pdf