This module is currently under development
We will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory, if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__compute_B_and_Bstagger_from_A__C = os.path.join(IGM_src_dir_path,"compute_B_and_Bstagger_from_A.C")
%%writefile $outfile_path__compute_B_and_Bstagger_from_A__C
#include "cctk.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <sys/time.h>
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "IllinoisGRMHD_headers.h"
#define LOOP_DEFINE_SIMPLE \
_Pragma("omp parallel for") \
for(int k=0;k<cctk_lsh[2];k++) \
for(int j=0;j<cctk_lsh[1];j++) \
for(int i=0;i<cctk_lsh[0];i++)
extern "C" void IllinoisGRMHD_compute_B_and_Bstagger_from_A(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_REAL dxi = 1.0/CCTK_DELTA_SPACE(0);
CCTK_REAL dyi = 1.0/CCTK_DELTA_SPACE(1);
CCTK_REAL dzi = 1.0/CCTK_DELTA_SPACE(2);
CCTK_REAL gridfunc_syms_Ax[3] = {-1, 1, Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Ax , gridfunc_syms_Ax,0,1,1);
CCTK_REAL gridfunc_syms_Ay[3] = { 1,-1, Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Ay , gridfunc_syms_Ay,1,0,1);
CCTK_REAL gridfunc_syms_Az[3] = { 1, 1,-Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Az , gridfunc_syms_Az,1,1,0);
CCTK_REAL gridfunc_syms_psi6phi[3] = { 1, 1, 1};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,psi6phi , gridfunc_syms_psi6phi,1,1,1);
LOOP_DEFINE_SIMPLE {
int index=CCTK_GFINDEX3D(cctkGH,i,j,k);
psi_bssn[index] = exp(phi_bssn[index]);
}
LOOP_DEFINE_SIMPLE {
// Look Mom, no if() statements!
int shiftedim1 = (i-1)*(i!=0); // This way, i=0 yields shiftedim1=0 and shiftedi=1, used below for our COPY boundary condition.
int shiftedi = shiftedim1+1;
int shiftedjm1 = (j-1)*(j!=0);
int shiftedj = shiftedjm1+1;
int shiftedkm1 = (k-1)*(k!=0);
int shiftedk = shiftedkm1+1;
int index,indexim1,indexjm1,indexkm1;
int actual_index = CCTK_GFINDEX3D(cctkGH,i,j,k);
CCTK_REAL Psi = psi_bssn[actual_index];
CCTK_REAL Psim3 = 1.0/(Psi*Psi*Psi);
// For the lower boundaries, the following applies a "copy"
// boundary condition on Bi_stagger where needed.
// E.g., Bx_stagger(i,jmin,k) = Bx_stagger(i,jmin+1,k)
// We find the copy BC works better than extrapolation.
// For the upper boundaries, we do the following copy:
// E.g., Psi(imax+1,j,k)=Psi(imax,j,k)
/**************/
/* Bx_stagger */
/**************/
index = CCTK_GFINDEX3D(cctkGH,i,shiftedj,shiftedk);
indexjm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedjm1,shiftedk);
indexkm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedj,shiftedkm1);
// Set Bx_stagger = \partial_y A_z - partial_z A_y
// "Grid" Ax(i,j,k) is actually Ax(i,j+1/2,k+1/2)
// "Grid" Ay(i,j,k) is actually Ay(i+1/2,j,k+1/2)
// "Grid" Az(i,j,k) is actually Az(i+1/2,j+1/2,k)
// Therefore, the 2nd order derivative \partial_z A_y at (i+1/2,j,k) is:
// ["Grid" Ay(i,j,k) - "Grid" Ay(i,j,k-1)]/dZ
Bx_stagger[actual_index] = (Az[index]-Az[indexjm1])*dyi - (Ay[index]-Ay[indexkm1])*dzi;
// Now multiply Bx and Bx_stagger by 1/sqrt(gamma(i+1/2,j,k)]) = 1/sqrt(1/2 [gamma + gamma_ip1]) = exp(-6 x 1/2 [phi + phi_ip1] )
int imax_minus_i = (cctk_lsh[0]-1)-i;
int indexip1jk = CCTK_GFINDEX3D(cctkGH,i + ( (imax_minus_i > 0) - (0 > imax_minus_i) ),j,k);
CCTK_REAL Psi_ip1 = psi_bssn[indexip1jk];
Bx_stagger[actual_index] *= Psim3/(Psi_ip1*Psi_ip1*Psi_ip1);
/**************/
/* By_stagger */
/**************/
index = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedk);
indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,shiftedk);
indexkm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedkm1);
// Set By_stagger = \partial_z A_x - \partial_x A_z
By_stagger[actual_index] = (Ax[index]-Ax[indexkm1])*dzi - (Az[index]-Az[indexim1])*dxi;
// Now multiply By and By_stagger by 1/sqrt(gamma(i,j+1/2,k)]) = 1/sqrt(1/2 [gamma + gamma_jp1]) = exp(-6 x 1/2 [phi + phi_jp1] )
int jmax_minus_j = (cctk_lsh[1]-1)-j;
int indexijp1k = CCTK_GFINDEX3D(cctkGH,i,j + ( (jmax_minus_j > 0) - (0 > jmax_minus_j) ),k);
CCTK_REAL Psi_jp1 = psi_bssn[indexijp1k];
By_stagger[actual_index] *= Psim3/(Psi_jp1*Psi_jp1*Psi_jp1);
/**************/
/* Bz_stagger */
/**************/
index = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedj,k);
indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,shiftedj,k);
indexjm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedjm1,k);
// Set Bz_stagger = \partial_x A_y - \partial_y A_x
Bz_stagger[actual_index] = (Ay[index]-Ay[indexim1])*dxi - (Ax[index]-Ax[indexjm1])*dyi;
// Now multiply Bz_stagger by 1/sqrt(gamma(i,j,k+1/2)]) = 1/sqrt(1/2 [gamma + gamma_kp1]) = exp(-6 x 1/2 [phi + phi_kp1] )
int kmax_minus_k = (cctk_lsh[2]-1)-k;
int indexijkp1 = CCTK_GFINDEX3D(cctkGH,i,j,k + ( (kmax_minus_k > 0) - (0 > kmax_minus_k) ));
CCTK_REAL Psi_kp1 = psi_bssn[indexijkp1];
Bz_stagger[actual_index] *= Psim3/(Psi_kp1*Psi_kp1*Psi_kp1);
}
LOOP_DEFINE_SIMPLE {
// Look Mom, no if() statements!
int shiftedim1 = (i-1)*(i!=0); // This way, i=0 yields shiftedim1=0 and shiftedi=1, used below for our COPY boundary condition.
int shiftedi = shiftedim1+1;
int shiftedjm1 = (j-1)*(j!=0);
int shiftedj = shiftedjm1+1;
int shiftedkm1 = (k-1)*(k!=0);
int shiftedk = shiftedkm1+1;
int index,indexim1,indexjm1,indexkm1;
int actual_index = CCTK_GFINDEX3D(cctkGH,i,j,k);
// For the lower boundaries, the following applies a "copy"
// boundary condition on Bi and Bi_stagger where needed.
// E.g., Bx(imin,j,k) = Bx(imin+1,j,k)
// We find the copy BC works better than extrapolation.
/******/
/* Bx */
/******/
index = CCTK_GFINDEX3D(cctkGH,shiftedi,j,k);
indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,k);
// Set Bx = 0.5 ( Bx_stagger + Bx_stagger_im1 )
// "Grid" Bx_stagger(i,j,k) is actually Bx_stagger(i+1/2,j,k)
Bx[actual_index] = 0.5 * ( Bx_stagger[index] + Bx_stagger[indexim1] );
/******/
/* By */
/******/
index = CCTK_GFINDEX3D(cctkGH,i,shiftedj,k);
indexjm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedjm1,k);
// Set By = 0.5 ( By_stagger + By_stagger_im1 )
// "Grid" By_stagger(i,j,k) is actually By_stagger(i,j+1/2,k)
By[actual_index] = 0.5 * ( By_stagger[index] + By_stagger[indexjm1] );
/******/
/* Bz */
/******/
index = CCTK_GFINDEX3D(cctkGH,i,j,shiftedk);
indexkm1 = CCTK_GFINDEX3D(cctkGH,i,j,shiftedkm1);
// Set Bz = 0.5 ( Bz_stagger + Bz_stagger_im1 )
// "Grid" Bz_stagger(i,j,k) is actually Bz_stagger(i,j+1/2,k)
Bz[actual_index] = 0.5 * ( Bz_stagger[index] + Bz_stagger[indexkm1] );
}
// Finish up by setting symmetry ghostzones on Bx, By, Bz, and their staggered variants.
CCTK_REAL gridfunc_syms_Bx[3] = {-1, 1,-Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bx , gridfunc_syms_Bx,0,0,0);
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bx_stagger, gridfunc_syms_Bx,1,0,0);
CCTK_REAL gridfunc_syms_By[3] = { 1,-1,-Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, By , gridfunc_syms_By,0,0,0);
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, By_stagger, gridfunc_syms_By,0,1,0);
CCTK_REAL gridfunc_syms_Bz[3] = { 1, 1, Sym_Bz};
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bz , gridfunc_syms_Bz,0,0,0);
IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bz_stagger, gridfunc_syms_Bz,0,0,1);
}
Writing ../src/compute_B_and_Bstagger_from_A.C
# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code
# First download the original IllinoisGRMHD source code
import urllib
from os import path
original_IGM_file_url = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/compute_B_and_Bstagger_from_A.C"
original_IGM_file_name = "compute_B_and_Bstagger_from_A-original.C"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)
# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
try:
original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
# If all else fails, hope wget does the job
!wget -O $original_IGM_file_path $original_IGM_file_url
# Perform validation
Validation__compute_B_and_Bstagger_from_A__C = !diff $original_IGM_file_path $outfile_path__compute_B_and_Bstagger_from_A__C
if Validation__compute_B_and_Bstagger_from_A__C == []:
# If the validation passes, we do not need to store the original IGM source code file
!rm $original_IGM_file_path
print("Validation test for compute_B_and_Bstagger_from_A.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for compute_B_and_Bstagger_from_A.C: FAILED!")
# We also print out the difference between the code generated
# in this tutorial module and the original IGM source code
print("Diff:")
for diff_line in Validation__compute_B_and_Bstagger_from_A__C:
print(diff_line)
Validation test for compute_B_and_Bstagger_from_A.C: FAILED! Diff: 56c56 < // For the lower boundaries, the following applies a "copy" --- > // For the lower boundaries, the following applies a "copy" 114c114 < --- > 133c133 < // For the lower boundaries, the following applies a "copy" --- > // For the lower boundaries, the following applies a "copy" 175a176 >
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-IllinoisGRMHD__compute_B_and_Bstagger_from_A.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means).
latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex
!rm -f Tut*.out Tut*.aux Tut*.log