  We describe a high-order numerical magnetohydrodynamics (MHD) solver built
upon a novel non-linear entropy stable numerical flux function that supports
eight travelling wave solutions. By construction the solver conserves mass,
momentum, and energy and is entropy stable. The method is designed to treat the
divergence-free constraint on the magnetic field in a similar fashion to a
hyperbolic divergence cleaning technique. The solver described herein is
especially well-suited for flows involving strong discontinuities. Furthermore,
we present a new formulation to guarantee positivity of the pressure. We
present the underlying theory and implementation of the new solver into the
multi-physics, multi-scale adaptive mesh refinement (AMR) simulation code
$\texttt{FLASH}$ (http://flash.uchicago.edu). The accuracy, robustness and
computational efficiency is demonstrated with a number of tests, including
comparisons to available MHD implementations in $\texttt{FLASH}$.
