  Three-dimensional direct numerical simulations are used to study the energy
cascade rate in isothermal compressible magnetohydrodynamic turbulence. Our
analysis is guided by a two-point exact law derived recently for this problem
in which flux, source, hybrid, and mixed terms are present. The relative
importance of each term is studied for different initial subsonic Mach numbers
$M_S$ and different magnetic guide fields ${\bf B}_0$. The dominant
contribution to the energy cascade rate comes from the compressible flux, which
depends weakly on the magnetic guide field ${\bf B}_0$, unlike the other terms
whose modulus increase significantly with $M_S$ and ${\bf B}_0$. In particular,
for strong ${\bf B}_0$ the source and hybrid terms are dominant at small scales
with almost the same amplitude but with a different sign. A statistical
analysis made with an isotropic decomposition based on the SO(3) rotation group
is shown to generate spurious results in presence of ${\bf B}_0$, when compared
with an axisymmetric decomposition better suited to the geometry of the
problem. Our numerical results are compared with previous analyses made with
in-situ measurements in the solar wind and the terrestrial magnetosheath.
