Exponential Integrators for MHD: Matrix-free Leja interpolation and efficient adaptive time stepping

We propose a novel algorithm for the temporal integration of the magnetohydrodynamics (MHD) equations. The approach is based on exponential Rosenbrock schemes in combination with Leja interpolation. It naturally preserves Gauss's law for magnetism and is unencumbered by the stability constraints observed for explicit methods. Remarkable progress has been achieved in designing exponential integrators and computing the required matrix functions efficiently. However, employing them in realistic MHD scenarios require matrix-free implementations that are competent on modern computer hardware. We show how an efficient algorithm based on Leja interpolation that only uses the right-hand side of the differential equation (i.e. matrix-free), can be constructed. We further demonstrate that it outperforms, in the context of magnetic reconnection and the Kelvin--Helmholtz instability, earlier work on Krylov-based exponential integrators as well as explicit methods. Furthermore, an adaptive step size strategy is employed that gives an excellent and predictable performance, particularly in the lenient to intermediate tolerance regime that is often of importance in practical applications.