This module is currently under development
This module is organized as follows
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__mhdflux__C = os.path.join(IGM_src_dir_path,"mhdflux.C")
In this tutorial module we will compute the flux terms for the conservative variables $U=\left\{\rho_{\star},\tilde{\tau},\tilde{S}_{i}\right\}$, which are defined in terms of the primitive variables as
$$ \left( \begin{matrix} \rho_{\star}\\ \tilde{\tau}\\ \tilde{S}_{i} \end{matrix} \right) = \left( \begin{matrix} \alpha\sqrt{\gamma}\rho_{b}u^{0}\\ \alpha^{2}\sqrt{\gamma}T^{00}-\rho_{\star}\\ \left(\rho_{\star}h + \alpha u^{0}\sqrt{\gamma}b^{2}\right)u_{i} - \alpha\sqrt{\gamma}b^{0}b_{i} \end{matrix} \right)\ . $$The flux terms for these conservative variables are
$$ \boldsymbol{F} = \left( \begin{matrix} F^{i}_{\rho_{\star}}\\ F^{i}_{\tilde{\tau}}\\ \left(F_{\tilde{S}}\right)^{j}_{\ i} \end{matrix} \right) = \left( \begin{matrix} \rho_{\star}v^{i}\\ \alpha^{2}\sqrt{\gamma}T^{0j} - \rho_{\star}v^{j}\\ \alpha\sqrt{\gamma}T^{j}_{\ i} \end{matrix} \right)\ . $$The MHD flux algorithm computes, for each of the fluxes above, the standard Harten-Lax-van Leer (HLL) flux
$$ \boxed{F^{\rm HLL} = \frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\left(U_{r} - U_{l}\right)}{c^{+} + c^{-}}}\ . $$As we go through the implementation, we will notice that will need a lot of the functions defined within the inlined_functions.C
code file of IllinoisGRMHD
. We will therefore explain the algorithm of each function in quite some detail, so that the reader can go through this tutorial without having to stop and read the inlined_functions.C
tutorial module. We do encourage the reader to go through the module as well, though, since we will be presenting the math behind each function, but not the code itself.
%%writefile $outfile_path__mhdflux__C
//-----------------------------------------------------------------------------
// Compute the flux for advecting rho_star, tau (Font's energy variable),
// and S_i .
//-----------------------------------------------------------------------------
static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th,
CCTK_REAL &cmax,CCTK_REAL &cmin,
CCTK_REAL &rho_star_flux,CCTK_REAL &tau_flux,CCTK_REAL &st_x_flux,CCTK_REAL &st_y_flux,CCTK_REAL &st_z_flux) {
CCTK_REAL psi4 = FACEVAL_LAPSE_PSI4[PSI4];
CCTK_REAL psi6 = FACEVAL_LAPSE_PSI4[PSI4]*FACEVAL_LAPSE_PSI4[PSI2];
CCTK_REAL psim4 = 1.0/(psi4);
CCTK_REAL alpha_sqrt_gamma = FACEVAL_LAPSE_PSI4[LAPSE]*psi6;
CCTK_REAL ONE_OVER_LAPSE = 1.0/FACEVAL_LAPSE_PSI4[LAPSE];
CCTK_REAL ONE_OVER_LAPSE_SQUARED=SQR(ONE_OVER_LAPSE);
Writing ../src/mhdflux.C
Next we compute the quantities $P_{\rm cold}$, $\epsilon_{\rm cold}$, $\frac{dP_{\rm cold}}{d\rho}$, $\epsilon_{\rm th}$, $h$, and $\Gamma_{\rm cold}$. We have written $\rho_{b}\equiv\rho$ so that the discussion contains a notation that is slightly less cluttered with indices.
The function compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()
is defined within the inlined_functions.C
code file of IllinoisGRMHD
. It currently implements 3 different cases:
Case 1: $\boldsymbol{\rho_{b} = 0}$
In this case, we set $\Gamma_{\rm cold} = \Gamma_{\rm tab}$, i.e. its tabulated input value, and all other quantities to zero: $P_{\rm cold}=\epsilon_{\rm cold}=\frac{dP_{\rm cold}}{d\rho}=h=\epsilon_{\rm th}=0$
Case 2: Single Polytrope EOS
In this case, we have the relations:
$$ \boxed{ \begin{align} P_{\rm cold} &= \kappa \rho^{\Gamma_{\rm cold}}\ ,\\ \epsilon_{\rm cold} &= \left.\left(\frac{P_{\rm cold}}{\rho}\right)\middle/\left(\Gamma_{\rm cold}-1\right)\right.\ ,\\ \frac{dP_{\rm cold}}{d\rho} &= \frac{\Gamma_{\rm cold}P_{\rm cold}}{\rho}\ ,\\ \epsilon_{\rm th} &= \left.\left(\frac{P-P_{\rm cold}}{\rho}\right)\middle/\left(\Gamma_{\rm cold}-1\right)\right.\ ,\\ h &= 1+\epsilon_{\rm cold}+\epsilon_{\rm th} + \frac{P}{\rho}\ ,\\ \Gamma_{\rm cold} &= \Gamma_{\rm tab}\ . \end{align}} $$Case 3: Piecewise Polytrope EOS
We now consider a set of dividing densities $\rho_{\min} < \rho_{1} < \cdots < \rho_{n-1} < \rho_{\max}$ such that the pressure and energy density are everywhere continuous. In each interval, the EOS is a single polytrope, with its own $K_{i}$ and $\left(\Gamma_{cold}\right)_{i}\equiv\Gamma_{i}$, so that
$$ \boxed{ P_{\rm cold} = \left\{ \begin{matrix} K_{0}\rho^{\Gamma_{0}} & , & \rho \leq \rho_{0}\\ K_{1}\rho^{\Gamma_{1}} & , & \rho_{0} \leq \rho \leq \rho_{1}\\ \vdots & & \vdots\\ K_{j}\rho^{\Gamma_{j}} & , & \rho_{j-1} \leq \rho \leq \rho_{j}\\ \vdots & & \vdots\\ K_{N-1}\rho^{\Gamma_{N-1}} & , & \rho_{N-2} \leq \rho \leq \rho_{N-1}\\ K_{N}\rho^{\Gamma_{N}} & , & \rho \geq \rho_{N-1} \end{matrix} \right. }\ . $$Then, for each set $\left\{P_{i},K_{i},\Gamma_{i}\right\}$ the desired quantities are computed using the same equations used in Case 2.
%%writefile -a $outfile_path__mhdflux__C
// First compute P_{cold}, \epsilon_{cold}, dP_{cold}/drho, \epsilon_{th}, h, and \Gamma_{cold},
// for right and left faces:
CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,Gamma_coldr;
compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ur,eos,Gamma_th,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,Gamma_coldr);
CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,Gamma_coldl;
compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ul,eos,Gamma_th,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,Gamma_coldl);
Appending to ../src/mhdflux.C
We now call upon the impose_speed_limit_output_u0()
function inside the inlined_functions.C
code file of IllinoisGRMHD
. The basic algorithm performed by this function is summarized here. We start by evaluating the quantity
where when going from line 1 to 2 and from line 3 to 4 we have used Eqs. (53) and (56) from Duez et al., respectively. Then we construct the "speed limit quantity"
$$ {\rm ONE\_MINUS\_ONE\_OVER\_GAMMA\_SPEED\_LIMIT\_SQUARED} \equiv B = 1-\frac{1}{\gamma^{2}_{\rm speed\ limit}}\ . $$If $A > B$, then we construct the correction factor $C\equiv A / B$, and adjust the velocities using
$$ \boxed{v^{i} \to \left(v^{i}+\beta^{i}\right)C - \beta^{i}}\ . $$Finally, since $A$ is evaluated using the first line above, namely
$$ A = \gamma_{ij}\left(\frac{v^{i}+\beta^{i}}{\alpha}\right)\left(\frac{v^{j}+\beta^{j}}{\alpha}\right)\ , $$but we know, from the last line how $A$ and $u^{0}$ are related, we compute
$$ \boxed{u^{0} = \frac{1}{\alpha\sqrt{1-A}}}\ . $$%%writefile -a $outfile_path__mhdflux__C
//Compute face velocities
// Begin by computing u0
output_stats stats; stats.failure_checker=0;
CCTK_REAL u0_r,u0_l;
impose_speed_limit_output_u0(FACEVAL,Ur,psi4,ONE_OVER_LAPSE,stats,u0_r);
impose_speed_limit_output_u0(FACEVAL,Ul,psi4,ONE_OVER_LAPSE,stats,u0_l);
Appending to ../src/mhdflux.C
Now we use the function compute_smallba_b2_and_u_i_over_u0_psi4()
from the inlined_functions.C
code file of IllinoisGRMHD
.
We will need the following identities
$$ \begin{align} v^{i} &= \frac{u^{i}}{u^{0}}\ ,\\ B^{0}_{(u)} &= \frac{u_{i}B^{i}}{\alpha}\ ,\\ B^{i}_{(u)} &= \frac{1}{u^{0}}\left(\frac{B^{i}}{\alpha} + u^{i}B^{0}_{(u)}\right)\ ,\\ b^{\mu} &= \frac{B^{\mu}_{(u)}}{\sqrt{4\pi}}\ . \end{align} $$We start by setting the relation
$$ b^{0} = \frac{u_{i}B^{i}}{\alpha\sqrt{4\pi}} \implies \boxed{\alpha\sqrt{4\pi}b^{0} = u_{i}B^{i}}\ . $$Then we compute
$$ \begin{align} b^{i} &= \frac{B^{i}_{(u)}}{\sqrt{4\pi}}\\ &= \frac{1}{u^{0}\sqrt{4\pi}}\left(\frac{B^{i}}{\alpha} + B^{0}_{(u)}u^{i}\right)\\ &= \frac{1}{u^{0}\sqrt{4\pi}}\left(\frac{B^{i}}{\alpha} + \sqrt{4\pi}b^{0}u^{i}\right)\\ &= \frac{1}{\alpha\sqrt{4\pi}}\left(\frac{B^{i}}{u^{0}} + \alpha\sqrt{4\pi}b^{0}\frac{u^{i}}{u^{0}}\right)\\ \implies &\boxed{b^{i} = \frac{1}{\alpha\sqrt{4\pi}}\left(\frac{B^{i}}{u^{0}} + \alpha\sqrt{4\pi}b^{0}v^{i}\right)}\ . \end{align} $$Finally, we compute
$$ \begin{align} b^{2} &= g_{\mu\nu}b^{\mu}b^{\nu}\\ &= g_{00}\left(b^{0}\right)^{2} + g_{ij}b^{i}b^{j} + 2g_{0i}b^{0}b^{i}\\ &= \left(-\alpha^{2} + \gamma_{ij}\beta^{i}\beta^{j}\right)\left(b^{0}\right)^{2} + \gamma_{ij}b^{i}b^{j} + 2b^{0}\gamma_{ij}\beta^{j}b^{i}\\ &= -\left(\alpha b^{0}\right)^{2} + \gamma_{ij}\left[b^{i}b^{j} + 2b^{0}b^{i}\beta^{j} + \left(b^{0}\right)^{2}\beta^{i}\beta^{j}\right]\\ \implies &\boxed{b^{2} = -\left(\alpha b^{0}\right)^{2} + \gamma_{ij}\left(b^{i} + b^{0}\beta^{i}\right)\left(b^{j} + b^{0}\beta^{j}\right)} \end{align} $$%%writefile -a $outfile_path__mhdflux__C
//Next compute b^{\mu}, the magnetic field measured in the comoving fluid frame:
CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI = ONE_OVER_LAPSE*ONE_OVER_SQRT_4PI;
/***********************************************************/
/********** RIGHT FACE ************/
// Note that smallbr[4] = b^a defined in Gammie's paper, on the right face.
CCTK_REAL u_x_over_u0_psi4r,u_y_over_u0_psi4r,u_z_over_u0_psi4r;
CCTK_REAL smallbr[NUMVARS_SMALLB];
// Compute b^{a}, b^2, and u_i over u^0
compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI,
u_x_over_u0_psi4r,u_y_over_u0_psi4r,u_z_over_u0_psi4r,smallbr);
// Then compute u_xr,u_yr, and u_zr. We need to set the zeroth component so we can specify U_LOWER{r,l}[{UX,UY,UZ}] (UX=1,UY=2,UZ=3).
CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4],
u_z_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4] };
/********** LEFT FACE ************/
// Note that smallbl[4] = b^a defined in Gammie's paper, on the left face.
CCTK_REAL u_x_over_u0_psi4l,u_y_over_u0_psi4l,u_z_over_u0_psi4l;
CCTK_REAL smallbl[NUMVARS_SMALLB];
// Compute b^{a}, b^2, and u_i over u^0
compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI,
u_x_over_u0_psi4l,u_y_over_u0_psi4l,u_z_over_u0_psi4l,smallbl);
// Then compute u_xr,u_yr, and u_zr. We need to set the zeroth component so we can specify U_LOWER{r,l}[{UX,UY,UZ}]
CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4],
u_z_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4] };
/***********************************************************/
Appending to ../src/mhdflux.C
We will now explain two different functions contained in the inlined_functions.C
code file of IllinoisGRMHD
. These functions are compute_v02()
and find_cp_cm
. These functions are called with the objective of computing the minimum ($-$) and maximum ($+$) characteristic speeds at each cell interface, $c_{\pm}^{r,l}$.
We approximate the general GRMHD dispersion relation (eq. 27 of Gammie & McKinney (2003)) by the simpler expression
$$ \omega_{\rm cm}^{2} = \left[v_{\rm A}^{2} + c_{\rm s}^{2}\left(1-v_{\rm A}^{2}\right)\right]k_{\rm cm}^{2}\ , $$where $\omega_{\rm cm}=-k_{\mu}u^{\mu}$ is the frequency and $k_{\rm cm}^{2} = K_{\mu}K^{\mu}$ the wavenumber of an MHD wave mode in the frame comoving with the fluid, where $K_{\mu}$ is defined as the projection of the wave vector $k^{\nu}$ onto the direction normal to $u^{\nu}$: $K_{\mu} = \left(g_{\mu\nu}+u_{\mu}u_{\nu}\right)k^{\nu}$. $c_{\rm s}$ is the sound speed, and $v_{\rm A}$ is the Alfvén speed, given by
$$ v_{\rm A} = \sqrt{\frac{b^{2}}{\rho_{b}h + b^{2}}}\ . $$With these definitions, we may then solve the approximate dispersion relation above along direction $i$, noting that in the comoving frame $k_{\mu} = \left(-\omega,k_{j}\delta^{j}_{\ i}\right)$ and the wave (phase) velocity is $c_{\pm} = \left.\omega\middle/\left(k_{j}\delta^{j}_{\ i}\right)\right.$. The dispersion can then be written as a quadratic equation for $c_{\pm}$:
$$ ac_{\pm}^{2} + bc_{\pm} + c = 0\ , $$with
$$ \boxed{ \begin{align} a &= \left(1-v_{0}^{2}\right)\left(u^{0}\right)^{2} - v_{0}^{2}g^{00}\ ,\\ b &= 2v_{0}^{2}g^{i0} - 2u^{i}u^{0}\left(1-v^{2}_{0}\right)\ ,\\ c &= \left(1-v_{0}^{2}\right)\left(u^{i}\right)^{2} - v_{0}^{2}g^{ii}\ ,\\ v_{0}^{2} &= v_{\rm A}^{2} + c_{\rm s}^{2}\left(1-v_{\rm A}^{2}\right)\ ,\\ c_{\rm s} &= \left.\left[\frac{dP_{\rm cold}}{d\rho_{b}} + \Gamma_{\rm th}\left(\Gamma_{\rm th}-1\right)\epsilon_{\rm th}\right]\middle/h\right.\ ,\\ c_{+} &= \max\left(\frac{-b \pm \sqrt{b^{2}-4ac}}{2a}\right)\ ,\\ c_{-} &= \min\left(\frac{-b \pm \sqrt{b^{2}-4ac}}{2a}\right)\ . \end{align} } $$Finally, after the maximum and minimum speeds $c_{\pm}$ have been computed at left and right facs, the standard Harten-Lax-van Leer (HLL), approximate Riemann solve is applied to compute fluxes for the three conservative variables $U = \left\{\rho_{\star},\tilde{\tau},\tilde{S}_{i}\right\}$:
$$ F^{\rm HLL} = \frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\left(U_{r} - U_{l}\right)}{c^{+} + c^{-}}\ , $$where
$$ \boxed{c^{\pm} = \pm\max\left(0,c^{r}_{\pm},c^{l}_{\pm}\right)}\ . $$%%writefile -a $outfile_path__mhdflux__C
// Compute v02 = v_A^2 + c_s^2*(1.0-v_A^2), where c_s = sound speed, and v_A = Alfven velocity
CCTK_REAL v02r,v02l;
// First right face
compute_v02(dPcold_drhor,Gamma_th,eps_thr,h_r,smallbr,Ur,v02r);
// Then left face.
compute_v02(dPcold_drhol,Gamma_th,eps_thl,h_l,smallbl,Ul,v02l);
int offset=flux_dirn-1;
// Now that we have computed v02 = v_A^2 + c_s^2*(1.0-v_A^2), we can
// next compute c_+ and c_- using a simplified dispersion relation.
// Note that, in solving the simplified disp. relation, we overestimate
// c_+ and c_- by around a factor of 2, making the MHD evolution more
// diffusive (and potentially more *stable*) than it could be.
CCTK_REAL cplusr,cminusr,cplusl,cminusl;
find_cp_cm(cplusr,cminusr,v02r,u0_r,
Ur[VX+offset],ONE_OVER_LAPSE_SQUARED,FACEVAL[SHIFTX+offset],psim4,FACEVAL[GUPXX+offset]);
find_cp_cm(cplusl,cminusl,v02l,u0_l,
Ul[VX+offset],ONE_OVER_LAPSE_SQUARED,FACEVAL[SHIFTX+offset],psim4,FACEVAL[GUPXX+offset]);
// Then compute cmax, cmin. This is required for the HLL flux.
CCTK_REAL cmaxL = MAX(0.0,MAX(cplusl,cplusr));
CCTK_REAL cminL = -MIN(0.0,MIN(cminusl,cminusr));
Appending to ../src/mhdflux.C
We now evaluate the flux for $U=\rho_{\star}$. First, remember that
$$ \rho_{\star} = \alpha\sqrt{\gamma}\rho_{b}u^{0}\ . $$$$ F^{i}_{\rho_{\star}} = \rho_{\star} v^{i}\ , $$where $i$ is the current flux direction. We can then evaluate the HLL flux in the $i$-direction,
$$ \boxed{F^{\rm HLL} = \frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\left(U_{r} - U_{l}\right)}{c^{+} + c^{-}}}\ . $$%%writefile -a $outfile_path__mhdflux__C
//*********************************************************************
// density flux = \rho_* v^m, where m is the current flux direction (the m index)
//*********************************************************************
CCTK_REAL rho_star_r = alpha_sqrt_gamma*Ur[RHOB]*u0_r;
CCTK_REAL rho_star_l = alpha_sqrt_gamma*Ul[RHOB]*u0_l;
CCTK_REAL Fr = rho_star_r*Ur[VX+offset]; // flux_dirn = 2, so offset = 1, implies Ur[VX] -> Ur[VY]
CCTK_REAL Fl = rho_star_l*Ul[VX+offset]; // flux_dirn = 2, so offset = 1, implies Ul[VX] -> Ul[VY]
// HLL step for rho_star:
rho_star_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(rho_star_r-rho_star_l) )/(cmaxL + cminL);
Appending to ../src/mhdflux.C
We then evaluate the flux for $U=\tilde{\tau}$. First, let us remember that
$$ \tilde{\tau} = \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star}\ . $$We then have the flux term
$$ F^{i}_{\tilde{\tau}} = \alpha^{2}\sqrt{\gamma}T^{0j} - \rho_{\star}v^{j}\ , $$where $i$ is the current flux direction. We can then evaluate the HLL flux in the $i$-direction,
$$ \boxed{F^{\rm HLL} = \frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\left(U_{r} - U_{l}\right)}{c^{+} + c^{-}}}\ . $$%%writefile -a $outfile_path__mhdflux__C
//*********************************************************************
// energy flux = \alpha^2 \sqrt{\gamma} T^{0m} - \rho_* v^m, where m is the current flux direction (the m index)
//*********************************************************************
// First compute some useful metric quantities:
CCTK_REAL alpha_squared_sqrt_gamma = FACEVAL_LAPSE_PSI4[LAPSE]*alpha_sqrt_gamma;
CCTK_REAL g4uptm = ONE_OVER_LAPSE_SQUARED*FACEVAL[SHIFTX+offset];
CCTK_REAL g4uptt = -ONE_OVER_LAPSE_SQUARED;
/********** RIGHT FACE ************/
// Compute a couple useful hydro quantities:
CCTK_REAL rho0_h_plus_b2_r = (Ur[RHOB]*h_r + smallbr[SMALLB2]);
CCTK_REAL P_plus_half_b2_r = (Ur[PRESSURE]+0.5*smallbr[SMALLB2]);
// Then compute T^{0m} and the flux:
CCTK_REAL TUP0m_r = rho0_h_plus_b2_r*SQR(u0_r)*Ur[VX+offset] + P_plus_half_b2_r*g4uptm - smallbr[SMALLBT]*smallbr[SMALLBX+offset];
Fr = alpha_squared_sqrt_gamma * TUP0m_r - rho_star_r * Ur[VX+offset];
// Finally compute tau
CCTK_REAL TUP00_r = rho0_h_plus_b2_r*u0_r*u0_r + P_plus_half_b2_r*g4uptt - smallbr[SMALLBT]*smallbr[SMALLBT];
CCTK_REAL tau_r = alpha_squared_sqrt_gamma * TUP00_r - rho_star_r;
/********** LEFT FACE *************/
// Compute a couple useful hydro quantities:
CCTK_REAL rho0_h_plus_b2_l = (Ul[RHOB]*h_l + smallbl[SMALLB2]);
CCTK_REAL P_plus_half_b2_l = (Ul[PRESSURE]+0.5*smallbl[SMALLB2]);
// Then compute T^{0m} and the flux:
CCTK_REAL TUP0m_l = rho0_h_plus_b2_l*SQR(u0_l)*Ul[VX+offset] + P_plus_half_b2_l*g4uptm - smallbl[SMALLBT]*smallbl[SMALLBX+offset];
Fl = alpha_squared_sqrt_gamma * TUP0m_l - rho_star_l * Ul[VX+offset];
// Finally compute tau
CCTK_REAL TUP00_l = rho0_h_plus_b2_l*u0_l*u0_l + P_plus_half_b2_l*g4uptt - smallbl[SMALLBT]*smallbl[SMALLBT];
CCTK_REAL tau_l = alpha_squared_sqrt_gamma * TUP00_l - rho_star_l;
// HLL step for tau:
tau_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(tau_r-tau_l) )/(cmaxL + cminL);
Appending to ../src/mhdflux.C
We then evaluate the flux for $U=\tilde{S}_{i}$. First, let us remember that
$$ \tilde{S}_{i} = \left(\rho_{\star} h + \alpha u^{0} \sqrt{\gamma}b^{2}\right)u_{i} - \alpha\sqrt{\gamma}b^{0}b_{i}\ . $$We then have the flux term
$$ \left(F_{\tilde{S}}\right)^{j}_{\ i} = \alpha\sqrt{\gamma}T^{j}_{\ i}\ , $$where $j$ is the current flux direction and $i$ correspond to the component of $\tilde{S}_{i}$ we are interested in. We can then evaluate the HLL flux in the $i$-direction,
$$ \boxed{F^{\rm HLL} = \frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\left(U_{r} - U_{l}\right)}{c^{+} + c^{-}}}\ . $$%%writefile -a $outfile_path__mhdflux__C
//*********************************************************************
// momentum flux = \alpha \sqrt{\gamma} T^m_j, where m is the current flux direction (the m index)
//*********************************************************************
// b_j = g_{ij} (b^i + b^t shift^i), g_{ij} = physical metric
//CCTK_REAL sbtr=0,sbtl=0;
CCTK_REAL smallb_lowerr[NUMVARS_SMALLB],smallb_lowerl[NUMVARS_SMALLB];
lower_4vector_output_spatial_part(psi4,FACEVAL,smallbr,smallb_lowerr);
lower_4vector_output_spatial_part(psi4,FACEVAL,smallbl,smallb_lowerl);
/********** Flux for S_x **********/
// [S_x flux] = \alpha \sqrt{\gamma} T^m_x, where m is the current flux direction (the m index)
// Again, offset = 0 for reconstruction in x direction, 1 for y, and 2 for z
// Note that kronecker_delta[flux_dirn][0] = { 1 if flux_dirn==1, 0 otherwise }.
Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX]
+ P_plus_half_b2_r*kronecker_delta[flux_dirn][0] - smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBX] );
Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX]
+ P_plus_half_b2_l*kronecker_delta[flux_dirn][0] - smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBX] );
// S_x =\alpha\sqrt{\gamma}( T^0_x )
CCTK_REAL st_x_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UX] - smallbr[SMALLBT]*smallb_lowerr[SMALLBX] );
CCTK_REAL st_x_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UX] - smallbl[SMALLBT]*smallb_lowerl[SMALLBX] );
// HLL step for Sx:
st_x_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_x_r-st_x_l) )/(cmaxL + cminL);
/********** Flux for S_y **********/
// [S_y flux] = \alpha \sqrt{\gamma} T^m_y, where m is the current flux direction (the m index)
// Again, offset = 1 for reconstruction in x direction, 2 for y, and 3 for z
// Note that kronecker_delta[flux_dirn][1] = { 1 if flux_dirn==2, 0 otherwise }.
Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1]
- smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBY] );
Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1]
- smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBY] );
// S_y =\alpha\sqrt{\gamma}( T^0_y )
CCTK_REAL st_y_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UY] - smallbr[SMALLBT]*smallb_lowerr[SMALLBY] );
CCTK_REAL st_y_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UY] - smallbl[SMALLBT]*smallb_lowerl[SMALLBY] );
// HLL step for Sy:
st_y_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_y_r-st_y_l) )/(cmaxL + cminL);
/********** Flux for S_z **********/
// [S_z flux] = \alpha \sqrt{\gamma} T^m_z, where m is the current flux direction (the m index)
// Again, offset = 1 for reconstruction in x direction, 2 for y, and 3 for z
// Note that kronecker_delta[flux_dirn][2] = { 1 if flux_dirn==3, 0 otherwise }.
Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2]
- smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBZ] );
Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2]
- smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBZ] );
// S_z =\alpha\sqrt{\gamma}( T^0_z )
CCTK_REAL st_z_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UZ] - smallbr[SMALLBT]*smallb_lowerr[SMALLBZ] );
CCTK_REAL st_z_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UZ] - smallbl[SMALLBT]*smallb_lowerl[SMALLBZ] );
// HLL step for Sz:
st_z_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_z_r-st_z_l) )/(cmaxL + cminL);
cmax = cmaxL;
cmin = cminL;
}
Appending to ../src/mhdflux.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/mhdflux.C"
original_IGM_file_name = "mhdflux-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__mhdflux__C = !diff $original_IGM_file_path $outfile_path__mhdflux__C
if Validation__mhdflux__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 mhdflux.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for mhdflux.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__mhdflux__C:
print(diff_line)
Validation test for mhdflux.C: FAILED! Diff: 2c2 < // Compute the flux for advecting rho_star, tau (Font's energy variable), --- > // Compute the flux for advecting rho_star, tau (Font's energy variable), 5c5 < static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, --- > static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th, 9d8 < 17a17 > 20,23c20,24 < CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,gamma_coldr; < compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(Ur,eos,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,gamma_coldr); < CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,gamma_coldl; < compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(Ul,eos,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,gamma_coldl); --- > CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,Gamma_coldr; > compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ur,eos,Gamma_th,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,Gamma_coldr); > CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,Gamma_coldl; > compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ul,eos,Gamma_th,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,Gamma_coldl); > 31a33 > 40c42 < compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI, --- > compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI, 43c45 < CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], --- > CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], 50c52 < compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI, --- > compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI, 53c55 < CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], --- > CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], 56a59 > 60c63 < compute_v02(dPcold_drhor,eos.gamma_th,eps_thr,h_r,smallbr,Ur,v02r); --- > compute_v02(dPcold_drhor,Gamma_th,eps_thr,h_r,smallbr,Ur,v02r); 62c65 < compute_v02(dPcold_drhol,eos.gamma_th,eps_thl,h_l,smallbl,Ul,v02l); --- > compute_v02(dPcold_drhol,Gamma_th,eps_thl,h_l,smallbl,Ul,v02l); 67,68c70,71 < // next compute c_+ and c_- using a simplified dispersion relation. < // Note that, in solving the simplified disp. relation, we overestimate --- > // next compute c_+ and c_- using a simplified dispersion relation. > // Note that, in solving the simplified disp. relation, we overestimate 80a84 > 91a96 > 122a128 > 136c142 < Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX] --- > Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX] 138c144 < Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX] --- > Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX] 152c158 < Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1] --- > Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1] 154c160 < Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1] --- > Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1] 168c174 < Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2] --- > Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2] 170c176 < Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2] --- > Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2] 182a189 >
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__mhdflux.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__mhdflux.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex
!rm -f Tut*.out Tut*.aux Tut*.log