A highly accurate drag solver for multi-fluid dust and gas hydrodynamics on GPUs

TL;DR

This work introduces a GPU-optimized, high-precision drag solver for multi-fluid dust–gas hydrodynamics, enabling accurate friction calculations via a matrix-exponential treatment of the drag subsystem. The key idea is to solve $\partial \tilde{\bm{U}}/\partial t = \bm{\Omega}\tilde{\bm{U}}$ with $\tilde{\bm{U}}(t)=e^{t\bm{\Omega}}\tilde{\bm{U}}(t_0)$ and to compute the matrix exponential using a scaling-and-squaring approach built on Taylor/Pade expansions, achieving machine-precision accuracy while keeping the drag step efficient. The solver is implemented in two GPU codes, Dyablo (Kokkos) and Shamrock (Sycl), and validated with Dustybox, Dustywave, and Dustyshock tests, demonstrating stable performance with up to 16 dust species. Performance analyses on A100 GPUs show the drag step costs remain manageable relative to the hydro step, supporting future inclusion of growth and fragmentation and enabling broader applications to radiative transfer or chemistry. The approach and codes are publicly available, highlighting its practical impact for exascale astrophysical simulations.

Abstract

Exascale supercomputing unleashes the potential for simulations of astrophysical systems with unprecedented resolution. Taking full advantage of this computing power requires the development of new algorithms and numerical methods that are GPU friendly and scalable. In the context of multi-fluid dust-gas dynamics, we propose a highly accurate algorithm that is specifically designed for GPUs. We developed a multi-fluid gas-dust algorithm capable of computing friction terms on GPU architectures to machine precision, with the constraint for the drag-time step to remain a fraction of the global hydrodynamic time step for computational efficiency in practice. We present a scaling-and-squaring algorithm tailored to modern architectures for computing the exponential of the drag matrix, enabling high accuracy in friction calculations across relevant astrophysical regimes. The algorithm was validated through the Dustybox Dustywave and Dustyshock tests. The algorithm was implemented and tested in two multi-GPU codes with different architectures and GPU programming models: Dyablo, an adaptive mesh refinement code based on the Kokkos library, and Shamrock, a multi-method code based on Sycl. On current architectures, the friction computation remains acceptable for both codes (below the typical hydro time step) up to 16 species, enabling a further implementation of growth and fragmentation. This algorithm might be applied to other physical processes, such as radiative transfer or chemistry.

Figures (7)

• Figure 1: Collisions test. Evolution of the gas and dust velocities (from top to bottom) for tests B and C (left and right columns, respectively). The dot-dashed black line represents the analytical solution. The other lines compare the numerical solutions obtained with the EXP integrator (orange) with other integrators, namely the first-order implicit RK1 (blue), the second-order fully implicit RK2 (purple), the MDIRK method with $\gamma = 1 - 1/\sqrt{2}$ (green), and the MDIRK method with $\gamma = 1 + 1/\sqrt{2}$ (red)
• Figure 2: Collisions tests B (left) and C (right). The evolution of error $Er_1$ is defined by Eq. (\ref{['eq:Error1_eq']}) for the EXP (orange) and other integrators: RK1 (blue), RK2 (purple), MKDIRK- (green), and MKDIRK+ (red).
• Figure 3: Evolution of the gas and dust velocities for a Collisions test with 20 dust fluids obtained with the EXP solver and compared to the analytical solution given by Eq. (\ref{['eq:V_analytic']}) (dashed black lines). For the first five dust species, their stopping time is short, such that they are strongly coupled to the gas. We therefore chose to plot the first and the fifth species. In addition to the $6$th dust species, we chose to skip every second fluid from the $10$th to the $20$th dust species.
• Figure 4: Dustywave test. Evolution of normalised densities (upper panels) and velocities (lower panels) for one (left panel) and four dust fluids (right panel). The dashed black lines correspond to the analytical solution. The other markers and colours represent the numerical results obtained with the EXP algorithm. We additionally plot the L1 error given by Eq. (\ref{['eq:Error2_eq']}) for the one-dust fluid case as an inset in the bottom left panel using a fixed time step of $\Delta t = 10^{-4}$ for EXP and RK1.
• Figure 5: Dustyshock test. Normalised velocities (first and third panels for one- and three-dust fluid(s), respectively) and density (second and fourth panels for one- and three-dust fluid(s), respectively). The dashed lines show the analytical solutions, and the circle markers show the numerical results we obtained with the EXP algorithm. With the EXP solver, the shock is captured within one to two cells.

Theorems & Definitions (1)

• theorem 1