Multidimensional finite volume methods with constrained transport for MHD simulations of the ISM

Author: Mervin Yap

The interstellar medium (ISM) is composed of complex physics and turbulent features. In order to study the dynamics of the ISM, computational simulations are often needed, where numerical methods are needed to solve the physical equations of the fluid. For the case of the ISM, we assume the ideal magnetohydrodynamical (MHD) equations. In this assumption, the compressible Euler equations, along with the induction equation, are solved. Other cases can involve solving the relativistic hydrodynamics equations for cosmological simulations, non-ideal MHD equations for plasma flows, or shallow water equations for atmospheric/oceanic simulations.

In the finite volume (FV) method, the solutions to the compressible Euler equations are computed using the cell averages and the fluxes are then used to evolve the hydrodynamical components. The cell averages are often expanded with higher-order approximations and then evolved with a higher-order time discretisation to achieve better spatial and temporal accuracy respectively (see Monthly Highlight in October 2025). In addition, the upwind constrained transport (UCT) method (Londrillo & Del Zanna, 2004) evolves the magnetic fields by computing the (averaged) solutions to the induction equation.

To preserve the higher-order accuracy (> 2nd-order) of the overall simulation, all components in the numerical scheme are required to have this high-order accuracy. The central WENO (CWENO) scheme (Levy et al., 1999) and the modified WENO (WENO-Z) scheme (Borges et al., 2008) are some examples of spatial reconstructions. Examples of temporal discretisation include the strong stability-preserving Runge-Kutta (SSP-RK) methods (Gottlieb et al., 2001), the diagonally split Runge-Kutta (DSRK) method (Horváth, 1998), and Patankar-type methods (Patankar, 1980). Even the UCT method requires inverse spatial reconstructions for its averaged values to maintain the accuracy.

Figure 1 shows the density of a two-dimensional MHD FV-UCT simulation with (Ornstein-Uhlenbeck) turbulence driving performed with the astrea code to test the stability of the numerical methods. The simulation involves the 4th-order WENO-Z spatial reconstruction, the 4th-order SSP-RK (5,4) temporal discretisation, and the local Lax-Friedrich approximate Riemann solver on a fixed, uniform (staggered) Cartesian grid with 512^2 grid cells. In the simulation, the numerical scheme was able to handle large Mach numbers, sometimes up to M ~ 200, without the simulation crashing (stable). Furthermore, the simulation does not seem to produce NaN values (robust), which is usually an indication of negative densities/pressures since they create infinite or complex wavespeeds. The UCT method was also able to preserve the numerical accuracy of the magnetic field within the simulation to machine precision (prevent the creation of magnetic monopoles). More validation of the accuracy can be performed with dissipation scale plots for the turbulence cascade or conservation plots for the variables.

At present, we are incorporating these methods into the DISPATCH code (Nordlund et al., 2018). The aim would be to have a GPU-accelerated code with these robust and stable numerical methods that can produce highly accurate simulations at larger scales.