This notebook documents the function from the original GiRaFFE that calculates the flux for $\tilde{S}_i$ according to the method of Harten, Lax, and von Leer and Einfeldt (HLLE), assuming that we have calculated the values of the velocity and magnetic field on the cell faces according to the piecewise-parabolic method (PPM) of Colella and Woodward (1984), modified for the case of GRFFE.
Notebook Status: Validated
Validation Notes: This code has been validated to round-off level agreement with the corresponding code in the original GiRaFFE
The differential equations that GiRaFFE evolves are written in conservation form, and thus have two different terms that contribute to the time evolution of any conserved quantity $\mathcal{C}$, the flux term $\mathcal{F}$ and the source term $\mathcal{S}$:
We will use a high-resolution shock capturing scheme, which solves the Riemann problem using the HLL approximate Riemann solver to implement the above equation, to ensure that sharp features in our GRFFE fields are properly modeled.
In GRFFE, the evolution equation for the Poynting (i.e., electromagnetic momentum) flux $\tilde{S}_i$ is given as $$ \boxed{\partial_t \tilde{S}_i + \underbrace{ \partial_j \left( \alpha \sqrt{\gamma} T^j_{{\rm EM} i} \right)}_{\rm Flux\ term} = \underbrace{\frac{1}{2} \alpha \sqrt{\gamma} T^{\mu \nu}_{\rm EM} \partial_i g_{\mu \nu}}_{\rm Source\ term}.} $$ We can then see that, if we rewrite this, the right-hand side (RHS) describing the time evolution $\partial_t \tilde{S}_i$ consists of two terms: the flux term and the source term. The flux term in particular can be tricky, as it may be discontinuous due to sharp features that may appear e.g., at current sheets or inside black hole horizons. If we were to simply approximate the term using finite-difference derivatives, such sharp features would lead to a Gibbs-like phenomenon. So, we implement a different algorithm to take the derivative.
The flux term itself is, as written above,
$$ \alpha \sqrt{\gamma} T^i_{{\rm EM} j} = \alpha \sqrt{\gamma} g_{j \mu} T^{\mu i}_{\rm EM}, $$where
$$ T^{\mu \nu}_{\rm EM} = b^2 u^\mu u^\nu + \frac{1}{2} b^2 g^{\mu \nu} - b^\mu b^\nu. $$The C functions implemented in this notebook will compute the flux at a given point so that we can easily take its derivative later. Having reconstructed the values of $v^i_{(n)}$ and $B^i$ on the cell faces, we then compute the value of the flux of $\tilde{S}_i$ on each face using the Harten, Lax, and von Leer and Einfeldt (hereafter HLLE) approximate Riemann solver. In particular, for each component of $\tilde{S}_i$ in each direction, we compute the HLL flux as $$ F^{\rm HLL} = \frac{c_{\rm min} f_{\rm R} + c_{\rm max} f_{\rm L} - c_{\rm min} c_{\rm max} (U_{\rm R}-U_{\rm L})}{c_{\rm min} + c_{\rm max}}, $$ where $$ f = \alpha \sqrt{\gamma} T^j_{{\rm EM} i} $$ and $$ U = \tilde{S}_j. $$ Here, $i$ is direction in which we are computing the Poynting flux, and $j$ is the component of the momentum we are computing it for.
These two quantities are computed on both the left and right sides of each cell face. We will be able to draw heavily on the GRFFE module (Tutorial) and the GRHD module (Tutorial) to compute $u^0$, $u^i$, and $b^\mu$, as well as the index-lowered forms of those vectors. Critically, these quantities depend on the Valencia 3-velocity $v^i_{(n)}$ and magnetic field $B^i$. We will not be using the usual gridfunctions for these, but rather the ones that we have previously calculated on the left and right sides of the cell faces using the Piecewise Parabolic Method.
The speeds $c_\min$ and $c_\max$ reflect the characteristic speeds of the plasma waves. In GRFFE, the expressions defining them reduce to a function of only the metric quantities:
\begin{align} c_\min &= - \min(c_-,0) \\ c_\max &= \max(c_+,0), \end{align}where $$ c_\pm = \left. \left(-b \pm \sqrt{b^2-4ac}\right)\middle/ \left(2a\right) \right., $$ and $$a = 1/\alpha^2,$$ $$b = 2 \beta^i / \alpha^2$$ and $$c = g^{ii} - (\beta^i)^2/\alpha^2.$$
These are coded in Tutorial-GiRaFFE_NRPy-Characteristic_Speeds.
Note that $c_\pm$ must be computed on each cell face, meaning that all the above metric quantities must be interpolated to cell faces. Metric face values are computed as described in this notebook.
The algorithm for finite-volume methods in general is as follows:
This notebook is organized as follows
# Step 0: 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)
from outputC import outCfunction, lhrh, add_to_Cfunction_dict, outC_function_dict # NRPy+: Core C code output module
import finite_difference as fin # NRPy+: Finite difference C code generation module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
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 GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds as chsp # GRFFE: the characteristic speeds
thismodule = "GiRaFFE_NRPy-Stilde-flux"
Finally, we can compute the flux in each direction. This momentum flux in the $m$ direction is defined as $\alpha \sqrt{\gamma} T^m_{\ \ j}$, based on the input flux_dirn. We have already defined $\alpha \sqrt{\gamma}$, so all we need to do is calculate $T^m_{\ \ j}$, where $T^{\mu \nu}_{\rm EM} = b^2 u^\mu u^\nu + \frac{1}{2} b^2 g^{\mu \nu} - b^\mu b^\nu$. In doing this index-lowering operation, recall that $g^{\mu \nu} g_{\nu \alpha} = \delta^\mu_\alpha$. We will do so in accordance with the method published by Harten, Lax, and von Leer and Einfeldt (hereafter HLLE) to solve the Riemann problem. So, we define $f(u) = T^m_{\ \ j}$ on each face as
$$
f = \alpha \sqrt{\gamma} \left( (\rho+b^2)(u^0 v^m) u_j + (P+\frac{1}{2}b^2) \delta^m_j - b^m b_j \right);
$$
Because $\rho = P = 0$ in GRFFE and $u^0 v^m = u^m$ in general (since $v^m$ is the drift velocity here), this simplifies to
$$
f = \alpha \sqrt{\gamma} \left( b^2 u^m u_j + \frac{1}{2}b^2 \delta^m_j - b^m b_j \right).
$$
We use $j$ to correspond to the component of the flux we are calculating; that is, $j=0$ corresponds to $x$, and so forth (however, remember that in a NRPy+ 3-vector, the numbers will still run from 0 to 2). $\delta^i_j$ is the standard Kronecker delta. We also define U_{\rm R} and U_{\rm L}:
$$
U = \alpha \sqrt{\gamma} \left( (\rho+b^2) u^0 u_j - b^0 b_j \right),
$$
which, in GRFFE, simplifies to
$$
U = \alpha \sqrt{\gamma} \left( b^2 u^0 u_j - b^0 b_j \right).
$$
In NRPy+, we'll let the GRHD and GRFFE modules handle these.
and combine based on eq. 3.15 in the HLLE paper, $$ F^{\rm HLL} = \frac{c_{\rm min} f_{\rm R} + c_{\rm max} f_{\rm L} - c_{\rm min} c_{\rm max} (U_{\rm R}-U_{\rm L})}{c_{\rm min} + c_{\rm max}}, $$
We'll write the HLLE step as a function so that we can loop over flux_dirn and mom_comp and write each version needed as we need it.
# We'll rewrite this assuming that we've passed the entire reconstructed
# gridfunctions. You could also do this with only one point, but then you'd
# need to declare everything as a Cparam in NRPy+
import GRHD.equations as GRHD
import GRFFE.equations as GRFFE
def calculate_GRFFE_Tmunu_and_contractions(flux_dirn, mom_comp, gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi):
GRHD.compute_sqrtgammaDET(gammaDD)
GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)
GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
GRFFE.compute_TEM4UU(gammaDD, betaU, alpha, GRFFE.smallb4U, GRFFE.smallbsquared, GRHD.u4U_ito_ValenciavU)
GRFFE.compute_TEM4UD(gammaDD, betaU, alpha, GRFFE.TEM4UU)
# Compute conservative variables in terms of primitive variables
GRHD.compute_S_tildeD(alpha, GRHD.sqrtgammaDET, GRFFE.TEM4UD)
global U,F
# Flux F = alpha*sqrt{gamma}*T^i_j
F = alpha*GRHD.sqrtgammaDET*GRFFE.TEM4UD[flux_dirn+1][mom_comp+1]
# U = alpha*sqrt{gamma}*T^0_j = Stilde_j
U = GRHD.S_tildeD[mom_comp]
def HLLE_solver(cmax, cmin, Fr, Fl, Ur, Ul):
# This solves the Riemann problem for the mom_comp component of the momentum
# flux StildeD in the flux_dirn direction.
# st_j_flux = (c_\min f_R + c_\max f_L - c_\min c_\max ( st_j_r - st_j_l )) / (c_\min + c_\max)
return (cmin*Fr + cmax*Fl - cmin*cmax*(Ur-Ul) )/(cmax + cmin)
Finally, we write the function that computes the actual flux. We take the parameter flux_dirn as input, so we can eventually create one C file for each flux direction. In each file, we will include the math to calculate each momentum-flux component mom_comp in that direction by looping over mom_comp.
We have written the function HLLE_solver() so that we can easily compute the flux as specified by those two indices.
The characteristic wave speeds are calculated in the module Tutorial-GiRaFFE_NRPy-Characteristic_Speeds; we will simply import that here.
def calculate_Stilde_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi):
chsp.find_cmax_cmin(flux_dirn,gamma_faceDD,beta_faceU,alpha_face)
global Stilde_fluxD
Stilde_fluxD = ixp.zerorank3()
for mom_comp in range(3):
calculate_GRFFE_Tmunu_and_contractions(flux_dirn, mom_comp, gamma_faceDD,beta_faceU,alpha_face,\
Valenciav_rU,B_rU,sqrt4pi)
Fr = F
Ur = U
calculate_GRFFE_Tmunu_and_contractions(flux_dirn, mom_comp, gamma_faceDD,beta_faceU,alpha_face,\
Valenciav_lU,B_lU,sqrt4pi)
Fl = F
Ul = U
Stilde_fluxD[mom_comp] = HLLE_solver(chsp.cmax, chsp.cmin, Fr, Fl, Ur, Ul)
There is some additional complexity to consider in generating the C code for these expressions. The flux term we need to add to the RHSs is a finite difference of the fluxes we have calculated so far, so these cannot be simple pointwise operations. However, we also cannot use NRPy+'s build-in finite-differencing tools because of how we store the reconstructed quantities (that is, quantities reconstructed on the $i-1/2$ face is stored at point $i$ in memory), which makes the FD template we need look just like a forward finite-differencing template, which NRPy+ cannot do. So, we must write the code to read and write data from and to memory ourselves. We will do so algorithmically with python string operations, so as to reduce the possibility from human error.
We will loop over the interior, but use string replacement to include an extra point in the $+x,y,z$ ghostzones. The finite differences here will be done backwards from what we're used to; while we would normally read from several points and then write to a single point, we will read from a single point and then write to two points. In this particular case, this will reduce the number of times we recalculate things.
def generate_C_code_for_Stilde_flux(out_dir,inputs_provided = False, alpha_face=None, gamma_faceDD=None, beta_faceU=None,
Valenciav_rU=None, B_rU=None, Valenciav_lU=None, B_lU=None, sqrt4pi=None,
outCparams = "outCverbose=False,CSE_sorting=none", write_cmax_cmin=False):
if not inputs_provided:
# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")
# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
# We'll also need to store the results of the HLLE step between functions.
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")
input_params_for_Stilde_flux = "const paramstruct *params,REAL *auxevol_gfs,REAL *rhs_gfs"
if write_cmax_cmin:
name_suffixes = ["_x","_y","_z"]
for flux_dirn in range(3):
calculate_Stilde_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi)
Stilde_flux_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED0"),rhs=Stilde_fluxD[0]),
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED1"),rhs=Stilde_fluxD[1]),
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED2"),rhs=Stilde_fluxD[2])
]
if write_cmax_cmin:
Stilde_flux_to_print = Stilde_flux_to_print \
+[
lhrh(lhs=gri.gfaccess("out_gfs","cmax"+name_suffixes[flux_dirn]),rhs=chsp.cmax),
lhrh(lhs=gri.gfaccess("out_gfs","cmin"+name_suffixes[flux_dirn]),rhs=chsp.cmin)
]
desc = "Compute the flux term of all 3 components of tilde{S}_i on the left face in the " + str(flux_dirn) + "direction for all components."
name = "calculate_Stilde_flux_D" + str(flux_dirn)
Ccode_function = outCfunction(
outfile = "returnstring", desc=desc, name=name,
params = input_params_for_Stilde_flux,
body = fin.FD_outputC("returnstring",Stilde_flux_to_print,params=outCparams),
loopopts ="InteriorPoints",
rel_path_to_Cparams=os.path.join("../")).replace("NGHOSTS+Nxx0","NGHOSTS+Nxx0+1").replace("NGHOSTS+Nxx1","NGHOSTS+Nxx1+1").replace("NGHOSTS+Nxx2","NGHOSTS+Nxx2+1")
with open(os.path.join(out_dir,name+".h"),"w") as file:
file.write(Ccode_function)
pre_body = """// Notice in the loop below that we go from 3 to cctk_lsh-3 for i, j, AND k, even though
// we are only computing the flux in one direction. This is because in the end,
// we only need the rhs's from 3 to cctk_lsh-3 for i, j, and k.
const REAL invdxi[4] = {1e100,invdx0,invdx1,invdx2};
const REAL invdx = invdxi[flux_dirn];"""
FD_body = """const int index = IDX3S(i0,i1,i2);
const int indexp1 = IDX3S(i0+kronecker_delta[flux_dirn][0],i1+kronecker_delta[flux_dirn][1],i2+kronecker_delta[flux_dirn][2]);
rhs_gfs[IDX4ptS(STILDED0GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED0GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED0GF,indexp1)] ) * invdx;
rhs_gfs[IDX4ptS(STILDED1GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED1GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED1GF,indexp1)] ) * invdx;
rhs_gfs[IDX4ptS(STILDED2GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED2GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED2GF,indexp1)] ) * invdx;"""
desc = "Compute the difference in the flux of StildeD on the opposite faces in flux_dirn for all components."
name = "calculate_Stilde_rhsD"
outCfunction(
outfile = os.path.join(out_dir,name+".h"), desc=desc, name=name,
params = "const int flux_dirn,const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs",
preloop = pre_body,
body = FD_body,
loopopts = "InteriorPoints",
rel_path_to_Cparams=os.path.join("../")
)
def add_to_Cfunction_dict__Stilde_flux(includes=None, rel_path_to_Cparams=os.path.join("../"),
path_from_rootsrcdir_to_this_Cfunc=os.path.join("RHSs/"),
inputs_provided = False, alpha_face=None, gamma_faceDD=None, beta_faceU=None,
Valenciav_rU=None, B_rU=None, Valenciav_lU=None, B_lU=None, sqrt4pi=None,
outCparams = "outCverbose=False,CSE_sorting=none", write_cmax_cmin=False):
if not inputs_provided:
# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")
# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
# We'll also need to store the results of the HLLE step between functions.
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")
input_params_for_Stilde_flux = "const paramstruct *params,REAL *auxevol_gfs,REAL *rhs_gfs"
if write_cmax_cmin:
name_suffixes = ["_x","_y","_z"]
for flux_dirn in range(3):
calculate_Stilde_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi)
Stilde_flux_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED0"),rhs=Stilde_fluxD[0]),
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED1"),rhs=Stilde_fluxD[1]),
lhrh(lhs=gri.gfaccess("out_gfs","Stilde_flux_HLLED2"),rhs=Stilde_fluxD[2])
]
if write_cmax_cmin:
Stilde_flux_to_print = Stilde_flux_to_print \
+[
lhrh(lhs=gri.gfaccess("out_gfs","cmax"+name_suffixes[flux_dirn]),rhs=chsp.cmax),
lhrh(lhs=gri.gfaccess("out_gfs","cmin"+name_suffixes[flux_dirn]),rhs=chsp.cmin)
]
desc = "Compute the flux term of all 3 components of tilde{S}_i on the left face in the " + str(flux_dirn) + "direction for all components."
name = "calculate_Stilde_flux_D" + str(flux_dirn)
body = fin.FD_outputC("returnstring",Stilde_flux_to_print,params=outCparams)
loopopts ="InteriorPoints"
add_to_Cfunction_dict(
includes=includes,
desc=desc,
name=name, params=input_params_for_Stilde_flux,
body=body, loopopts=loopopts,
path_from_rootsrcdir_to_this_Cfunc = path_from_rootsrcdir_to_this_Cfunc,
rel_path_to_Cparams=rel_path_to_Cparams)
outC_function_dict[name] = outC_function_dict[name].replace("NGHOSTS+Nxx0","NGHOSTS+Nxx0+1").replace("NGHOSTS+Nxx1","NGHOSTS+Nxx1+1").replace("NGHOSTS+Nxx2","NGHOSTS+Nxx2+1")
pre_body = """// Notice in the loop below that we go from 3 to cctk_lsh-3 for i, j, AND k, even though
// we are only computing the flux in one direction. This is because in the end,
// we only need the rhs's from 3 to cctk_lsh-3 for i, j, and k.
const REAL invdxi[4] = {1e100,invdx0,invdx1,invdx2};
const REAL invdx = invdxi[flux_dirn];"""
FD_body = """const int index = IDX3S(i0,i1,i2);
const int indexp1 = IDX3S(i0+kronecker_delta[flux_dirn][0],i1+kronecker_delta[flux_dirn][1],i2+kronecker_delta[flux_dirn][2]);
rhs_gfs[IDX4ptS(STILDED0GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED0GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED0GF,indexp1)] ) * invdx;
rhs_gfs[IDX4ptS(STILDED1GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED1GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED1GF,indexp1)] ) * invdx;
rhs_gfs[IDX4ptS(STILDED2GF,index)] += (auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED2GF,index)] - auxevol_gfs[IDX4ptS(STILDE_FLUX_HLLED2GF,indexp1)] ) * invdx;"""
desc = "Compute the difference in the flux of StildeD on the opposite faces in flux_dirn for all components."
name = "calculate_Stilde_rhsD"
params = "const int flux_dirn,const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs"
preloop = pre_body
body = FD_body
loopopts = "InteriorPoints"
add_to_Cfunction_dict(
includes=includes,
desc=desc,
name=name, params=params,
preloop=pre_body, body=body, loopopts=loopopts,
path_from_rootsrcdir_to_this_Cfunc = path_from_rootsrcdir_to_this_Cfunc,
rel_path_to_Cparams=rel_path_to_Cparams)
GiRaFFE_NRPy.Stilde_flux NRPy+ Module [Back to top]¶Here, as a code validation check, we verify agreement in the SymPy expressions for the GiRaFFE evolution equations and auxiliary quantities we intend to use between
This first validation directly compares the sympy expressions. This is generally quicker and more reliable, but might overlook some complexities in implementing the C code.
all_passed=True
def comp_func(expr1,expr2,basename,prefixname2="Sf."):
if str(expr1-expr2)!="0":
print(basename+" - "+prefixname2+basename+" = "+ str(expr1-expr2))
all_passed=False
def gfnm(basename,idx1,idx2=None,idx3=None):
if idx2 is None:
return basename+"["+str(idx1)+"]"
if idx3 is None:
return basename+"["+str(idx1)+"]["+str(idx2)+"]"
return basename+"["+str(idx1)+"]["+str(idx2)+"]["+str(idx3)+"]"
# These are the standard gridfunctions we've used before.
#ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU",DIM=3)
#gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
#betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU")
#alpha = gri.register_gridfunctions("AUXEVOL",["alpha"])
#AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD",DIM=3)
#BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU",DIM=3)
# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")
# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
# We'll also need to store the results of the HLLE step between functions.
Stilde_flux_HLLED = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")
# ...and some more for the fluxes we calculate here. These three gridfunctions will each store
# the momentum flux of one component of StildeD in one direction; we'll be able to reuse them
# as we loop over each direction, reducing our memory costs.
Stilde_fluxD = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_fluxD",DIM=3)
import GiRaFFE_NRPy.Stilde_flux as Sf
for flux_dirn in range(3):
expr_list = []
exprcheck_list = []
namecheck_list = []
print("Checking the flux in direction "+str(flux_dirn))
calculate_Stilde_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi)
Sf.calculate_Stilde_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi)
for mom_comp in range(3):
namecheck_list.extend([gfnm("Stilde_fluxD",mom_comp)])
exprcheck_list.extend([Sf.Stilde_fluxD[mom_comp]])
expr_list.extend([Stilde_fluxD[mom_comp]])
for mom_comp in range(len(expr_list)):
comp_func(expr_list[mom_comp],exprcheck_list[mom_comp],namecheck_list[mom_comp])
if all_passed:
print("ALL TESTS PASSED!")
else:
print("ERROR: AT LEAST ONE TEST DID NOT PASS")
sys.exit(1)
Checking the flux in direction 0 Checking the flux in direction 1 Checking the flux in direction 2 ALL TESTS PASSED!
Our next test will generate the C files and compare the output directly. Unfortunately, we will need keep CSE sorting on for this test if we don't want false negatives, which is very slow for these functions.
import difflib
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
subdir = os.path.join("RHSs")
out_dir = os.path.join("GiRaFFE_standalone_Ccodes")
cmd.mkdir(out_dir)
cmd.mkdir(os.path.join(out_dir,subdir))
valdir = os.path.join("GiRaFFE_Ccodes_validation")
cmd.mkdir(valdir)
cmd.mkdir(os.path.join(valdir,subdir))
# These must be registered if we want to test write_cmax_cmin=True
gri.register_gridfunctions("AUXEVOL","cmax_x")
gri.register_gridfunctions("AUXEVOL","cmin_x")
gri.register_gridfunctions("AUXEVOL","cmax_y")
gri.register_gridfunctions("AUXEVOL","cmin_y")
gri.register_gridfunctions("AUXEVOL","cmax_z")
gri.register_gridfunctions("AUXEVOL","cmin_z")
generate_C_code_for_Stilde_flux(out_dir,True, alpha_face,gamma_faceDD,beta_faceU,
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi,write_cmax_cmin=True)
Sf.generate_C_code_for_Stilde_flux(valdir,True, alpha_face,gamma_faceDD,beta_faceU,
Valenciav_rU,B_rU,Valenciav_lU,B_lU,sqrt4pi,write_cmax_cmin=True)
print("Printing difference between original C code and this code...")
# Open the files to compare
files = ["calculate_Stilde_rhsD.h",
"calculate_Stilde_flux_D0.h",
"calculate_Stilde_flux_D1.h",
"calculate_Stilde_flux_D2.h"]
for file in files:
print("Checking file " + file)
with open(os.path.join(valdir,file)) as file1, open(os.path.join(out_dir,file)) as file2:
# Read the lines of each file
file1_lines = file1.readlines()
file2_lines = file2.readlines()
num_diffs = 0
for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir,file), tofile=os.path.join(out_dir,file)):
sys.stdout.writelines(line)
num_diffs = num_diffs + 1
if num_diffs == 0:
print("No difference. TEST PASSED!")
else:
print("ERROR: Disagreement found with .py file. See differences above.")
sys.exit(1)
Output C function calculate_Stilde_rhsD() to file GiRaFFE_standalone_Ccodes/calculate_Stilde_rhsD.h Output C function calculate_Stilde_rhsD() to file GiRaFFE_Ccodes_validation/calculate_Stilde_rhsD.h Printing difference between original C code and this code... Checking file calculate_Stilde_rhsD.h No difference. TEST PASSED! Checking file calculate_Stilde_flux_D0.h No difference. TEST PASSED! Checking file calculate_Stilde_flux_D1.h No difference. TEST PASSED! Checking file calculate_Stilde_flux_D2.h No difference. TEST PASSED!
This computes phase speeds in the direction given by flux_dirn. Note that we replace the full dispersion relation with a simpler one, which overestimates the maximum speeds by a factor of ~2. See full discussion around Eqs. 49 and 50 in Duez, et al.. In summary, we solve the dispersion relation (in, e.g., the $x$-direction) with a wave vector of $k_\mu = (-\omega,k_x,0,0)$. So, we solve the approximate dispersion relation $\omega_{\rm cm}^2 = [v_A^2 + c_s^2 (1-v_A^2)]k_{\rm cm}^2$ for the wave speed $\omega/k_x$, where the sound speed $c_s = \sqrt{\Gamma P/(h \rho_0)}$, the Alfvén speed $v_A = 1$ (in GRFFE), $\omega_{\rm cm} = -k_\mu k^\mu$ is the frequency in the comoving frame, $k_{\rm cm}^2 = K_\mu K^\mu$ is the wavenumber squared in the comoving frame, and $K_\mu = (g_{\mu\nu} + u_\mu u_\nu)k^\nu$ is the part of the wave vector normal to the four-velocity $u^\mu$. See below for a complete derivation.
What follows is a complete derivation of the quadratic we solve. We start from the following relations: \begin{align} w_{\rm cm} &= (-k_0 u^0 - k_x u^x) \\ k_{\rm cm}^2 &= K_{\mu} K^{\mu}, \\ K_{\mu} K^{\mu} &= (g_{\mu a} + u_{\mu} u_a) k^a g^{\mu b} [ (g_{c b} + u_c u_b) k^c ] \\ \end{align} The last term of the above can be written as follow: $$ (g_{c b} + u_{c} u_{b}) k^c = (\delta^{\mu}_c + u_c u^{\mu} ) k^c $$
Then, \begin{align} K_{\mu} K^{\mu} &= (g_{\mu a} + u_{\mu} u_a) k^a g^{\mu b} [ (g_{c b} + u_c u_b) k^c ] \\ &= (g_{\mu a} + u_{\mu} u_a) k^a (\delta^{\mu}_c + u_c u^{\mu} ) k^c \\ &=[(g_{\mu a} + u_{\mu} u_a) \delta^{\mu}_c + (g_{\mu a} + u_{\mu} u_a) u_c u^{\mu} ] k^c k^a \\ &=[(g_{c a} + u_c u_a) + (u_c u_a - u_a u_c] k^c k^a \\ &=(g_{c a} + u_c u_a) k^c k^a \\ &= k_a k^a + u^c u^a k_c k_a \\ k^a = g^{\mu a} k_{\mu} &= g^{0 a} k_0 + g^{x a} k_x \\ k_a k^a &= k_0 g^{0 0} k_0 + k_x k_0 g^{0 x} + g^{x 0} k_0 k_x + g^{x x} k_x k_x \\ &= g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2 \\ u^c u^a k_c k_a &= (u^0 k_0 + u^x k_x) (u^0 k_0 + u^x k_x) = (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 \\ (k_0 u^0)^2 + 2 k_x u^x k_0 u^0 + (k_x u^x)^2 &= v_0^2 [ (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 + g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2] \\ (1-v_0^2) (u^0 k_0 + u^x k_x)^2 &= v_0^2 (g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2) \\ (1-v_0^2) (u^0 k_0/k_x + u^x)^2 &= v_0^2 (g^{00} (k_0/k_x)^2 + 2 g^{x0} k_0/k_x + g^{xx}) \\ (1-v_0^2) (u^0 X + u^x)^2 &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\ (1-v_0^2) ((u^0)^2 X^2 + 2 u^x (u^0) X + (u^x)^2) &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\ 0 &= X^2 ( (1-v_0^2) (u^0)^2 - v_0^2 g^{00}) + X (2 u^x u^0 (1-v_0^2) - 2 v_0^2 g^{x0}) + (1-v_0^2) (u^x)^2 - v_0^2 g^{xx} \\ a &= (1-v_0^2) (u^0)^2 - v_0^2 g^{00} = (1-v_0^2) (u^0)^2 + v_0^2/\alpha^2 \leftarrow {\rm VERIFIED} \\ b &= 2 u^x u^0 (1-v_0^2) - 2 v_0^2 \beta^x/\alpha^2 \leftarrow {\rm VERIFIED,\ } X\rightarrow -X, {\rm because\ } X = -w/k_1, {\rm \ and\ we\ are\ solving\ for} -X. \\ c &= (1-v_0^2) (u^x)^2 - v_0^2 (\gamma^{xx}\psi^{-4} - (\beta^x/\alpha)^2) \leftarrow {\rm VERIFIED} \\ v_0^2 &= v_A^2 + c_s^2 (1 - v_A^2) \\ \end{align}
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-GiRaFFE_NRPy-Stilde-flux.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-GiRaFFE_NRPy-Stilde-flux",location_of_template_file=os.path.join(".."))
Created Tutorial-GiRaFFE_NRPy-Stilde-flux.tex, and compiled LaTeX file to
PDF file Tutorial-GiRaFFE_NRPy-Stilde-flux.pdf