  The magnetohydrodynamics (MHD) equations are continuum models used in the
study of a wide range of plasma physics systems, including the evolution of
complex plasma dynamics in tokamak disruptions. However, efficient numerical
solution methods for MHD are extremely challenging due to disparate time and
length scales, strong hyperbolic phenomena, and nonlinearity. Therefore the
development of scalable, implicit MHD algorithms and high-resolution adaptive
mesh refinement strategies is of considerable importance. In this work, we
develop a high-order stabilized finite-element algorithm for the reduced
visco-resistive MHD equations based on the MFEM finite element library
(mfem.org). The scheme is fully implicit, solved with the Jacobian-free
Newton-Krylov (JFNK) method with a physics-based preconditioning strategy. Our
preconditioning strategy is a generalization of the physics-based
preconditioning methods in [Chacon, et al, JCP 2002] to adaptive, stabilized
finite elements. Algebraic multigrid methods are used to invert sub-block
operators to achieve scalability. A parallel adaptive mesh refinement scheme
with dynamic load-balancing is implemented to efficiently resolve the
multi-scale spatial features of the system. Our implementation uses the MFEM
framework, which provides arbitrary-order polynomials and flexible adaptive
conforming and non-conforming meshes capabilities. Results demonstrate the
accuracy, efficiency, and scalability of the implicit scheme in the presence of
large scale disparity. The potential of the AMR approach is demonstrated on an
island coalescence problem in the high Lundquist-number regime ($\ge 10^7$)
with the successful resolution of plasmoid instabilities and thin current
sheets.
