Early numerical analysis of MHD began with the seminal paper [1]}, which proves existence and convergence of a discrete solution to the stationary, incompressible MHD equations in convex domains with \(C^{1,1}\)  boundary. The paper considers \(\inf \) -\(\sup \)  stable mixed elements for the kinetic variables, and \(H^1\) -conforming approximations of the magnetic unknowns. A more appropriate functional setting for the magnetic variables in more general domains is that of \(H(\operatornamewithlimits{curl};\Omega )\)  . Such a setting is provided by [2]}, which proposes a mixed finite element approximation shown to converge optimally in non-convex domains under a small data assumption. Convergence of a stabilised finite element method (FEM) for linearised MHD is shown in [3]}, where the nonlinear problem is approximated by a coupled Picard algorithm which requires resolving the linear case at each iteration. A convergent method for the non-stationary problem is provided by [4]}. The early work of [5]} provides analysis and numerical theory of MHD with non-ideal boundary, which resolves some of the many difficulties of modelling realistic magnetic flow problems. A review of the various numerical methods at the time for the discretisation of three-dimensional incompressible MHD is given in [6]}. A more detailed description of the known results for MHD is provided by [7]}. Due to the saddle-point structure of both the fluid and magnetic equations, many of the aforementioned conforming methods require elements satisfying an \(\inf \) -\(\sup \)  condition. This becomes problematic as the Reynolds and magnetic Reynolds number become too large. A stabilised finite element discretisation is proposed in [8]} and shown to converge unconditionally to the physical solution – even in the presence of singular solutions. One of the earliest works focused on conserving the divergence-free constraints on the magnetic field is found in [9]}, which provides a stabilised mixed method for the non-stationary problem. However, only the magnetic equations are discretised, and the fluid velocity is assumed to be prescribed. More recently, a stabilised FEM for stationary MHD is designed in [10]} which preserves the divergence free constraints of both the velocity and magnetic fields at the discrete level. We also make note of the recent publications [11]}, [12]}, [13]}. A non-conforming approximation of the linearised model is proposed in [14]} using a mixed discontinuous Galerkin (DG) approach. Following this work, a DG approximation of the non-linear model is provided by [15]} which conserves the divergence free conditions. The DG approach has the drawback of leading to a large number of degrees of freedom. This issue can be alleviated, as shown in [15]}, through the use of a hybridisable discontinuous Galerikin (HDG) approximation, for which element degrees of freedom can be locally eliminated through static condensation. The recent article [17]} devises a HDG method which locally conserves the divergence free constraints of both the fluid and magnetic field, and satisfies a global momentum balance. For standard diffusion models, HDG and Hybrid-High Order methods share many similarities [18]}; however, the HDG method of [15]}, [17]} use hybrid spaces also for the fluid and magnetic pressure, which is not the usual choice in HHO for Navier–Stokes, and leads to many more degrees of freedom – as expected to achieve improved conservation of divergence (see discussion in Remark  below.
