MHDuet : a high-order General Relativistic Radiation MHD code for CPU and GPU architectures

TL;DR

MHDuet addresses the challenge of simulating fully relativistic magnetohydrodynamics with neutrino transport in compact objects by coupling the CCZ4 formulation of Einstein gravity to a magnetized fluid and an M1 neutrino scheme. The code employs Simflowny to auto-generate highly optimized AMReX- or SAMRAI-based implementations, enabling high-order HRSC methods, LES, IMEX time integration, and sub-cycling on CPU and GPU architectures. Validation on cold and hot neutron-star tests, plus binary neutron-star mergers with tabulated EoS, demonstrates convergence faster than second order, accurate reproduction of quasi-normal modes, and stable gravitational-wave signals. Performance benchmarks reveal dramatic GPU speedups (over an order of magnitude in some configurations) and excellent strong/weak scaling, confirming MHDuet as a powerful tool for multi-messenger astrophysics and large-scale surveys of compact objects.

Abstract

We present MHDuet, an open source evolution code for general relativistic magnetohydrodynamics with neutrino transport. The code solves the full set of Einstein equations coupled to a relativistic, magnetized fluid with an M1 neutrino radiation scheme using advanced techniques, including adaptive mesh and large eddy simulation techniques, to achieve high accuracy. The Simflowny platform generates the code from a high-level specification of the computational system, producing code that runs with either the SAMRAI or AMReX infrastructure. The choice of AMReX enables compilation and execution on GPUs, running an order of magnitude faster than on CPUs at the node level. We validate the code against benchmark tests, reproducing previous results obtained with the SAMRAI infrastructure, and demonstrate its capabilities with simulations of neutron stars employing realistic tabulated equations of state. Resolution studies clearly demonstrate convergence faster than second order in the grid spacing. Scaling tests reveal excellent strong and weak scaling performance when running on GPUs. The goal of the code is to provide a powerful tool for studying the dynamics of compact objects within multi-messenger astrophysics.

Figures (12)

• Figure 1: Synchronization between two levels with no sub-cycling (left) and with any of the Berger-Oliger schemes (right), using a Runge-Kutta with three substeps. Prolongation operations between levels is denoted with red arrows, while that restriction is marked with green ones.
• Figure 2: Single NS with the LS-K220 EoS.Top panel: The central rest-mass density, defined as $\rho_{c}(t)\equiv \rho(r\!=\!0,t)$, is normalized by its initial value and plotted as a function of time. The different curves represent solutions obtained with versions of the code running on SAMRAI, AMReX-CPU, and AMReX-GPU infrastructures. The agreement demonstrates that the SAMRAI and AMReX versions of the code generate the same results. Bottom panel: The power spectral density of the central density for the low and high resolutions evolved with AMReX-GPU. Vertical lines indicate frequencies obtained by previous codes as discussed in Section \ref{['sec:coldstar']}.
• Figure 3: Rotating, hot, magnetized star using the DD2 tabulated EoS. Snapshots on a meridional plane: from left to right, density, magnetic field and temperature, at the initial (top) and final time (bottom) of our simulation. The agreement within the stellar interior for the two times indicates a stationary solution is achieved as discussed in Section \ref{['sec:hotstar']}.
• Figure 4: Rotating, hot magnetized NS using the DD2 tabulated EoS. Stellar properties at the center: from top to bottom, density, $z$-component of the magnetic field, and temperature. The evolution oscillates about a unique, stationary solution with an amplitude slowly decreasing with time. Tiny differences between the high and low resolution runs become apparent only at late times. The density and magnetic field converge between second and third order, whereas the temperature converges even faster.
• Figure 5: Rotating, hot magnetized NS using the DD2 tabulated EoS. Estimates of the accuracy of the solution for the high-resolution simulation. The fractional change in baryonic mass of the star $\Delta M/ M_0$ (i.e., which ideally should remain constant during the evolution) and the L2-norm of the energy and momentum constraint violations. Also displayed are the time integrals of the energy constraint residual, $\theta$, and the solenoidal constraint violation, properly normalized (i.e., $\phi/|B|$).