IllinoisGRMHD: A Compact, Dynamic-Spacetime General Relativistic Magnetohydrodynamics Code for Easy User Adoption

Original author: Zachariah B. Etienne.

Date: 2015-10-12 12:00:00 -0600 (Mon, 12 Oct 2015)

Abstract

IllinoisGRMHD solves the equations of General Relativistic MagnetoHydroDynamics (GRMHD) using a high-resolution shock capturing scheme. It is a rewrite of the Illinois Numerical Relativity (ILNR) group’s GRMHD code, and generates results that agree to roundoff error with that original code. Its feature set coincides with the features of the ILNR group’s recent code (ca. 2009–2014), which was used in their modeling of the following systems:

  1. Magnetized circumbinary disk accretion onto binary black holes
  2. Magnetized black hole–neutron star mergers
  3. Magnetized Bondi flow, Bondi-Hoyle-Littleton accretion
  4. White dwarf–neutron star mergers

IllinoisGRMHD is particularly good at modeling GRMHD flows into black holes without the need for excision. Its HARM-based conservative-to-primitive solver has also been modified to check the physicality of conservative variables prior to primitive inversion, and move them into the physical range if they become unphysical.

1 Introduction

Currently IllinoisGRMHD consists of

  1. the Piecewise Parabolic Method (PPM) for reconstruction,
  2. the Harten, Lax, van Leer (HLL/HLLE) approximate Riemann solver, and
  3. a modified HARM Conservative-to-Primitive solver (see REQUIRED CITATION #2 below).

IllinoisGRMHD evolves the vector potential Aμ (on staggered grids) instead of the magnetic fields (Bi) directly, to guarantee that the magnetic fields will remain divergenceless even at AMR boundaries. On uniform resolution grids, this vector potential formulation produces results equivalent to those generated using the standard, staggered flux-CT scheme. This scheme is based on that of Del Zanna (2003, see below OPTIONAL CITATION #1).

For further information about motivations, basic equations, how IllinoisGRMHD works, as well as basic code test results, please see the IllinoisGRMHD code announcement paper [1]. If you use IllinoisGRMHD for your research, you are asked to include the REQUIRED CITATIONS listed below in your citations.

For a quick “Guide to Getting Started”, please visit the IllinoisGRMHD web page:
http://math.wvu.edu/~zetienne/ILGRMHD/

2 Citations

2.1 Required citations

  1. IllinoisGRMHD code announcement paper [1]
  2. Harm con2prim paper [2]

2.2 Optional citations

  1. MHD paper [3]

3 Acknowledgements

Note that IllinoisGRMHD is based on the GRMHD code of the Illinois Numerical Relativity group (ca. 2014), 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.

References

[1]   Z. B. Etienne, V. Paschalidis, R. Haas, P. Mösta and S. L. Shapiro, “IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes,” Class. Quant. Grav. 32, 175009 (2015) doi:10.1088/0264-9381/32/17/175009 [arXiv:1501.07276 [astro-ph.HE]].

[2]   S. C. Noble, C. F. Gammie, J. C. McKinney and L. Del Zanna, “Primitive variable solvers for conservative general relativistic magnetohydrodynamics,” Astrophys. J. 641, 626-637 (2006) doi:10.1086/500349 [arXiv:astro-ph/0512420 [astro-ph]].

[3]   L. Del Zanna, N. Bucciantini and P. Londrillo, “An efficient shock-capturing central-type scheme for multidimensional relativistic flows. 2. Magnetohydrodynamics,” Astron. Astrophys. 400, 397-414 (2003) doi:10.1051/0004-6361:20021641 [arXiv:astro-ph/0210618 [astro-ph]].

4 Parameters




damp_lorenz
Scope: private  REAL



Description: Damping factor for the generalized Lorenz gauge. Has units of 1/length = 1/M. Typically set this parameter to 1.5/(maximum Delta t on AMR grids).



Range   Default: 0.0
*:*
any real






conserv_to_prims_debug
Scope: restricted  INT



Description: 0: no, 1: yes



Range   Default: (none)
0:1
zero (no) or one (yes)






em_bc
Scope: restricted  KEYWORD



Description: EM field boundary condition



Range   Default: copy
copy
Copy data from nearest boundary point
frozen
Frozen boundaries






gamma_speed_limit
Scope: restricted  REAL



Description: Maximum relativistic gamma factor.



Range   Default: 10.0
1:*
Positive > 1, though you’ll likely have troubles far above 10.






gamma_th
Scope: restricted  REAL



Description: thermal gamma parameter



Range   Default: -1
0:*
Physical values
-1
forbidden value to make sure it is explicitly set in the parfile






k_poly
Scope: restricted  REAL



Description: initial polytropic constant



Range   Default: 1.0
0:*
Positive






matter_bc
Scope: restricted  KEYWORD



Description: Chosen Matter boundary condition



Range   Default: outflow
outflow
Outflow boundary conditions
frozen
Frozen boundaries






neos
Scope: restricted  INT



Description: number of parameters in EOS table. If you want to increase from the default max value, you MUST also set eos_params_arrays1 and eos_params_arrays2 in interface.ccl to be consistent!



Range   Default: 1
1:10
Any integer between 1 and 10






psi6threshold
Scope: restricted  REAL



Description: Where Psi6̂ > Psi6threshold, we assume we’re inside the horizon in the primitives solver, and certain limits are relaxed or imposed



Range   Default: 1e100
*:*
Can set to anything






rho_b_atm
Scope: restricted  REAL



Description: Floor value on the baryonic rest mass density rho_b (atmosphere). Given the variety of systems this code may encounter, there *is no reasonable default*. Your run will die unless you override this default value in your initial data thorn.



Range   Default: 1e200
*:*
Allow for negative values. This enables us to debug the code and verify if rho_b_atm is properly set.






rho_b_max
Scope: restricted  REAL



Description: Ceiling value on the baryonic rest mass density rho_b. The enormous value effectively disables this ceiling by default. It can be quite useful after a black hole has accreted a lot of mass, leading to enormous densities inside the BH. To enable this trick, set rho_b_max in your initial data thorn! You are welcome to change this parameter mid-run (after restarting from a checkpoint).



Range   Default: 1e300
0:*
Note that you will have problems unless rho_b_atm<rho_b_max






sym_bz
Scope: restricted  REAL



Description: In-progress equatorial symmetry support: Symmetry parameter across z axis for magnetic fields = +/- 1



Range   Default: 1.0
-1.0:1.0
Set to +1 or -1.






symmetry
Scope: restricted  KEYWORD



Description: Currently only no symmetry supported, though work has begun in adding equatorial-symmetry support. FIXME: Extend ET symmetry interface to support symmetries on staggered gridfunctions



Range   Default: none
none
no symmetry, full 3d domain






tau_atm
Scope: restricted  REAL



Description: Floor value on the energy variable tau (cf. tau_stildefix_enable). Given the variety of systems this code may encounter, there *is no reasonable default*. Effectively the current (enormous) value should disable the tau_atm floor. Please set this in your initial data thorn, and reset at will during evolutions.



Range   Default: 1e100
0:*
Positive






update_tmunu
Scope: restricted  BOOLEAN



Description: Update Tmunu, for RHS of Einstein’s equations?



  Default: yes






verbose
Scope: restricted  KEYWORD



Description: Determines how much evolution information is output



Range   Default: essential+iteration output
no
Complete silence
essential
”Essential health monitoring of the GRMHD evolution: Information about conservative-to-prim itive fixes, etc.”
see [1] below
Outputs health monitoring information, plus a record of which RK iteration. Very useful for backtracing a crash.



[1]

essential+iteration output




lapse_timelevels
Scope: shared from ADMBASE INT






metric_timelevels
Scope: shared from ADMBASE INT






shift_timelevels
Scope: shared from ADMBASE INT



5 Interfaces

General

Implements:

illinoisgrmhd

Inherits:

admbase

boundary

spacemask

tmunubase

hydrobase

grid

Grid Variables

5.0.1 PRIVATE GROUPS




  Group Names     Variable Names    Details   




diagnostic_gfs   compact0
failure_checker   descriptionGridfunction to track conservative-to-primitives solver fixes. Beware that this gridfunction is overwritten at each RK substep.
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




grmhd_primitives_reconstructed_temps   compact0
ftilde_gf   descriptionTemporary variables used for primitives reconstruction
temporary   dimensions3
rho_br   distributionDEFAULT
Pr   group typeGF
vxr   tagsprolongation=”none” Checkpoint=”no”
vyr   timelevels1
vzr  variable typeREAL




grmhd_conservatives_rhs   compact0
rho_star_rhs   descriptionStorage for the right-hand side of the partial_t rho_star
    descriptionpartial_t tau
rho_star_rhs   descriptionand partial_t tilde{S}_i equations. Needed for MoL timestepping.
tau_rhs   dimensions3
st_x_rhs   distributionDEFAULT
st_y_rhs   group typeGF
st_z_rhs   tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




em_ax_rhs   compact0
Ax_rhs   descriptionStorage for the right-hand side of the partial_t A_x equation. Needed for MoL timestepping.
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




em_ay_rhs   compact0
Ay_rhs   descriptionStorage for the right-hand side of the partial_t A_y equation. Needed for MoL timestepping.
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




em_az_rhs   compact0
Az_rhs   descriptionStorage for the right-hand side of the partial_t A_z equation. Needed for MoL timestepping.
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL








  Group Names     Variable Names    Details   




em_psi6phi_rhs   compact0
psi6phi_rhs   descriptionStorage for the right-hand side of the partial_t (psi6̂ Phi) equation. Needed for MoL timestepping.
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




grmhd_cmin_cmax_temps   compact0
cmin_x   descriptionStore min and max characteristic speeds in all three directions.
cmax_x   dimensions3
cmin_y   distributionDEFAULT
cmax_y   group typeGF
cmin_z   tagsprolongation=”none” Checkpoint=”no”
cmax_z   timelevels1
 variable typeREAL




grmhd_flux_temps   compact0
rho_star_flux   descriptionTemporary variables for storing the flux terms of tilde{S}_i.
tau_flux   dimensions3
st_x_flux   distributionDEFAULT
st_y_flux   group typeGF
st_z_flux   tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




tupmunu   compact0
TUPtt   descriptionT{̂mu nu}
    descriptionstored to avoid expensive recomputation
TUPtx   dimensions3
TUPty   distributionDEFAULT
TUPtz   group typeGF
TUPxx   tagsprolongation=”none” Checkpoint=”no”
TUPxy   timelevels1
TUPxz  variable typeREAL




5.0.2 PUBLIC GROUPS




  Group Names     Variable Names    Details   




grmhd_conservatives   compact0
rho_star   descriptionEvolved mhd variables
tau   dimensions3
mhd_st_x   distributionDEFAULT
mhd_st_y   group typeGF
mhd_st_z   timelevels3
 variable typeREAL




em_ax   compact0
Ax   descriptionx-component of the vector potential
    descriptionevolved when constrained_transport_scheme==3
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsProlongation=”STAGGER011”
  timelevels3
 variable typeREAL




em_ay   compact0
Ay   descriptiony-component of the vector potential
    descriptionevolved when constrained_transport_scheme==3
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsProlongation=”STAGGER101”
  timelevels3
 variable typeREAL




em_az   compact0
Az   descriptionz-component of the vector potential
    descriptionevolved when constrained_transport_scheme==3
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsProlongation=”STAGGER110”
  timelevels3
 variable typeREAL




em_psi6phi   compact0
psi6phi   descriptionsqrt{gamma} Phi
    descriptionwhere Phi is the em scalar potential
  dimensions3
  distributionDEFAULT
  group typeGF
  tagsProlongation=”STAGGER111”
  timelevels3
 variable typeREAL




grmhd_primitives_allbutbi   compact0
rho_b   descriptionPrimitive variables density
    descriptionpressure
rho_b   descriptionand components of three velocity vî. Note that vî is defined in terms of 4-velocity as: vî = uî/u0̂. Note that this definition differs from the Valencia formalism.
P   dimensions3
vx   distributionDEFAULT
vy   group typeGF
vz   tagsInterpNumTimelevels=1 prolongation=”none”
  timelevels1
 variable typeREAL








  Group Names     Variable Names    Details   




grmhd_primitives_bi   compact0
Bx   descriptionB-field components defined at vertices.
By   dimensions3
Bz   distributionDEFAULT
  group typeGF
  tagsInterpNumTimelevels=1 prolongation=”none”
  timelevels1
 variable typeREAL




grmhd_primitives_bi_stagger   compact0
Bx_stagger   descriptionB-field components defined at staggered points [Bx_stagger at (i+1/2
    descriptionj
Bx_stagger   descriptionk)
Bx_stagger   descriptionBy_stagger at (i
Bx_stagger   descriptionj+1/2
Bx_stagger   descriptionk)
Bx_stagger   descriptionBz_stagger at (i
Bx_stagger   descriptionj
Bx_stagger   descriptionk+1/2)].
By_stagger   dimensions3
Bz_stagger   distributionDEFAULT
  group typeGF
  tagsInterpNumTimelevels=1 prolongation=”none”
  timelevels1
 variable typeREAL




eos_params_arrays1   compact0
rho_tab   descriptionEquation of state (EOS) parameters
P_tab   dimensions1
eps_tab   distributionCONSTANT
  group typeARRAY
  size10
  timelevels1
 variable typeREAL




eos_params_arrays2   compact0
k_tab   descriptionEquation of state (EOS) parameters
gamma_tab   dimensions1
  distributionCONSTANT
  group typeARRAY
  size11
  timelevels1
 variable typeREAL




eos_params_scalar   compact0
n_poly   descriptionpolytropic index
  dimensions0
  distributionCONSTANT
  group typeSCALAR
  timelevels1
 variable typeREAL




bssn_quantities   compact0
gtxx   descriptionBSSN quantities
    descriptioncomputed from ADM quantities
gtxy   dimensions3
gtxz   distributionDEFAULT
gtyy   group typeGF
gtyz   tagsprolongation=”none” Checkpoint=”no”
gtzz   timelevels1
gtupxx  variable typeREAL




Adds header:

IllinoisGRMHD_headers.h

Uses header:

Symmetry.h

6 Schedule

This section lists all the variables which are assigned storage by thorn WVUThorns/IllinoisGRMHD. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

 

Always:  
ADMBase::metric[metric_timelevels] ADMBase::curv[metric_timelevels] ADMBase::lapse[lapse_timelevels] ADMBase::shift[shift_timelevels]  
IllinoisGRMHD::BSSN_quantities  
grmhd_conservatives[3] em_Ax[3] em_Ay[3] em_Az[3] em_psi6phi[3]  
grmhd_primitives_allbutBi grmhd_primitives_Bi grmhd_primitives_Bi_stagger grmhd_primitives_reconstructed_temps grmhd_conservatives_rhs em_Ax_rhs em_Ay_rhs em_Az_rhs em_psi6phi_rhs grmhd_cmin_cmax_temps grmhd_flux_temps TUPmunu diagnostic_gfs 
eos_params_arrays1 eos_params_arrays2 eos_params_scalar  
   

Scheduled Functions

MoL_Register

  illinoisgrmhd_registervars

  register evolved, rhs variables in illinoisgrmhd for mol

 

 After: bssn_registervars
   lapse_registervars
 Language:c
 Options: meta
 Type: function

CCTK_BASEGRID

  illinoisgrmhd_initsymbound

  schedule symmetries

 

 After: lapse_initsymbound
 Language:c
 Type: function

HydroBase_Boundaries

  illinoisgrmhd_compute_b_and_bstagger_from_a

  compute b and b_stagger from a, sync: grmhd_primitives_bi,grmhd_primitives_bi_stagger

 

 After: outer_boundaries_on_a_mu
  Language:c
 Sync: grmhd_primitives_bi
   grmhd_primitives_bi_stagger
 Type: function

AddToTmunu

  illinoisgrmhd_conserv_to_prims

  compute primitive variables from conservatives. this is non-trivial, requiring a newton-raphson root-finder.

 

 After: compute_b_and_bstagger_from_a
 Language:c
 Type: function

AddToTmunu

  illinoisgrmhd_outer_boundaries_on_p_rho_b_vx_vy_vz

  apply outflow-only, flat bcs on {p,rho_b,vx,vy,vz}. outflow only == velocities pointed inward from outer boundary are set to zero.

 

 After: illinoisgrmhd_conserv_to_prims
 Language:c
 Sync: grmhd_primitives_allbutbi
 Type: function

CCTK_POSTPOSTINITIAL

  illinoisgrmhd_postpostinitial

  hydrobase_con2prim in cctk_postpostinitial sets conserv to prim then outer boundaries (obs, which are technically disabled). the post ob syncs actually reprolongate the conservative variables, making cons and prims inconsistent. so here we redo the con2prim, avoiding the sync afterward, then copy the result to other timelevels

 

 After: hydrobase_con2prim
 Before:mol_poststep
 Type: group

IllinoisGRMHD_PostPostInitial

  illinoisgrmhd_initsymbound

  schedule symmetries – actually just a placeholder function to ensure prolongations / processor syncs are done before outer boundaries are updated.

 

 Before: compute_b
 Language:c
 Sync: grmhd_conservatives
   em_ax
   em_ay
   em_az
   em_psi6phi
 Type: function

IllinoisGRMHD_PostPostInitial

  illinoisgrmhd_compute_b_and_bstagger_from_a

  compute b and b_stagger from a sync: grmhd_primitives_bi,grmhd_primitives_bi_stagger

 

 After: postid
   empostid
   lapsepostid
 Language:c
 Sync: grmhd_primitives_bi
   grmhd_primitives_bi_stagger
 Type: function

IllinoisGRMHD_PostPostInitial

  illinoisgrmhd_conserv_to_prims

  compute primitive variables from conservatives. this is non-trivial, requiring a newton-raphson root-finder.

 

 After: compute_b
 Language:c
 Type: function

IllinoisGRMHD_PostPostInitial

  illinoisgrmhd_postpostinitial_set_symmetries__copy_timelevels

  compute post-initialdata quantities

 

 After: compute_b
   p2c
 Language:c
 Type: function

MoL_CalcRHS

  illinoisgrmhd_driver_evaluate_mhd_rhs

  evaluate rhss of gr hydro & grmhd equations

 

 After: bssn_rhs
    shift_rhs
 Language:c
 Type: function

HydroBase_Boundaries

  illinoisgrmhd_initsymbound

  schedule symmetries – actually just a placeholder function to ensure prolongations / processor syncs are done before outer boundaries are updated.

 

 Before: compute_b_postrestrict
 Language:c
 Sync: grmhd_conservatives
   em_ax
   em_ay
   em_az
   em_psi6phi
 Type: function

HydroBase_Boundaries

  illinoisgrmhd_outer_boundaries_on_a_mu

  apply linear extrapolation bcs on a_{mu}, so that bcs are flat on bî.

 

 After: boundary_syncs
 Before: mhd_conserv2prims_postrestrict
 Language:c
 Type: function

Aliased Functions

 

Alias Name:         Function Name:
IllinoisGRMHD_InitSymBound Boundary_SYNCs
IllinoisGRMHD_PostPostInitial_Set_Symmetries__Copy_Timelevelsmhdpostid
IllinoisGRMHD_compute_B_and_Bstagger_from_A compute_b
IllinoisGRMHD_driver_evaluate_MHD_rhs IllinoisGRMHD_RHS_eval