Notebook

Tutorial-IllinoisGRMHD: `driver_evaluate_MHD_rhs.C`¶

Authors: Leo Werneck & Zach Etienne¶

This module is currently under development

In this tutorial module we explain the driver functions that compute the right-hand side (RHS) of the MHD equations 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: Load up necessary ETK and IllinoisGRMHD header files
Step 3: The IllinoisGRMHD_driver_evaluate_MHD_rhs() function
1. Step 3.1: Equation of state (EOS) parameters
2. Step 3.2: Set pointers to GRMHD gridfunctions
3. Step 3.3: Convert ADM variables to BSSN variables
4. Step 3.4: Setting up pointers to the metric and stress-energy tensor gridfunctions
5. Step 3.5: Initialization of the RHS variables
6. Step 3.6: Compute extrinsic curvature terms from the RHS of $\partial_{t}\tilde\tau$ and $T^{\mu\nu}$
7. Step 3.7: Computing ${\rm ftilde}$
8. Step 3.8: The RHSs of $\rho_{\star}$ , $\tilde\tau$ , $\tilde{S}_{i}$ , and $A_{i}$
  1. Step 3.8.1: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $x$ -direction
  2. Step 3.8.2: Evaluating $\partial_{x}\boldsymbol{F}$
  3. Step 3.8.3: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $y$ -direction
  4. Step 3.8.4: Evaluating $\partial_{y}\boldsymbol{F}$
  5. Step 3.8.5: Evaluating $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$
  6. Step 3.8.6: Multiple reconstructions
  7. Step 3.8.7: Evaluating $\partial_{z}\boldsymbol{F}$
  8. Step 3.8.8: Evaluating $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$
  9. Step 3.8.9: Evaluating $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$
  10. Step 3.8.10: Evaluating $\partial_{t}\left[\psi^{6}\Phi\right]$ and $\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms}$
Step 4: The driver_evaluate_MHD_rhs.h header file
Step 5: Code validation
1. Step 5.a: driver_evaluate_MHD_rhs.C
2. Step 5.b: driver_evaluate_MHD_rhs.h
Step 6: 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__driver_evaluate_MHD_rhs__C = os.path.join(IGM_src_dir_path,"driver_evaluate_MHD_rhs.C")
outfile_path__driver_evaluate_MHD_rhs__h = os.path.join(IGM_src_dir_path,"driver_evaluate_MHD_rhs.h")

Step 1: Introduction [Back to top]¶

$\label{introduction}$

We will start by creating the file driver_evaluate_MHD_rhs.C and writing dwn the preamble of the file, which contains useful references and information to the user.

We remind the reader of the "generalized Lorenz gauge condition",

$\nabla_{\mu}\mathcal{A^{\mu}} = \xi n_{\mu}\mathcal{A^{\mu}}\ ,$

where $n_{\mu} = \left(\alpha,0,0,0\right)$ is the normal unit vector, $\mathcal{A}_{\mu}$ is the magnetic 4-vector potential, and $xi$ is a parameter with dimensions 1/Length, just like the $\eta$ parameter in the gamma-driving shift condition, so its value is chosen so that the CFL condition remains satisfied.

In [2]:

%%writefile $outfile_path__driver_evaluate_MHD_rhs__C
/*********************************************
 * Evaluate RHS of GRMHD & induction equations
 * (vector potential prescription), using the
 * generalized Lorenz gauge condition for the
 * EM gauge.
 *
 * Based originally on the Illinois GRMHD code,
 * written by Matt Duez, Yuk Tung Liu, and Branson
 * Stephens (original version), and then developed
 * primarily by Zachariah Etienne, Yuk Tung Liu,
 * and Vasileios Paschalidis.
 *
 * Rewritten for public release in 2013
 *      by Zachariah B. Etienne
 *
 * References:
 * Original unigrid GRMHD evolution prescription:
 *    http://arxiv.org/abs/astro-ph/0503420
 * Vector potential formulation in full GR:
 *    http://arxiv.org/abs/1007.2848
 * Improved EM gauge conditions for AMR grids:
 *    http://arxiv.org/abs/1110.4633
 * Generalized Lorenz gauge prescription:
 *    http://arxiv.org/abs/1207.3354
 *
 * Note that the Generalized Lorenz gauge strength
 *  parameter has units of 1/M, just like the \eta
 *  parameter in the gamma-driving shift condition,
 *  so setting it too large will result in violation
 *  of the CFL condition.
 *
 * This version of PPM implements the standard
 * Colella & Woodward PPM, though modified as in GRHydro
 * to have 3 ghostzones instead of 4.
 *********************************************/

Writing ../src/driver_evaluate_MHD_rhs.C

Step 2: Load up necessary ETK and `IllinoisGRMHD` header files [Back to top]¶

$\label{header_files}$

Here we load all necessary ETK and IllinoisGRMHD file, as well as some standard C++ libraries.

In [3]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


#include "cctk.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <sys/time.h>
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "IllinoisGRMHD_headers.h" /* Generic #define's and function prototypes */
#include "driver_evaluate_MHD_rhs.h" /* Function prototypes for this file only */
#include "IllinoisGRMHD_EoS_lowlevel_functs.C"
#include "inlined_functions.C"

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3: The `IllinoisGRMHD_driver_evaluate_MHD_rhs()` function [Back to top]¶

$\label{driver_mhd_rhs_function}$

This is the basic function declaration. We set up basic ETK parameters and verify double precision is being used.

In [4]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


extern "C" void IllinoisGRMHD_driver_evaluate_MHD_rhs(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  int levelnumber = GetRefinementLevel(cctkGH);

  if(CCTK_Equals(verbose, "essential+iteration output")) {
    CCTK_VInfo(CCTK_THORNSTRING,"***** Iter. # %d, Lev: %d, Integrating to time: %e *****",cctk_iteration,levelnumber,cctk_delta_time/cctk_levfac[0]+cctk_time);
  }

  if( sizeof(CCTK_REAL) < 8 ) CCTK_VError(VERR_DEF_PARAMS,"Error: IllinoisGRMHD assumes that CCTK_REAL is a double precision number. Setting otherwise will likely cause havoc with the conserv_to_prims solver.");

  if(cctk_nghostzones[0]<3 || cctk_nghostzones[1]<3 || cctk_nghostzones[2]<3) { CCTK_VError(VERR_DEF_PARAMS,"ERROR. Need at least 3 ghostzones for IllinoisGRMHD evolutions."); }

  CCTK_REAL dX[3] = { CCTK_DELTA_SPACE(0), CCTK_DELTA_SPACE(1), CCTK_DELTA_SPACE(2) };

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.1: Equation of state (EOS) parameters [Back to top]¶

$\label{eos_parameters}$

Next we set up the EOS struct, which is defined in the IllinoisGRMHD_headers.h header file. We set the following parameters:

$\rm neos$ : number of EOS (currently set to 1, which is a $\Gamma$ -law EOS)
$\rm K\_poly$ : this is the constant $\kappa$ from the polytropic EOS $P = \kappa \rho_{0}^{\Gamma}$
$\rm rho\_poly$ : $\rho_{0}$ , fluid rest-mass
$\rm P\_poly$ : $P$ , pressure
$\rm gamma\_th$ : $\Gamma_{\rm th}$ , the constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks
$\rm eps\_poly$ : $\epsilon$ , specific internal energy
$\rm k\_poly$ : $\kappa$ , polytropic EOS constant
$\rm gamma\_poly$ : $\Gamma$ , the $\Gamma$ -law EOS constant

In [5]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  /**********************************
   * Piecewise Polytropic EOS Patch *
   *   Setting up the EOS struct    *
   **********************************/
  /*
   * The short piece of code below takes care
   * of initializing the EOS parameters.
   * Please refer to the "inlined_functions.C"
   * source file for the documentation on the
   * function.
   */
  eos_struct eos;
  initialize_EOS_struct_from_input(eos);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.2: Set pointers to GRMHD gridfunctions [Back to top]¶

$\label{set_pointers_grmhd_gfs}$

In [6]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // in_prims,out_prims_r, and out_prims_l are arrays of pointers to the actual gridfunctions.
  gf_and_gz_struct in_prims[MAXNUMVARS],out_prims_r[MAXNUMVARS],out_prims_l[MAXNUMVARS];
  int which_prims_to_reconstruct[MAXNUMVARS],num_prims_to_reconstruct;

  /* SET POINTERS TO GRMHD GRIDFUNCTIONS */
  // The order here MATTERS, and must be consistent with the global variable declarations in
  //   evaluate_MHD_rhs_headers.h (look for RHOB=0, etc.)
  //   For example, in_prims[0] _must_ be rho_b.
  int ww=0;
  in_prims[ww].gf=rho_b; out_prims_r[ww].gf=rho_br; out_prims_l[ww].gf=rho_bl; ww++;
  in_prims[ww].gf=P;     out_prims_r[ww].gf=Pr;     out_prims_l[ww].gf=Pl;     ww++;
  in_prims[ww].gf=vx;    out_prims_r[ww].gf=vxr;    out_prims_l[ww].gf=vxl;    ww++;
  in_prims[ww].gf=vy;    out_prims_r[ww].gf=vyr;    out_prims_l[ww].gf=vyl;    ww++;
  in_prims[ww].gf=vz;    out_prims_r[ww].gf=vzr;    out_prims_l[ww].gf=vzl;    ww++;
  in_prims[ww].gf=Bx;    out_prims_r[ww].gf=Bxr;    out_prims_l[ww].gf=Bxl;    ww++;
  in_prims[ww].gf=By;    out_prims_r[ww].gf=Byr;    out_prims_l[ww].gf=Byl;    ww++;
  in_prims[ww].gf=Bz;    out_prims_r[ww].gf=Bzr;    out_prims_l[ww].gf=Bzl;    ww++;
  in_prims[ww].gf=Bx_stagger; out_prims_r[ww].gf=Bx_staggerr; out_prims_l[ww].gf=Bx_staggerl; ww++;
  in_prims[ww].gf=By_stagger; out_prims_r[ww].gf=By_staggerr; out_prims_l[ww].gf=By_staggerl; ww++;
  in_prims[ww].gf=Bz_stagger; out_prims_r[ww].gf=Bz_staggerr; out_prims_l[ww].gf=Bz_staggerl; ww++;
  in_prims[ww].gf=vxr;    out_prims_r[ww].gf=vxrr;    out_prims_l[ww].gf=vxrl;    ww++;
  in_prims[ww].gf=vyr;    out_prims_r[ww].gf=vyrr;    out_prims_l[ww].gf=vyrl;    ww++;
  in_prims[ww].gf=vzr;    out_prims_r[ww].gf=vzrr;    out_prims_l[ww].gf=vzrl;    ww++;
  in_prims[ww].gf=vxl;    out_prims_r[ww].gf=vxlr;    out_prims_l[ww].gf=vxll;    ww++;
  in_prims[ww].gf=vyl;    out_prims_r[ww].gf=vylr;    out_prims_l[ww].gf=vyll;    ww++;
  in_prims[ww].gf=vzl;    out_prims_r[ww].gf=vzlr;    out_prims_l[ww].gf=vzll;    ww++;

  // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:
  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }
  // Left/right variables are not yet defined, yet we set the # of gz's to zero by default:
  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.3: Convert ADM variables to BSSN variables [Back to top]¶

$\label{admbase_to_bssnbase}$

We summarize here the algorithm of the IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij() function, which is explained in detail in this tutorial module (Link not available yet - TODO):

First, $\gamma\equiv \det\left(\gamma_{ij}\right)$ , where $\gamma_{ij}$ is the physical spatial metric, is evaluated.
Then, $\phi$ , the conformal factor, is computed via the relation $\phi = \frac{1}{12}\log\gamma$ .
Next, we compute $e^{-4\phi}$ .
Then, the conformal metric, $\bar{\gamma}_{ij} = e^{-4\phi}\gamma_{ij}$ , is computed.
Next, the condition $\bar\gamma = 1$ is enforced, by first computing $\gamma$ then performing $\bar\gamma_{ij}\to\left(\frac{1}{\bar\gamma}\right)^{1/3}\bar\gamma_{ij}$ .
Then, $\gamma_{ij}$ is computed from $\bar\gamma_{ij}$ after the condition $\bar\gamma = 1$ is enforced via the inverse relation $\gamma_{ij} = e^{4\phi}\bar\gamma_{ij}$ .
Finally, we compute the inverse conformal metric $\bar\gamma^{ij}$ .

A note on notation: in the C code, we have the following identifications between the quantities described above and the C variables:

Input and temporary variables:

$\gamma_{ij} := {\rm gij\_physL}$
$\det(\gamma_{ij}) := {\rm gijdet}$
$\phi := {\rm phiL}$
$\psi \equiv e^{\phi} := {\rm psiL}$
$\bar\gamma_{ij} := {\rm gtijL}$ (the "t" stands for the notation where the conformal metric is written as $\tilde\gamma_{ij}$ instead of $\bar\gamma_{ij}$ . In our discussion we use the latter to keep our notation consistent with other NRPy notebooks).
$\det(\bar\gamma_{ij}) := {\rm gtijdet}$
$\left(\frac{1}{\bar\gamma}\right)^{1/3} := {\rm gtijdet\_Fm1o3}$

Output/gridfunction variables:

$\gamma_{ij} := {\rm gij}$ (Physical metric)
$\bar\gamma_{ij} := {\rm gtij}$ (Conformal metric)
$\bar\gamma^{ij} := {\rm gtupij}$ (Inverse conformal metric)
$\phi := {\rm phi}$
$\psi := {\rm psi}$

In [7]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // Convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.
  IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh,  gxx,gxy,gxz,gyy,gyz,gzz,alp,
                                                                gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
                                                                gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
                                                                phi_bssn,psi_bssn,lapm1);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.4: Setting up pointers to the metric and stress-energy tensor gridfunctions [Back to top]¶

$\label{pointers_metric_tmunu_gfs}$

In [8]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  /* SET POINTERS TO METRIC GRIDFUNCTIONS */
  CCTK_REAL *metric[NUMVARS_FOR_METRIC_FACEVALS]; // "metric" here is array of pointers to the actual gridfunctions.
  ww=0;
  metric[ww]=phi_bssn;ww++;
  metric[ww]=psi_bssn;ww++;
  metric[ww]=gtxx;    ww++;
  metric[ww]=gtxy;    ww++;
  metric[ww]=gtxz;    ww++;
  metric[ww]=gtyy;    ww++;
  metric[ww]=gtyz;    ww++;
  metric[ww]=gtzz;    ww++;
  metric[ww]=lapm1;   ww++;
  metric[ww]=betax;   ww++;
  metric[ww]=betay;   ww++;
  metric[ww]=betaz;   ww++;
  metric[ww]=gtupxx;  ww++;
  metric[ww]=gtupyy;  ww++;
  metric[ww]=gtupzz;  ww++;

  /* SET POINTERS TO STRESS-ENERGY TENSOR GRIDFUNCTIONS */
  CCTK_REAL *TUPmunu[10];// "TUPmunu" here is array of pointers to the actual gridfunctions.
  ww=0;
  TUPmunu[ww]=TUPtt; ww++;
  TUPmunu[ww]=TUPtx; ww++;
  TUPmunu[ww]=TUPty; ww++;
  TUPmunu[ww]=TUPtz; ww++;
  TUPmunu[ww]=TUPxx; ww++;
  TUPmunu[ww]=TUPxy; ww++;
  TUPmunu[ww]=TUPxz; ww++;
  TUPmunu[ww]=TUPyy; ww++;
  TUPmunu[ww]=TUPyz; ww++;
  TUPmunu[ww]=TUPzz; ww++;

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.5: Initialization of the RHS variables [Back to top]¶

$\label{initialize_rhss}$

In [9]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // 1) First initialize {rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs} to zero
#pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);
        Ax_rhs[index]=0.0;
        Ay_rhs[index]=0.0;
        Az_rhs[index]=0.0;
        psi6phi_rhs[index]=0.0;

        tau_rhs[index]=0.0;
        rho_star_rhs[index]=0.0;
        st_x_rhs[index]=0.0;
        st_y_rhs[index]=0.0;
        st_z_rhs[index]=0.0;

        //if(i==17 && j==19 && k==26) CCTK_VInfo(CCTK_THORNSTRING,"CONSSS: %.15e %.15e %.15e %.15e %.15e | %.15e",rho_star[index],mhd_st_x[index],mhd_st_y[index],mhd_st_z[index],tau[index],P[index]);
      }

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.6: Compute extrinsic curvature terms from the RHS of $\partial_{t}\tilde\tau$ and $T^{\mu\nu}$ [Back to top]¶

$\label{tau_rhs_ext_curv_and_tupmunu}$

In [10]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  // Here, we:
  // 1) Compute tau_rhs extrinsic curvature terms, and
  // 2) Compute TUPmunu.
  // This function is housed in the file: "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C"
  compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,
                                                        eos, Gamma_th,
                                                        gtupxy,gtupxz,gtupyz,
                                                        kxx,kxy,kxz,kyy,kyz,kzz,
                                                        tau_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.7: Computing ${\rm ftilde}$ [Back to top]¶

$\label{computing_ftilde}$

This is part of the flattening scheme of the PPM algorithm. The main reference to look at is Colella & Woodward (1983). The equations implemented can be found in Appendix A (particularly eqs. (A.1) and (A.2)), while the flattening method is introduced and discussed in section 4. More will follow when we talk about the reconstruct_set_of_prims_PPM.C file of IllinoisGRMHD.

In [11]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  int flux_dirn;
  flux_dirn=1;
  // First compute ftilde, which is used for flattening left and right face values
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8: The RHSs of $\rho_{\star}$ , $\tilde\tau$ , $\tilde{S}_{i}$ , and $A_{i}$ [Back to top]¶

$\label{rhs_mhd_and_a_i}$

This part of the code evaluates the RHSs of $\rho_{\star}$ , $\tilde\tau$ , and $\tilde{S}_{i}$ , i.e.

$\partial_{t} \begin{bmatrix} \rho_{\star}\\ \tilde\tau\\ \tilde{S}_{i} \end{bmatrix} = -\partial_{j} \underbrace{\begin{bmatrix} \rho_{\star}v^{j}\\ \alpha^{2}\sqrt{\gamma}T^{0j} - \rho_{\star}v^{j}\\ \alpha\sqrt{\gamma}T^{j}_{\ i} \end{bmatrix}}_{\rm Flux\ terms} + \underbrace{\begin{bmatrix} 0\\ s\\ \frac{1}{2}\alpha\sqrt{\gamma}T^{\alpha\beta}\partial_{i}g_{\alpha\beta} \end{bmatrix}}_{\rm Source\ terms}\ .$

At the same time, we are also interested in evaluating the RHS of the evolution equation for $A_{i}$ , namely

$\partial_{t}A_{i} = \epsilon_{ijk}v^{j}\tilde{B}^{k} - \partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right) = \psi^{6}\epsilon_{ijk}v^{j}B^{k} - \underbrace{\partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right)}_{\rm Gauge\ terms}$

The following summary greatly oversimplifies what the code below does, but it is enough for the user to understand the purpose of the algorithm:

Compute $\partial_{x}\boldsymbol{F}$ , then $\partial_{y}\boldsymbol{F}$ , and finally $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$
Compute $\partial_{y}\boldsymbol{F}$ , then $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$
Compute $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$
Add gauge terms to $\partial_{t}A_{i}$

Now, in between every step of the summary above, care must be taken to evaluate the gridfunctions at the appropriate gridpoints (see table below for the location of each variable in the computational grid).

Variable(s)	Gridpoint location in the computational grid
Metric terms, $\vec{P}$ , $\rho_*$ , $\tilde{S}_i$ , $\tilde{\tau}$	$(i,j,k)$
$B^x$ , $\tilde{B}^x$	$(i+\frac{1}{2},j,k)$
$B^y$ , $\tilde{B}^y$	$(i,j+\frac{1}{2},k)$
$B^z$ , $\tilde{B}^z$	$(i,j,k+\frac{1}{2})$
$A_x$	$(i,j+\frac{1}{2},k+\frac{1}{2})$
$A_y$	$(i+\frac{1}{2},j,k+\frac{1}{2})$
$A_z$	$(i+\frac{1}{2},j+\frac{1}{2},k)$
$\sqrt{\gamma}\Phi$	$(i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2})$

$\label{table_staggerings}$

For example, we know that

$\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms} \equiv \psi^{6}\left(v^{x}B^{y} - v^{y}B^{x}\right)\ .$

But once we evaluate $v^{x}$ and $v^{y}$ , we know them at the point $(i,j,k)$ . Similarly, the gridfunction $B^{x}$ is known at $(i+\frac{1}{2},j,k)$ , while $B^{y}$ is known at $(i,j+\frac{1}{2},k)$ . This means that we are not able to immediately evaluate the equation above, since determining $A_{z}$ at $(i+\frac{1}{2},j+\frac{1}{2},k)$ requires knowning $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ at $(i+\frac{1}{2},j+\frac{1}{2},k)$ as well. To this end, we reconstruct the variables $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ using the PPM method at the desired staggered point. An analogous procedure is required in order to determine the RHS of $\partial_{t}A_{x}$ and $\partial_{t}A_{y}$ .

Step 3.8.1: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $x$ -direction [Back to top]¶

$\label{reconstructing_vx_vy_by_along_x}$

We want to evaluate $\partial_{x}\boldsymbol{F}$ . It is important that we keep in the back of our minds our intention of evaluating the RHS of $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ as well, since then we can reconstruct $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ cleverly, as we need them at the same gridpoint as $A_{z}$ (see the table in the end of step 3.8).

We start by reconstructing $\left\{\rho_{0},P,v^{i},B^{i}, B^{y}_{\rm stagger}\right\}$ in the $x$ -direction, keeping in mind that after the reconstruction we will know:

The flux variables at $\left(i-\frac{1}{2},j,k\right)$ , so that we can evaluate $\partial_{x}\boldsymbol{F}_{i,j,k}=dx^{-1}\left(\boldsymbol{F}_{i+1/2,j,k}-\boldsymbol{F}_{i-1/2,j,k}\right)$
The velocities $\left\{v^{x},v^{y}\right\}$ at $\left(i-\frac{1}{2},j,k\right)$
The staggered value $B^{y}_{\rm stagger}$ at $\left(i-\frac{1}{2},j+\frac{1}{2},k\right)$

In [12]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C



  /* There are two stories going on here:
   * 1) Computation of \partial_x on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}},
   *    via PPM reconstruction onto (i-1/2,j,k), so that
   *    \partial_x F = [ F(i+1/2,j,k) - F(i-1/2,j,k) ] / dx
   * 2) Computation of \partial_t A_i, where A_i are *staggered* gridfunctions,
   *    where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.
   *    Ai_rhs = \partial_t A_i = \epsilon_{ijk} \psi^{6} v^j B^k,
   *    where \epsilon_{ijk} is the flat-space antisymmetric operator.
   * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},
   *     so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these
   *     staggered points. For example:
   * 2Aa) vx and vy are at (i,j,k), and we reconstruct them to (i-1/2,j,k) below. After
   *      this, we'll reconstruct again in the y-dir'n to get {vx,vy} at (i-1/2,j-1/2,k)
   * 2Ab) By_stagger is at (i,j+1/2,k), and we reconstruct below to (i-1/2,j+1/2,k). */
  ww=0;
  which_prims_to_reconstruct[ww]=RHOB;      ww++;
  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;
  which_prims_to_reconstruct[ww]=VX;        ww++;
  which_prims_to_reconstruct[ww]=VY;        ww++;
  which_prims_to_reconstruct[ww]=VZ;        ww++;
  //which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BY_STAGGER;ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.2: Evaluating $\partial_{x}\boldsymbol{F}$ [Back to top]¶

$\label{fluxes_x_dirn}$

Next we set the face values of $B^{x}$ (which are needed for the computation of MHD flux terms) by making them consistent with $B^{x}_{\rm stagger}$ .

After that, we evaluate $\partial_{x}\boldsymbol{F}$ and add it to the RHS of $\partial_{t}\left(\rho_{\star},\tilde\tau,\tilde{S}_{i}\right)$ . It is important to notice that, as we mentioned, $A_{z}$ is defined at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$ , but other functions, like $v^{x}$ and $v^{y}$ , are now known only at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$ . The function add_fluxes_and_source_terms_to_hydro_rhss() below takes care of this, and we will study the process in more detail when we look at the add_fluxes_and_source_terms_to_hydro_rhss.C file of IllinoisGRMHD.

In [13]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).
  //   Instead of reconstructing, we simply set B^x face values to be consistent with BX_STAGGER.
#pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
        int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexim1=CCTK_GFINDEX3D(cctkGH,i-1+(i==0),j,k); /* indexim1=0 when i=0 */
        out_prims_r[BX_CENTER].gf[index]=out_prims_l[BX_CENTER].gf[index]=in_prims[BX_STAGGER].gf[indexim1]; }
  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:
  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,
                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,
                                            cmax_x,cmin_x,
                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,
                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.3: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $y$ -direction [Back to top]¶

$\label{reconstructing_vx_vy_by_along_y}$

We want to evaluate $\partial_{y}\boldsymbol{F}$ . At this point we must remember that $v^{x}$ and $v^{y}$ have already been reconstruct along the $x$ -direction and are now known at $\left(i-\frac{1}{2},j,k\right)$ . Our goal is to reconstruct these quantities at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$ .

We then reconstruct $\left\{\rho_{0},P,v^{i},B^{i}, B^{i}_{\rm stagger}\right\}$ in the $y$ -direction, keeping in mind that after the reconstruction we will know:

The flux variables at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$ , so that we can evaluate $\partial_{y}\boldsymbol{F}_{i,j,k}=dy^{-1}\left(\boldsymbol{F}_{i,j+1/2,k}-\boldsymbol{F}_{i,j-1/2,k}\right)$
The velocities $\left\{v^{x},v^{y}\right\}$ at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$
The staggered value $B^{x}_{\rm stagger}$ at $\left(i+\frac{1}{2},j-\frac{1}{2},k\right)$
The staggered value $B^{y}_{\rm stagger}$ at $\left(i-\frac{1}{2},j+\frac{1}{2},k\right)$
The staggered value $B^{z}_{\rm stagger}$ at $\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)$

In [14]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // Note that we have already reconstructed vx and vy along the x-direction,
  //   at (i-1/2,j,k). That result is stored in v{x,y}{r,l}.  Bx_stagger data
  //   are defined at (i+1/2,j,k).
  // Next goal: reconstruct Bx, vx and vy at (i+1/2,j+1/2,k).
  flux_dirn=2;
  // First compute ftilde, which is used for flattening left and right face values
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);

  // in_prims[{VXR,VXL,VYR,VYL}].gz_{lo,hi} ghostzones are set to all zeros, which
  //    is incorrect. We fix this below.
  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]
  in_prims[VXR]=out_prims_r[VX];
  in_prims[VXL]=out_prims_l[VX];
  in_prims[VYR]=out_prims_r[VY];
  in_prims[VYL]=out_prims_l[VY];

  /* There are two stories going on here:
   * 1) Computation of \partial_y on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}},
   *    via PPM reconstruction onto (i,j-1/2,k), so that
   *    \partial_y F = [ F(i,j+1/2,k) - F(i,j-1/2,k) ] / dy
   * 2) Computation of \partial_t A_i, where A_i are *staggered* gridfunctions,
   *    where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.
   *    Ai_rhs = \partial_t A_i = \epsilon_{ijk} \psi^{6} v^j B^k,
   *    where \epsilon_{ijk} is the flat-space antisymmetric operator.
   * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},
   *     so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these
   *     staggered points. For example:
   * 2Aa) VXR = [right-face of vx reconstructed along x-direction above] is at (i-1/2,j,k),
   *      and we reconstruct it to (i-1/2,j-1/2,k) below. Similarly for {VXL,VYR,VYL}
   * 2Ab) Bx_stagger is at (i+1/2,j,k), and we reconstruct to (i+1/2,j-1/2,k) below
   * 2Ac) By_stagger is at (i-1/2,j+1/2,k) already for Az_rhs, from the previous step.
   * 2B) Ax_rhs is defined at (i,j+1/2,k+1/2), and it depends on {By,Bz,vy,vz}.
   *     Again the trick is to reconstruct these onto these staggered points.
   * 2Ba) Bz_stagger is at (i,j,k+1/2), and we reconstruct to (i,j-1/2,k+1/2) below */
  ww=0;
  // NOTE! The order of variable reconstruction is important here,
  //   as we don't want to overwrite {vxr,vxl,vyr,vyl}!
  which_prims_to_reconstruct[ww]=VXR;       ww++;
  which_prims_to_reconstruct[ww]=VYR;       ww++;
  which_prims_to_reconstruct[ww]=VXL;       ww++;
  which_prims_to_reconstruct[ww]=VYL;       ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);
  ww=0;
  // Reconstruct other primitives last!
  which_prims_to_reconstruct[ww]=RHOB;      ww++;
  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;
  which_prims_to_reconstruct[ww]=VX;        ww++;
  which_prims_to_reconstruct[ww]=VY;        ww++;
  which_prims_to_reconstruct[ww]=VZ;        ww++;
  which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
  //which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BX_STAGGER;ww++;
  which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.4: Evaluating $\partial_{y}\boldsymbol{F}$ [Back to top]¶

$\label{fluxes_y_dirn}$

In [15]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).
  //   Instead of reconstructing, we simply set B^y face values to be consistent with BY_STAGGER.
#pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
        int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexjm1=CCTK_GFINDEX3D(cctkGH,i,j-1+(j==0),k); /* indexjm1=0 when j=0 */
        out_prims_r[BY_CENTER].gf[index]=out_prims_l[BY_CENTER].gf[index]=in_prims[BY_STAGGER].gf[indexjm1]; }
  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:
  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,
                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,
                                            cmax_y,cmin_y,
                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,
                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.5: Evaluating $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

$\label{rhs_az_no_gauge_terms}$

As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:

Staggered and unstaggered variables	Gridpoint location at which the variable is known
$\left(v^{x}\right)_{r,l},\left(v^{y}\right)_{r,l}$	$(i-\frac{1}{2},j-\frac{1}{2},k)$
$\left(B^{z}_{\rm stagger}\right)_{r,l}$	$(i+\frac{1}{2},j-\frac{1}{2},k)$
$\left(B^{y}_{\rm stagger}\right)_{r,l}$	$(i-\frac{1}{2},j+\frac{1}{2},k)$
$\phi$	$(i,j,k)$

We start by interpolating $\phi$ to $\left(i+\frac{1}{2},j,k\right)$ , followed by a second interpolation so that $\phi$ is known at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$ .

In [16]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  /*****************************************
   * COMPUTING RHS OF A_z, BOOKKEEPING NOTE:
   * We want to compute
   * \partial_t A_z - [gauge terms] = \psi^{6} (v^x B^y - v^y B^x).
   * A_z is defined at (i+1/2,j+1/2,k).
   * ==========================
   * Where defined  | Variables
   * (i-1/2,j-1/2,k)| {vxrr,vxrl,vxlr,vxll,vyrr,vyrl,vylr,vyll}
   * (i+1/2,j-1/2,k)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i-1/2,j+1/2,k)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i,j,k)        | {phi}
   * ==========================
   ******************************************/
  // Interpolates to i+1/2
#define IPH(METRICm1,METRICp0,METRICp1,METRICp2) (-0.0625*((METRICm1) + (METRICp2)) + 0.5625*((METRICp0) + (METRICp1)))
  // Next compute phi at (i+1/2,j+1/2,k):
#pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=1;j<cctk_lsh[1]-2;j++) for(int i=1;i<cctk_lsh[0]-2;i++) {
        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j-1,k)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j  ,k)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j+1,k)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j+2,k)]));
      }

Appending to ../src/driver_evaluate_MHD_rhs.C

Then we update the RHS of $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ . Keep in mind that the function A_i_rhs_no_gauge_terms() takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$ . We will look at it in more detail when we see the A_i_rhs_no_gauge_terms.C file from IllinoisGRMHD.

In [17]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  int A_directionz=3;
  A_i_rhs_no_gauge_terms(A_directionz,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_x,cmin_x,cmax_y,cmin_y, Az_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.6: Multiple reconstructions [Back to top]¶

$\label{multiple_reconstructions}$

$\left\{\rho_{\star},P,v^{i}, B^{i}\right\}$ from $\left(i,j,k\right)$ to $\left(i,j,k-\frac{1}{2}\right)$
[∂tAx]no gauge terms is defined at (i,j+12,k+12)
1. $\left(v^{y}\right)_{r,l}$ and $\left(v^{z}\right)_{r,l}$ are at $\left(i,j-\frac{1}{2},k\right)$ , so we reconstruct them to $\left(i,j-\frac{1}{2},k-\frac{1}{2}\right)$
2. $\left(B^{z}_{\rm stagger}\right)_{r,l}$ is already known at $\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)$
3. $B^{y}_{\rm stagger}$ is at $\left(i,j+\frac{1}{2},k\right)$ , so we reconstruct it to $\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)$
[∂tAy]no gauge terms is defined at (i+12,j,k+12)
1. $v^{x}$ and $v^{z}$ are reconstructed to $\left(i,j,k-\frac{1}{2}\right)$ . We'll reconstruct them to $\left(i,j-\frac{1}{2},k-\frac{1}{2}\right)$ later
2. $B^{z}_{\rm stagger}$ is already known at $\left(i,j,k+\frac{1}{2}\right)$ . We'll reconstruct them to $\left(i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}\right)$ later
3. $B^{x}_{\rm stagger}$ is at $\left(i,j+\frac{1}{2},k\right)$ , so we reconstruct it to $\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)$

In [18]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix
  //    this below.
  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]
  in_prims[VYR]=out_prims_r[VY];
  in_prims[VYL]=out_prims_l[VY];
  in_prims[VZR]=out_prims_r[VZ];
  in_prims[VZL]=out_prims_l[VZ];

  flux_dirn=3;
  // First compute ftilde, which is used for flattening left and right face values
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);

  /* There are two stories going on here:
   * 1) Single reconstruction to (i,j,k-1/2) for {rho,P,vx,vy,vz,Bx,By,Bz} to compute
   *    z-dir'n advection terms in \partial_t {rho_star,tau,mhd_st_{x,y,z}} at (i,j,k)
   * 2) Multiple reconstructions for *staggered* gridfunctions A_i:
   *    Ai_rhs = \partial_t A_i = \epsilon_{ijk} \psi^{6} v^j B^k,
   *    where \epsilon_{ijk} is the flat-space antisymmetric operator.
   * 2A) Ax_rhs is defined at (i,j+1/2,k+1/2), depends on v{y,z} and B{y,z}
   * 2Aa) v{y,z}{r,l} are at (i,j-1/2,k), so we reconstruct here to (i,j-1/2,k-1/2)
   * 2Ab) Bz_stagger{r,l} are at (i,j-1/2,k+1/2) already.
   * 2Ac) By_stagger is at (i,j+1/2,k), and below we reconstruct its value at (i,j+1/2,k-1/2)
   * 2B) Ay_rhs is defined at (i+1/2,j,k+1/2), depends on v{z,x} and B{z,x}.
   * 2Ba) v{x,z} are reconstructed to (i,j,k-1/2). Later we'll reconstruct again to (i-1/2,j,k-1/2).
   * 2Bb) Bz_stagger is at (i,j,k+1/2). Later we will reconstruct to (i-1/2,j,k+1/2).
   * 2Bc) Bx_stagger is at (i+1/2,j,k), and below we reconstruct its value at (i+1/2,j,k-1/2)
   */
  ww=0;
  // NOTE! The order of variable reconstruction is important here,
  //   as we don't want to overwrite {vxr,vxl,vyr,vyl}!
  which_prims_to_reconstruct[ww]=VYR;       ww++;
  which_prims_to_reconstruct[ww]=VZR;       ww++;
  which_prims_to_reconstruct[ww]=VYL;       ww++;
  which_prims_to_reconstruct[ww]=VZL;       ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);
  // Reconstruct other primitives last!
  ww=0;
  which_prims_to_reconstruct[ww]=RHOB;      ww++;
  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;
  which_prims_to_reconstruct[ww]=VX;        ww++;
  which_prims_to_reconstruct[ww]=VY;        ww++;
  which_prims_to_reconstruct[ww]=VZ;        ww++;
  which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
  //which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
  which_prims_to_reconstruct[ww]=BX_STAGGER; ww++;
  which_prims_to_reconstruct[ww]=BY_STAGGER; ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);
  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).
  //   Instead of reconstructing, we simply set B^z face values to be consistent with BZ_STAGGER.
#pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
        int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexkm1=CCTK_GFINDEX3D(cctkGH,i,j,k-1+(k==0)); /* indexkm1=0 when k=0 */
        out_prims_r[BZ_CENTER].gf[index]=out_prims_l[BZ_CENTER].gf[index]=in_prims[BZ_STAGGER].gf[indexkm1]; }

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.7: Evaluating $\partial_{z}\boldsymbol{F}$ [Back to top]¶

$\label{fluxes_z_dirn}$

In [19]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:
  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,
                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,
                                            cmax_z,cmin_z,
                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,
                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);

  // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not set correcty.
  //    We fix this below.
  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]
  in_prims[VXR]=out_prims_r[VX];
  in_prims[VZR]=out_prims_r[VZ];
  in_prims[VXL]=out_prims_l[VX];
  in_prims[VZL]=out_prims_l[VZ];
  // FIXME: lines above seem to be inconsistent with lines below.... Possible bug, not major enough to affect evolutions though.
  in_prims[VZR].gz_lo[1]=in_prims[VZR].gz_hi[1]=0;
  in_prims[VXR].gz_lo[1]=in_prims[VXR].gz_hi[1]=0;
  in_prims[VZL].gz_lo[1]=in_prims[VZL].gz_hi[1]=0;
  in_prims[VXL].gz_lo[1]=in_prims[VXL].gz_hi[1]=0;

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.8: Evaluating $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

$\label{rhs_ax_no_gauge_terms}$

As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:

Staggered and unstaggered variables	Gridpoint location at which the variable is known
$\left(v^{y}\right)_{r,l},\left(v^{z}\right)_{r,l}$	$(i,j-\frac{1}{2},k-\frac{1}{2})$
$\left(B^{y}_{\rm stagger}\right)_{r,l}$	$(i,j+\frac{1}{2},k-\frac{1}{2})$
$\left(B^{z}_{\rm stagger}\right)_{r,l}$	$(i,j-\frac{1}{2},k+\frac{1}{2})$
$\phi$	$(i,j,k)$

We start by interpolating $\phi$ to $\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)$ .

In [20]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C

  /*****************************************
   * COMPUTING RHS OF A_x, BOOKKEEPING NOTE:
   * We want to compute
   * \partial_t A_x - [gauge terms] = \psi^{6} (v^y B^z - v^z B^y).
   * A_x is defined at (i,j+1/2,k+1/2).
   * ==========================
   * Where defined  | Variables
   * (i,j-1/2,k-1/2)| {vyrr,vyrl,vylr,vyll,vzrr,vzrl,vzlr,vzll}
   * (i,j+1/2,k-1/2)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i,j-1/2,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i,j,k)        | {phi}
   * ==========================
   ******************************************/
  // Next compute phi at (i,j+1/2,k+1/2):
#pragma omp parallel for
  for(int k=1;k<cctk_lsh[2]-2;k++) for(int j=1;j<cctk_lsh[1]-2;j++) for(int i=0;i<cctk_lsh[0];i++) {
        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k-1)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k  )]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k+1)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k+2)]));
      }

Appending to ../src/driver_evaluate_MHD_rhs.C

Then we update the RHS of $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$ . Keep in mind that the function A_i_rhs_no_gauge_terms() takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)$ . We will look at it in more detail when we see the A_i_rhs_no_gauge_terms.C file from IllinoisGRMHD.

In [21]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  int A_directionx=1;
  A_i_rhs_no_gauge_terms(A_directionx,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_y,cmin_y,cmax_z,cmin_z, Ax_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.9: Evaluating $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

$\label{rhs_ay_no_gauge_terms}$

As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:

Staggered and unstaggered variables	Gridpoint location at which the variable is known
$\left(v^{x}\right)_{r,l},\left(v^{z}\right)_{r,l}$	$(i-\frac{1}{2},j,k-\frac{1}{2})$
$\left(B^{x}_{\rm stagger}\right)_{r,l}$	$(i+\frac{1}{2},j,k-\frac{1}{2})$
$\left(B^{z}_{\rm stagger}\right)_{r,l}$	$(i-\frac{1}{2},j,k+\frac{1}{2})$
$\phi$	$(i,j,k)$

We start by interpolating $\phi$ to $\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)$ .

In [22]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // We reprise flux_dirn=1 to finish up computations of Ai_rhs's!
  flux_dirn=1;
  // First compute ftilde, which is used for flattening left and right face values
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);

  ww=0;
  // NOTE! The order of variable reconstruction is important here,
  //   as we don't want to overwrite {vxr,vxl,vyr,vyl}!
  which_prims_to_reconstruct[ww]=VXR;       ww++;
  which_prims_to_reconstruct[ww]=VZR;       ww++;
  which_prims_to_reconstruct[ww]=VXL;       ww++;
  which_prims_to_reconstruct[ww]=VZL;       ww++;
  which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"
  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);

  /*****************************************
   * COMPUTING RHS OF A_y, BOOKKEEPING NOTE:
   * We want to compute
   * \partial_t A_y - [gauge terms] = \psi^{6} (v^z B^x - v^x B^z).
   * A_y is defined at (i+1/2,j,k+1/2).
   * ==========================
   * Where defined  | Variables
   * (i-1/2,j,k-1/2)| {vxrr,vxrl,vxlr,vxll,vzrr,vzrl,vzlr,vzll}
   * (i+1/2,j,k-1/2)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i-1/2,j,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)
   * (i,j,k)        | {phi}
   * ==========================
   ******************************************/
  // Next compute phi at (i+1/2,j,k+1/2):
#pragma omp parallel for
  for(int k=1;k<cctk_lsh[2]-2;k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=1;i<cctk_lsh[0]-2;i++) {
        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k-1)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k  )]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k+1)]),
              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k+2)]));
      }

Appending to ../src/driver_evaluate_MHD_rhs.C

Then we update the RHS of $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$ . Keep in mind that the function A_i_rhs_no_gauge_terms() takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)$ . We will look at it in more detail when we see the A_i_rhs_no_gauge_terms.C file from IllinoisGRMHD.

In [23]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  int A_directiony=2;
  A_i_rhs_no_gauge_terms(A_directiony,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_z,cmin_z,cmax_x,cmin_x, Ay_rhs);

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 3.8.10: Evaluating $\partial_{t}\left[\psi^{6}\Phi\right]$ and $\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms}$ [Back to top]¶

$\label{rhs_psi6phi_and_ai_gauge_terms}$

Finally, we compute

$\partial_{t}\left[\sqrt{\gamma}\Phi\right] = \partial_{t}\left[\psi^{6}\Phi\right] = -\partial_{j}\left(\alpha\sqrt{\gamma}A^{j} - \beta^{j}\left[\sqrt{\gamma}\Phi\right]\right) -\xi\alpha\left[\sqrt{\gamma}\Phi\right]\ ,$

and

$\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms} = -\partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right)\ .$

Notice that we will need $A^{i}$ to compute $\partial_{t}\left[\psi^{6}\Phi\right]$ , but we only have $A_{i}$ , so we need to determine $\bar\gamma^{ij}$ ( ${\rm gtupij}$ ).

In [24]:

%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C


  // Next compute psi6phi_rhs, and add gauge terms to A_i_rhs terms!
  //   Note that in the following function, we don't bother with reconstruction, instead interpolating.
  // We need A^i, but only have A_i. So we add gtupij to the list of input variables.
  CCTK_REAL *interp_vars[MAXNUMINTERP];
  ww=0;
  interp_vars[ww]=betax;   ww++;
  interp_vars[ww]=betay;   ww++;
  interp_vars[ww]=betaz;   ww++;
  interp_vars[ww]=gtupxx;  ww++;
  interp_vars[ww]=gtupxy;  ww++;
  interp_vars[ww]=gtupxz;  ww++;
  interp_vars[ww]=gtupyy;  ww++;
  interp_vars[ww]=gtupyz;  ww++;
  interp_vars[ww]=gtupzz;  ww++;
  interp_vars[ww]=psi_bssn;ww++;
  interp_vars[ww]=lapm1;   ww++;
  interp_vars[ww]=Ax;      ww++;
  interp_vars[ww]=Ay;      ww++;
  interp_vars[ww]=Az;      ww++;
  int max_num_interp_variables=ww;
  if(max_num_interp_variables>MAXNUMINTERP) {CCTK_VError(VERR_DEF_PARAMS,"Error: Didn't allocate enough space for interp_vars[]."); }
  // We are FINISHED with v{x,y,z}{r,l} and P{r,l} so we use these 8 gridfunctions' worth of space as temp storage.
  Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi,
                                                 vxr,vyr,vzr,vxl,vyl,vzl,Pr,Pl,
                                                 psi6phi_rhs,Ax_rhs,Ay_rhs,Az_rhs);


  return;
  /*
  // FUN DEBUGGING TOOL (trust me!):
  #pragma omp parallel for
  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
  int index=CCTK_GFINDEX3D(cctkGH,i,j,k);
  //st_x_rhs[index]=0.0;
  //st_y_rhs[index]=0.0;
  //st_z_rhs[index]=0.0;
  //rho_star_rhs[index]=0.0;
  //tau_rhs[index]=0.0;

  psi6phi_rhs[index] = 0.0;
  Ax_rhs[index] = 0.0;
  Ay_rhs[index] = 0.0;
  Az_rhs[index] = 0.0;
  }
  */
}

// We add #include's here instead of compiling these separately to help ensure that functions are properly inlined.
//    These files only include about 800 lines of code in total (~1200 lines in total), but it's arguably more
//    convenient to edit a 600 line file than an 1800 line file, so I'd prefer to leave this unconventional structure
//    alone.
#include "reconstruct_set_of_prims_PPM.C"
#include "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C"
#include "add_fluxes_and_source_terms_to_hydro_rhss.C"
#include "mhdflux.C"
#include "A_i_rhs_no_gauge_terms.C"
#include "Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C"

Appending to ../src/driver_evaluate_MHD_rhs.C

Step 4: The `driver_evaluate_MHD_rhs.h` header file [Back to top]¶

$\label{driver_evaluate_MHD_rhs__h}$

Now we generate the header file for the driver_evaluate_MHD_rhs.C file.

In [25]:

%%writefile $outfile_path__driver_evaluate_MHD_rhs__h
#ifndef DRIVER_EVALUATE_MHD_RHS_H_
#define DRIVER_EVALUATE_MHD_RHS_H_

/* PRIVATE FUNCTIONS, Called within driver_evaluate_MHD_rhs.C ONLY */
static void ftilde_gf_compute(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,gf_and_gz_struct *input,CCTK_REAL *ftilde_gf);
static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,
                                         eos_struct &eosi,gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,
                                         CCTK_REAL *ftilde_gf,CCTK_REAL *temporary);

static void compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu
(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **metric,gf_and_gz_struct *prims,
 CCTK_REAL **TUPmunu,eos_struct &eos, CCTK_REAL Gamma_th,
 CCTK_REAL *gupxy,CCTK_REAL *gupxz,CCTK_REAL *gupyz,
 CCTK_REAL *kxx,CCTK_REAL *kxy,CCTK_REAL *kxz,CCTK_REAL *kyy,CCTK_REAL *kyz,CCTK_REAL *kzz,
 CCTK_REAL *tau_rhs);
static void A_i_rhs_no_gauge_terms(const int A_dirn, const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,
                                   CCTK_REAL *phi_interped,CCTK_REAL *cmax_1,CCTK_REAL *cmin_1,CCTK_REAL *cmax_2,CCTK_REAL *cmin_2, CCTK_REAL *A3_rhs);

static void Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **interp_vars,CCTK_REAL *psi6phi,
                                                           CCTK_REAL *shiftx_iphjphkph,CCTK_REAL *shifty_iphjphkph,CCTK_REAL *shiftz_iphjphkph,
                                                           CCTK_REAL *alpha_iphjphkph,CCTK_REAL *alpha_Phi_minus_betaj_A_j_iphjphkph,CCTK_REAL *alpha_sqrtg_Ax_interp,
                                                           CCTK_REAL *alpha_sqrtg_Ay_interp,CCTK_REAL *alpha_sqrtg_Az_interp,
                                                           CCTK_REAL *psi6phi_rhs,CCTK_REAL *Ax_rhs,CCTK_REAL *Ay_rhs,CCTK_REAL *Az_rhs);

static void add_fluxes_and_source_terms_to_hydro_rhss(const int flux_dirn,const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,
                                                      CCTK_REAL **metric,gf_and_gz_struct *in_prims,CCTK_REAL **TUPmunu,
                                                      int numvars_reconstructed,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,eos_struct &eos,
                                                      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 *rho_star_rhs,CCTK_REAL *tau_rhs,CCTK_REAL *st_x_rhs,CCTK_REAL *st_y_rhs,CCTK_REAL *st_z_rhs);

#include "harm_primitives_headers.h"

#endif /* DRIVER_EVALUATE_MHD_RHS_H_ */

Writing ../src/driver_evaluate_MHD_rhs.h

Step 5: Code validation [Back to top]¶

$\label{code_validation}$

Step 5.a: `driver_evaluate_MHD_rhs.C` [Back to top]¶

$\label{code_validation_driver_evaluate_MHD_rhs__c}$

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

In [26]:

# 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/driver_evaluate_MHD_rhs.C"
original_IGM_file_name = "driver_evaluate_MHD_rhs-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__driver_evaluate_MHD_rhs__C  = !diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__C

if Validation__driver_evaluate_MHD_rhs__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 driver_evaluate_MHD_rhs.C: PASSED!")
else:
    # If the validation fails, we keep the original IGM source code file
    print("Validation test for driver_evaluate_MHD_rhs.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__driver_evaluate_MHD_rhs__C:
        print(diff_line)

Validation test for driver_evaluate_MHD_rhs.C: FAILED!
Diff:
3,4c3,4
<  * (vector potential prescription), using the 
<  * generalized Lorenz gauge condition for the 
---
>  * (vector potential prescription), using the
>  * generalized Lorenz gauge condition for the
7c7
<  * Based originally on the Illinois GRMHD code, 
---
>  * Based originally on the Illinois GRMHD code,
10,11c10,11
<  * primarily by Zachariah Etienne, Yuk Tung Liu, 
<  * and Vasileios Paschalidis. 
---
>  * primarily by Zachariah Etienne, Yuk Tung Liu,
>  * and Vasileios Paschalidis.
13c13
<  * Rewritten for public release in 2013 
---
>  * Rewritten for public release in 2013
17c17
<  * Original unigrid GRMHD evolution prescription: 
---
>  * Original unigrid GRMHD evolution prescription:
32c32
<  * This version of PPM implements the standard 
---
>  * This version of PPM implements the standard
34c34
<  * to have 3 ghostzones instead of 4. 
---
>  * to have 3 ghostzones instead of 4.
36a37
> 
45a47
> #include "IllinoisGRMHD_EoS_lowlevel_functs.C"
47a50
> 
56c59
<   
---
> 
63c66,77
<   // FIXME: only for single gamma-law EOS.
---
> 
>   /**********************************
>    * Piecewise Polytropic EOS Patch *
>    *   Setting up the EOS struct    *
>    **********************************/
>   /*
>    * The short piece of code below takes care
>    * of initializing the EOS parameters.
>    * Please refer to the "inlined_functions.C"
>    * source file for the documentation on the
>    * function.
>    */
65,72c79,80
<   eos.neos=neos;
<   eos.K_poly=K_poly;
<   eos.rho_tab[0]=rho_tab[0];
<   eos.P_tab[0]=P_tab[0];
<   eos.gamma_th=gamma_th;
<   eos.eps_tab[0]=eps_tab[0];
<   eos.k_tab[0]=k_tab[0];   eos.k_tab[1]=k_tab[1];
<   eos.gamma_tab[0]=gamma_tab[0]; eos.gamma_tab[1]=gamma_tab[1];
---
>   initialize_EOS_struct_from_input(eos);
> 
106a115
> 
112a122
> 
145a156
> 
162a174
> 
167c179,180
<   compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,eos,
---
>   compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,
>                                                         eos, Gamma_th,
170a184
> 
177a192
> 
179,180c194,195
<    * 1) Computation of \partial_x on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}}, 
<    *    via PPM reconstruction onto (i-1/2,j,k), so that 
---
>    * 1) Computation of \partial_x on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}},
>    *    via PPM reconstruction onto (i-1/2,j,k), so that
205a221
> 
220c236,237
<   // Note that we have already reconstructed vx and vy along the x-direction, 
---
> 
>   // Note that we have already reconstructed vx and vy along the x-direction,
237,239c254,256
<   /* There are two stories going on here: 
<    * 1) Computation of \partial_y on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}}, 
<    *    via PPM reconstruction onto (i,j-1/2,k), so that 
---
>   /* There are two stories going on here:
>    * 1) Computation of \partial_y on RHS of \partial_t {rho_star,tau,mhd_st_{x,y,z}},
>    *    via PPM reconstruction onto (i,j-1/2,k), so that
264c281
<   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, 
---
>   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
280c297
<   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, 
---
>   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
281a299
> 
295a314
> 
298c317
<    * We want to compute 
---
>    * We want to compute
314c333
<         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= 
---
>         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
320a340
> 
324c344,345
<   // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix 
---
> 
>   // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix
337c358
<   /* There are two stories going on here: 
---
>   /* There are two stories going on here:
361c382
<   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, 
---
>   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
377c398
<   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, 
---
>   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
384a406
> 
396,397c418,419
<   in_prims[VXR]=out_prims_r[VX]; 
<   in_prims[VZR]=out_prims_r[VZ]; 
---
>   in_prims[VXR]=out_prims_r[VX];
>   in_prims[VZR]=out_prims_r[VZ];
404a427
> 
407c430
<    * We want to compute 
---
>    * We want to compute
421c444
<         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= 
---
>         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
427a451
> 
430a455
> 
447c472
<   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, 
---
>   reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,
452c477
<    * We want to compute 
---
>    * We want to compute
466c491
<         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= 
---
>         temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=
472c497,498
<   
---
> 
> 
475a502
> 
498c525
<   Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi, 
---
>   Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi,
505c532
<   // FUN DEBUGGING TOOL (trust me!): 
---
>   // FUN DEBUGGING TOOL (trust me!):
532a560
>

Step 5.b: `driver_evaluate_MHD_rhs.h` [Back to top]¶

$\label{code_validation_driver_evaluate_MHD_rhs__h}$

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

In [27]:

# 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/driver_evaluate_MHD_rhs.h"
original_IGM_file_name = "driver_evaluate_MHD_rhs-original.h"
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__driver_evaluate_MHD_rhs__h  = !diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__h

if Validation__driver_evaluate_MHD_rhs__h == []:
    # 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 driver_evaluate_MHD_rhs.h: PASSED!")
else:
    # If the validation fails, we keep the original IGM source code file
    print("Validation test for driver_evaluate_MHD_rhs.h: 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__driver_evaluate_MHD_rhs__h:
        print(diff_line)

Validation test for driver_evaluate_MHD_rhs.h: FAILED!
Diff:
6c6
< static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct, 
---
> static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,
12c12
<  CCTK_REAL **TUPmunu,eos_struct &eos,
---
>  CCTK_REAL **TUPmunu,eos_struct &eos, CCTK_REAL Gamma_th,
34a35
>

Step 6: 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__driver_evaluate_MHD_rhs.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means).

In [28]:

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__driver_evaluate_MHD_rhs.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
!rm -f Tut*.out Tut*.aux Tut*.log

Tutorial-IllinoisGRMHD: driver_evaluate_MHD_rhs.C¶

Authors: Leo Werneck & Zach Etienne¶

In this tutorial module we explain the driver functions that compute the right-hand side (RHS) of the MHD equations 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: Load up necessary ETK and IllinoisGRMHD header files [Back to top]¶

Step 3: The IllinoisGRMHD_driver_evaluate_MHD_rhs() function [Back to top]¶

Step 3.1: Equation of state (EOS) parameters [Back to top]¶

Step 3.2: Set pointers to GRMHD gridfunctions [Back to top]¶

Step 3.3: Convert ADM variables to BSSN variables [Back to top]¶

Step 3.4: Setting up pointers to the metric and stress-energy tensor gridfunctions [Back to top]¶

Step 3.5: Initialization of the RHS variables [Back to top]¶

Step 3.6: Compute extrinsic curvature terms from the RHS of ∂t˜τ\partial_{t}\tilde\tau and TμνT^{\mu\nu} [Back to top]¶

Step 3.7: Computing ftilde{\rm ftilde} [Back to top]¶

Step 3.8: The RHSs of ρ⋆\rho_{\star}, ˜τ\tilde\tau, ˜Si\tilde{S}_{i}, and AiA_{i} [Back to top]¶

Step 3.8.1: Reconstructing {vx,vy,Bystagger}\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\} along the xx-direction [Back to top]¶

Step 3.8.2: Evaluating ∂xF\partial_{x}\boldsymbol{F} [Back to top]¶

Step 3.8.3: Reconstructing {vx,vy,Bystagger}\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\} along the yy-direction [Back to top]¶

Step 3.8.4: Evaluating ∂yF\partial_{y}\boldsymbol{F} [Back to top]¶

Step 3.8.5: Evaluating [∂tAz]no gauge terms\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms} [Back to top]¶

Step 3.8.6: Multiple reconstructions [Back to top]¶

Step 3.8.7: Evaluating ∂zF\partial_{z}\boldsymbol{F} [Back to top]¶

Step 3.8.8: Evaluating [∂tAx]no gauge terms\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms} [Back to top]¶

Step 3.8.9: Evaluating [∂tAy]no gauge terms\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms} [Back to top]¶

Step 3.8.10: Evaluating ∂t[ψ6Φ]\partial_{t}\left[\psi^{6}\Phi\right] and [∂tAi]gauge terms\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms} [Back to top]¶

Step 4: The driver_evaluate_MHD_rhs.h header file [Back to top]¶

Step 5: Code validation [Back to top]¶

Step 5.a: driver_evaluate_MHD_rhs.C [Back to top]¶

Step 5.b: driver_evaluate_MHD_rhs.h [Back to top]¶

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

Tutorial-IllinoisGRMHD: `driver_evaluate_MHD_rhs.C`¶

In this tutorial module we explain the driver functions that compute the right-hand side (RHS) of the MHD equations within `IllinoisGRMHD`¶

Step 2: Load up necessary ETK and `IllinoisGRMHD` header files [Back to top]¶

Step 3: The `IllinoisGRMHD_driver_evaluate_MHD_rhs()` function [Back to top]¶

Step 3.6: Compute extrinsic curvature terms from the RHS of $\partial_{t}\tilde\tau$ and $T^{\mu\nu}$ [Back to top]¶

Step 3.7: Computing ${\rm ftilde}$ [Back to top]¶

Step 3.8: The RHSs of $\rho_{\star}$ , $\tilde\tau$ , $\tilde{S}_{i}$ , and $A_{i}$ [Back to top]¶

Step 3.8.1: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $x$ -direction [Back to top]¶

Step 3.8.2: Evaluating $\partial_{x}\boldsymbol{F}$ [Back to top]¶

Step 3.8.3: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $y$ -direction [Back to top]¶

Step 3.8.4: Evaluating $\partial_{y}\boldsymbol{F}$ [Back to top]¶

Step 3.8.5: Evaluating $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

Step 3.8.7: Evaluating $\partial_{z}\boldsymbol{F}$ [Back to top]¶

Step 3.8.8: Evaluating $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

Step 3.8.9: Evaluating $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$ [Back to top]¶

Step 3.8.10: Evaluating $\partial_{t}\left[\psi^{6}\Phi\right]$ and $\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms}$ [Back to top]¶

Step 4: The `driver_evaluate_MHD_rhs.h` header file [Back to top]¶

Step 5.a: `driver_evaluate_MHD_rhs.C` [Back to top]¶

Step 5.b: `driver_evaluate_MHD_rhs.h` [Back to top]¶

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