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, at
http://arxiv.org/abs/1501.07276. 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/

===================  
REQUIRED CITATIONS:  
===================

  1. IllinoisGRMHD code announcement paper: Class. Quantum Grav. 32 (2015) 175009,
    (http://arxiv.org/abs/1501.07276)
  2. Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L.  2006, Astrophysical Journal, 641, 626.
===================  
OPTIONAL CITATIONS:  
===================

  1. Del Zanna, Bucciantini & Londrillo A&A 400, 397 (2003)

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

3 Parameters




damp_lorenz
Scope: private  REAL



Description: Damping factor for Lorenz gauge. Has units of 1/length = 1/M. Try using the same value as used for the eta parameter, or smaller! Brian tried 0.1 with ADM mass of 1.



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






illinoisgrmhd_maxnumevolvedvars
Scope: restricted  INT



Description: The maximum number of evolved variables used by BSSN



Range   Default: 9
9:9
”Just 9: rho_star,tau,mhd_st_ x,mhd_st_y,mhd_st_z, Ax,Ay,Az, and psi6phi”






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: 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: Symmetry: 0=NOSYMM, 1=EQUATORIAL. FIXME: Use ET symmetry interface. PROBLEM: ET does not support staggered gridfunctions



Range   Default: none
none
no symmetry, full 3d domain
equatorial
NOT READY FOR PRIME TIME: reflection across z=0






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






mol_num_evolved_vars
Scope: shared from METHODOFLINES INT



4 Interfaces

General

Implements:

illinoisgrmhd

Inherits:

admbase

boundary

spacemask

tmunubase

hydrobase

grid

Grid Variables

4.0.1 PRIVATE GROUPS




  Group Names     Variable Names    Details   




diagnostic_gfs   compact0
failure_checker   dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




grmhd_primitives_reconstructed_temps   compact0
ftilde_gf   dimensions3
temporary   distributionDEFAULT
rho_br   group typeGF
Pr   tagsprolongation=”none” Checkpoint=”no”
vxr   timelevels1
vyr  variable typeREAL




grmhd_conservatives_rhs   compact0
rho_star_rhs   dimensions3
tau_rhs   distributionDEFAULT
st_x_rhs   group typeGF
st_y_rhs   tagsprolongation=”none” Checkpoint=”no”
st_z_rhs   timelevels1
 variable typeREAL




em_ax_rhs   compact0
Ax_rhs   dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




em_ay_rhs   compact0
Ay_rhs   dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




em_az_rhs   compact0
Az_rhs   dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL








  Group Names     Variable Names    Details   




em_psi6phi_rhs   compact0
psi6phi_rhs   dimensions3
  distributionDEFAULT
  group typeGF
  tagsprolongation=”none” Checkpoint=”no”
  timelevels1
 variable typeREAL




grmhd_cmin_cmax_temps   compact0
cmin_x   dimensions3
cmax_x   distributionDEFAULT
cmin_y   group typeGF
cmax_y   tagsprolongation=”none” Checkpoint=”no”
cmin_z   timelevels1
cmax_z  variable typeREAL




grmhd_flux_temps   compact0
rho_star_flux   dimensions3
tau_flux   distributionDEFAULT
st_x_flux   group typeGF
st_y_flux   tagsprolongation=”none” Checkpoint=”no”
st_z_flux   timelevels1
 variable typeREAL




tupmunu   compact0
TUPtt   dimensions3
TUPtx   distributionDEFAULT
TUPty   group typeGF
TUPtz   tagsprolongation=”none” Checkpoint=”no”
TUPxx   timelevels1
TUPxy  variable typeREAL




4.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   dimensions3
P   distributionDEFAULT
vx   group typeGF
vy   tagsInterpNumTimelevels=1 prolongation=”none”
vz   timelevels1
 variable typeREAL








  Group Names     Variable Names    Details   




grmhd_primitives_bi   compact0
Bx   dimensions3
By   distributionDEFAULT
Bz   group typeGF
  tagsInterpNumTimelevels=1 prolongation=”none”
  timelevels1
 variable typeREAL




grmhd_primitives_bi_stagger   compact0
Bx_stagger   dimensions3
By_stagger   distributionDEFAULT
Bz_stagger   group typeGF
  tagsInterpNumTimelevels=1 prolongation=”none”
  timelevels1
 variable typeREAL




eos_params_arrays1   compact0
rho_tab   descriptioneos_parameters
P_tab   dimensions1
eps_tab   distributionCONSTANT
  group typeARRAY
  size10
  timelevels1
 variable typeREAL




eos_params_arrays2   compact0
k_tab   descriptioneos_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   dimensions3
gtxy   distributionDEFAULT
gtxz   group typeGF
gtyy   tagsprolongation=”none” Checkpoint=”no”
gtyz   timelevels1
gtzz  variable typeREAL




Adds header:

IllinoisGRMHD_headers.h

Uses header:

Symmetry.h

5 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: 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