  We present a divergence-free and $H(div)$-conforming hybridized discontinuous
Galerkin (HDG) method and a computationally efficient variant called
embedded-HDG (E-HDG) for solving stationary incompressible viso-resistive
magnetohydrodynamic (MHD) equations. The proposed E-HDG approach uses
continuous facet unknowns for the vector-valued solutions (velocity and
magnetic fields) while it uses discontinuous facet unknowns for the scalar
variable (pressure and magnetic pressure). This choice of function spaces makes
E-HDG computationally far more advantageous, due to the much smaller number of
degrees of freedom, compared to the HDG counterpart. The benefit is even more
significant for three-dimensional/high-order/fine mesh scenarios. On simplicial
meshes, the proposed methods with a specific choice of approximation spaces are
well-posed for linear(ized) MHD equations. For nonlinear MHD problems, we
present a simple approach exploiting the proposed linear discretizations by
using a Picard iteration. The beauty of this approach is that the
divergence-free and $H(div)$-conforming properties of the velocity and magnetic
fields are automatically carried over for nonlinear MHD equations. We study the
accuracy and convergence of our E-HDG method for both linear and nonlinear MHD
cases through various numerical experiments, including two- and
three-dimensional problems with smooth and singular solutions. The numerical
examples show that the proposed methods are pressure robust, and the divergence
of the resulting velocity and magnetic fields is machine zero for both smooth
and singular problems.
