General Implicit Runge-Kutta Integrators for Multifluid Gas-Dust Aerodynamic Drag

TL;DR

This work addresses the computational challenge of integrating aerodynamic drag in multifluid gas–dust systems with many dust species. It introduces the General Implicit Runge-Kutta (GIRK) integrator, designed to work with Strang splitting and general-purpose hydrodynamics codes, while preserving linear complexity in the number of dust bins $N_ ext{d}$ by building on the analytical drag framework of Krapp_2020/2024. Through two Strang-splitting schemes and a tunable set of parameters, GIRK achieves high-order accuracy and strong asymptotic stability across stiff and non-stiff regimes, and its performance scales linearly with $N_ ext{d}$ as demonstrated on standard benchmarks (DUSTYBOX, DUSTYWAVE, DUSTYSHOCK) and a steady-state 1D shearing box. Compared to the MDIRK method, GIRK offers greater implementational flexibility in general hydrodynamics codes while maintaining competitive accuracy, enabling efficient simulations with large numbers of dust sizes for astrophysical applications such as protoplanetary disks.

Abstract

The integration of aerodynamic drag is a fundamental step in simulating dust dynamics in hydrodynamical simulations. We propose a novel integration scheme, designed to be compatible with Strang splitting techniques, which allows for the straightforward integration of external forces and hydrodynamic fluxes in general-purpose hydrodynamic simulation codes. Moreover, this solver leverages an analytical solution to the problem of drag acceleration, ensuring linear complexity even in cases with multiple dust grain sizes, as opposed to the cubic scaling of methods that require a matrix inversion step. This new General Implicit Runge-Kutta integrator (GIRK) is evaluated using standard benchmarks for dust dynamics such as DUSTYBOX, DUSTYWAVE, and DUSTYSHOCK. The results demonstrate not only the accuracy of the method but also the expected scalings in terms of accuracy, convergence to equilibrium, and execution time. GIRK can be easily implemented in hydrodynamical simulations alongside hydrodynamical steps and external forces, and is especially useful in simulations with a large number of dust grain sizes.

Figures (8)

• Figure 1: Simulation outputs and analytical solutions for the three DUSTYBOX benchmark tests described in Huang_2022. The system is composed of a gas fluid and two dust fluids. The simulation outputs are shown as dots, while the analytical solution is drawn as a continuous line. For the non-stiff test A, the timestep is chosen according to the CFL condition Eq. \ref{['eq:CFL']}, while for the stiff tests B and C, the timestep is, respectively, $\Delta t = 0.005$ and $\Delta t = 0.05$. The integrator used (GIRK with $\mathcal{D}_{\Delta t/2}\mathcal{H}_{\Delta t}\mathcal{D}_{\Delta t/2}$ Strang splitting) is able to correctly capture the evolution of the mixture for all stiffness regimes. For Test C, we also show how much the integrator undershoots the analytical solution for the gas velocity, before converging.
• Figure 2: Convergence orders for the DIRK and the GIRK integrations for different Strang operator splitting schemes on Test A of Huang_2022. The error is computed using Eq. \ref{['eq:rel_error']}. All integrators follow their expected scalings, given the properties imposed as in Section \ref{['sec:impose_properties']} and \ref{['sec:parameters_GIRK']}.
• Figure 3: External force damping benchmark test. The dots represent the simulation output, while the solid line shows the analytical solution. The simulation was carried out using the GIRK drag integrator with the $\mathcal{D}_{\Delta t/2}\mathcal{H}_{\Delta t}\mathcal{D}_{\Delta t/2}$ Strang splitting scheme. As expected, the simulation follows the expected solution and converges to the correct asymptotic value.
• Figure 4: Relative error of the final simulation output compared to the analytical asymptotic solution, computed using Eq. \ref{['eq:rel_error']}. The MDIRK method performs best, while the DIRK integrator coupled with Strang splitting converges only at first order in the stiff regime. The GIRK integrators, however, show significantly improved accuracy and faster convergence across both splitting schemes.
• Figure 5: Simulation outputs and analytical solution of the DUSTYWAVE benchmark test for both a single dust specie (left panel) and four dust species (central panel). The plots show the time evolution of the normalized density perturbation. The simulation was carried out using the $\mathcal{D}_{\Delta t/2}\,\mathcal{H}_{\Delta t}\,\mathcal{D}_{\Delta t/2}$ Strang splitting scheme. In the right panel, we show the scaling of the error with the number of cells, for both Strang splitting schemes, computed on the test shown in the central panel.