Notebook

Tutorial-IllinoisGRMHD: mhdflux.C¶

Authors: Leo Werneck & Zach Etienne¶

This module is currently under development

In this tutorial module we explain how the MHD fluxes are evaluated within IllinoisGRMHD.¶

Required and recommended citations:¶

(Required) Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. (arxiv:1501.07276).
(Required) Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) (astro-ph/0512420).
(Recommended) Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 (astro-ph/0210618).

Table of Contents¶

$\label{toc}$

This module is organized as follows

Step 0: Source directory creation
Step 1: Introduction
Step 2: Basic quantities
Step 3: Equation of State (EOS) quantities
Step 4: Impose a speed limit and determine $u^{0}$
Step 5: The magnetic field measured in the comoving fluid frame, $b^{\mu}$
Step 6: The minimum and maximum characteristic speeds
Step 7: MHD flux: $\rho_{\star}$
Step 8: MHD flux: $\tilde{\tau}$
Step 9: MHD flux: $\tilde{S}_{i}$
Step 10: Code validation
Step 11: Output this notebook to $\LaTeX$ -formatted PDF file

Step 0: Source directory creation [Back to top]¶

$\label{src_dir}$

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.

In [1]:

# 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")

Step 1: Introduction [Back to top]¶

$\label{introduction}$

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.

Step 2: Basic quantities [Back to top]¶

$\label{basic_quantities}$

We start by defining useful quantities, such as $\psi^{\pm4}$ , $\psi^{6}$ , $\alpha\sqrt{\gamma}$ , $\alpha^{-1}$ , and $\alpha^{-2}$ .

In [2]:

%%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

Step 3: Equation of State (EOS) quantities [Back to top]¶

$\label{eos_quantities}$

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 constains 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 implementes 3 different cases:

Case 1: $\boldsymbol{\rho_{b} = 0}$

In this case, we set $\Gamma_{\rm cold} = \Gamma_{\rm tab}$ , i.e. it's 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 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.

In [3]:

%%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

Step 4: Impose a speed limit and determine $u^{0}$ [Back to top]¶

$\label{impose_speed_limit_output_u0}$

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

$\begin{align} {\rm one\_minus\_one\_over\_alpha\_u0\_squared} \equiv A &= \gamma_{ij}\left(\frac{v^{i}+\beta^{i}}{\alpha}\right)\left(\frac{v^{j}+\beta^{j}}{\alpha}\right)\\ &= \frac{\gamma_{ij}}{\alpha^{2}}\left[\frac{\gamma^{ik}u_{k}}{u^{0}} - \beta^{i} + \beta^{i}\right]\left[\frac{\gamma^{j\ell}u_{\ell}}{u^{0}} - \beta^{j} + \beta^{j}\right]\\ &=\frac{\gamma_{ij}u^{i}u^{j}}{\left(\alpha u^{0}\right)^{2}}\\ &=\frac{\left(\alpha u^{0}\right)^{2}-1}{\left(\alpha u^{0}\right)^{2}}\\ &=1 - \frac{1}{\left(\alpha u^{0}\right)^{2}}\ \\ \implies \boxed{A = 1 - \frac{1}{\left(\alpha u^{0}\right)^{2}}}\ , \end{align}$

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}}}\ .$

In [4]:

%%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

Step 5: The magnetic field measured in the comoving fluid frame, $b^{\mu}$ [Back to top]¶

$\label{magnetic_field_comoving_frame}$

Now we use th 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}$

In [5]:

%%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

Step 6: The minimum and maximum characteristic speeds [Back to top]¶

$\label{min_max_characteristic_speeds}$

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)}\ .$

In [6]:

%%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

Step 7: MHD flux: $\rho_{\star}$ [Back to top]¶

$\label{rho_star_flux}$

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^{-}}}\ .$

In [7]:

%%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

Step 8: MHD flux: $\tilde{\tau}$ [Back to top]¶

$\label{energy_flux}$

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^{-}}}\ .$

In [8]:

%%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

Step 9: MHD flux: $\tilde{S}_{i}$ [Back to top]¶

$\label{momentum_flux}$

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^{-}}}\ .$

In [9]:

%%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

Step 10: Code validation [Back to top]¶

$\label{code_validation}$

First we download the original IllinoisGRMHD source code and then compare it to the source code generated by this tutorial notebook.

In [10]:

# 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
>

Step 11: Output this notebook to $\LaTeX$ -formatted PDF file [Back to top]¶

$\label{latex_pdf_output}$

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).

In [11]:

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

Tutorial-IllinoisGRMHD: mhdflux.C¶

Authors: Leo Werneck & Zach Etienne¶

In this tutorial module we explain how the MHD fluxes are evaluated within IllinoisGRMHD.¶

Required and recommended citations:¶

Table of Contents¶

Step 0: Source directory creation [Back to top]¶

Step 1: Introduction [Back to top]¶

Step 2: Basic quantities [Back to top]¶

Step 3: Equation of State (EOS) quantities [Back to top]¶

Step 4: Impose a speed limit and determine u0u^{0} [Back to top]¶

Step 5: The magnetic field measured in the comoving fluid frame, bμb^{\mu} [Back to top]¶

Step 6: The minimum and maximum characteristic speeds [Back to top]¶

Step 7: MHD flux: ρ⋆\rho_{\star} [Back to top]¶

Step 8: MHD flux: ˜τ\tilde{\tau} [Back to top]¶

Step 9: MHD flux: ˜Si\tilde{S}_{i} [Back to top]¶

Step 10: Code validation [Back to top]¶

Step 11: Output this notebook to LATEX\LaTeX-formatted PDF file [Back to top]¶

Step 4: Impose a speed limit and determine $u^{0}$ [Back to top]¶

Step 5: The magnetic field measured in the comoving fluid frame, $b^{\mu}$ [Back to top]¶

Step 7: MHD flux: $\rho_{\star}$ [Back to top]¶

Step 8: MHD flux: $\tilde{\tau}$ [Back to top]¶

Step 9: MHD flux: $\tilde{S}_{i}$ [Back to top]¶

Step 11: Output this notebook to $\LaTeX$ -formatted PDF file [Back to top]¶