  The magnetohydrodynamics (MHD) equations are generally known to be difficult
to solve numerically, due to their highly nonlinear structure and the strong
coupling between the electromagnetic and hydrodynamic variables, especially for
high Reynolds and coupling numbers. In this work, we present a scalable
augmented Lagrangian preconditioner for a finite element discretization of the
$\mathbf{B}$-$\mathbf{E}$ formulation of the incompressible viscoresistive MHD
equations. For stationary problems, our solver achieves robust performance with
respect to the Reynolds and coupling numbers in two dimensions and good results
in three dimensions. We extend our method to fully implicit methods for
time-dependent problems which we solve robustly in both two and three
dimensions. Our approach relies on specialized parameter-robust multigrid
methods for the hydrodynamic and electromagnetic blocks. The scheme ensures
exactly divergence-free approximations of both the velocity and the magnetic
field up to solver tolerances. We confirm the robustness of our solver by
numerical experiments in which we consider fluid and magnetic Reynolds numbers
and coupling numbers up to 10,000 for stationary problems and up to 100,000 for
transient problems in two and three dimensions.
