  We propose a new Harten-Lax-van Leer discontinuities (HLLD) approximate
Riemann solver to improve the stability of shocks and the accuracy of low-speed
flows in multidimensional magnetohydrodynamic (MHD) simulations. Stringent
benchmark tests verify that the new solver is more robust against a numerical
shock instability and is more accurate for low-speed, nearly incompressible
flows than the original solver, whereas additional computational costs are
quite low. The novel ability of the new solver enables us to tackle MHD
systems, including both high and low Mach number flows.
