This module is currently under development
This module is organized as follows
set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables.C
param.ccl
interface.ccl
schedule.ccl
make.code.defn
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: Load up cmdline_helper and create the directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
import cmdline_helper as cmd
CtH_dir_path = os.path.join("..","Convert_to_HydroBase")
cmd.mkdir(CtH_dir_path)
CtH_src_dir_path = os.path.join(CtH_dir_path,"src")
cmd.mkdir(CtH_src_dir_path)
# Step 0b: Create the output file path
outfile_path__Convert_to_HydroBase__source = os.path.join(CtH_src_dir_path,"Convert_to_HydroBase.C")
outfile_path__Convert_to_HydroBase__make = os.path.join(CtH_src_dir_path,"make.code.defn")
outfile_path__Convert_to_HydroBase__param = os.path.join(CtH_dir_path,"param.ccl")
outfile_path__Convert_to_HydroBase__interface = os.path.join(CtH_dir_path,"interface.ccl")
outfile_path__Convert_to_HydroBase__schedule = os.path.join(CtH_dir_path,"schedule.ccl")
%%writefile $outfile_path__Convert_to_HydroBase__source
#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"
void Convert_to_HydroBase(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
// Generally, we only need the HydroBase variables for diagnostic purposes, so we run the below loop only at iterations in which diagnostics are run.
if(Convert_to_HydroBase_every==0 || cctk_iteration%Convert_to_HydroBase_every!=0) return;
/***************
* PPEOS Patch *
***************
* We will need to set up our EOS in
* order to be able to compute eps below
*/
eos_struct eos;
initialize_EOS_struct_from_input(eos);
#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++) {
int index = CCTK_GFINDEX3D(cctkGH,i,j,k);
/* Note that we currently do not set Abar, Y_e, temperature, entropy, Avec[3], Aphi, Avec_stag[3], Aphi_stag */
CCTK_REAL PRIMS[MAXNUMVARS];
int ww=0;
PRIMS[ww] = rho_b[index]; ww++;
PRIMS[ww] = P[index]; ww++;
PRIMS[ww] = vx[index]; ww++;
PRIMS[ww] = vy[index]; ww++;
PRIMS[ww] = vz[index]; ww++;
PRIMS[ww] = Bx[index]; ww++;
PRIMS[ww] = By[index]; ww++;
PRIMS[ww] = Bz[index]; ww++;
rho[index] = PRIMS[RHOB];
press[index] = PRIMS[PRESSURE];
/***************
* PPEOS Patch *
***************
* For our hybrid piecewise polytropic EOS,
* we have
* .------------------------------------------------------.
* | eps = eps_cold + (P - P_cold)/( rho*(Gamma_th - 1) ) |
* .------------------------------------------------------.
*/
/* Compute P_cold and eps_cold */
CCTK_REAL P_cold, eps_cold;
compute_P_cold__eps_cold(eos,PRIMS[RHOB], P_cold,eps_cold); /* <- This function is defined in inlined_functions.C */
/* Compute eps as described above */
eps[index] = (PRIMS[PRESSURE]-P_cold)/PRIMS[RHOB]/(Gamma_th-1.0);
// IllinoisGRMHD defines v^i = u^i/u^0.
// Meanwhile, the ET/HydroBase formalism, called the Valencia
// formalism, splits the 4 velocity into a purely spatial part
// and a part that is normal to the spatial hypersurface:
// u^a = G (n^a + U^a), (Eq. 14 of arXiv:1304.5544; G=W, U^a=v^a)
// where n^a is the unit normal vector to the spatial hypersurface,
// n_a = {-\alpha,0,0,0}, and U^a is the purely spatial part, which
// is defined in HydroBase as the vel[] vector gridfunction.
// Then u^a n_a = - \alpha u^0 = G n^a n_a = -G, and
// of course \alpha u^0 = 1/sqrt(1+γ^ij u_i u_j) = \Gamma,
// the standard Lorentz factor.
// Note that n^i = - \beta^i / \alpha, so
// u^a = \Gamma (n^a + U^a)
// -> u^i = \Gamma ( U^i - \beta^i / \alpha )
// which implies
// v^i = u^i/u^0
// = \Gamma/u^0 ( U^i - \beta^i / \alpha ) <- \Gamma = \alpha u^0
// = \alpha ( U^i - \beta^i / \alpha )
// = \alpha U^i - \beta^i
CCTK_REAL lapseL=alp[index];
CCTK_REAL lapseL_inv=1.0/lapseL;
vel[CCTK_GFINDEX4D(cctkGH,i,j,k,0)] = (PRIMS[VX] + betax[index])*lapseL_inv;
vel[CCTK_GFINDEX4D(cctkGH,i,j,k,1)] = (PRIMS[VY] + betay[index])*lapseL_inv;
vel[CCTK_GFINDEX4D(cctkGH,i,j,k,2)] = (PRIMS[VZ] + betaz[index])*lapseL_inv;
// \alpha u^0 = 1/sqrt(1+γ^ij u_i u_j) = \Gamma = w_lorentz
// First compute u^0:
// Derivation of first equation:
// \gamma_{ij} (v^i + \beta^i)(v^j + \beta^j)/(\alpha)^2
// = \gamma_{ij} 1/(u^0)^2 ( \gamma^{ik} u_k \gamma^{jl} u_l /(\alpha)^2 <- Using Eq. 53 of arXiv:astro-ph/0503420
// = 1/(u^0 \alpha)^2 u_j u_l \gamma^{jl} <- Since \gamma_{ij} \gamma^{ik} = \delta^k_j
// = 1/(u^0 \alpha)^2 ( (u^0 \alpha)^2 - 1 ) <- Using Eq. 56 of arXiv:astro-ph/0503420
// = 1 - 1/(u^0 \alpha)^2 <= 1
CCTK_REAL shiftxL = betax[index];
CCTK_REAL shiftyL = betay[index];
CCTK_REAL shiftzL = betaz[index];
CCTK_REAL gxxL = gxx[index];
CCTK_REAL gxyL = gxy[index];
CCTK_REAL gxzL = gxz[index];
CCTK_REAL gyyL = gyy[index];
CCTK_REAL gyzL = gyz[index];
CCTK_REAL gzzL = gzz[index];
CCTK_REAL one_minus_one_over_alpha_u0_squared = (gxxL* SQR(PRIMS[VX] + shiftxL) +
2.0*gxyL*(PRIMS[VX] + shiftxL)*(PRIMS[VY] + shiftyL) +
2.0*gxzL*(PRIMS[VX] + shiftxL)*(PRIMS[VZ] + shiftzL) +
gyyL* SQR(PRIMS[VY] + shiftyL) +
2.0*gyzL*(PRIMS[VY] + shiftyL)*(PRIMS[VZ] + shiftzL) +
gzzL* SQR(PRIMS[VZ] + shiftzL) )*SQR(lapseL_inv);
/*** Check for superluminal velocity ***/
//FIXME: Instead of >1.0, should be one_minus_one_over_alpha_u0_squared > ONE_MINUS_ONE_OVER_GAMMA_SPEED_LIMIT_SQUARED, for consistency with conserv_to_prims routines
if(one_minus_one_over_alpha_u0_squared > 1.0) {
CCTK_VInfo(CCTK_THORNSTRING,"Convert_to_HydroBase WARNING: Found superluminal velocity. This should have been caught by IllinoisGRMHD.");
}
// A = 1.0-one_minus_one_over_alpha_u0_squared = 1-(1-1/(al u0)^2) = 1/(al u0)^2
// 1/sqrt(A) = al u0
CCTK_REAL alpha_u0 = 1.0/sqrt(1.0-one_minus_one_over_alpha_u0_squared);
if(std::isnan(alpha_u0*lapseL_inv)) printf("BAD FOUND NAN ALPHAU0 CALC: %.15e %.15e %.15e\n",alpha_u0,lapseL_inv,one_minus_one_over_alpha_u0_squared);
w_lorentz[index] = alpha_u0;
Bvec[CCTK_GFINDEX4D(cctkGH,i,j,k,0)] = PRIMS[BX_CENTER];
Bvec[CCTK_GFINDEX4D(cctkGH,i,j,k,1)] = PRIMS[BY_CENTER];
Bvec[CCTK_GFINDEX4D(cctkGH,i,j,k,2)] = PRIMS[BZ_CENTER];
}
}
Overwriting ../Convert_to_HydroBase/src/Convert_to_HydroBase.C
%%writefile $outfile_path__Convert_to_HydroBase__param
# Parameter definitions for thorn convert_to_HydroBase
# $Header:$
#############################################################################
### import HydroBase & ADMBase parameters
shares: HydroBase
USES CCTK_INT timelevels
shares: ADMBase
USES CCTK_INT lapse_timelevels
USES CCTK_INT shift_timelevels
USES CCTK_INT metric_timelevels
shares: IllinoisGRMHD
USES CCTK_INT neos
USES CCTK_REAL Gamma_th
USES CCTK_REAL K_ppoly_tab0
USES CCTK_REAL rho_ppoly_tab_in[10]
USES CCTK_REAL Gamma_ppoly_tab_in[10]
#############################################################################
private:
INT Convert_to_HydroBase_every "How often to convert IllinoisGRMHD primitive variables to HydroBase (Valencia formulation) primitive variables? Needed for some ET-based diagnostics. NOT needed for pure IllinoisGRMHD runs."
{
0:* :: "zero (disable) or positive (every N iterations)"
} 0
Overwriting ../Convert_to_HydroBase/param.ccl
%%writefile $outfile_path__Convert_to_HydroBase__interface
# Interface definition for thorn Convert_to_HydroBase
# $Header:$
implements: Convert_to_HydroBase
inherits: grid HydroBase ADMBase IllinoisGRMHD
uses include header: IllinoisGRMHD_headers.h
Overwriting ../Convert_to_HydroBase/interface.ccl
%%writefile $outfile_path__Convert_to_HydroBase__schedule
# Schedule definitions for thorn Convert_to_HydroBase
# $Header:$
SCHEDULE Convert_to_HydroBase AT CCTK_INITIAL AFTER SetTmunu
{
LANG: C
} "Convert IllinoisGRMHD-native variables to HydroBase"
SCHEDULE Convert_to_HydroBase AT CCTK_ANALYSIS BEFORE compute_bi_b2_Poyn_fluxET BEFORE particle_tracerET BEFORE VolumeIntegralGroup BEFORE convert_to_MHD_3velocity AFTER ML_BSSN_evolCalcGroup
{
OPTIONS: GLOBAL-EARLY,LOOP-LOCAL
LANG: C
} "Convert IllinoisGRMHD-native variables to HydroBase"
Overwriting ../Convert_to_HydroBase/schedule.ccl
%%writefile $outfile_path__Convert_to_HydroBase__make
# Main make.code.defn file for thorn Convert_to_HydroBase
# $Header:$
# Source files in this directory
SRCS = Convert_to_HydroBase.C
# Subdirectories containing source files
SUBDIRS =
Overwriting ../Convert_to_HydroBase/src/make.code.defn
# # 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/A_i_rhs_no_gauge_terms.C"
# original_IGM_file_name = "A_i_rhs_no_gauge_terms-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__A_i_rhs_no_gauge_terms__C = !diff $original_IGM_file_path $outfile_path__A_i_rhs_no_gauge_terms__C
# if Validation__A_i_rhs_no_gauge_terms__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 A_i_rhs_no_gauge_terms.C: PASSED!")
# else:
# # If the validation fails, we keep the original IGM source code file
# print("Validation test for A_i_rhs_no_gauge_terms.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__A_i_rhs_no_gauge_terms__C:
# print(diff_line)
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__A_i_rhs_no_gauge_terms.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__A_i_rhs_no_gauge_terms.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex
!rm -f Tut*.out Tut*.aux Tut*.log