PaperlayerView exact revision on arXiv

Source-linked AI summary

An Extension of the Athena++ Code Framework for GRMHD Based on Advanced Riemann Solvers and Staggered-Mesh Constrained Transport

Christopher J. White, James M. Stone, Charles F. Gammie

arXiv:1511.00943v2astro-ph.HEastro-ph.IM

TL;DR

Existing GRMHD codes commonly use relatively diffusive Riemann solvers, while combining advanced solvers with staggered constrained transport in arbitrary stationary spacetimes remains limited. This paper extends Athena++ with HLLC/HLLD solvers through local frame transformations and staggered-mesh constrained transport, then demonstrates reliable tests, black-hole simulations, and performance scalability. The framework’s scope is bounded by ongoing development of higher-order reconstruction, higher-order time integration, and radiative-transfer capabilities.

  • Problem

    Many GRMHD codes use simple, more diffusive Riemann solvers, while advanced solvers and staggered constrained transport had not been combined in GRMHD.

  • Method

    The paper develops an Athena++ GRMHD algorithm for arbitrary stationary spacetimes using local frame transformations for HLLC/HLLD solvers and staggered-mesh constrained transport.

  • Results

    The code passes general-relativistic tests, simulates equilibrium and turbulent accretion flows around spinning black holes, and demonstrates performance and scalability.

  • Takeaways & Limitations

    The framework supports accurate GRMHD simulations with advanced solvers and divergence-free magnetic-field evolution across stationary curved spacetimes.

  • Takeaways & Limitations

    Higher-order reconstruction, higher-order temporal integration, and GR-compatible radiative transfer were still under development.

Abstract

from arXiv · show

We present a new general relativistic magnetohydrodynamics (GRMHD) code integrated into the Athena++ framework. Improving upon the techniques used in most GRMHD codes, ours allows the use of advanced, less diffusive Riemann solvers, in particular HLLC and HLLD. We also employ a staggered-mesh constrained transport algorithm suited for curvilinear coordinate systems in order to maintain the divergence-free constraint of the magnetic field. Our code is designed to work with arbitrary stationary spacetimes in one, two, or three dimensions, and we demonstrate its reliability in a number of tests. We also report on its promising performance and scalability.

1. INTRODUCTION

GRMHD simulations require algorithms that evolve fluids accurately in curved spacetimes. The paper combines advanced Riemann solvers with staggered constrained transport in Athena++ to address accuracy and magnetic-field divergence control.

  • Motivation: GRMHD simulations are needed for astrophysical phenomena occurring in sufficiently strong gravitational fields.Examples include black-hole accretion, collapsars, and merging neutron-star binaries.
  • Godunov framework: Finite-volume Godunov methods evolve cell-averaged conserved quantities and obtain interface fluxes from Riemann problems.Their shock-capturing accuracy and stability motivate their use in compressible-fluid simulations.
  • Algorithmic challenges: Riemann-solver accuracy is critical, while simple solvers such as local Lax-Friedrichs and HLLE can be overly diffusive for subsonic flows.The paper therefore targets more advanced approximate solvers for GRMHD.
  • Magnetic-field constraint: The magnetic-field divergence-free constraint is enforced through constrained transport, which maintains a discretized version of the constraint using staggered fields.Alternative approaches include divergence cleaning and vector-potential evolution.
  • Paper contribution: The authors extend the Gardiner–Stone staggered constrained-transport algorithm to arbitrary stationary coordinate systems within Athena++.The implementation supports magnetic fields at cell interfaces and is part of a GRMHD algorithm designed for multiple dynamical regimes.

2. EQUATIONS

The paper formulates GRMHD in flux-conservative form on stationary curved-coordinate grids. Hydrodynamic variables are volume averaged, while magnetic fields use area averages and are updated from edge-centered electric fields.

  • Governing equations: The GRMHD equations are presented in a form that makes the paper’s discretization and notation explicit.The formulation rewrites covariant equations using partial derivatives plus connection-related source terms.
  • Physical variables: The ideal-MHD variables include fluid density, pressure, enthalpy, four-velocity, and magnetic-field quantities related through the stress-energy and electromagnetic tensors.In ideal MHD, electric fields are inferred from magnetic fields and fluid velocities.
  • Finite-volume discretization: Stationary metrics permit time-independent, coordinate-dependent cell volumes and face areas for finite-volume updates.These geometric quantities define the grid averages used by the discretized equations.
  • Finite-volume discretization: The timestep update advances conserved quantities using fluxes and source terms evaluated at intermediate integration times.The intermediate time level n+1/2 depends on the integration substep.
  • Magnetic-field discretization: Magnetic fields are stored as area averages, and edge-centered electric fields update them through the constrained-transport discretization.This choice is central to the Evans–Hawley constrained-transport scheme.

3. NUMERICAL ALGORITHM

The numerical algorithm addresses nonlinear primitive-variable recovery, higher-order interface reconstruction, and magnetic-field handling within a finite-volume GRMHD update. It uses iterative inversion, slope-limited reconstruction, and staggered magnetic variables.

  • 3.1. Variable Inversion: Variable inversion recovers primitive variables from conserved variables despite nonlinear relations and physically admissible-state restrictions.The problem includes possible nonconvergence and states with negative densities, pressures, or superluminal velocities.
  • 3.1. Variable Inversion: In GRMHD, primitive recovery uses projected spatial velocities and Newton–Raphson iteration on the relativistic enthalpy W = γ^2ρh.The projected variables describe subluminal motion for all values and include ũ^0 = 0.
  • 3.1. Variable Inversion: Special-relativistic MHD recovery likewise uses Newton–Raphson iteration, while special-relativistic hydrodynamics can instead solve a quartic equation for |v|.The quartic has coefficients a4 through a0.
  • 3.1. Variable Inversion: Recovered densities, pressures, and adiabatic indices are bounded by admissibility floors and limits, with failed iterations assigned floor values.The conserved variables are reevaluated after any such modification.
  • 3.2. Reconstruction: Reconstruction infers left and right interface states from a larger stencil while preserving discontinuities and avoiding spurious extrema.The modified van Leer limiter uses four neighboring cell values and is TVD for orthogonal coordinate systems.
  • 3.2. Reconstruction: The longitudinal magnetic field is already defined at the interface, so reconstruction is applied only to transverse magnetic components.The same procedure is used on the other coordinate surfaces.

3.3. Frame Transformation

The method transforms interface states from arbitrary stationary coordinates into a locally Minkowski orthonormal frame, applies special-relativistic Riemann solvers, and transforms fluxes back to the coordinate basis.

  • Frame construction: Left and right interface states are explicitly transformed into a locally Minkowski frame so special-relativistic Riemann algorithms can be reused.Scalar density and pressure transform trivially, while vector quantities require a metric-orthonormal basis.
  • Frame construction: The basis is orthonormal, normalized, normal to constant-time surfaces, and spatially aligned with the interface direction.These properties make the transformed problem compatible with one-dimensional Riemann fluxes.
  • Flux interpretation: The basis conditions ensure that the Riemann solver’s fluxes correspond to evolution in the coordinate-time direction and allow transverse fluxes to be omitted.Only the x-fluxes are needed to infer the flux across the selected coordinate interface.
  • Flux transformation: Transformation matrices map tensors between the orthonormal and coordinate bases, including stress-energy fluxes and four-velocity components.The same construction applies to interfaces normal to the other coordinate directions by cyclic permutation.
  • Moving interfaces: A constant-coordinate interface generally moves in the orthonormal frame when the shift has a component normal to that interface.The interface velocity is related to the shift and lapse through v̂^x = β^1/(αγ^11).
  • Scope and caveat: The normal-observer transformation remains applicable inside ergospheres and horizon-penetrating event horizons, but can become singular at coordinate singularities.Polar-axis singularities therefore require careful treatment.

3.4. Riemann solver

The code supports HLL-family Riemann solvers in a locally Minkowski frame, including less diffusive HLLC and HLLD variants that resolve internal wave structure. Fluxes are selected from the appropriate wavefan region, including for moving interfaces.

  • HLL solvers use left and right states plus leftgoing and rightgoing signal speeds, while guaranteeing physically admissible intermediate states for admissible inputs.Signal speeds are chosen from sound speeds in hydrodynamics or fast magnetosonic speeds in MHD.
  • LLF and HLLE use a single intermediate state and are computationally simple, but can be diffusive for flows slower than the extremal wavespeed.HLLE differs by allowing distinct leftgoing and rightgoing signal speeds; the two solvers usually produce similar results.
  • HLLC resolves the contact discontinuity, while HLLD additionally resolves Alfvén waves to lessen numerical diffusion at the cost of nonlinear root finds.HLLC is implemented for relativistic hydrodynamics, whereas HLLD is used for relativistic MHD.
  • Figure 1 distinguishes extremal waves, contact discontinuities, and Alfvén waves across LLF, HLLC, and HLLD wavefans.Fluxes for a stationary interface are evaluated along x = 0, which may cross 3, 4, or 6 wavefan regions.
  • The solver constructs a wavefan with internal wave speeds, conserved states, and normal fluxes, then selects the region relevant to the interface.For moving interfaces, the relevant region is determined by comparing interface velocity with the wavespeeds.

3.5. Constrained Transport (CT)

The constrained-transport scheme evolves face-centered magnetic fields to preserve the divergence-free constraint without divergence cleaning. It constructs edge-centered electric fields from face fields, cell-centered fields, gradients, and upwinding in arbitrary coordinates.

  • The implementation extends the Gardiner–Stone constrained-transport algorithm to relativity and uses a face-centered magnetic-field update to maintain divergence freedom without divergence cleaning.It differs from the original formulation by using cell-centered rather than face-centered velocities.
  • The five-step procedure obtains face-centered electric fields, computes cell-centered fields and gradients, upwinds gradients, and combines them into an edge-centered field.For E3, the final combination uses four face-centered fields and four upwinded gradients.
  • Figure 2 depicts how face-centered fields, cell-centered fields, gradients, and upwinding produce the edge-centered electric field E3.The diagram is a four-cell timeslice with the x3 direction out of the plane.
  • Face-centered electric fields are taken from Riemann-solver fluxes, while cell-centered electric fields are computed from velocities and interpolated magnetic fields.The fluxes are transformed back to global coordinates before supplying the face-centered fields.
  • Electric-field gradients are computed between cell centers and faces, then shifted onto faces by upwinding according to the sign of the mass flux.The procedure requires cell widths, which can approximately cancel from the final expression when metric and grid spacing vary smoothly.

3.6. Source Terms

The method accounts for geometric source terms in non-Cartesian coordinates while using a second-order van Leer update for conserved variables. Stationarity removes geometric energy source terms, and metric symmetries reduce storage requirements.

  • Non-Cartesian updates include geometric source terms for the energy-momentum equations, while the continuity and magnetic-field equations have none.Additional physics source terms, such as heating and cooling, can be added alongside geometric terms.
  • Because the metric is stationary, geometric source terms in the energy equation vanish.The remaining momentum source terms are obtained by contracting stress-energy components with connection coefficients.
  • The code uses a temporally second-order van Leer integrator that computes primitives, solves Riemann problems, evaluates source terms, and advances conserved variables in two stages.The van Leer integrator is TVD and does not introduce spurious extrema.
  • Metric quantities can be precomputed for stationary spacetimes, and symmetries plus separable formulas avoid full three-dimensional geometric storage.The required geometric data include cell volumes, interface areas, edge lengths, cell widths, metric components, and connection coefficients.
  • For Kerr–Schild coordinates, the code stores 22N_r + 32N_θ + 9N_φ + 20N_rN_θ metric values and requires no more than 2D storage for discussed metrics.This is substantially smaller than storing all evolving variables, which requires 20–50 N_1N_2N_3 values.

4. TESTS

The tests assess linear-wave convergence, nonlinear shock-tube behavior, coordinate-system consistency, and curved-spacetime convergence. Advanced solvers reduce errors in selected waves while the code retains second-order convergence and reproduces expected multidimensional and stationary solutions.

  • Linear Waves: Both hydrodynamic and MHD linear-wave tests exhibit second-order convergence with the expected second-order reconstruction and time integration.
  • Linear Waves: HLLC reduces the general-relativistic hydrodynamic entropy-wave error by a factor of 2.8 relative to HLLE.
  • Linear Waves: HLLD reduces general-relativistic MHD errors relative to HLLE by factors of 4.8, 1.8, and 1.4 for entropy, slow, and Alfvén waves.
  • Shock Tubes: The HLLD shock-tube calculations reproduce expected nonlinear behavior while retaining some known artifacts and reducing selected overshoots at the same resolution.
  • Coordinate Tests: Blast-wave solutions remain qualitatively consistent across Minkowski and transformed snake coordinates, with differences generally confined to cell-level shock-front variations and instability fingering.
  • Curved-Spacetime Tests: The curved-spacetime stationary tests achieve full second-order spatial convergence for hydrodynamics with HLLC and MHD with HLLD.
  • Curved-Spacetime Tests: A highly nontrivial-spacetime test also finds second-order convergence using a frame-transforming HLLE solver, supporting correct equation solving and frame transformations.

5. TORUS SIMULATIONS

The torus simulations test equilibrium preservation, convergence, and long-term behavior in hydrodynamic and magnetized Fishbone–Moncrief tori. The results show second-order convergence, improved equilibrium errors with HLLC, and broadly comparable but resolution-dependent turbulent profiles.

  • Hydrodynamic torus: Second-order convergence is observed for density errors in 2D Fishbone–Moncrief hydrodynamic tori.The test uses Nr × Nr/2 grids and excludes edge effects through a density-based region T.
  • Hydrodynamic torus: At t = 430, the equilibrium error is 0.030 with HLLE and 0.021 with HLLC on the 64 × 32 grid.This corresponds to three orbital periods at the pressure maximum.
  • Magnetized torus: The 3D magnetized torus is evolved to t = 12,000 at effective resolution 288 × 256 × 192, allowing MRI-driven turbulence to disrupt the initial torus.The magnetic field is initialized from a purely azimuthal vector potential and normalized to a maximum gas-to-magnetic pressure ratio of 100.
  • Magnetized torus: Electron temperature broadly matches the reference magnitude and radial trend, while the simulated disk is hotter inside r = 2.04420 and shows a kink near r = 5–6.The gas-to-magnetic pressure profile has a similar shape but reaches about 200 in outer regions, compared with about 60 in the reference cases.
  • Magnetized torus: The simulated disks have larger azimuthal correlation lengths than the reference, especially at large radius, but correlation lengths decrease as resolution increases.The authors interpret the larger outer-region lengths as evidence that the fully saturated turbulent state has not yet been reached there.

6. PERFORMANCE

The performance evaluation measures Athena++ across geometries, fluid models, and Riemann solvers, while describing optimization and scaling strategies. MHD is slower than hydrodynamics, realistic GR metrics add cost, and scaling remains strong to 6144 cores.

  • Optimization: Precomputation, memory-layout-aware array sweeps, and innermost-loop vectorization are used to improve Athena++ performance.The implementation minimizes expensive metric evaluations and cache misses while addressing vectorization barriers in complex solver functions.
  • Single-core performance: Realistic Kerr-Schild geometry incurs the full GR cost because nontrivial geometric factors enter most calculations.The SR-to-GR slowdown in Minkowski coordinates reflects the cost of GR variable inversion and wavespeed formulas.
  • Single-core performance: MHD problems run at 1/4 to 3/4 the speed of pure hydrodynamics for fixed geometry and Riemann solver.HLLE runs at 93 ± 6% the speed of corresponding LLF tests, with similar accuracies.
  • Scaling: At 6144 cores, hydrodynamics retains over 97% of its 24-core performance and MHD is indistinguishable from 100%.The main observed penalty is 20% per core when moving from one core to a full 24-core node; adding nodes has negligible cost.

7. SUMMARY

The Athena++ GRMHD implementation combines advanced HLL solvers with staggered-mesh constrained transport and supports arbitrary stationary spacetimes. Tests establish reliability, while performance and scalability support higher-accuracy and higher-resolution simulations.

  • The code supports advanced HLLC and HLLD Riemann solvers through local frame transformations and staggered-mesh constrained transport for arbitrary stationary spacetimes.These techniques had not previously been combined in GRMHD according to the authors.
  • The framework’s modularity supports static mesh refinement, while adaptive refinement, higher-order reconstruction, higher-order time integration, and GR-compatible radiative transfer remain under development.These capabilities define the current development boundary of the framework.
  • Tests cover standard MHD problems, general-relativistic coordinate transformations, Schwarzschild Bondi accretion, and fluids around spinning black holes.The spinning-black-hole problem exercises hydrodynamics, magnetic fields, metric terms, frame transformations, geometrical source terms, and horizon-penetrating coordinates together.
  • Performance and scalability are intended to enable numerical analyses using both more accurate methods and improved resolution.The authors identify magnetically arrested disks and tilted disks as example applications.
  • Figure 17 reports performance per core for a three-dimensional GRMHD simulation on a cluster.

A. FRAME TRANSFORMATION

The frame transformation constructs an orthonormal basis adapted to an interface in a stationary curved spacetime. Its geometric constraints simplify flux conversion by eliminating the need for transverse fluxes from the special-relativistic Riemann solver.

  • Basis properties: The transformed basis is required to be mutually orthogonal, normalized, future-directed in time, and adapted to constant-coordinate hypersurfaces.The timelike basis vector is orthogonal to constant-x0 surfaces, while the spatial vectors satisfy the stated interface geometry.
  • Interface construction: The construction is developed for an interface at constant x1, with the other directions obtained by cyclic permutation.
  • Spatial basis choice: The x-direction is chosen to enforce the projection condition rather than simply aligning with the coordinate basis or making it orthogonal to the constant-x1 hypersurface.This choice avoids complicating the conversion of fluxes back to the global frame.
  • Normalization: The basis vectors are normalized using the metric, with the timelike vector selected to point forward in time and identified with the unit normal of a 3+1 foliation.The construction is checked against the Minkowski metric for sensible orientation and normalization.
  • Flux simplification: The y- and z-basis vectors lie in the constant-x0 plane, and the resulting geometry means the special-relativistic Riemann solver needs no y- or z-fluxes.The chosen orientation completes the transformation matrix and its inverse.
  • Matrix form: The transformation matrices are lower-diagonal, reflecting the imposed basis properties and the selected spatial rotation.
Loading 1511.00943v2…