  We develop a new numerical scheme for ideal magnetohydrodynamic (MHD)
simulations, which is robust against one- and multi-dimensional shocks, and is
accurate for low Mach number flows and discontinuities. The scheme belongs to a
family of the advection upstream splitting method employed in computational
aerodynamics, and it splits the inviscid flux in MHD equations into advection,
pressure, and magnetic tension parts, and then individually evaluates mass,
pressure, and magnetic tension fluxes at the interface of a computational cell.
The mass flux is designed to avoid numerical shock instability in
multidimension, while preserving contact discontinuity. The pressure flux
possesses a proper scaling for low Mach number flows, allowing reliable
simulations of nearly incompressible flows. The magnetic tension flux is built
to be consistent with the HLLD approximate Riemann solver to preserve
rotational discontinuity. We demonstrate various benchmark tests to verify the
novel performance of the scheme. Our results indicate that the scheme must be a
promising tool to tackle astrophysical systems that include both low and high
Mach number flows, as well as magnetic field inhomogeneities.
