Analysis of an HDG Method for Linearized Incompressible Resistive MHD Equations

TL;DR

This paper introduces a hybridized discontinuous Galerkin method for stationary linearized incompressible MHD equations, providing error estimates and numerical verification of optimal convergence in multiple dimensions.

Contribution

It develops a novel HDG flux formulation for MHD equations enabling hybridization and derives rigorous error estimates with numerical validation.

Findings

Optimal convergence for fluid velocity and magnetic variables

Quasi-optimal convergence for other quantities

Numerical results confirm theoretical error estimates

Abstract

We present a hybridized discontinuous Galerkin (HDG) method for stationary linearized incompressible magnetohydrodynamics (MHD) equations. At the heart of the paper is the introduction of an HDG flux of the dual saddle-point form of the MHD equations that facilitates the hybridization of discontinuous Galerkin (DG) method. We carry out the $\textit{a priori}$ error estimates for the proposed HDG method on simplicial meshes in both two- and three-dimensions. The analysis provides optimal convergence for the fluid velocity and the magnetic variables, and quasi-optimal convergence for the remaining quantities. Numerical examples are presented to verify the theoretical findings.

Equations

$${\left(\nabla \mathbf{u}\right)}_{ij}=\frac{\partial u_i}{\partial x_j},{\left(\nabla \cdot \mathbf{L}\right)}_i=\nabla \cdot \mathbf{L}\left(i,:\right)=\sum _{j=1}^3\frac{\partial {\mathbf{L}}_{ij}}{\partial x_j},{\left(\mathbf{a}\otimes \mathbf{b}\right)}_{ij}=a_i{b}_j={\left(\mathbf{a}{\mathbf{b}}^T\right)}_{ij}.$$

$${\mathcal{P}}_k\left({\Omega }_h\right)=\left\{u\in L^2(\Omega ):u{|}_K\in {\mathcal{P}}_k\left(K\right)\forall K\in {\Omega }_h\right\}.$$

$$-\frac{1}{\text{Re}}\Delta \mathbf{u}+\nabla p+(\mathbf{w}\cdot \nabla )\mathbf{u}+\kappa \mathbf{d}\times (\nabla \times \mathbf{b})$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\kappa }{\text{Rm}}\nabla \times (\nabla \times \mathbf{b})+\nabla r-\kappa \nabla \times (\mathbf{u}\times \mathbf{d})$$

$$\nabla \cdot \mathbf{b}$$

$$\text{Re}\mathbf{L}-\nabla \mathbf{u}$$

$$-\nabla \cdot \mathbf{L}+\nabla p+(\mathbf{w}\cdot \nabla )\mathbf{u}+\kappa \mathbf{d}\times (\nabla \times \mathbf{b})$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\text{Rm}}{\kappa }\mathbf{J}-\nabla \times \mathbf{b}$$

$$\nabla \times \mathbf{J}+\nabla r-\kappa \nabla \times (\mathbf{u}\times \mathbf{d})$$

$$\nabla \cdot \mathbf{b}$$

$$\mathbf{u}$$

$${\mathbf{b}}^t:=-\mathbf{n}\times (\mathbf{n}\times \mathbf{b})$$

$$r$$

$${⟨{\mathbf{u}}_D\cdot \mathbf{n},1⟩}_{\partial \Omega }=0.$$

$$(p,1{)}_{\Omega }=0.$$

$$\left[\begin{array}[]{c}{\hat{\boldsymbol{F}}}^{1}\cdot\boldsymbol{n}\\ {\hat{\boldsymbol{F}}}^{2}\cdot\boldsymbol{n}\\ {\hat{\boldsymbol{F}}}^{3}\cdot\boldsymbol{n}\\ {\hat{\boldsymbol{F}}}^{4}\cdot\boldsymbol{n}\\ {\hat{\boldsymbol{F}}}^{5}\cdot\boldsymbol{n}\\ {\hat{\boldsymbol{F}}}^{6}\cdot\boldsymbol{n}\end{array}\right]=\left[\begin{array}[]{c}-{\hat{\boldsymbol{u}}}\otimes\boldsymbol{n}\\ -\boldsymbol{L}\boldsymbol{n}+m\boldsymbol{u}+{p}\boldsymbol{n}+\frac{1}{2}\kappa\boldsymbol{{d}}\times\left(\boldsymbol{n}\times\left(\boldsymbol{b}^{t}+\hat{\boldsymbol{b}}^{t}\right)\right)+\alpha_{1}\left(\boldsymbol{u}-{\hat{\boldsymbol{u}}}\right)\\ {\hat{\boldsymbol{u}}}\cdot\boldsymbol{n}\\ -\boldsymbol{n}\times\hat{\boldsymbol{b}}^{t}\\ \boldsymbol{n}\times\boldsymbol{J}+\hat{r}\boldsymbol{n}-\frac{1}{2}\kappa\boldsymbol{n}\times\left(\left(\boldsymbol{u}+{\hat{\boldsymbol{u}}}\right)\times\boldsymbol{{d}}\right)+\alpha_{2}\left(\boldsymbol{b}^{t}-\hat{\boldsymbol{b}}^{t}\right)\\ \boldsymbol{b}\cdot\boldsymbol{n}+\alpha_{3}\left(r-\hat{r}\right)\end{array}\right]$$

$${\mathbf{G}}_h:={\left[{\mathcal{P}}_k({\Omega }_h)\right]}^{d\times d},{\mathbf{V}}_h:={\left[{\mathcal{P}}_k({\Omega }_h)\right]}^d,{\mathbf{Q}}_h:={\mathcal{P}}_k({\Omega }_h),$$

$${\mathbf{H}}_h:={\left[{\mathcal{P}}_k({\Omega }_h)\right]}^{\tilde{d}},{\mathbf{C}}_h:={\left[{\mathcal{P}}_k({\Omega }_h)\right]}^d,{\mathbf{S}}_h:={\mathcal{P}}_k({\Omega }_h),{\mathbf{M}}_h:={\left[{\mathcal{P}}_k({\mathcal{E}}_h)\right]}^d,$$

$${\mathbf{\Lambda }}_h^t:=\left\{\mathbf{\lambda }\in {\left[{\mathcal{P}}_k({\mathcal{E}}_h)\right]}^d:\mathbf{\lambda }\cdot {\mathbf{n}}_e=0\forall e\in {\mathcal{E}}_h\right\},{\mathbf{\Gamma }}_h:={\left[{\mathcal{P}}_k({\mathcal{E}}_h)\right]}^d,$$

$${\left(\mathbf{u},\mathbf{d}\times \left(\nabla \times \mathbf{b}\right)\right)}_K$$

$$\left[\mathbf{d}\times \left(\mathbf{n}\times \mathbf{b}\right)\right]\cdot \mathbf{u}$$

$$\text{Re}{\left({\mathbf{L}}_h,\mathbf{G}\right)}_K+{\left({\mathbf{u}}_h,\nabla \cdot \mathbf{G}\right)}_K+{⟨{\hat{\mathbf{F}}}_h^1\cdot \mathbf{n},\mathbf{G}⟩}_{\partial K}$$

$${\left({\mathbf{L}}_h,\nabla \mathbf{v}\right)}_K-{\left(p_h,\nabla \cdot \mathbf{v}\right)}_K-{\left({\mathbf{u}}_h\otimes \mathbf{w},\nabla \mathbf{v}\right)}_K$$

$$+\kappa {\left({\mathbf{b}}_h,\nabla \times \left(\mathbf{v}\times \mathbf{d}\right)\right)}_K+{⟨{\hat{\mathbf{F}}}_h^2\cdot \mathbf{n},\mathbf{v}⟩}_{\partial K}$$

$$-{\left({\mathbf{u}}_h,\nabla q\right)}_K+{⟨{\hat{\mathbf{F}}}_h^3\cdot \mathbf{n},q⟩}_{\partial K}$$

$$\frac{\text{Rm}}{\kappa }{\left({\mathbf{J}}_h,\mathbf{H}\right)}_K-{\left({\mathbf{b}}_h,\nabla \times \mathbf{H}\right)}_K+{⟨{\hat{\mathbf{F}}}_h^4\cdot \mathbf{n},\mathbf{H}⟩}_{\partial K}$$

$${\left({\mathbf{J}}_h,\nabla \times \mathbf{c}\right)}_K-{\left(r_h,\nabla \cdot \mathbf{c}\right)}_K-\kappa {\left({\mathbf{u}}_h,\mathbf{d}\times \left(\nabla \times \mathbf{c}\right)\right)}_K$$

$$+{⟨{\hat{\mathbf{F}}}_h^5\cdot \mathbf{n},\mathbf{c}⟩}_{\partial K}$$

$$-{\left({\mathbf{b}}_h,\nabla s\right)}_K+{⟨{\hat{\mathbf{F}}}_h^6\cdot \mathbf{n},s⟩}_{\partial K}$$

$${⟨[[{\hat{\mathbf{F}}}_h^2\cdot \mathbf{n}]],\mathbf{\mu }⟩}_e=0,{⟨[[{\hat{\mathbf{F}}}_h^5\cdot \mathbf{n}]],{\mathbf{\lambda }}^t⟩}_e=0,{⟨[[{\hat{\mathbf{F}}}_h^6\cdot \mathbf{n}]],\gamma ⟩}_e=0$$

$${⟨{\hat{\mathbf{u}}}_h,\mathbf{\mu }⟩}_e={⟨{\mathbf{u}}_D,\mathbf{\mu }⟩}_e,{⟨{\hat{\mathbf{b}}}_h^t,{\mathbf{\lambda }}^t⟩}_e={⟨{\mathbf{h}}_D,{\mathbf{\lambda }}^t⟩}_e,{⟨{\hat{r}}_h,\gamma ⟩}_e=0$$

$$(p_h,1)$$

$$\text{Re}{\left({\mathbf{L}}_h,{\mathbf{L}}_h\right)}_K+{\left({\mathbf{u}}_h,\nabla \cdot {\mathbf{L}}_h\right)}_K+{⟨-{\hat{\mathbf{u}}}_h\otimes \mathbf{n},{\mathbf{L}}_h⟩}_{\partial K}$$

$${\left({\mathbf{L}}_h,\nabla {\mathbf{u}}_h\right)}_K-{\left(p_h,\nabla \cdot {\mathbf{u}}_h\right)}_K-{\left({\mathbf{u}}_h\otimes \mathbf{w},\nabla {\mathbf{u}}_h\right)}_K+\kappa {\left({\mathbf{b}}_h,\nabla \times \left({\mathbf{u}}_h\times \mathbf{d}\right)\right)}_K$$

$$+{⟨-{\mathbf{L}}_h\mathbf{n}+m{\mathbf{u}}_h+p_h\mathbf{n}+\frac{1}{2}\kappa \mathbf{d}\times \left(\mathbf{n}\times \left({\mathbf{b}}_h^t+{\hat{\mathbf{b}}}_h^t\right)\right)+{\alpha }_1\left({\mathbf{u}}_h-{\hat{\mathbf{u}}}_h\right),{\mathbf{u}}_h⟩}_{\partial K}$$

$$-{\left({\mathbf{u}}_h,\nabla q\right)}_K+{⟨{\hat{\mathbf{u}}}_h\cdot \mathbf{n},q⟩}_{\partial K}$$

$$\frac{\text{Rm}}{\kappa }{\left({\mathbf{J}}_h,{\mathbf{J}}_h\right)}_K-{\left({\mathbf{b}}_h,\nabla \times {\mathbf{J}}_h\right)}_K+{⟨-\mathbf{n}\times {\hat{\mathbf{b}}}_h^t,{\mathbf{J}}_h⟩}_{\partial K}$$

$${\left({\mathbf{J}}_h,\nabla \times {\mathbf{b}}_h\right)}_K-{\left(r_h,\nabla \cdot {\mathbf{b}}_h\right)}_K-\kappa {\left({\mathbf{u}}_h,\mathbf{d}\times \left(\nabla \times {\mathbf{b}}_h\right)\right)}_K$$

$$+{⟨\mathbf{n}\times {\mathbf{J}}_h+{\hat{r}}_h\mathbf{n}-\frac{1}{2}\kappa \mathbf{n}\times \left(\left({\mathbf{u}}_h+{\hat{\mathbf{u}}}_h\right)\times \mathbf{d}\right)+{\alpha }_2\left({\mathbf{b}}_h^t-{\hat{\mathbf{b}}}_h^t\right),{\mathbf{b}}_h⟩}_{\partial K}$$

$$-{\left({\mathbf{b}}_h,\nabla r_h\right)}_K+{⟨{\mathbf{b}}_h\cdot \mathbf{n}+{\alpha }_3\left(r_h-{\hat{r}}_h\right),r_h⟩}_{\partial K}$$

$$\text{Re}{\Vert {\mathbf{L}}_h\Vert }_0^2+\frac{\text{Rm}}{\kappa }{\Vert {\mathbf{J}}_h\Vert }_0^2-⟨{\hat{\mathbf{u}}}_h\otimes \mathbf{n},{\mathbf{L}}_h⟩+⟨\frac{m}{2}{\mathbf{u}}_h,{\mathbf{u}}_h⟩+⟨{\alpha }_1({\mathbf{u}}_h-{\hat{\mathbf{u}}}_h),{\mathbf{u}}_h⟩$$

$$+⟨\frac{1}{2}\kappa \mathbf{d}\times \left(\mathbf{n}\times {\hat{\mathbf{b}}}_h^t\right),{\mathbf{u}}_h⟩+⟨{\hat{\mathbf{u}}}_h\cdot \mathbf{n},p_h⟩-⟨\mathbf{n}\times {\hat{\mathbf{b}}}_h^t,{\mathbf{J}}_h⟩+⟨{\hat{r}}_h\mathbf{n},{\mathbf{b}}_h⟩$$

Full text

{AMS}

65N30, 65N12, 65N15, 76W05

1 Introduction

An important base-level representation for continuum approximation of the dynamics of charged fluids in the presence of electromagnetic fields is the resistive magnetohydrodynamics (MHD) model. MHD models describe important physical phenomena in astrophysical systems (e.g. solar flares, and planetary magnetic field generation) and in critical science and technological applications (e.g., magnetically confined fusion energy devices) [28], for example. The single fluid resistive MHD model involes the partial differential equations (PDEs) describing conservation of mass, momentum, and energy, coupled with the low-frequency Maxwell’s equations. This multiphysics PDE system is highly nonlinear and characterized by multiple interacting physical phenomena spanning a wide range of length- and time-scales. These characteristics make the task of developing scalable, robust, accurate, and efficient computational methods extremely challenging.

The most common computational solution strategies for MHD have been the use of explicit and partially implicit time integration methods. Notably are implicit-explicit [2, 43, 35], semi-implicit [47, 30, 50], and operator-splitting [33, 45] techniques that include some use of implicit operators in the formulation. The implicitness of these approaches is used to enhance efficiency by removing stringent explicit time-scale constraints in the problem, either from diffusion or from fast-wave phenomena [9, 36].

In addition to the challenges associated with designing robust and efficient time integrators, there are a number of spatial discretization issues including the dual saddle-point structure of the velocity-pressure $\left(\boldsymbol{u},p\right)$ and the enforcement of the solenoidal involution/constraint on the magnetic induction $(\nabla\cdot\boldsymbol{b}=0)$. This adds considerable complexity to the numerical approximation of resistive MHD system. In the context of finite volume and finite element methods, there are four popular approaches to deal with these difficulties: 1) physics-compatible discretizations that directly enforce key mathematical properties of the continuous problem (see e.g. [39, 34, 3]); 2) methods that transform to potential-based formulations to eliminate one or both saddle-point sub-systems [10, 43, 37, 48]; 3) exact and weighted-exact penalty formulations [29, 25, 19, 20]; and 4) and stabilization methods that regularize the dual saddle-point structure [46, 18, 49].

In this paper we propose a hybridized discontinuous Galerkin (HDG) formulation for a linearized version of the resistive MHD system. The hybridization technique and post-processing have been proposed to reduce computational costs of saddle-point problems and to improve the accuracy of numerical solutions [1]. HDG methods were developed by Cockburn, coauthors, and others to mitigate the computational costs of classical discontinuous Galerkin (DG) methods. They have been proposed for various types of PDEs including, but not limited to, Poisson-type equations [13, 15, 40, 22], the Stokes equation [12, 41], the Oseen equations [8], and the incompressible Navier-Stokes equations [42].

In HDG discretizations, the coupled unknowns are single-valued traces introduced on the mesh skeleton, i.e., the faces, and for high order implicit systems the resulting matrix is substantially smaller and sparser compared to standard DG approaches. Once they are solved for, the volume DG unknowns can be recovered in an element-by-element fashion, completely independent of one another. Therefore HDG methods have an intrinsic structure for parallel computing which is essential for large scale applications. Nevertheless, devising an HDG method for coupled PDE systems is challenging because construction of a consistent and robust HDG flux is nontrivial. We adopt the upwind HDG framework proposed in [5, 7, 6] since it provides a systematic construction of HDG methods for a large class of PDEs.

Our work starts with section 2 where notations and conventions are introduced to enable the construction of HDG method in section 3. Specifically, the proposed HDG method is introduced directly on the dual saddle-point structure of the MHD system and its well-posedness is analyzed using an energy approach. This is followed by the a priori error estimation in section 4 where we combine an energy analysis, specially designed projections, and a duality argument to provide convergence rates for all variables. Our development can serve as a standalone high-order solver for linearized MHD equations, or can be used as the fast-time scale solver in an implicit/explicit time integration method, and/or the solver for a sub-step in a fixed-point nonlinear solver. Various numerical results will be presented in section 5 to verify our theoretical findings. Section 6 concludes the paper with future work. This is followed by four appendices in which we detail the definition and analysis of projection operators, state some auxiliary results, discuss the well-posedness of the adjoint equation, and present a postprocessing procedure to enforce the solenoidal constraints.

2 Notation

In this section we introduce common notations and conventions to be used in the rest of the paper. Let $\Omega\subset{\mathbb{R}}^{d}$, $d=2,3$, be a bounded domain such that it is simply-connected and its boundary ${\partial\Omega}$ is a Lipschitz manifold with only one component. Suppose that we have a triangulation of $\Omega$, i.e., a partition of $\Omega$ into a finite number of nonoverlapping $d$-dimensional simplices. We assume that the triangulation is shape-regular, i.e., for all $d$-dimensional simplices in the triangulation, the ratio of the diameter of the simplex and the radius of an inscribed $d$-dimensional ball is uniformly bounded. We will use ${\Omega_{h}}$ and ${\mathcal{E}_{h}}$ to denote the sets of $d$- and ($d-1$)-dimensional simplices of the triangulation, and call ${\mathcal{E}_{h}}$ the mesh skeleton of the triangulation. The boundary and interior mesh skeletons are defined by ${\mathcal{E}_{h}^{\partial}}:=\{e\in{\mathcal{E}_{h}}\,:\,e\subset\partial\Omega\}$ and ${\mathcal{E}_{h}^{o}}:={\mathcal{E}_{h}}\setminus{\mathcal{E}_{h}^{\partial}}$. We also define ${{\partial\Omega}_{h}}:=\left\{{\partial K}:K\in{\Omega_{h}}\right\}$. The mesh size of triangulations is $h:=\max_{K\in{\Omega_{h}}}\text{diam}(K)$.

We use $\left(\cdot,\cdot\right)_{D}$ (respectively $\left<\cdot,\cdot\right>_{D}$) to denote the $L^{2}$-inner product on $D$ if $D$ is a $d$- (respectively $(d-1)$-) dimensional domain. The standard notation $W^{s,p}(D)$, $s\geq 0$, $1\leq p\leq\infty$, is used for the Sobolev space on $D$ based on $L^{p}$-norm with differentiability $s$ (see, e.g., [23]) and $\left\|\cdot\right\|_{W^{s,p}(D)}$ denotes the associated norm. In particular, if $p=2$, we use $H^{s}(D):=W^{s,p}(D)$ and $\left\|\cdot\right\|_{s,D}$. $W^{s,p}({\Omega_{h}})$ denotes the space of functions whose restrictions on $K$ reside in $W^{s,p}(K)$ for each $K\in{\Omega_{h}}$ and its norm is $\left\|u\right\|_{W^{s,p}({\Omega_{h}})}^{p}:=\sum_{K\in{\Omega_{h}}}\left\|u|_{K}\right\|_{W^{s,p}(K)}^{p}$ if $1\leq p<\infty$ and $\left\|u\right\|_{W^{s,\infty}({\Omega_{h}})}:=\max_{K\in{\Omega_{h}}}\left\|u|_{K}\right\|_{W^{s,\infty}(K)}$. For simplicity, we use $\left(\cdot,\cdot\right)$, $\left<\cdot,\cdot\right>$, $\left\|\cdot\right\|_{s}$, $\left\|\cdot\right\|_{{{\partial\Omega}_{h}}}$, and $\left\|\cdot\right\|_{W^{s,\infty}}$ for $\left(\cdot,\cdot\right)_{\Omega}$, $\left<\cdot,\cdot\right>_{{{\partial\Omega}_{h}}}$, $\left\|\cdot\right\|_{s,\Omega}$, $\left\|\cdot\right\|_{0,{{\partial\Omega}_{h}}}$, and $\left\|\cdot\right\|_{W^{s,\infty}({\Omega_{h}})}$, respectively. We define $\|u,v\|:=\|u\|+\|v\|$. Furthermore, we denote by $A\lesssim B$ the inequality $A\leq\lambda B$ with a constant $\lambda>0$ independent of the mesh size, and by $A\sim B$ the combination of $A\lesssim B$ and $B\lesssim A$.

For vector- or matrix-valued functions these notations are naturally extended with a component-wise inner product. We define similar spaces (respectively inner products and norms) on a single element and a single skeleton face/edge by replacing ${\Omega_{h}}$ with $K$ and ${\mathcal{E}_{h}}$ with $e$. We define the gradient of a vector, the divergence of a matrix, and the outer product symbol $\otimes$ as:

[TABLE]

In this paper $\boldsymbol{n}$ denotes a unit outward normal vector field on faces/edges. If $\partial K^{-}\cap\partial K^{+}\in{\mathcal{E}_{h}}$ for two distinct simplices $K^{-},K^{+}$, then $\boldsymbol{n}^{-}$ and $\boldsymbol{n}^{+}$ denote the outward unit normal vector fields on $\partial K^{-}$ and $\partial K^{+}$, respectively, and $\boldsymbol{n}^{-}=-\boldsymbol{n}^{+}$ on $\partial K^{-}\cap\partial K^{+}$. We simply use $\boldsymbol{n}$ to denote either $\boldsymbol{n}^{-}$ or $\boldsymbol{n}^{+}$ in an expression that is valid for both cases, and this convention is also used for other quantities (restricted) on a face/edge $e\in{\mathcal{E}_{h}}$. For a scalar quantity $u$ which is double-valued on $e:=\partial K^{-}\cap\partial K^{+}$, the jump term on $e$ is defined by $[\![{u\boldsymbol{n}}]\!]|_{e}=u^{+}\boldsymbol{n}^{+}+u^{-}\boldsymbol{n}^{-}$ where $u^{+}$ and $u^{-}$ are the traces of $u$ from $K^{+}$- and $K^{-}$-sides, respectively. For double-valued vector quantity $\boldsymbol{u}$ and matrix quantity $\boldsymbol{L}$, jump terms are $[\![{\boldsymbol{u}\cdot\boldsymbol{n}}]\!]|_{e}=\boldsymbol{u}^{+}\cdot\boldsymbol{n}^{+}+\boldsymbol{u}^{-}\cdot\boldsymbol{n}^{-}$ and $[\![{\boldsymbol{L}\boldsymbol{n}}]\!]|_{e}=\boldsymbol{L}^{+}\boldsymbol{n}^{+}+\boldsymbol{L}^{-}\boldsymbol{n}^{-}$ where $\boldsymbol{L}\boldsymbol{n}$ denotes the matrix-vector product.

We define ${\mathcal{P}}_{k}\left(K\right)$ as the space of polynomials of degree at most $k$ on a domain $K$, and

[TABLE]

The space of polynomials on the mesh skeleton ${\mathcal{P}}_{k}\left({\mathcal{E}_{h}}\right)$ is similarly defined, and their extensions to vector- or matrix-valued polynomials $\left[{\mathcal{P}}_{k}({\Omega_{h}})\right]^{d}$, $\left[{\mathcal{P}}_{k}({\Omega_{h}})\right]^{d\times d}$, $\left[{\mathcal{P}}_{k}({\mathcal{E}_{h}})\right]^{d}$, etc, are straightforward.

3 HDG Formulation

We consider the following nondimensional linearized incompressible MHD system [32]

[TABLE]

where $\boldsymbol{u}$ is velocity of the fluid (plasma or liquid metal), $\boldsymbol{b}$ the magnetic field, $p$ the fluid pressure, and $r$ a scalar potential. The following are constant parameters: a fluid Reynolds number $\textup{Re}>0$, a magnetic Reynolds number $\textup{Rm}>0$, and a coupling parameter $\kappa=\textup{Ha}^{2}/(\textup{Re}\textup{Rm})$, with the Hartmann number $\textup{Ha}>0$. Here, $\boldsymbol{{d}}\in\left[W^{1,\infty}\left(\Omega\right)\right]^{d}$ is a prescribed magnetic field and $\boldsymbol{w}\in\left[W^{1,\infty}\left({\Omega_{h}}\right)\right]^{d}\cap H(div,\Omega)$ with $\nabla\cdot\boldsymbol{w}=0$ is a prescribed velocity.

By introducing auxiliary variables $\boldsymbol{L}$ and $\boldsymbol{J}$, we cast (1) into a first order hyperbolic system:

[TABLE]

In this paper, we consider the following (Dirichlet) boundary conditions for the MHD system (2)

[TABLE]

where we have defined the tangent component of a vector $\boldsymbol{a}$ as $\boldsymbol{a}^{t}:=-\boldsymbol{n}\times(\boldsymbol{n}\times\boldsymbol{a})$. Additionally, we require the compatibility condition for $\boldsymbol{u}_{D}$:

[TABLE]

For the uniqueness of the pressure, $p$, we require that the pressure has zero mean, i.e.,

[TABLE]

Following the upwind HDG framework in [5] we define the HDG flux as

[TABLE]

where $\hat{\boldsymbol{b}}^{t}$, ${\hat{\boldsymbol{u}}}$, and $\hat{r}$ are the single-valued trace quantities residing on the mesh skeleton ${\mathcal{E}_{h}}$. They will be new unknowns in the discretizations, which will be described later, to hybridize the DG method. Here, $m:=\boldsymbol{w}\cdot\boldsymbol{n}$, and $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{3}$ are constant parameters. As will be shown later, the conditions $\alpha_{1}>\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}}$, $\alpha_{2}>0$, and $\alpha_{3}>0$ are sufficient for the well-posedness of our HDG formulation. Note that all 6 components of the HDG flux, ${\hat{\boldsymbol{F}}}$, for simplicity are denoted in the same fashion (by a bold italic symbol). It is, however, clear from (2) that ${\hat{\boldsymbol{F}}}^{1}$ is a third order tensor, ${\hat{\boldsymbol{F}}}^{2}$ is a second order tensor, ${\hat{\boldsymbol{F}}}^{3}$ is a vector, etc, and that the normal HDG flux components, ${\hat{\boldsymbol{F}}}^{i}\cdot\boldsymbol{n}$, defined in (6), are tensors of one order lower.

For discretization we introduce the discontinuous piecewise polynomial spaces

[TABLE]

where $\tilde{d}=3$ if $d=3$, and $\tilde{d}=1$ if $d=2$.

Let us introduce two identities which are useful throughout the paper:

[TABLE]

These identities follow from integration by parts and vector product identities.

Next, we multiply (2a) through (2f) by test functions $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)$, integrate by parts all terms, and introduce the HDG flux (6) in the boundary terms. This results in a local discrete weak formulation, which we shall call the local solver of the HDG method, for the MHD system (2):

[TABLE]

for all $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)\in\mathbf{G}_{h}\left(K\right)\times\mathbf{V}_{h}\left(K\right)\times\mathbf{Q}_{h}\left(K\right)\times\mathbf{H}_{h}\left(K\right)\times\mathbf{C}_{h}\left(K\right)\times\mathbf{S}_{h}\left(K\right)$ and for all $K\in{\Omega_{h}}$, where $\boldsymbol{u}_{h}$, $\boldsymbol{L}_{h}$, …, are the discrete counterparts of $\boldsymbol{u}$, $\boldsymbol{L}$, …, and ${\hat{\boldsymbol{F}}}_{h}^{i}$ is the discrete counterpart of ${\hat{\boldsymbol{F}}}^{i}$ in (6) by replacing the unknowns $\boldsymbol{u}$, $\boldsymbol{L}$, …, with their discrete counterparts.

Since $\hat{\boldsymbol{b}}^{t}_{h}$, ${\hat{\boldsymbol{u}}}_{h}$, and $\hat{r}_{h}$ are the new unknowns, we need to equip extra equations to close the system (8). To that end, we observe that an element $K$ communicates with its neighbors only through the trace unknowns. For the HDG method to be conservative, we (weakly) enforce the continuity of the HDG flux (6) across each interior edge, i.e., $\left<[\![{{\hat{\boldsymbol{F}}}_{h}\cdot\boldsymbol{n}}]\!],\boldsymbol{\delta}\right>_{e}=0,\forall e\in{\mathcal{E}_{h}^{o}}$. Since ${\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h}$ and $\hat{r}_{h}$ are single-valued on ${\mathcal{E}_{h}}$, $[\![{{\hat{\boldsymbol{F}}}_{h}^{1}\cdot\boldsymbol{n}}]\!]=\boldsymbol{0}$, $[\![{{\hat{\boldsymbol{F}}}_{h}^{3}\cdot\boldsymbol{n}}]\!]=0$, and $[\![{{\hat{\boldsymbol{F}}}_{h}^{4}\cdot\boldsymbol{n}}]\!]=\boldsymbol{0}$. The conservation constraints to be enforced reduce to

[TABLE]

for all $(\boldsymbol{\mu},\boldsymbol{\lambda}^{t},\gamma)\in\mathbf{M}_{h}\left(e\right)\times\mathbf{\Lambda}^{t}_{h}\left(e\right)\times\mathbf{\Gamma}_{h}\left(e\right)$, and for all $e$ in ${\mathcal{E}_{h}^{o}}$. Finally, we enforce the Dirichlet boundary conditions through the trace unknowns:

[TABLE]

for all $(\boldsymbol{\mu},\boldsymbol{\lambda}^{t},\gamma)\in\boldsymbol{M}_{h}\left(e\right)\times\boldsymbol{\Lambda}^{t}_{h}\left(e\right)\times\Gamma_{h}\left(e\right)$ for all $e$ in ${\mathcal{E}_{h}^{\partial}}$.

In (8), (9), and (10), we seek $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})\in\mathbf{G}_{h}\times\mathbf{V}_{h}\times\mathbf{Q}_{h}\times\mathbf{H}_{h}\times\mathbf{C}_{h}\times\mathbf{S}_{h}$ and $({\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h})\in\mathbf{M}_{h}\times\mathbf{\Lambda}^{t}_{h}\times\mathbf{\Gamma}_{h}$. From this point forward, for simplicity in writing, we will not state explicitly that equations hold for all test functions, for all elements, or for all edges.

We will refer to $\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},$ and $r_{h}$ as the local variables, and to equation (8) on each element as the local solver. This reflects the fact that we can solve for local variables element-by-element as function of ${\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},$ and $\hat{r}_{h}$. On the other hand, we will refer to ${\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},$ and $\hat{r}_{h}$ as the global variables, which are governed by equations (9) and (10) on the mesh skeleton. Finally, for the uniqueness of the discrete pressure $p_{h}$, we enforce the discrete counterpart of (5):

[TABLE]

3.1 Well-posedness of the HDG formulation

In this subsection we discuss well-posedness of the system (8)–(10).

Theorem 3.1**.**

*The HDG system (8)–(11) is well-posed. *

Proof 3.2**.**

Since the problem is a system of linear equations with the same number of equations and unknowns, without loss of generality, we only need show that $\left(\boldsymbol{g},\boldsymbol{f},\boldsymbol{u}_{D},\boldsymbol{h}_{D}\right)=\boldsymbol{0}$ implies $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h},{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h})=\boldsymbol{0}$. To begin, we take $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)=(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$ and get

[TABLE]

If we integrate by parts the first four terms of the second equation and the first term of the fifth equation, sum the resulting equations, and sum over all elements, then we arrive at

[TABLE]

Here, we have used

[TABLE]

which is obtained from $\nabla\cdot\boldsymbol{w}=0$.

Next, we set $({\boldsymbol{\mu},\boldsymbol{\lambda}^{t},\gamma})=({{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h}})$, and sum (9) over all interior edges to obtain

[TABLE]

where we have used the continuity of $\boldsymbol{{d}}$ to eliminate $\langle{\boldsymbol{{d}}\times({\boldsymbol{n}\times\hat{\boldsymbol{b}}^{t}_{h}}),{\hat{\boldsymbol{u}}}_{h}}\rangle_{{{\partial\Omega}_{h}}\setminus{\partial\Omega}}$ and $\langle{\boldsymbol{n}\times\left({\hat{\boldsymbol{u}}}_{h}\times\boldsymbol{{d}}\right),\hat{\boldsymbol{b}}^{t}_{h}}\rangle_{{{\partial\Omega}_{h}}\setminus{\partial\Omega}}$.

Since $\boldsymbol{u}_{D}=\boldsymbol{0}$ and $\boldsymbol{h}_{D}=\boldsymbol{0}$ by assumption, we conclude from the boundary conditions (10) that ${\hat{\boldsymbol{u}}}_{h}=\boldsymbol{0}$, $\hat{\boldsymbol{b}}^{t}_{h}=\boldsymbol{0}$, and $\hat{r}_{h}=0$ on ${\partial\Omega}$. The integrals in (14) can then be written over ${{\partial\Omega}_{h}}$ since the contribution on the domain boundary, ${\partial\Omega}$, is zero. Subtracting (14) from (13) we arrive at

[TABLE]

Finally, using the fact that $\boldsymbol{w}\in H(div,\Omega)$ and ${\hat{\boldsymbol{u}}}_{h}=\boldsymbol{0}$ on ${\partial\Omega}$, we can freely add $0=\left<\frac{m}{2}{\hat{\boldsymbol{u}}}_{h},{\hat{\boldsymbol{u}}}_{h}\right>$ to rewrite (15) as

[TABLE]

Choosing $\alpha_{1}>\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}}$, $\alpha_{2}>0$, and $\alpha_{3}>0$, we can conclude that $\boldsymbol{L}_{h}=\boldsymbol{0}$ and $\boldsymbol{J}_{h}=\boldsymbol{0}$, that $\boldsymbol{u}_{h}={\hat{\boldsymbol{u}}}_{h}$, $\boldsymbol{b}_{h}^{t}=\hat{\boldsymbol{b}}^{t}_{h}$, and $r_{h}=\hat{r}_{h}$ on ${\mathcal{E}_{h}^{o}}$, and that $\boldsymbol{u}_{h}=\boldsymbol{0}$, $\boldsymbol{b}_{h}^{t}=\boldsymbol{0}$, and $r_{h}=0$ on ${\partial\Omega}$.

Now, we integrate (8a) by parts to obtain $\nabla\boldsymbol{u}_{h}=\boldsymbol{0}$ in $K$, which implies that $\boldsymbol{u}_{h}$ must be elementwise constant. The fact that $\boldsymbol{u}_{h}={\hat{\boldsymbol{u}}}_{h}$ on ${\mathcal{E}_{h}^{o}}$ means $\boldsymbol{u}_{h}$ is also continuous on ${\mathcal{E}_{h}}$, and since $\boldsymbol{u}_{h}=\boldsymbol{0}$ on ${\partial\Omega}$ we conclude that $\boldsymbol{u}_{h}=\boldsymbol{0}$, and therefore ${\hat{\boldsymbol{u}}}_{h}=\boldsymbol{0}$.

Note that since $\boldsymbol{b}_{h}^{t}=\hat{\boldsymbol{b}}^{t}_{h}$ on ${\mathcal{E}_{h}^{o}}$, $\boldsymbol{b}_{h}^{t}$ is continuous on $\Omega$. Furthermore, the third conservation constraint in (9) implies that $\boldsymbol{b}_{h}\cdot\boldsymbol{n}$ is continuous on $\Omega$. Integrating both (8d) and (8f) by parts, we can conclude that $\nabla\times\boldsymbol{b}_{h}=\boldsymbol{0}$ and $\nabla\cdot\boldsymbol{b}_{h}=0$ on $\Omega$. When $\boldsymbol{b}_{h}\in H(div,\Omega)\cap H(curl,\Omega)$ and $\boldsymbol{b}_{h}^{t}=\boldsymbol{0}$ on ${\partial\Omega}$, and when $\Omega$ is simply-connected with one component to the boundary, there is a constant $C>0$ such that $\left\lVert\boldsymbol{b}_{h}\right\rVert_{0}\leq C(\left\lVert\nabla\cdot\boldsymbol{b}_{h}\right\rVert_{0}+\left\lVert\nabla\times\boldsymbol{b}_{h}\right\rVert_{0})$ (see, e.g., [26, Lemma 3.4]). This implies that $\boldsymbol{b}_{h}=\boldsymbol{0}$, and hence $\hat{\boldsymbol{b}}^{t}_{h}=\boldsymbol{0}$.

*Considering the vanishing unknowns above, integrating by parts reduces (8b) and (8e) to $(\nabla p_{h},\boldsymbol{v})_{K}=0$ and $(\nabla r_{h},\boldsymbol{{c}})_{K}=0$, respectively. Thus, $p_{h}$ and $r_{h}$ are elementwise constants. Since $r_{h}=\hat{r}_{h}$ on ${\mathcal{E}_{h}^{o}}$, then $r_{h}$ is continuous on $\Omega$, and since $r_{h}=0$ on ${\partial\Omega}$, we can conclude that $r_{h}=0$, and hence $\hat{r}_{h}=0$. Finally, we use the first conservation constraint in (9) to conclude $p_{h}$ is continuous and hence constant on $\Omega$. Using the zero-average condition (11) yields $p_{h}=0$. *

3.2 Well-posedness of the local solver

The key design of the HDG method is that it allows us to separate the computation of the volume (DG) unknowns $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$ and the trace unknowns $({{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h}})$. In practice, we first solve (8) for local unknowns $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$ as a function of $({{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h}})$. These are then substituted into the conservative algebraic equation (9) on the mesh skeleton to solve for the unknown $({{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h}})$. Finally, the local unknowns $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},\allowbreak p_{h},\allowbreak\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$ are computed, as in the first step, using $({{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h}})$ from the second step. It is therefore important to study the well-posedness of the local solver.

Similar to HDG methods for Stokes equation [41, 16, 5], it turns out that the local solver is not well-posed unless extra conditions are imposed on the pressure. Two methods for achieving the well-posedness of the local solver for the Stokes equations are proposed in [41]. One is a pseudotransient approach, and the other involves introducing the element average edge pressure as global unknowns. These methods are both suitable for our setting here. Here, we present a new approach in which we introduce the elementwise pressure integral as a global unknown and require their sum to vanish. Toward this goal, we introduce the space of elementwise constants, $\mathbf{X}_{h}:={\mathcal{P}}_{0}({\Omega_{h}})$. Next we augment (8c) to read

[TABLE]

with $\rho_{h}\in\mathbf{X}_{h}$, $\bar{q}|_{K}:=|K|^{-1}\left(q,1\right)_{K}$ the average of $q$ in $K$, and $|K|$ the volume of element $K$. Next we augment the global solver with

[TABLE]

for all $\xi$ in $\mathbf{X}_{h}$, and remove the constraint (11), which will be automatically satisfied by this construction.

To justify (17) and (18) we make the following observations. First, summing (18) over all elements and using the compatibility condition on $\boldsymbol{u}_{D}$, (4), we conclude

[TABLE]

and

[TABLE]

Next, setting $q=1$ on $K$ in (17) and using (20), we can conclude that

[TABLE]

and therefore that (8c) holds for each $K$. Additionally, (21) and (19) imply that (11) holds. Finally, we note that we have added the same number of new unknowns $\rho_{h}$ as the number of equations in (18).

For this modified HDG scheme we claim well-posedness of the local solver.

Theorem 3.3**.**

*The local solver given by (8) such that (8c) is replaced by (17), is well-posed. In other words, given $({\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h},\boldsymbol{g},\boldsymbol{f},\rho_{h})$, there exists a unique solution

$(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$ to the system. *

Proof 3.4**.**

We show that $({\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h},\boldsymbol{g},\boldsymbol{f},\rho_{h})=\boldsymbol{0}$ implies $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})=\boldsymbol{0}$. To begin, set $({\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h},\boldsymbol{g},\boldsymbol{f},\rho_{h})=\boldsymbol{0}$. Then (17) reduces to $-\left(\boldsymbol{u}_{h},\nabla q\right)_{K}+\left(p_{h},\bar{q}\right)_{K}=0$, and taking $q$ as constant gives $\left(p_{h},\bar{q}\right)_{K}=0$ and hence $-\left(\boldsymbol{u}_{h},\nabla q\right)_{K}=0$.

Next, choose $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)=(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h})$, integrate by parts the first four terms in (8b) and the first term in (8e), and sum the resulting equations in the local solver to conclude

[TABLE]

Recalling we have set $\alpha_{1}>\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}}$, $\alpha_{2}>0$, and $\alpha_{3}>0$, we can conclude that

[TABLE]

*Using an argument similar to that in Section 3.1 we can conclude that $\boldsymbol{u}_{h}=\boldsymbol{b}_{h}=\boldsymbol{0}$ in $K$. From (8b) and (8e) we can conclude $\left(\nabla p_{h},\boldsymbol{v}\right)_{K}=0$ and $\left(\nabla r_{h},\boldsymbol{{c}}\right)_{K}=0$, respectively. Thus, $p_{h}$ and $r_{h}$ must be constant, and since $r_{h}=0$ on ${\partial K}$, $r_{h}$ is identically zero in $K$. Now since $(p_{h},\bar{q})_{K}=0$, we have $p_{h}=0$ in $K$. *

Remark 3.5**.**

*Note that introducing $\rho_{h}$ and equations (17) and (18) does not alter the solution of the original HDG scheme. Indeed, if $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h},\allowbreak{\hat{\boldsymbol{u}}}_{h},\allowbreak\hat{\boldsymbol{b}}^{t}_{h},\allowbreak\hat{r}_{h},\allowbreak\rho_{h})$ is a solution of the modified scheme, it is also a solution of the original one because the modified scheme contains all the equations of the original one, except for (8c) and (11). But we have already shown that (17), (18), and (4) imply (8c), and that (21) and (19) imply (11). Conversely, reversing these arguments, if $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h},\allowbreak{\hat{\boldsymbol{u}}}_{h},\allowbreak\hat{\boldsymbol{b}}^{t}_{h},\allowbreak\hat{r}_{h})$ is a solution of the original scheme, we can define $\rho_{h}$ as in (21) and add (21) to (8c) to recover (17). Then, taking $q$ as constant in (8c) implies (20). Also, (21) and (11) imply (19). Finally, adding (20) and (19) we recover (18), showing that $({\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h},{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h},\rho_{h}})$ is a solution of the modified one. Since the original HDG scheme is well-posed, so is the modified one. *

4 Error analysis

For an unknown $\sigma$ we use ${\varepsilon}_{\sigma}$ to denote the error between the exact solution $\sigma$ and its finite element approximation $\sigma_{h}$. For example, ${\varepsilon}_{\boldsymbol{L}}:=\boldsymbol{L}-\boldsymbol{L}_{h}$ and ${\varepsilon}_{{\hat{\boldsymbol{u}}}}:={\hat{\boldsymbol{u}}}-{\hat{\boldsymbol{u}}}_{h}$, where ${\hat{\boldsymbol{u}}}$ is the trace of the exact solution $\boldsymbol{u}$ on the mesh skeleton. We use $\Pi\sigma$ to denote some interpolation (or projection) of the unknown $\sigma$ into its associated finite element space, and decompose ${\varepsilon}_{\sigma}$ into ${\varepsilon}_{\sigma}^{I}+{\varepsilon}_{\sigma}^{h}$ where

[TABLE]

Here the superscript $I$ of ${\varepsilon}^{I}$ denotes the ‘I’nterpolation (in fact projection) error, and the superscript $h$ of $e^{h}$ indicates the difference between the interpolation of the exact solution and the finite element approximation. We will see that $\Pi\sigma$ may not depend only on $\sigma$. In particular, we define a collective projection $\boldsymbol{\Pi}(\boldsymbol{L},\boldsymbol{u},{p},\boldsymbol{J},\boldsymbol{b},r,{\hat{\boldsymbol{u}}},\hat{\boldsymbol{b}}^{t},\hat{r})$ in Appendix B for our HDG scheme. Each component of $\boldsymbol{\Pi}$ may depend on the others. Nonetheless, for simplicity of presentation we use $\Pi\boldsymbol{L}$ to denote the $\boldsymbol{L}$-component of $\boldsymbol{\Pi}$ for example. The analysis and the properties of the proposed projection can be referred to Appendix B. This projection simplifies the error equation substantially as we now see.

Lemma 4.1**.**

*Assume that the exact solution $(\boldsymbol{L},\boldsymbol{u},{p},\boldsymbol{J},\boldsymbol{b},r)$ of (2)-(3) is sufficiently regular. Then the exact solution satisfies (8)–(10). That is, if we replace $(\boldsymbol{L}_{h},\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\boldsymbol{b}_{h},r_{h},{\hat{\boldsymbol{u}}}_{h},\hat{\boldsymbol{b}}^{t}_{h},\hat{r}_{h})$ with $(\boldsymbol{L},\boldsymbol{u},{p},\boldsymbol{J},\boldsymbol{b},r,\boldsymbol{u},\boldsymbol{b}^{t},r)$ in (8)–(10), then (8)–(10) hold true for all $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s,\boldsymbol{\mu},\boldsymbol{\lambda}^{t},\gamma)\in\mathbf{G}_{h}\left(K\right)\times\mathbf{V}_{h}\left(K\right)\times\mathbf{Q}_{h}\left(K\right)\times\mathbf{H}_{h}\left(K\right)\times\mathbf{C}_{h}\left(K\right)\times\mathbf{S}_{h}\left(K\right)\times\mathbf{M}_{h}(e)\times\mathbf{\Lambda}^{t}_{h}(e)\times\mathbf{\Gamma}_{h}(e)$. *

Proof 4.2**.**

*Multiply (2) by $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)\in\mathbf{G}_{h}\left(K\right)\times\mathbf{V}_{h}\left(K\right)\times\mathbf{Q}_{h}\left(K\right)\times\mathbf{H}_{h}\left(K\right)\times\mathbf{C}_{h}\left(K\right)\times\mathbf{S}_{h}\left(K\right)$, integrate over $\Omega$, and integrate by parts. By the regularity assumptions on the solution, the solution components $((\boldsymbol{L}\boldsymbol{n})\otimes\boldsymbol{n},\boldsymbol{u},{p},\boldsymbol{J}^{t},\boldsymbol{b},r)$ are single valued on $\mathcal{E}$. and the exact solution satisfies (8). With the additional fact that $(\boldsymbol{w}\cdot\boldsymbol{n})\boldsymbol{n}$ and $\boldsymbol{{d}}$ are also single-valued on $\mathcal{E}_{h}$, we have that the exact solution satisfies (9). Finally, the boundary conditions (3) trivially imply that the exact solution satisfies (10). *

Lemma 4.3** (Error equation).**

The discretization errors satisfy

[TABLE]

Proof 4.4**.**

Using the fact that the numerical solution and exact solution both satisfy (8) (see Lemma 4.1), the linearity of the operators lead to the following error equations:

[TABLE]

Next, we split the error terms into their interpolation and approximation components as in (23) using the projections $\Pi$ defined in Appendix B. Due to the cancellation properties of $\Pi$ we obtain reduced error equations:

[TABLE]

More details of the cancellation properties of $\Pi$ used in the above formula are:

[TABLE]

Notice that (26) looks like (8), but with the approximation error replacing the finite element solution, and with some nonzero right hand side terms. Since the approximation error is in the finite element spaces, we can choose the test functions to be the approximation error terms. Similar to the procedure to arrive at (13), we take $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)=({\varepsilon}_{\boldsymbol{L}}^{h},{\varepsilon}_{\boldsymbol{u}}^{h},{\varepsilon}_{p}^{h},{\varepsilon}_{\boldsymbol{J}}^{h},{\varepsilon}_{\boldsymbol{b}}^{h},{\varepsilon}_{r}^{h})$, integrate by parts the first four terms of (26) and the first term of (26e), and sum the resulting equations in (26) to arrive at

[TABLE]

For the boundary conditions and conservation conditions, since the exact solution satisfies (9)–(10), we have

[TABLE]

We split the error terms into their interpolation and approximation components as before, and use the projections defined in Appendix B to cancel terms. More detailed roles for cancellations are

[TABLE]

and recall that the definitions of $\Pi$ for ${\hat{\boldsymbol{u}}}$, $\hat{\boldsymbol{b}}$, $\hat{r}$ are $L^{2}$ projections. Then we have

[TABLE]

In (29b), we can also erase ${\varepsilon}_{\hat{\boldsymbol{u}}}^{I}$ but we did not do it on purpose.

Equations (29d)–(29f) imply that on ${\partial\Omega}$, ${\varepsilon}_{\hat{\boldsymbol{u}}}^{h}=\boldsymbol{0}$, ${\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{h}=\boldsymbol{0}$, and ${\varepsilon}_{\hat{r}}^{h}=0$. In consideration of this zero contribution on $\partial\Omega$, summing of the formulae (29a)–(29c) with $({\boldsymbol{\mu},\boldsymbol{\lambda}^{t},\gamma})=({{\varepsilon}_{\hat{\boldsymbol{u}}}^{h},{\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{h},{\varepsilon}_{\hat{r}}^{h}})$ including $\partial\Omega$ gives

[TABLE]

Subtracting (4.4) from (27), we arrive at

[TABLE]

*Following the same procedure to get (16) from (15), we obtain the conclusion. *

Lemma 4.5**.**

There holds:

[TABLE]

Proof 4.6**.**

It is clear that bounding the energy is the same as bounding the right hand side of (24). The estimate of $\textup{Re}\left({\varepsilon}_{\boldsymbol{L}}^{I},{\varepsilon}_{\boldsymbol{L}}^{h}\right)$ is straightforward by Cauchy-Schwarz inequality. To estimate $\kappa\left({\varepsilon}_{\boldsymbol{b}}^{I},\nabla\times\left({\varepsilon}_{\boldsymbol{u}}^{h}\times\boldsymbol{{d}}\right)\right)$, note first that an algebraic computation gives

[TABLE]

The boundedness of the left hand side can be obtained by

[TABLE]

where $\mathbb{P}_{0}\boldsymbol{{d}}$ is the $L^{2}$ projection of $\boldsymbol{{d}}$ to the piecewise constant space on $K$, and here we used (57f), Hölder inequality $\|f_{1}f_{2}f_{3}\|_{L^{1}}\leq\|f_{1}\|_{L^{2}}\|f_{2}\|_{L^{\infty}}\|f_{3}\|_{L^{2}}$, the Bramble–Hilbert lemma (see e.g., [4]) and the inverse estimate in the last two inequalities.

For an estimate of $\kappa\left({\varepsilon}_{\boldsymbol{u}}^{I},\boldsymbol{{d}}\times\left(\nabla\times{\varepsilon}_{\boldsymbol{b}}^{h}\right)\right)$, we first note that

[TABLE]

due to (57j). A similar argument as above gives

[TABLE]

*Finally, we simply use the Cauchy-Schwarz/Hölder inequality for the last term on the right hand side of (24). *

Corollary 4.7** (Energy estimate).**

There holds:

[TABLE]

Proof 4.8**.**

*Apply Young’s inequality to each of the terms on the right side of (31) involving $\left\|{\varepsilon}_{\boldsymbol{L}}^{h}\right\|_{0}$ and $\|{{\varepsilon}_{\boldsymbol{b}^{t}}^{h}-{\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{h}}\|_{{\partial\Omega}_{h}}$. Note also that $\Pi{\hat{\boldsymbol{u}}}$ is the best approximation of $\boldsymbol{u}$ on ${{\partial\Omega}_{h}}$, so $\left\|{\varepsilon}_{\hat{\boldsymbol{u}}}^{I}\right\|_{{\partial\Omega}_{h}}$ is bounded by $\left\|{\varepsilon}_{\boldsymbol{u}}^{I}\right\|_{{\partial\Omega}_{h}}$. *

In the energy estimate (33), we do not have direct control on $\|{\varepsilon}_{\boldsymbol{u}}^{h}\|_{0}$ and $\|{\varepsilon}_{\boldsymbol{b}}^{h}\|_{0}$. In the following, we employ an indirect approach to control these quantities via a duality argument. A similar approach for the Oseen equation has been conducted in [8]. It is more complicated for our MHD system due to the coupling between fluids and electromagnetics. In particular, ${\varepsilon}_{\boldsymbol{u}}^{h}$ and ${\varepsilon}_{\boldsymbol{b}}^{h}$ are coupled and have to be simultaneously analyzed. To begin, we define a dual problem of the MHD system (2) as

[TABLE]

with homogeneous boundary conditions. Here, $\boldsymbol{\theta}$ and $\boldsymbol{\sigma}$ are two given functions in $L^{2}\left(\Omega\right)$, and the superscript “*” is used to denote the corresponding unknowns in the adjoint equation. We assume the following regularity estimate holds for the adjoint problem (34)

[TABLE]

The well-posedness of (34) and the conditions under which the regularity estimate (35) holds are discussed in Appendix C.

We use the interpolation operators $\Pi^{*}$ defined in Appendix A (77), (78) below and ${\varepsilon}_{\boldsymbol{L}^{*}}^{I}$, ${\varepsilon}_{p^{*}}^{I}$, … will denote $\boldsymbol{L}^{*}-\Pi^{*}\boldsymbol{L}^{*}$, ${p}^{*}-\Pi^{*}{p}^{*}$, etc. Testing (34b) with ${\varepsilon}_{\boldsymbol{u}}^{h}$ and (34e) with ${\varepsilon}_{\boldsymbol{b}}^{h}$ we have

[TABLE]

where we have used integration by parts in the second equality, the properties of the $\Pi^{*}$ operators (77) and (78) in the third equality, and integration by parts again in the last equality.

Lemma 4.9**.**

The following identities hold true:

[TABLE]

Proof 4.10**.**

In (26), choose test functions $(\boldsymbol{G},\boldsymbol{v},q,\boldsymbol{H},\boldsymbol{{c}},s)$ to be the adjoint projections $(\Pi^{*}\boldsymbol{L}^{*},\Pi^{*}\boldsymbol{u}^{*},\Pi^{*}{p}^{*},\Pi^{*}\boldsymbol{J}^{*},\Pi^{*}\boldsymbol{b}^{*},\Pi^{*}r^{*})$. Noting the identities

[TABLE]

*all the results follow immediately except the fifth one. In the fifth identity, note that the first equality is trivial by the definitions of $\Pi^{*}\boldsymbol{b}^{*}$ and ${\varepsilon}_{\boldsymbol{b}^{*}}^{I}$, and the second equality comes from (26) with the orthogonality of ${\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{I}$ and ${\varepsilon}_{\hat{r}}^{I}$. *

Substituting the result of Lemma 4.9 into (36) we obtain

[TABLE]

Note first that the last term in the second line and the first term in the last line can be simplified as $\kappa\left({\varepsilon}_{\boldsymbol{b}}^{h},\nabla\times({\varepsilon}_{\boldsymbol{u}^{*}}^{I}\times\boldsymbol{{d}})\right)$. Note also that we can use (34a), (34c), (34d), (34f) to get

[TABLE]

We can subtract the $L^{2}$ projections on the mesh skeleton, $\mathbb{P}_{e}\boldsymbol{u}^{*}$, $\mathbb{P}_{e}\boldsymbol{b}^{*}$, $\mathbb{P}_{e}r^{*}$ of $\boldsymbol{u}^{*}$, $\boldsymbol{b}^{*}$, $r^{*}$, in flux terms using the flux conservation of interpolation operators and numerical solutions, i.e., (29a), (29b), (29c). Moreover, we can subtract them on the domain boundary since ${\varepsilon}_{\hat{r}}^{h}$, $\boldsymbol{u}^{*}$, $(\boldsymbol{b}^{*})^{t}$, and $r^{*}$ are zero there. If we use the above three observations, i.e., the simplification of the two terms, substitution of volume integral terms using (37), and the subtraction of $L^{2}$ projections on the mesh skeleton, then we have (some terms are colored for readers’ convenience)

[TABLE]

Note that the integration by parts and properties of ${\varepsilon}_{\boldsymbol{u}^{*}}^{I}$ and $\mathbb{P}_{e}\boldsymbol{u}^{*}$ yield

[TABLE]

These two identities simplify the second and third lines of the above formula and yield

[TABLE]

To simplify it further, we note the identities

[TABLE]

where cancellations come from the properties of ${\varepsilon}_{\boldsymbol{b}^{*}}^{I}$, $\mathbb{P}_{e}\boldsymbol{b}^{*}$, $\mathbb{P}_{e}r^{*}$ and the orthogonalities of ${\varepsilon}_{\hat{r}}^{I}$ and ${\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{I}$ to polynomials. We use these identities to the previous formula, more precisely, to the last two terms in the third line, to the fourth line, to the last terms in the last two lines. Then we obtain (here some terms are colored only for readers’ convenience later)

[TABLE]

Here we used the facts $\boldsymbol{n}\times{\varepsilon}_{\boldsymbol{b}}^{h}=\boldsymbol{n}\times{\varepsilon}_{\boldsymbol{b}^{t}}^{h}$ and $\left<{\varepsilon}_{\boldsymbol{b}}^{h},\boldsymbol{n}\times{\varepsilon}_{\boldsymbol{J}^{*}}^{I}\right>=\left<{\varepsilon}_{\boldsymbol{b}^{t}}^{h},\boldsymbol{n}\times{\varepsilon}_{\boldsymbol{J}^{*}}^{I}\right>$. Algebraic manipulations with (7) yield

[TABLE]

We reduce this further, particularly from the third to ninth lines. For the third line, note that ${\varepsilon}_{\hat{\boldsymbol{u}}}^{h}$, $\boldsymbol{{d}}$, and the exact solution are single-valued, and $\left<{\varepsilon}_{\hat{\boldsymbol{u}}}^{h},\Pi^{*}\boldsymbol{u}^{*}-\mathbb{P}_{e}\boldsymbol{u}^{*}\right>=\left<{\varepsilon}_{\hat{\boldsymbol{u}}}^{h},{\varepsilon}_{\boldsymbol{u}^{*}}^{I}\right>$. Note also that $\Pi^{*}\boldsymbol{b}^{*}-\mathbb{P}_{e}\boldsymbol{b}^{*}=-{\varepsilon}_{\boldsymbol{b}^{*}}^{I}+{\varepsilon}_{\hat{\boldsymbol{b}}^{*}}^{I}$ holds. Then the third line is

[TABLE]

where the vanishing equality comes from (78c). For the fourth line, note first that

[TABLE]

where the second equality is valid because $\alpha_{1}$ is constant on each facet. Then the fourth line is

[TABLE]

where we use (78c) again. Note that

[TABLE]

From these the fifth and sixth lines are rewritten as

[TABLE]

Further, the seventh, eighth, ninth lines are vanishing, i.e.,

[TABLE]

where we used (77c) for the first two identities, and the fact $\boldsymbol{n}\times{\varepsilon}_{\boldsymbol{b}^{*}}^{I}=\boldsymbol{n}\times{\varepsilon}^{I}_{(\boldsymbol{b}^{*})^{t}}$ for the third identity. Finally, \left({\varepsilon}_{\boldsymbol{J}}^{h},-{\varepsilon}_{\boldsymbol{J}^{*}}^{I}\right)={\color[rgb]{1,0,0}0} by the definition of ${\varepsilon}_{\boldsymbol{J}^{*}}^{I}$. As a consequence, we have a reduced formula

[TABLE]

Estimation for $I_{1}$: Combining the estimate for ${\varepsilon}_{\boldsymbol{L}^{*}}^{I}$ and the regularity estimate (see (80)) gives

[TABLE]

Estimation for $I_{2}$: Using (34a), (57i), (57j), Lemma A.3, and the regularity of the adjoint solutions, we have

[TABLE]

Estimation for $I_{4}$: By the identity (32), it suffices to estimate

[TABLE]

By the triangle inequality, the inverse estimate, and (80), we have

[TABLE]

and we also have

[TABLE]

thus

[TABLE]

Estimation for $I_{5}$: By an argument similar to the estimate of $\left\|\nabla\Pi^{*}\boldsymbol{u}^{*}\right\|_{0}$ above, $\left\|\nabla\times\Pi^{*}\boldsymbol{b}^{*}\right\|_{0}\lesssim\left\|\boldsymbol{\theta},\boldsymbol{\sigma}\right\|_{0}$. Since $I_{5}=\left({\varepsilon}_{\boldsymbol{u}}^{I},(\boldsymbol{{d}}-\mathbb{P}_{0}\boldsymbol{{d}})\times(\nabla\times\Pi^{*}\boldsymbol{b}^{*})\right)$,

[TABLE]

Estimation for $I_{6}$ and $I_{7}$: Integrating $I_{6}$ by parts (see (7)) we have

[TABLE]

Now we can write $-\kappa I_{6}+I_{7}$ as

[TABLE]

For the first term, as in the estimate of $I_{4}$, it suffices to estimate

[TABLE]

Invoking the Hölder’s inequality and an inverse estimate we can bound the upper bounds of the first term as

[TABLE]

For the second term, we first observe that

[TABLE]

By the Hölder inequality,

[TABLE]

where we used the fact that $\Pi^{*}\hat{\boldsymbol{b}^{*}}$ and $\Pi^{*}\hat{\boldsymbol{u}^{*}}$ are the best approximations on ${{\partial\Omega}_{h}}$. By Lemma A.2 this can be estimated by $\left(h\kappa\left\|\boldsymbol{{d}}\right\|_{W^{1,\infty}}+h^{2}\left\|\boldsymbol{w}\right\|_{W^{1,\infty}}\right)\left\|{\varepsilon}_{\boldsymbol{u}}^{h}\right\|_{0}\left\|\boldsymbol{\theta},\boldsymbol{\sigma}\right\|_{0}$. As a consequence,

[TABLE]

Estimation for $I_{3}$, $I_{8}$, and $I_{9}$: Integrating $I_{3}$ by parts (see (7)) gives

[TABLE]

Some algebraic manipulations give

[TABLE]

The first term is easily estimated by

[TABLE]

For the second term we have

[TABLE]

where we used Lemma A.2 and the discrete trace inequality. Using the Cauchy–Schwarz inequality, (80), and Lemma A.2, the third term is bounded by

[TABLE]

Combining the above estimates we conclude

[TABLE]

Estimation for $I_{10}$: Using the approximation capability of the projector $\boldsymbol{\Pi}$ (see Appendix B) we have

[TABLE]

At this point we are ready to estimate the discretization errors for $\boldsymbol{L}$, $\boldsymbol{J}$, $\boldsymbol{u}$, and $\boldsymbol{b}$. For readability let us absorb $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$, Re, Rm, $\kappa$, and the norms on $\boldsymbol{{d}}$ and $\boldsymbol{w}$ into the implicit constants.

Theorem 4.11**.**

Suppose that $\alpha_{1}-\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}}$, $\alpha_{2}$, and $\alpha_{3}$ are chosen to be positive constants independent of $h$, Re, Rm, $\kappa$, $\boldsymbol{{d}}$, and $\boldsymbol{w}$. Suppose also that $h$ is sufficiently small, i.e.,

[TABLE]

with $C$ (depending on the coefficients in estimates (43) and (44)). Then it holds that

[TABLE]

and the following error estimates hold:

[TABLE]

Proof 4.12**.**

We proceed by taking $\boldsymbol{\theta}={\varepsilon}_{\boldsymbol{u}}^{h}$, $\boldsymbol{\sigma}={\varepsilon}_{\boldsymbol{b}}^{h}$ in (34). If we use (38), the estimates (39)–(45), and Young’s inequality, we can obtain

[TABLE]

which can be simplified to become

[TABLE]

if the constants that multiply $\left\|{\varepsilon}_{\boldsymbol{u}}^{h}\right\|_{0}$ and $\left\|{\varepsilon}_{\boldsymbol{b}}^{h}\right\|_{0}$ in (43) and (44) are sufficiently small, which is true under the assumptions we have made on $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{3}$ together with (46). The discretization error terms on the right hand side of (51) (i.e., terms with superscript “$h$”) are bounded by $E_{h}$ (see definition of $E_{h}$ in (31)). This implies

[TABLE]

Applying Young’s inequality to the right side of (33) for $\left\|{\varepsilon}_{\boldsymbol{u}}^{h}\right\|_{0}$, $\left\|{\varepsilon}_{\boldsymbol{b}}^{h}\right\|_{0}$, and using (52), we get

[TABLE]

Then (47) follows from the approximation properties of ${\varepsilon}_{\boldsymbol{J}}^{I},{\varepsilon}_{\boldsymbol{u}}^{I},{\varepsilon}_{\boldsymbol{b}}^{I},{\varepsilon}_{\boldsymbol{L}}^{I}$ with the trace inequality Lemma A.2 discussed in Appendix B and A, respectively. Further, it gives

[TABLE]

*where the first inequality is from the definition of $E_{h}$ in (24). Then (48) follows from the triangle inequality. Finally, the above estimate with (52) and the triangle ineqality give (49). *

What remains is to estimate $\left\|{\varepsilon}_{p}^{h}\right\|_{0}$ and $\left\|{\varepsilon}_{r}^{h}\right\|_{0}$, and to that end we extend the argument in [8] for our MHD system.

Theorem 4.13**.**

There holds:

[TABLE]

Proof 4.14**.**

We consider a projection operator $\tilde{\Pi}:\left[H^{1}(\Omega)\right]^{d}\rightarrow\left[{\mathcal{P}}_{k}({\Omega_{h}})\right]^{d}$ defined by

[TABLE]

for $\boldsymbol{\vartheta}\in\left[H^{1}(\Omega)\right]^{d}$ and $K\in{\Omega_{h}}$, where ${\mathcal{P}}_{k}^{\perp}(K)$ is the subspace of ${\mathcal{P}}_{k}(K)$ which is orthogonal to ${\mathcal{P}}_{k-1}(K)$ in $L^{2}(K)$. Its well-posedness is based on the orthogonal decomposition (see [17, Lemma 4.1])

[TABLE]

and it has optimal approximation property by the Bramble–Hilbert lemma.

Since $\Pi{p}$ is the $L^{2}$-projection, from (11) we have

[TABLE]

that is, ${\varepsilon}_{p}^{h}$ has zero mean. It is known [27] that there exists a function $\boldsymbol{\vartheta}\in\left[H_{0}^{1}(\Omega)\right]^{d}$ that $\nabla\cdot\boldsymbol{\vartheta}={\varepsilon}_{p}^{h}$ and $\left\|\boldsymbol{\vartheta}\right\|_{1}\lesssim\left\|{\varepsilon}_{p}^{h}\right\|_{0}$. For such a $\boldsymbol{\vartheta}$, we have

[TABLE]

where we have performed integration by parts twice and used the definition of the projector $\tilde{\Pi}$. Since the exact solution $\left(\boldsymbol{L},{p},\boldsymbol{u},\boldsymbol{b}\right)$ and its trace also satisfy the HDG local (sub-) equation (8b), we can add and subtract the corresponding projections in (57) to obtain

[TABLE]

where we have taken $\tilde{\Pi}\boldsymbol{\vartheta}$ as the test function in (8b).

Combining (53) and (4.14) yields

[TABLE]

which can be further simplified using two facts: first, integrating by parts twice and using the definition of $\tilde{\Pi}$ give

[TABLE]

and second, we combine the first equation in (9) and (57k) to have

[TABLE]

In particular, we obtain

[TABLE]

By the triangle and Hölder inequalities,

[TABLE]

*where we have used Lemma A.3, the approximation capability of $\tilde{\Pi}$ and the $L^{2}$-projection, definition of $E_{h}$ in (24), the property of $\boldsymbol{\vartheta}$, and we absorb all mesh independent parameters into the implicit constant in the final inequality. As a consequence, we have $\left\|{\varepsilon}_{p}^{h}\right\|_{0}\lesssim\left\|{\varepsilon}_{\boldsymbol{L}}^{h},{\varepsilon}_{\boldsymbol{u}}^{h},{\varepsilon}_{\boldsymbol{b}}^{h},{\varepsilon}_{\boldsymbol{b}}^{I}\right\|_{0}+h^{\frac{1}{2}}E_{h}$. Then the conclusion follows from the triangle inequality and the estimates of $\left\|{\varepsilon}_{\boldsymbol{L}}^{h},{\varepsilon}_{\boldsymbol{u}}^{h},{\varepsilon}_{\boldsymbol{b}}^{h},{\varepsilon}_{\boldsymbol{b}}^{I}\right\|_{0}$ and $E_{h}$. *

For an analogous result for $\left\|{\varepsilon}_{r}^{h}\right\|_{0}$, we make use of the following Lemma.

Lemma 4.15**.**

*There exists some $\left[\boldsymbol{\vartheta}\in H^{1}(\Omega)\right]^{d}$ such that $\nabla\cdot\boldsymbol{\vartheta}={\varepsilon}_{r}^{h}$, $\boldsymbol{n}\times\boldsymbol{\vartheta}=\boldsymbol{0}$ on ${\partial\Omega}$, and $\left\|\boldsymbol{\vartheta}\right\|_{1}\lesssim\left\|{\varepsilon}_{r}^{h}\right\|_{0}$. *

Proof 4.16**.**

Let $\overline{{\varepsilon}_{r}^{h}}$ be the mean-value of ${\varepsilon}_{r}^{h}$ in $\Omega$ and $1_{\Omega}$ be the indicator function on $\Omega$. It is well known that there exists $\boldsymbol{\vartheta}_{0}$ in $\left[H_{0}^{1}(\Omega)\right]^{d}$ such that $\nabla\cdot\boldsymbol{\vartheta}_{0}={\varepsilon}_{r}^{h}-\overline{{\varepsilon}_{r}^{h}}1_{\Omega}$, $\left\|\boldsymbol{\vartheta}_{0}\right\|_{1}\lesssim\|{\varepsilon}_{r}^{h}-\overline{{\varepsilon}_{r}^{h}}1_{\Omega}\|_{0}$. Since $\Omega$ is a bounded domain with polyhedral boundary, $\partial\Omega$ consists of finitely many piecewise smooth $(d-1)$-dimensional components, so there exists a nonnegative smooth function $\phi_{0}\not=0$ on $\partial\Omega$ supported on the interior of the smooth components. Let $\phi$ be the renormalized function of $\phi_{0}$ satisfying $\int_{\partial\Omega}\phi=\overline{{\varepsilon}_{r}^{h}}|\Omega|$. Then $\phi\boldsymbol{n}\in[H^{1/2}(\partial\Omega)]^{d}$ holds because $\phi$ is a smooth function on $\partial\Omega$, and there is $\Phi\in[H^{1}(\Omega)]^{d}$ such that $\Phi|_{\partial\Omega}=\phi\boldsymbol{n}$ and (cf. [24, Theorem II.4.3])

[TABLE]

According to a result in [24, p.176], there exists $\boldsymbol{\vartheta}_{1}$ in $\left[H^{1}(\Omega)\right]^{d}$ such that $\nabla\cdot\boldsymbol{\vartheta}_{1}=\overline{{\varepsilon}_{r}^{h}}1_{\Omega}$, $\boldsymbol{\vartheta}_{1}=\phi\boldsymbol{n}$ on $\partial\Omega$, and also satisfies

[TABLE]

*The conclusion follows by taking $\boldsymbol{\vartheta}=\boldsymbol{\vartheta}_{0}+\boldsymbol{\vartheta}_{1}$. *

Theorem 4.17**.**

There holds:

[TABLE]

Proof 4.18**.**

*First choose a function $\boldsymbol{\vartheta}$ that satisfies the statements of Lemma 4.15, and consider *

[TABLE]

From (26e) with $\boldsymbol{{c}}=\tilde{\Pi}\boldsymbol{\vartheta}$, we have

[TABLE]

If we use this to the previous identity, and also use the flux condition (28b) with $\mathbb{P}_{e}\boldsymbol{\vartheta}$ as the test function, we have

[TABLE]

where the surface integrals on the domain boundary are zero due to the fact that $\boldsymbol{n}\times\boldsymbol{\vartheta}=\boldsymbol{0}$ on ${\partial\Omega}$ and (28f). Note that

[TABLE]

With these identities and (28b) we have

[TABLE]

By the triangle, Hölder, and trace inequalities,

[TABLE]

Since $\left\|\nabla\tilde{\Pi}\boldsymbol{\vartheta}\right\|_{0}\lesssim\left\|\boldsymbol{\vartheta}\right\|_{1}$ and $\left\|\boldsymbol{\vartheta}\right\|_{1}\lesssim\left\|{\varepsilon}_{r}^{h}\right\|_{0}$, we have

[TABLE]

*The conclusion follows from the triangle inequality and the estimates of $\left\|{\varepsilon}_{\boldsymbol{J}}^{h},{\varepsilon}_{\boldsymbol{u}}^{I},{\varepsilon}_{\boldsymbol{u}}^{h}\right\|_{0}$ and $E_{h}$. *

5 Numerical Results

In this section we apply the proposed HDG scheme for 2D MHD problems. The application for large-scale 3D problems is ongoing and will be presented elsewhere. The first problem we consider is the Hartmann flow whose analytical solution exists and is one dimensional in nature. The second problem is posed on a non-convex domain to demonstrate the approximation capability of our proposed HDG scheme though the non-convexity is not covered by our theory. To further challenge the proposed HDG approach, we will consider a singular problem as the third example.

Before presenting the results for each example, we make some general observations to differentiate the proposed HDG scheme from the DG method of [32]. First, recalling the primary motivation of HDG schemes, we have a reduced global system to solve for, which offers computational savings, especially for high-order solution of large-scale problems. Second, the HDG scheme allows the convenience of equal polynomial order approximations of all unknowns with direct approximations of $\nabla\boldsymbol{u}$ and $\nabla\times\boldsymbol{b}$. On the other hand, the DG method of [32] employs polynomials of order $k-1$ and $k+1$ for $p_{h}$ and $r_{h}$, respectively, where $k$ is the order of approximation of $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$, and $\nabla\boldsymbol{u}$, $\nabla\times\boldsymbol{b}$ are approximated by the derivatives of $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$. These differences make a direct comparison of the two methods difficult. The DG scheme is proven to be optimal (converging with $\mathcal{O}({h^{k}})$) in the DG energy norms for sufficiently smooth solutions. By the definitions of these norms in [32], $\boldsymbol{u}_{h}$, $\nabla\boldsymbol{u}_{h}$, $\boldsymbol{b}_{h}$, and $\nabla\times\boldsymbol{b}_{h}$ are proven to converge at worst with $\mathcal{O}(h^{k})$ in the $L^{2}$ norm. Optimal $L^{2}$ convergence of $\mathcal{O}(h^{k+1})$ for $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$ is not proven in [32], but is demonstrated in all numerical examples with smooth solutions. On the other hand, the HDG scheme is proven to be $L^{2}$ optimal ($\mathcal{O}({h^{k+1}})$) for $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$ and quasi-optimal ($\mathcal{O}({h^{k+\frac{1}{2}}})$) for the remaining local quantities. These rates are demonstrated (and sometimes exceeded) in the numerical examples involving smooth solutions. For the numerical example involving a singular solution, the HDG scheme gives similar $L^{2}$ error magnitudes and convergence rates as in the DG scheme for $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$ for the same polynomial order $k$.

We note that the numerical solutions $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$ in the numerical results below satisfy the divergence-free constraint in the weak sense. However, the pointwise divergence-free constraint can be fulfilled by a post-processing procedure outlined in Appendix D, which shows that convergence rates of the post-processed solutions to the exact solutions are as good as those of $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$. For that reason, we do not use the post-processed solution in the error computation.

5.1 Hartmann Flow

In this numerical study, we consider a conducting incompressible fluid (liquid metal, for example) in a domain $[-\infty,\infty]\times[-l_{0},l_{0}]\times[-\infty,\infty]$ (bounded by infinite parallel plates in the $x_{2}$ direction [32, 49]). The fluid is subject to a uniform pressure gradient $G:=-\frac{\partial p}{\partial x_{1}}$ in the $x_{1}$ direction, and a uniform external magnetic field $b_{0}$ in the $x_{2}$ direction. If we consider no-slip boundary conditions on the $x_{2}$ boundaries, the resulting flow pattern is known as Hartmann flow which admits an analytical solution which is one dimensional in nature. We assume that the infinite parallel plates are perfectly insulating. Here, we consider the Hartmann flow in a 2D domain $\Omega=[0,0.025]\times[-1,1]$. If we define the characteristic velocity as $u_{0}:=\sqrt{Gl_{0}/\rho}$, and consider the driving pressure gradient $G$ as a forcing term (incorporated in $\boldsymbol{g}$), the non-dimensionalized solution with $\boldsymbol{g}=(1,0)$, $\boldsymbol{f}=(0,0)$ reads

[TABLE]

where $\textup{Ha}:=\kappa\textup{Re}\textup{Rm}$, and $p_{0}$ is a constant that enables $p$ to satisfy the zero average pressure condition (5). We set $\boldsymbol{w}=\boldsymbol{u}$ and $\boldsymbol{{d}}=\boldsymbol{b}$, and we enforce the boundary conditions on ${\partial\Omega}$ using the exact solution, i.e., $\boldsymbol{u}_{D}=\boldsymbol{u}$, $\boldsymbol{h}_{D}=\boldsymbol{b}^{t}$, and $r_{D}=0$.

At refinement level $l$, the domain is divided into $l\times 80l$ squares, each of which is divided into two triangles from top right to bottom left. In Figure 1 are the convergence plots with $\textup{Re}=\textup{Rm}=7.07$ and $\kappa=200$. The convergence rates for $\boldsymbol{L}_{h}$, $\boldsymbol{u}_{h}$, $p_{h}$, $\boldsymbol{J}_{h}$, $\boldsymbol{b}_{h}$, and $r_{h}$ are observed to be approximately $k+\frac{1}{2}$, $k+1$, $k+\frac{1}{2}$, $k+1$, $k+1$, and $k+1$, respectively. These observed rates approximately match or exceed their respective theoretical rates of $k+\frac{1}{2}$, $k+1$, $k+\frac{1}{2}$, $k+\frac{1}{2}$, $k+1$, and $k+\frac{1}{2}$ which were proven in Section 4.

5.2 Non-convex Domain

This example illustrates the convergence of the HDG scheme applied to a problem posed on the non-convex domain $\Omega=(-1,1)\times(-1,1)\setminus[0,1)\times(-1,0]$ in Figure 2 (similar to Section 5.1.1 in [32]). We take $\textup{Re}=\textup{Rm}=\kappa=1$, $\boldsymbol{w}=(2,1)$, and $\boldsymbol{{d}}=(x_{1},-x_{2})$. We set $\boldsymbol{g}$ and $\boldsymbol{f}$ such that the manufactured solution for (1) is the following

[TABLE]

where $p_{0}$ is the constant that enables $p$ to satisfy the zero average pressure condition (5). We use the exact solution to enforce the boundary conditions ${\partial\Omega}$, i.e., $\boldsymbol{u}_{D}=\boldsymbol{u}$, $\boldsymbol{h}_{D}=\boldsymbol{b}^{t}$, and $r_{D}=r$.

At refinement level $l$, each quadrant of the domain (see Figure 2 for an example with $l=4$) is subdivided into $l\times l$ squares, each of which is divided into two triangles from top right to bottom left. In Figure 3 are the convergence plots. For this problem, we observe the optimal convergence rates of $k+1$ for all of the local variables, which matches or exceeds the rates proven in Section 4.

5.3 Singular Solution

Although we do not discuss the implications of singular solutions on the theoretical convergence rates of the HDG scheme, applying the scheme to such a problem is instructive in assessing its robustness. This example illustrates the convergence of the HDG scheme using a manufactured solution with a singularity (similar to the example in Section 5.2 of [32]). In particular, we consider the same non-convex domain and mesh refinement as in the previous example (see Figure 2). We take $\textup{Re}=\textup{Rm}=\kappa=1$, $\boldsymbol{w}=\boldsymbol{0}$, and $\boldsymbol{{d}}=(-1,1)$. We choose $\boldsymbol{g}$ and $\boldsymbol{f}$ such that the analytical solution of (1) has the form

[TABLE]

where

[TABLE]

On ${\partial\Omega}$ we use the exact solution to set the boundary condition, i.e., $\boldsymbol{u}_{D}=\boldsymbol{u}$, $\boldsymbol{h}_{D}=\boldsymbol{b}^{t}$, and $r_{D}=r$. For this problem, it is known that $\boldsymbol{u}\in\left[H^{1+\lambda}(\Omega)\right]^{2}$, ${p}\in H^{\lambda}(\Omega)$, and $\boldsymbol{b}\in\left[H^{2/3}(\Omega)\right]^{2}$, and that the solution contains magnetic and hydrodynamic singularities that are among the strongest singularities [32].

Convergence results for this problem are shown in Figure 4. For the fluid variables $\boldsymbol{L}_{h}$, $\boldsymbol{u}_{h}$, and $p_{h}$, we observe convergence rates of approximately $\lambda$, $2\lambda$, and $\lambda$, respectively. For the magnetic variables $\boldsymbol{J}_{h}$, $\boldsymbol{b}_{h}$, and $r_{h}$, we observe convergence rates of approximately $1/2$, $2/3$, and $1/3$, respectively.

5.4 3D numerical experiments on cubical meshes

We show numerical results for a three dimensional problem with our HDG method adapted to hexahedral meshes with tensor product polynomial spaces. Our theoretical analysis is only on the method on tetrahedral meshes, so it does not support this method on cubical meshes. Nonetheless we present this numerical result here in order to demonstrate that the HDG method can be applied to 3D problems and can be implemented using hexahedral meshes.

We set $\Omega=[0,1]^{3}$, $\boldsymbol{w}=(1,2,-4)$, $\boldsymbol{{d}}=(-3,1,5)$, and set the forcing functions and boundary conditions to solve for the manufactured solution

[TABLE]

with Dirichlet boundary conditions applied on $\partial\Omega$ for $\boldsymbol{u}$, $r$, and the tangential components of $\boldsymbol{b}$.

Convergence rates of $L^{2}$-errors with respect to uniform mesh refinements are given in Figure 5. The convergence rates of $p$ and $r$ are optimal with order $k+1$, and the convergence rate of $\boldsymbol{u}$ is suboptimal with $k+\frac{1}{2}$. On the contrary, the errors of $\boldsymbol{L}$, $\boldsymbol{J}$, $\boldsymbol{b}$ seem to have slightly lower convergence rates between $k$ and $k+\frac{1}{2}$.

6 Conclusions and future work

In this paper we have constructed an HDG method for a linearization of the incompressible resistive magnetohydrodynamics equations. We have carried out the a priori error analysis using elaborate interpolation operators, a duality argument with elliptic regularity assumptions, and an energy approach. Specifically, this allows us to prove optimal convergence for the velocity variable, $\boldsymbol{u}$, and the magnetic field variable, $\boldsymbol{b}$, and quasi-optimal convergence for the remaining quantities $\boldsymbol{L}$, ${p}$, $\boldsymbol{J}$, and $r$. Numerical performances of the method are tested on three examples: the Hartmann flow, a manufactured solution over a non-convex domain, and singular solution on a non-convex domain. The numerical results show that the theoretical convergence rates of all the unknowns are obtained for smooth solutions even when the elliptic regularity assumption of domain fails to hold. Ongoing work includes 3D computation on parallel computers for large-scale problems and extensions of our HDG method to nonlinear time-dependent magnetohydrodynamics equations.

Appendix A Auxiliary results

In this appendix we collect some technical results that are useful for our analysis.

Lemma A.1** (Inverse Inequality. [44, Lemma 1.44]).**

For $v\in\mathcal{P}_{k}(K)$ with $K\in{\Omega_{h}}$, there exists $C>0$ independent of $h$ such that

[TABLE]

Lemma A.2** (Trace inequality. [44, Lemma 1.49]).**

For $v\in H^{1}({\Omega_{h}})$ and for $K\in{\Omega_{h}}$ with $e\subset{\partial K}$, there exists $C>0$ independent of $h$ such that

[TABLE]

Applying the arithmetic-geometric mean inequality to the right side, we can derive

[TABLE]

If $v\in H^{1}({\Omega_{h}})$ is in piecewise polynomial spaces, we can derive the following inequality from Lemma A.2 and the inverse inequality (Lemma A.1):

[TABLE]

Lemma A.3**.**

Suppose that $\Pi:H^{1}(K)\rightarrow\mathcal{P}_{k}(K)$ is a bounded interpolation which is a projection on $\mathcal{P}_{k}(K)$. Then

[TABLE]

Proof A.4**.**

For any constant $c$, $\nabla\Pi\boldsymbol{v}=\nabla\left(\Pi\boldsymbol{v}-c\right)$, so

[TABLE]

*where we have used the Poincaré inequality and the approximation property of $\Pi$ in the last inequality. *

Appendix B Definition of projections and their properties

In this section we use $\mathcal{P}_{k}$, $\boldsymbol{\mathcal{P}}_{k}$, $\tilde{\boldsymbol{\mathcal{P}}}_{k}$ for spaces of scalar, $d$-dimensional vector, $d\times d$ matrix-valued polynomials. By $\mathcal{P}_{k}^{\perp}$, $\boldsymbol{\mathcal{P}}_{k}^{\perp}$, $\tilde{\boldsymbol{\mathcal{P}}}_{k}^{\perp}$ we denote the spaces of polynomials of order at most $k$ orthogonal to all polynomials of order at most $\left(k-1\right)$. $\boldsymbol{\mathcal{P}}_{k}^{t}(e)$ contains the tangential component of all polynomials in $\boldsymbol{\mathcal{P}}_{k}(e)$. We desire to have error equations conform to the original equations to facilitate the error analysis. To begin, we define a collective interpolation operator $\boldsymbol{\Pi}\left(\boldsymbol{L},\boldsymbol{u},{p},\boldsymbol{J},\boldsymbol{b},r,{\hat{\boldsymbol{u}}},\hat{\boldsymbol{b}}^{t},\hat{r}\right)$ implicitly through the interpolation errors ${\varepsilon}_{\boldsymbol{u}}^{I}=\boldsymbol{u}-\Pi\boldsymbol{u}$, ${\varepsilon}_{\boldsymbol{b}}^{I}=\boldsymbol{b}-\Pi\boldsymbol{b}$, etc, where $\Pi\boldsymbol{u}$, $\Pi\boldsymbol{b}$, etc, are components of the collective interpolator $\boldsymbol{\Pi}$ on $\boldsymbol{u},\boldsymbol{b}$, etc. Specifically:

• •

$L^{2}$ projections on $e\in{\mathcal{E}_{h}}$ or on $K\in{\Omega_{h}}$ are defined as:

[TABLE]

The well-definedness and optimality of the $L^{2}$-projections are clear. The coupled projector $\Pi\left(\boldsymbol{b},r\right):=\left(\Pi\boldsymbol{b},\Pi r\right)$ has been studied in [15], and in particular we have

[TABLE]

where, again, for simplicity we choose the same solution order $k$ for all the unknowns. Here, we assume that $\boldsymbol{b}$ and $r$ are sufficiently smooth, that is, $\boldsymbol{b}\in\left[H^{k+1}\left(\Omega\right)\right]^{d}$ and $r\in H^{k+1}\left(\Omega\right)$.

To understand the approximation capability of the coupled projector $\Pi\left(\boldsymbol{L},\boldsymbol{u}\right):=\left(\Pi\boldsymbol{L},\Pi\boldsymbol{u}\right)$ we need to recall a result in [15, Lemma A.1].

Lemma B.1**.**

*Suppose that $w\in\mathcal{P}_{k}^{\perp}(K)$. Then, for any $e\subset{\partial K}$, the map $w\mapsto w|_{e}\in\mathcal{P}_{k}(e)$ is an isomorphism and $\|w\|_{0,K}^{2}\sim h_{K}\|w\|_{0,e}^{2}$ holds with a constant independent of $h_{K}$. *

Lemma B.2** (Estimation for ${\varepsilon}_{\boldsymbol{u}}^{I}$).**

Suppose $\boldsymbol{u}\in\left[{H}^{k+1}\left(\Omega\right)\right]^{d}$, $\boldsymbol{L}\in\left[{H}^{k+1}\left(\Omega\right)\right]^{d\times d}$, $r\in H^{k+1}\left(\Omega\right)$, $\boldsymbol{b}\in\left[{H}^{k+1}\left(\Omega\right)\right]^{d}$, and ${p}\in H^{k+1}\left(\Omega\right)$. The projection $\Pi\boldsymbol{u}$ is well-defined and optimal, i.e.,

[TABLE]

*where $C(\alpha_{1},\boldsymbol{w})=1/(\alpha_{1}-\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}})$. *

Proof B.3**.**

We extend the proof of a result in [8]. To begin, we define $\boldsymbol{g}:=-\frac{1}{2}\kappa\boldsymbol{{d}}\times\left(\boldsymbol{n}\times\left({\varepsilon}_{\boldsymbol{b}^{t}}^{I}+{\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{I}\right)\right)$, take $\boldsymbol{\mu}=\boldsymbol{v}|_{\partial K}$ for some $\boldsymbol{v}\in\boldsymbol{\mathcal{P}}_{k}^{\perp}(K)$ which will be determined later, and rewrite (57k) as

[TABLE]

where we have used the integration by parts in the second equality, definitions of the projections $\Pi\boldsymbol{L}$, $\Pi\boldsymbol{u}$, and $\Pi{p}$, the orthogonality between $\nabla\cdot\left(\Pi\boldsymbol{L}\right),\nabla\left(\Pi{p}\right)$ and $\boldsymbol{v}\in\boldsymbol{\mathcal{P}}_{k}^{\perp}(K)$, and the orthogonality ${\varepsilon}_{p}^{I}\perp\nabla\cdot\boldsymbol{v}$ in the last equality. Now, let $\mathbb{P}_{k}\boldsymbol{u}$ be the $L^{2}$ projection of $\boldsymbol{u}$ and define $\delta_{\boldsymbol{u}}^{I}:=\boldsymbol{u}-\mathbb{P}_{k}\boldsymbol{u}$. By the triangle inequality, it suffices to estimate the approximation capability of $\delta_{\boldsymbol{u}}:=\mathbb{P}_{k}\boldsymbol{u}-\Pi\boldsymbol{u}$. From the above formula, we have

[TABLE]

Since $\delta_{\boldsymbol{u}}\in\boldsymbol{\mathcal{P}}_{k}^{\perp}(K)$, we can take $\boldsymbol{v}=\delta_{\boldsymbol{u}}$ to obtain

[TABLE]

where we have used $\nabla\cdot\boldsymbol{w}=0$ in the integration by parts

[TABLE]

By Lemma B.1 and the fact that $\alpha_{1}>\frac{1}{2}\left\|\boldsymbol{w}\right\|_{L^{\infty}}$, we have

[TABLE]

We now estimate $|F_{\boldsymbol{u}}(\delta_{\boldsymbol{u}})|$. Defining $\delta_{\boldsymbol{w}}=\boldsymbol{w}-\mathbb{P}_{0}\boldsymbol{w}$, note that $(\delta_{\boldsymbol{u}}^{I}\otimes\boldsymbol{w},\nabla\delta_{\boldsymbol{u}})_{K}=(\delta_{\boldsymbol{u}}^{I}\otimes\delta_{\boldsymbol{w}},\nabla\delta_{\boldsymbol{u}})_{K}$ holds due to the definition of $\delta_{\boldsymbol{u}}^{I}$. Thus, we have

[TABLE]

where we have used $\|\delta_{\boldsymbol{w}}\|_{L^{\infty}(K)}\lesssim h_{K}\|\boldsymbol{w}\|_{W^{1,\infty}(K)}$, the inverse inequality, and the continuous and discrete trace inequalies ((55) and (56), respectively) in the last step. Taking the approximation capability of $\mathbb{P}_{k}\boldsymbol{u}$ into account, we get

[TABLE]

The estimate of $|F_{\boldsymbol{L}}(\delta_{\boldsymbol{u}})|$ is straightforward since $\delta_{\boldsymbol{u}}\in\boldsymbol{\mathcal{P}}_{k}^{\perp}(K)$:

[TABLE]

For the estimate of $|F_{\boldsymbol{g}}(\delta_{\boldsymbol{u}})|$, note that

[TABLE]

Using the definition of $\Pi^{*}\hat{\boldsymbol{b}}^{t}$ as the $L^{2}$-projection, the continuous and discrete trace inequalies ((55) and (56), respectively), the above estimate of $\left\|\nabla{\varepsilon}_{\boldsymbol{b}}^{I}\right\|_{0,K}$ (63), and finally the estimate of $\left\|{\varepsilon}_{\boldsymbol{b}}^{I}\right\|_{0,K}$ in (58a), we have

[TABLE]

*Using (61), (62), and (64) in (60), and using the triangle inequality $\|{\varepsilon}_{\boldsymbol{u}}^{I}\|_{0,K}\lesssim\|\boldsymbol{u}-\mathbb{P}_{k}\boldsymbol{u}\|_{0,K}+\|\delta_{\boldsymbol{u}}\|_{0,K}$ ends the proof. *

To estimate $\|{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$, we need some auxiliary results. We first recall a result with a sketch of its proof.

Lemma B.4**.**

Let $e_{K}$ be a fixed face of the simplex $K$. For $\boldsymbol{R}\in[L^{2}(K)]^{d}$ and $g\in L^{2}({\partial K})$, we define $\Pi(\boldsymbol{R},g)\in\boldsymbol{\mathcal{P}}_{k}(K)$ as

[TABLE]

*Then, $\|\Pi(\boldsymbol{R},g)\|_{K}\lesssim\|\boldsymbol{R}\|_{K}+h_{K}^{1/2}\|g\|_{{\partial K}}$. *

Proof B.5**.**

We refer to [11] for the existence and uniqueness of $\Pi(\boldsymbol{R},g)$. Let $\boldsymbol{\sigma}_{1}=\Pi(\boldsymbol{R},0)$ and $\boldsymbol{\sigma}_{2}=\Pi(\boldsymbol{0},g)$. By the standard scaling argument,

[TABLE]

To estimate $\boldsymbol{\sigma}_{2}$, note that there exists $a_{e}\in{\mathbb{R}}$, $e\not=e_{K}$ such that $(1\;0\;0)^{T}=\sum_{e,e\not=e_{K}}a_{e}\boldsymbol{n}_{e}$, and the first component of $\boldsymbol{\sigma}_{2}$, say $\boldsymbol{\sigma}_{2}^{1}$, is

[TABLE]

*Since $\boldsymbol{\sigma}_{2}\cdot\boldsymbol{n}_{e}\perp\mathcal{P}_{k-1}(K)$ by the definition of $\boldsymbol{\sigma}_{2}$, $\|\boldsymbol{\sigma}_{2}\cdot\boldsymbol{n}_{e}\|_{0,K}\lesssim h_{K}^{1/2}\|\boldsymbol{\sigma}_{2}\cdot\boldsymbol{n}_{e}\|_{0,e}\lesssim h_{K}^{1/2}\|g\|_{0,e}$ by Lemma B.1. The estimate $\|\boldsymbol{\sigma}_{2}\|_{0,K}\lesssim h_{K}^{1/2}\|g\|_{0,{\partial K}}$ follows easily by using this inequality to each component of $\boldsymbol{\sigma}_{2}$. *

We now recall other known facts without proofs (cf. Lemma 4.8 in [14]).

Lemma B.6**.**

*For a face $e$ of $K$, let $\mathcal{B}_{e}$ be an orthogonal basis of the vectors orthogonal to $\boldsymbol{n}_{e}$, and let $\mathcal{B}=\left\{\mathbb{I}_{d}\right\}\cup\left\{\boldsymbol{t}\otimes\boldsymbol{n}_{e},\boldsymbol{t}\in\mathcal{B}_{e}\right\}$. This $\mathcal{B}$ is a basis of the space of $d\times d$ matrices. *

Lemma B.7** (Estimation for ${\varepsilon}_{\boldsymbol{L}}^{I}$).**

Assume $\boldsymbol{u}\in\left[H^{k+1}\left(\Omega\right)\right]^{d}$, $\boldsymbol{L}\in\left[{H}^{k+1}\left(\Omega\right)\right]^{d\times d}$, $r\in H^{k+1}\left(\Omega\right),\boldsymbol{b}\in\left[{H}^{k+1}\left(\Omega\right)\right]^{d}$, and ${p}\in H^{k+1}\left(\Omega\right)$. Furthermore, suppose the trace of the tensor $\boldsymbol{L}$ vanishes, i.e., $\operatorname{tr}\boldsymbol{L}=0$. There holds:

[TABLE]

Proof B.8**.**

We proceed in a manner similar to [8, Theorem 2.3] with adaptations corresponding to our more complicated projectors $\boldsymbol{\Pi}\left(\boldsymbol{L},\boldsymbol{u}\right)$.

The dual basis of $\mathcal{B}$ (see Lemma B.6) can be written as

[TABLE]

where $W_{e,\boldsymbol{t}}:(\boldsymbol{t}\otimes\boldsymbol{n}_{e})=1$ for the $e$ and $\boldsymbol{t}$ corresponding to the subscripts of $W$, and 0 otherwise. Any $d\times d$ matrix, $A$, can be written as

[TABLE]

so

[TABLE]

Since $W_{e,\boldsymbol{t}}$ is an element of $\mathcal{B}^{*}$ independent of mesh size, this identity reduces the estimate of $\|{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$ to the estimates of $\|{\varepsilon}_{\boldsymbol{L}}^{I}\boldsymbol{n}_{e}\cdot\boldsymbol{t}\|_{0,K}$ with $\boldsymbol{t}\in\mathcal{B}_{e}$ and $\|\operatorname{tr}{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$.

We first estimate $\|{\varepsilon}_{\boldsymbol{L}}^{I}\boldsymbol{n}\cdot\boldsymbol{t}\|_{0,K}$ with $\boldsymbol{n}=\boldsymbol{n}_{e}$ for some $e$. Let $e_{K}$ be a fixed face of $K$ and define $\Pi_{1}\boldsymbol{L},\Pi_{2}\boldsymbol{L}\in\tilde{\boldsymbol{\mathcal{P}}}_{1}(K)$ as

[TABLE]

and

[TABLE]

The existence and uniqueness of $\Pi_{1}\boldsymbol{L}$ and $\Pi_{2}\boldsymbol{L}$ follow from Lemma B.4. By the triangle inequality,

[TABLE]

Again by the triangle inequality, we bound the first term in (68) as

[TABLE]

where $\Pi^{RTN}$ is the row-wise canonical Raviart-Thomas-Nédélec (RTN) interpolation operator into the row-wise $(k+1)$-th order RTN element, which contains $\tilde{\boldsymbol{\mathcal{P}}}_{k}(K)$ for all $K\in{\Omega_{h}}$. From [11, Proposition 2.1 (vi)], we have that $\left\|\Pi^{RTN}\boldsymbol{L}-\Pi_{1}\boldsymbol{L}\right\|_{0,K}\lesssim h_{K}^{k+1}\left\|\mathbb{P}_{k}\nabla\cdot\boldsymbol{L}\right\|_{k}$, and from a well known property of the canonical RTN interpolation operator we have that $\left\|\boldsymbol{L}-\Pi^{RTN}\boldsymbol{L}\right\|_{0,K}\lesssim h_{K}^{k+1}\left\|\boldsymbol{L}\right\|_{k+1}$. Therefore, for the first term of (68) we have

[TABLE]

Note that $\Pi^{RTN}\boldsymbol{L}$ is not necessarily in $\tilde{\boldsymbol{\mathcal{P}}}_{k}({\Omega_{h}})$ but the above argument does not require $\Pi^{RTN}\boldsymbol{L}\in\tilde{\boldsymbol{\mathcal{P}}}_{k}({\Omega_{h}})$.

For the estimate of the second term in (68), note that the definitions of $\Pi_{1}$ and $\Pi_{2}$ give

[TABLE]

By Lemma B.4, we can estimate the second term in (68) as

[TABLE]

For the estimate of the third term in (68), recalling (57i), (57j), and (67a), we derive $\left(\Pi_{2}\boldsymbol{L}-\Pi\boldsymbol{L},\boldsymbol{G}\right)_{K}=0$ for all $\boldsymbol{G}\in\tilde{\boldsymbol{\mathcal{P}}}_{k-1}(K)$. Selecting $\boldsymbol{G}=\left(\boldsymbol{t}\otimes\boldsymbol{n}\right)q$ with $q\in\mathcal{P}_{k-1}(K)$, we have that $\left(\Pi_{2}-\Pi\right)\boldsymbol{L}\boldsymbol{n}\cdot\boldsymbol{t}\in\mathcal{P}_{k}^{\perp}(K)$, and by Lemma B.1,

[TABLE]

for any $e$ of ${\partial K}$. From (57k) and (67b), we have, for $e\neq e_{K}$

[TABLE]

with $\boldsymbol{g}:=-\frac{1}{2}\kappa\boldsymbol{{d}}\times\left(\boldsymbol{n}\times\left({\varepsilon}_{\boldsymbol{b}^{t}}^{I}+{\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{I}\right)\right)$. Choosing $\boldsymbol{\mu}=\left[\left(\Pi_{2}-\Pi\right)\boldsymbol{L}\boldsymbol{n}\cdot\boldsymbol{t}\right]\boldsymbol{t}$ and applying the Cauchy-Schwarz inequality to the above expression, we have

[TABLE]

where we have used the fact that $\left\|{\varepsilon}_{\hat{\boldsymbol{b}}^{t}}^{I}\right\|_{0,e}\leq\left\|{\varepsilon}_{\boldsymbol{b}^{t}}^{I}\right\|_{0,e}$, the continuous trace inequality (55), the bound on $\left\|\nabla{\varepsilon}_{\boldsymbol{b}}^{I}\right\|_{0,K}$ given by (63), and a similar bound for $\left\|\nabla{\varepsilon}_{\boldsymbol{u}}^{I}\right\|_{0,K}$ given by

[TABLE]

Combining the previous expressions, we have

[TABLE]

where we have used (58a) in the final step, but for simplicity in writing have not expanded $\left\|{\varepsilon}_{\boldsymbol{u}}^{I}\right\|_{0,K}$. Thus, from the three estimates (69), (71), (73) with (68), we have

[TABLE]

To complete the estimate of $\|{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$, we need to estimate $\|\operatorname{tr}{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$. First, by taking $\boldsymbol{G}=q\mathbb{I}$ in (57i) with $q\in\mathcal{P}_{k-1}(K)$, we get

[TABLE]

where $\mathbb{P}_{k-1}$ is the orthogonal $L^{2}$ projection into $\mathcal{P}_{k-1}(K)$. For a fixed $e_{K}\subset{\partial K}$, taking $\boldsymbol{\mu}=w\boldsymbol{n}_{e_{K}}$ with $w\in\mathcal{P}_{k}^{\perp}(K)$ in (57k) and using (65), we also get

[TABLE]

where $\zeta$ is a scalar function on $K$ defined by

[TABLE]

and $\mathbb{P}_{e_{K}}$ is the orthogonal $L^{2}$ projection into $\mathcal{P}_{k}(e_{K})$. We define $\Pi_{e_{K}}(f,g)\in\mathcal{P}_{k}(K)$ for $f\in L^{2}(K)$ and $g\in L^{2}(e_{K})$ as

[TABLE]

We refer to [11, Lemma 3.1] for well-posedness of this interpolation and optimal approximation property. By an argument similar to Lemma B.4, we have

[TABLE]

For simplicity, we will use $\Pi_{e_{K}}f$ if $g=f|_{e_{K}}$.

Note that $\operatorname{tr}{\varepsilon}_{\boldsymbol{L}}^{I}$ is an element-wise polynomial because $\operatorname{tr}\boldsymbol{L}=0$, so $\operatorname{tr}{\varepsilon}_{\boldsymbol{L}}^{I}=\Pi_{e_{K}}\operatorname{tr}{\varepsilon}_{\boldsymbol{L}}^{I}$. From this, the identities (75) and (76), and the above inequalities from scaling argument, we have

[TABLE]

The optimality of the orthogonal $L^{2}$ projection and the inverse trace inequality give

[TABLE]

and using this and the previous estimate, we can write

[TABLE]

*and here we used previous results on $\left\|{\varepsilon}_{\boldsymbol{L}}^{I}\boldsymbol{n}\right\|_{0,K}$ and $\left\|(\alpha_{1}+m){\varepsilon}_{\boldsymbol{u}}^{I}-\boldsymbol{g}\right\|_{0,K}$ in the final inequality. We completed the estimate of $\|{\varepsilon}_{\boldsymbol{L}}^{I}\|_{0,K}$. *

We now define the adjoint projection $\Pi^{*}(\boldsymbol{L}^{*},\boldsymbol{u}^{*},{p}^{*},\boldsymbol{J}^{*},\boldsymbol{b}^{*},r^{*},{\hat{\boldsymbol{u}}}^{*},(\hat{\boldsymbol{b}}^{*})^{t},\hat{r}^{*})$. As in the splitting of errors with $\Pi$, we define

[TABLE]

for an adjoint unknown $\sigma^{*}$. We first define $\Pi^{*}\boldsymbol{J}^{*}$, $\Pi^{*}{p}^{*}$, $\Pi^{*}{\hat{\boldsymbol{u}}}^{*}$, $\Pi^{*}(\hat{\boldsymbol{b}}^{*})^{t}$, $\Pi^{*}\hat{r}^{*}$ as $L^{2}$ projections into relevant polynomials spaces, and define $\Pi^{*}\boldsymbol{b}^{*}$, $\Pi^{*}r^{*}$ to satisfy

[TABLE]

We then choose $\Pi^{*}\boldsymbol{L}^{*}$, $\Pi^{*}\boldsymbol{u}^{*}$ to satisfy

[TABLE]

where

[TABLE]

Assuming that $(\boldsymbol{L}^{*},\boldsymbol{u}^{*},{p}^{*},\boldsymbol{J}^{*},\boldsymbol{b}^{*},r^{*},{\hat{\boldsymbol{u}}}^{*},(\hat{\boldsymbol{b}}^{*})^{t},\hat{r}^{*})$ are sufficiently regular, we can show that the interpolation $\boldsymbol{\Pi}^{*}$ is well-defined and provides optimal approximations. Due to the similarity between $(\Pi\boldsymbol{b},\Pi r)$ and $(\Pi^{*}\boldsymbol{b}^{*},\Pi^{*}r^{*})$, we can conclude

[TABLE]

It can also be shown that

[TABLE]

assuming $\operatorname{tr}\boldsymbol{L}^{*}=0$. The proofs are analogous to those for the $\boldsymbol{\Pi}$ projections, with the only differences resulting from the absence of $m$ from (78c). As a consequence, from the elliptic regularity assumption (35), we have

[TABLE]

and the implicit constant depends on $\boldsymbol{{d}}$, $\boldsymbol{w}$, $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ but not $h$.

Appendix C Well-posedness and regularity of the adjoint

problem

In this section, we discuss the well-posedness of the adjoint equation (34) and conditions under which the regularity result (35) holds. For the well-posedness of (34) we can adapt the approach in [32]. In particular, since (34b) and (34e) satisfy the following antisymmetry:

[TABLE]

we can invoke a similar argument as in [32] to conclude

[TABLE]

because $\boldsymbol{J}^{*}=\nabla\times\boldsymbol{b}^{*}$. We next assume that solutions of the Stokes and time-harmonic Maxwell equations (with [math] frequency) satisfy higher order regularities, i.e., when $\boldsymbol{w}=\boldsymbol{{d}}=0$ in (34), the solution $\boldsymbol{u}^{*}$, ${p}^{*}$, $\boldsymbol{b}^{*}$ satisfy $\left\|\boldsymbol{u}^{*}\right\|_{2}+\left\|{p}^{*}\right\|_{1}\lesssim\left\|\boldsymbol{\sigma}\right\|_{0}$ and $\left\|\boldsymbol{b}^{*}\right\|_{2}\lesssim\left\|\boldsymbol{\theta}\right\|_{0}$. Sufficient conditions for these assumptions are known but the details are beyond the scope of this paper, so we refer the readers to [21, 31, 38]. We now assume that $(\boldsymbol{u}^{*},\boldsymbol{L}^{*},{p}^{*},\boldsymbol{b}^{*},\boldsymbol{J}^{*},r^{*})$ is the solution of (34). From (34) and the regularity assumptions of the Stokes and Maxwell equations, we have

[TABLE]

The regularity of $\left\|r^{*}\right\|_{1}$ is already given in (81) and the regularity of $\left\|\boldsymbol{L}^{*}\right\|_{0}$ and $\left\|\boldsymbol{J}^{*}\right\|_{0}$ is obvious from the identities $\boldsymbol{L}^{*}=\nabla\boldsymbol{u}^{*}$ and $\boldsymbol{J}^{*}=\nabla\times\boldsymbol{b}^{*}$.

Appendix D Divergence-free post-processing of $\boldsymbol{u}_{h}$ and $\boldsymbol{b}_{h}$

In this appendix we adapt a postprocessing procedure in [14] to enforce the pointwise solenoidal constraint on $\boldsymbol{u}$ and $\boldsymbol{b}$. For simplicity the exposition is done only for $d=3$. For completeness, error estimates of the post-processed solutions will also be derived. To begin, let $\lambda_{0}$, $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ be the barycentric coordinates of a tetrahedron $K$, and $B_{K}$ be a symmetric matrix-valued bubble function defined by

[TABLE]

with the index $i$ counted modulo 4. Let $\mathcal{N}_{k}$ be the space of ${\mathbb{R}}^{3}$-valued polynomials

[TABLE]

where $\mathcal{N}_{k}^{\prime}(K)$ is the space of homogeneous ${\mathbb{R}}^{3}$-valued polynomials $\boldsymbol{v}$ of degree $k$ such that $\boldsymbol{v}\cdot\boldsymbol{x}=0$. We define $\mathcal{S}_{k}(K)$ by

[TABLE]

We recall an alternative characterization of the interpolation of the Brezzi–Douglas–Marini (BDM) element [14, Proposition A.1]:

[TABLE]

Our post-processed solution $\bar{\boldsymbol{u}}_{h}|_{K}\in\mathcal{P}_{k}(K)$ is defined as

[TABLE]

where $\boldsymbol{L}_{h}^{skw}=\left(\boldsymbol{L}_{h,32}-\boldsymbol{L}_{h,23},\boldsymbol{L}_{h,13}-\boldsymbol{L}_{h,31},\boldsymbol{L}_{h,21}-\boldsymbol{L}_{h,12}\right)$. In two dimensions, i.e. $d=2$, $\boldsymbol{n}\times\boldsymbol{a}$, $\nabla\times\boldsymbol{u}$, and $B_{K}$ are replaced by $n_{1}a_{2}-n_{2}a_{1}$, $\partial_{1}u_{2}-\partial_{2}u_{1}$, and the standard bubble function on $K$, respectively.

The post-processed solution $\bar{\boldsymbol{u}}_{h}$ is in the BDM space, so its divergence is well-defined. Further, it is divergence-free because, for any $q\in\mathcal{P}_{k-1}({\Omega_{h}})$,

[TABLE]

where the last equality is due to (8c). For error analysis, it is enough to estimate $\delta_{\boldsymbol{u}}:=\bar{\boldsymbol{u}}_{h}-\Pi^{BDM}\boldsymbol{u}$. Using $\nabla\times\boldsymbol{u}=\boldsymbol{L}^{skw}$ together with (82) and (83) gives

[TABLE]

Using a scaling argument yields

[TABLE]

On the other hand, by the triangle inequality we have

[TABLE]

Thus, the convergence rates of $\left\|\bar{\boldsymbol{u}}_{h}-\boldsymbol{u}\right\|_{0}$ and $\left\|\boldsymbol{u}_{h}-\boldsymbol{u}\right\|_{0}$ are the same.

To post-process $\boldsymbol{b}_{h}$, we define $\bar{\boldsymbol{b}}_{h}|_{K}\in\mathcal{P}_{k}(K)$ as

[TABLE]

Then the divergence of $\bar{\boldsymbol{b}}_{h}$ is well-defined and $\nabla\cdot\bar{\boldsymbol{b}}_{h}=0$ by (8f). By a completely analogous argument as above and the fact $\boldsymbol{J}=\nabla\times\boldsymbol{b}$, we can conclude

[TABLE]

hence the convergence rates of $\left\|\boldsymbol{b}-\bar{\boldsymbol{b}}_{h}\right\|_{0}$ and $\left\|\boldsymbol{b}-\boldsymbol{b}_{h}\right\|_{0}$ are the same.