  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 the first part of this work, we present a scalable augmented Lagrangian
preconditioner for a finite element discretisation 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. Our approach relies on specialised 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.
  In the second part, we focus on incompressible, resistive Hall MHD models and
derive structure-preserving finite element methods for these equations. We
present a variational formulation of Hall MHD that enforces the magnetic
Gauss's law precisely (up to solver tolerances) and prove the well-posedness of
a Picard linearisation. For the transient problem, we present time
discretisations that preserve the energy and magnetic and hybrid helicity
precisely in the ideal limit for two types of boundary conditions.
  In the third part, we investigate anisothermal MHD models. We start by
performing a bifurcation analysis for a magnetic Rayleigh--B\'enard problem at
a high coupling number $S=1{,}000$ by choosing the Rayleigh number in the range
between 0 and $100{,}000$ as the bifurcation parameter. We study the effect of
the coupling number on the bifurcation diagram and outline how we create
initial guesses to obtain complex solution patterns and disconnected branches
for high coupling numbers.
