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

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:

Magnetized circumbinary disk accretion onto binary black holes
Magnetized black hole–neutron star mergers
Magnetized Bondi flow, Bondi-Hoyle-Littleton accretion
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

IllinoisGRMHD evolves the vector potential

A_{μ}

(on staggered grids) instead of the magnetic fields (

B^{i}

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

2 Citations

2.1 Required citations

2.2 Optional citations

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.


lapse_timelevels	Scope: shared from ADMBASE	INT


metric_timelevels	Scope: shared from ADMBASE	INT


shift_timelevels	Scope: shared from ADMBASE	INT

5 Interfaces

General

Grid Variables

5.0.1 PRIVATE GROUPS


Group Names	Variable Names	Details

diagnostic_gfs		compact	0
	failure_checker	description	Gridfunction to track conservative-to-primitives solver fixes. Beware that this gridfunction is overwritten at each RK substep.
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

grmhd_primitives_reconstructed_temps		compact	0
	ftilde_gf	description	Temporary variables used for primitives reconstruction
	temporary	dimensions	3
	rho_br	distribution	DEFAULT
	Pr	group type	GF
	vxr	tags	prolongation=”none” Checkpoint=”no”
	vyr	timelevels	1
	vzr	variable type	REAL

grmhd_conservatives_rhs		compact	0
	rho_star_rhs	description	Storage for the right-hand side of the partial_t rho_star
		description	partial_t tau
	rho_star_rhs	description	and partial_t tilde{S}_i equations. Needed for MoL timestepping.
	tau_rhs	dimensions	3
	st_x_rhs	distribution	DEFAULT
	st_y_rhs	group type	GF
	st_z_rhs	tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

em_ax_rhs		compact	0
	Ax_rhs	description	Storage for the right-hand side of the partial_t A_x equation. Needed for MoL timestepping.
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

em_ay_rhs		compact	0
	Ay_rhs	description	Storage for the right-hand side of the partial_t A_y equation. Needed for MoL timestepping.
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

em_az_rhs		compact	0
	Az_rhs	description	Storage for the right-hand side of the partial_t A_z equation. Needed for MoL timestepping.
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL


Group Names	Variable Names	Details

em_psi6phi_rhs		compact	0
	psi6phi_rhs	description	Storage for the right-hand side of the partial_t (psi6^ Phi) equation. Needed for MoL timestepping.
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

grmhd_cmin_cmax_temps		compact	0
	cmin_x	description	Store min and max characteristic speeds in all three directions.
	cmax_x	dimensions	3
	cmin_y	distribution	DEFAULT
	cmax_y	group type	GF
	cmin_z	tags	prolongation=”none” Checkpoint=”no”
	cmax_z	timelevels	1
		variable type	REAL

grmhd_flux_temps		compact	0
	rho_star_flux	description	Temporary variables for storing the flux terms of tilde{S}_i.
	tau_flux	dimensions	3
	st_x_flux	distribution	DEFAULT
	st_y_flux	group type	GF
	st_z_flux	tags	prolongation=”none” Checkpoint=”no”
		timelevels	1
		variable type	REAL

tupmunu		compact	0
	TUPtt	description	T{^mu nu}
		description	stored to avoid expensive recomputation
	TUPtx	dimensions	3
	TUPty	distribution	DEFAULT
	TUPtz	group type	GF
	TUPxx	tags	prolongation=”none” Checkpoint=”no”
	TUPxy	timelevels	1
	TUPxz	variable type	REAL

5.0.2 PUBLIC GROUPS


Group Names	Variable Names	Details

grmhd_conservatives		compact	0
	rho_star	description	Evolved mhd variables
	tau	dimensions	3
	mhd_st_x	distribution	DEFAULT
	mhd_st_y	group type	GF
	mhd_st_z	timelevels	3
		variable type	REAL

em_ax		compact	0
	Ax	description	x-component of the vector potential
		description	evolved when constrained_transport_scheme==3
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	Prolongation=”STAGGER011”
		timelevels	3
		variable type	REAL

em_ay		compact	0
	Ay	description	y-component of the vector potential
		description	evolved when constrained_transport_scheme==3
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	Prolongation=”STAGGER101”
		timelevels	3
		variable type	REAL

em_az		compact	0
	Az	description	z-component of the vector potential
		description	evolved when constrained_transport_scheme==3
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	Prolongation=”STAGGER110”
		timelevels	3
		variable type	REAL

em_psi6phi		compact	0
	psi6phi	description	sqrt{gamma} Phi
		description	where Phi is the em scalar potential
		dimensions	3
		distribution	DEFAULT
		group type	GF
		tags	Prolongation=”STAGGER111”
		timelevels	3
		variable type	REAL

grmhd_primitives_allbutbi		compact	0
	rho_b	description	Primitive variables density
		description	pressure
	rho_b	description	and 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	dimensions	3
	vx	distribution	DEFAULT
	vy	group type	GF
	vz	tags	InterpNumTimelevels=1 prolongation=”none”
		timelevels	1
		variable type	REAL


Group Names	Variable Names	Details

grmhd_primitives_bi		compact	0
	Bx	description	B-field components defined at vertices.
	By	dimensions	3
	Bz	distribution	DEFAULT
		group type	GF
		tags	InterpNumTimelevels=1 prolongation=”none”
		timelevels	1
		variable type	REAL

grmhd_primitives_bi_stagger		compact	0
	Bx_stagger	description	B-field components defined at staggered points [Bx_stagger at (i+1/2
		description	j
	Bx_stagger	description	k)
	Bx_stagger	description	By_stagger at (i
	Bx_stagger	description	j+1/2
	Bx_stagger	description	k)
	Bx_stagger	description	Bz_stagger at (i
	Bx_stagger	description	j
	Bx_stagger	description	k+1/2)].
	By_stagger	dimensions	3
	Bz_stagger	distribution	DEFAULT
		group type	GF
		tags	InterpNumTimelevels=1 prolongation=”none”
		timelevels	1
		variable type	REAL

eos_params_arrays1		compact	0
	rho_tab	description	Equation of state (EOS) parameters
	P_tab	dimensions	1
	eps_tab	distribution	CONSTANT
		group type	ARRAY
		size	10
		timelevels	1
		variable type	REAL

eos_params_arrays2		compact	0
	k_tab	description	Equation of state (EOS) parameters
	gamma_tab	dimensions	1
		distribution	CONSTANT
		group type	ARRAY
		size	11
		timelevels	1
		variable type	REAL

eos_params_scalar		compact	0
	n_poly	description	polytropic index
		dimensions	0
		distribution	CONSTANT
		group type	SCALAR
		timelevels	1
		variable type	REAL

bssn_quantities		compact	0
	gtxx	description	BSSN quantities
		description	computed from ADM quantities
	gtxy	dimensions	3
	gtxz	distribution	DEFAULT
	gtyy	group type	GF
	gtyz	tags	prolongation=”none” Checkpoint=”no”
	gtzz	timelevels	1
	gtupxx	variable type	REAL

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

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

	After:	lapse_initsymbound
	Language:	c
	Type:	function

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

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

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

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

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

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

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

	After:	compute_b
	Language:	c
	Type:	function

	After:	compute_b
		p2c
	Language:	c
	Type:	function

	After:	bssn_rhs
		shift_rhs
	Language:	c
	Type:	function

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

Alias Name:	Function Name:
IllinoisGRMHD_InitSymBound	Boundary_SYNCs
IllinoisGRMHD_PostPostInitial_Set_Symmetries__Copy_Timelevels	mhdpostid
IllinoisGRMHD_compute_B_and_Bstagger_from_A	compute_b
IllinoisGRMHD_driver_evaluate_MHD_rhs	IllinoisGRMHD_RHS_eval

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