Iharm3D: Vectorized General Relativistic Magnetohydrodynamics

Iharm3D is an open-source C code for simulating black hole accretion systems in arbitrary stationary spacetimes using ideal general-relativistic magnetohydrodynamics (GRMHD). It is an implementation of the HARM ("High Accuracy Relativistic Magnetohydrodynamics") algorithm outlined in Gammie et al. (2003) with updates as outlined in McKinney&Gammie (2004) and Noble et al. (2006). The code is most directly derived from Ryan et al. (2015) but with radiative transfer portions removed. HARM is a conservative finite-volume scheme for solving the equations of ideal GRMHD, a hyperbolic system of partial differential equations, on a logically Cartesian mesh in arbitrary coordinates.

iharm3D Functionality and Purpose iharm3D 1 is an open-source C code for simulating black hole accretion systems in arbitrary stationary spacetimes using ideal general-relativistic magnetohydrodynamics (GRMHD). It is an implementation of the HARM ("High Accuracy Relativistic Magnetohydrodynamics") algorithm outlined in Gammie et al. (2003) with updates as outlined in McKinney & Gammie (2004) and Noble et al. (2006). The code is most directly derived from Ryan et al. (2015) but with radiative transfer portions removed. HARM is a conservative finite-volume scheme for solving the equations of ideal GRMHD, a hyperbolic system of partial differential equations, on a logically Cartesian mesh in arbitrary coordinates.

Statement of Need
Numerical simulations are crucial in modeling observations of active galactic nuclei, such as the recent horizon-scale results from the Event Horizon Telescope and GRAVITY collaborations. The computational simplicity of ideal GRMHD enables the generation of long, high-resolution simulations and broad parameter-exploration studies that can be compared to observations for parameter inference.
Multiple codes already exist for solving the ideal GRMHD equations on regular Eulerian meshes in 3D, including: As the length of this list illustrates, the field of GRMHD simulation is now well established, and many codes now exist to serve different needs. These codes can be distinguished by the tradeoffs they make in prioritizing speed, simplicity, and generality, with the latter encompassing, e.g., support for dynamical spacetimes, adaptive mesh refinement, or higher-order integration schemes.
In particular, iharm3D aims to be a simple and fast code capable of simulating the original systems of interest when designing HARM, even at the cost of features aimed at more general applicability. It provides a fast and scalable update to HARM, but maintains the conventions and structure of the original described in Gammie et al. (2003). The result is a code that is relatively easy to understand and modify, yet capable of running simulations at state-of-the-art scale.

Implementation Notes
In MHD, uncorrected discretization errors inevitably lead to violations of the no-monopoles condition ∇ · B = 0. As in the original HARM implementation, iharm3D uses the "Flux-CT" scheme for cell-centered constrained transport outlined in Tóth (2000).
iharm3D also retains numerical evaluation of all metric-dependent quantities, allowing trivial modification of the coordinate system or background spacetime so long as the line element is available in analytic form. This can be used as a form of static mesh refinement, since the coordinates can be adapted to place resolution in areas of interest (e.g., near the accretion disk midplane).
In GRMHD, "conserved" variables (energy and momentum densities) are complicated analytic functions of "primitive" variables (density, pressure, and velocity). Conserved variables are stepped forward in time and then inversion to primitives is done numerically. iharm3D uses the "1D W " scheme outlined in Noble et al. (2006).
As the equations of ideal GRMHD are rescalable, any consistent set of units may be chosen to evolve them in the code. For numerical stability, when simulating accretion systems we choose units in which GM = c = 1 with M the mass of the central object, and scale the density of the initial conditions such that the maximum value of ρ is 1.
To model a collisionless plasma, iharm3D implements an optional means of tracking and partitioning dissipation into ions and electrons (Ressler et al. (2015)). Currently the code implements five different heating models, described in Howes (2010) To avoid catastrophic failures caused by discretization error, especially in low density regions, fluid variables are bounded at the end of each step. Typically, the bounds in black hole accretion problems are enforced as follows: • Density ρ > 10 −6 k, for k ≡ 1 r 2 (1+r/10) , with r the radial coordinate, • Internal energy density u > 10 −8 k γ where γ ≡ adiabatic index, • ρ and u are incremented until σ ≡ 2P b ρ < 400 and β ≡ Pgas P b > 2.5 × 10 −5 where P b ≡ b 2 2 is the magnetic pressure, • ρ is incremented until u ρ < 100, • When evolving electron temperatures, u is decremented until

Tests
The convergence properties of HARM are well-studied in Gammie et al. (2003). iharm3D implements most of the tests presented in that paper as integration and regression tests. Figure  1 shows convergence results for linear modes and for un-magnetized Bondi flow. Fast wave 2 4 2 6 Alfvén wave Figure 1: Results of convergence tests with iharm3D's main branch, plotting L1 norm of the difference between the computed solution and the analytic or stable result with increasing domain size. Wave solutions were performed on a 3D cubic grid N zones to one side, the Bondi accretion problem was performed on a logically Cartesian 2D square grid N zones on one side.
iharm3D implements three additional tests which check that fluid evolution is identical for different domain decompositions: one which initializes a new fluid state, one which restarts from a checkpoint file, and one comparing the initialized state to an equivalent checkpoint file.

Scaling
Key iharm3D routines are written for effective compiler vectorization and prioritize simple memory access patterns to make good use of high memory bandwidth. Originally developed for Intel Knights Landing (KNL) chips on the Stampede2 supercomputer at Texas Advanced Computing Center (TACC), iharm3D also runs efficiently on TACC Frontera, which uses traditional Cascade Lake (CLX) CPUs. Figure 2 presents scaling results for iharm3D on both Stampede2 and Frontera.  Research projects using iharm3D iharm3D is one of several GRMHD codes used by the EHT Collaboration to produce its library of fluid simulations. Images produced from this library were used for validation tests in Event