Error estimates for structurepreserving discretization of the incompressible MHD system ^{1}^{1}1This material is based upon work supported in part by the US Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DESC0014400 and by Beijing International Center for Mathematical Research of Peking University, China.
Abstract.
In this paper, we carry out the error analysis for the structurepreserving discretization of the incompressible MHD system. This system, as a coupled system of NavierStokes equations and Maxwell’s equations, is nonlinear. We use its energy estimate and the underlying physical structure to facilitate the error analysis. Under certain CFL conditions, we prove the optimal order of convergence. To support the theoretical results, we also present numerical tests.
28
1. Introduction
An incompressible magnetohydrodynamic (MHD) system is a coupled partial differential equation system resulting from the incompressible NavierStokes equations and the (reduced) Maxwell’s equations. Assuming is a simply connected openbounded domain with Lipschitz boundary, the model problem we consider is
(1.1) 
The coefficients in this system are the Reynolds number , the magnetic Reynolds number , the coupling number , the relative electric conductivity , and the relative magnetic permeability . The initial conditions for this set of equations are
and the boundary conditions are
As discussed in the literature, the variables , and , once known, uniquely determine and . There are many different numerical methods to discretize MHD. We now briefly examine some existing literature on some of the numerical methods and their error analysis for two types of MHD systems: the stationary MHD system [12, 28] and the evolutionary MHD system [21, 26, 14].
For the stationary MHD system, Gunzburger, Meir and Peterson [12] propose a formulation with finite element discretization for the magnetic field, and analyze its wellposedness and convergent behavior. Schötzau [28], who also works on the stationary MHD system, proposes a new formulation with discretization for the magnetic field, and proves its wellposedness and the optimal order of convergence. There are also many other methods for stationary problems, for example, [10, 11, 29].
For the evolutionary MHD system, Prohl [26] studies the coupled and decoupled schemes based on conforming discretization of the magnetic field. He proves that the discrete solution converges to the weak solution under a strong CourantFriedrichsLewy (CFL) condition; that is, ( stands for the time step size, and for the mesh size). And He [14] studies the MHD system on a regular domain with conforming discretization of the magnetic field. He proves an unconditional optimal order of convergence.
In this paper, we study the convergence property of a structurepreserving discretization presented in [21]. This method is based on the mixed formulation [3], which comes from the idea of FEEC (finite element exterior calculus) [1, 2] and DEC (discrete exterior calculus) [5]. and conforming finite element discretization are used for the electric field and the magnetic field respectively. The advantage of this approach is that the important Gauss’s law for magnetic field is preserved exactly on the discrete level. Moreover, the incompressible MHD system we focus on is a timedependent nonlinear problem. Therefore, to conduct the error analysis, we work on an evolutionary nonlinear saddle point problem. Before approaching the detailed analysis, we briefly review the existing literatures for the error estimates of the (evolutionary) saddle point problems and nonlinear problems.
Abstract error estimates exist for standard (linear, nonevolutionary) saddle point systems [3]. Optimal order of convergence is ensured by the wellposedness of the discretization system and the approximation property of the finite element space. For the evolutionary saddle point problem, Boffi and Gastaldi [4] build a general framework for the semidiscretization of the evolutionary (linear) saddle point problem and provide sufficient conditions for a good approximation in the natural functional spaces.
For nonlinear saddle point problems, no abstract error estimate framework can be found in the literature. But various techniques have been developed for specific problems. For example, Temann [30] discusses the theory and numerical methods for NS equations. Heywood and Rannacher [17, 18, 19, 20] discuss the stability and error estimates of both semidiscretization and full discretization schemes for the NS systems. He [13] study linearized implicitexplicit schemes for this model.
General error estimates exist for nonlinear parabolic and elliptic problems. Thomée et al. [27, 24, 22, 23, 9, 31] investigate the error estimates of nonlinear parabolic problems intensively. Xu [32] uses the priori estimate to derive the error estimates of a general nonlinear elliptic problem. Brezzi, Rappaz and Raviart [6, 7, 8] build an abstract theory for finite element approximation of nonlinear problems.
Due to the nonlinearity and the loss of coercivity of the MHD model, the error estimate becomes difficult. To estimate the error of nonlinear problems, we usually need to prove that the norm (or a stronger norm) of the numerical solution is bounded. Generally, there are two ways to obtain this bound, one is using the mathematical induction method [13, 15], the other is introducing a semidiscrete problem [B.Li;W.Sun2013a, B.Li;W.Sun2013b]. Moreover, due to the loss of coercivity, we cannot use Cea’s Lemma to derive the error estimates directly.
In our analysis, we take advantage of the energy estimate of the structurepreserving discretization instead of estimating the norm of the numerical solution. We prove the unconditional error estimates for the velocity , the magnetic field , and the volume current density . And under certain constrains on the timestep size, we derive error estimates for the electric field and the pressure . Numerical tests support the theoretical results.
2. Magnetohydrodynamics model
In this section, we introduce some notation which follows mostly [21]. We first define the usual inner product
and the norm
For the sake of simplicity, we write both and as .
Given a linear operator , we define
and
where is the trace operator defined by
We define
When , we typically write instead of , and instead of . In the analysis, we also use the spaces and with norms
Another useful norm in the analysis is the discrete norm
where is the time step size (as is the in the following context). For example, stands for the discrete norm and for the discrete norm.
Next, we introduce some useful function spaces in the discretization.
We use (, , or ) to denote the dual space of , and to denote the finite element space of . The divergencefree subspace of is defined as
For the sake of simplicity, we assume that in the analysis. The Hilbert spaces and are equipped with norms and , which are defined as
Here,
We also use Sobolev space and , which are defined as
The corresponding norms are denoted by and , which are defined as
To facilitate the analysis, we also introduce a trilinear form of , namely,
(2.1) 
Based on the above notations, the variational formulation for system (1.1) is: find and such that for any and ,
(2.2) 
where , and .
3. Finite element discretization
In this section, we briefly go over the finite element discretization of (1.1). For the temporal discretization, we use the backward Euler method. For the spacial discretization, we recall the formulation of both nonlinear and linearized discretization here. These discretization formulations are reasonable in the sense that they inherit the energy estimate from the continuous level. At the end of this section, we go over their energy estimates.
Algorithm 3.1.
Find and such that for any and ,
(3.1) 
where , , and .
The above formulation uses linearization as a discretization scheme. In fact, we can discretize the nonlinear system directly and solve the nonlinear equation by Picard or Newton iteration.
Algorithm 3.2.
Find and such that for any and ,
(3.2) 
where .
As mentioned before, these above formulations admit desirable energy estimates. For the sake of completeness, we cite some of these estimates, which are established in [21].
Theorem 3.1 (Energy estimates).
Here, and depend on the initial data. A similar energy estimate holds for (3.1).
4. Error estimates
Before starting the detailed analysis, first we recall Gronwall’s inequality [20], which is an important tool in our analysis.
Theorem 4.1 (Gronwall’s inequality [20]).
Let , , and , , , , for integers , be nonnegative numbers such that
Suppose that (for all ), and set , then
We choose to be the th order polynomial space, the th order RaviartThomas elements, the th order Nédélec element, and the th order polynomial space.
Define , , and
(4.1)  
(4.2) 
Therefore, we can write the MHD system (2.2) in the form of a saddle point problem. That is, find such that for any ,
(4.3) 
where , and . Additionally, we can write (3.2) as: find such that for any ,
(4.4) 
Before giving the detailed error estimates, we define the projections of first. Assume that is the Stokes projection of . That is, for any given , we define such that satisfies
We choose as the projection of , and as the canonical interpolation of .
Define
By the definitions of and , we have
Similarly, we define
It follows that . Moreover, we define
For simplicity, we denote
where Noticing that
we can rewrite the error equation as
(4.5) 
where
4.1. Main results
We summarize main results of this paper for error estimates of (3.2) in the following theorem.
Theorem 4.2.
For any fixed time step such that , if is the solution to (4.3), and is the solution to (3.2), the following estimates hold:

There exists a constant , which only depends on , , and the computation domain, such that
(4.6) when the time step size is sufficiently small. And there also holds
(4.7) 
There exists a constant , which only depends on , , , and the computation domain, such that
(4.8) 
There exists a constant , which only depends on , , , and the computation domain, such that
(4.9)
By similar arguments similar to Theorem 4.2, we can get the error estimates of the Picard linearization scheme (3.1) as follows.
Theorem 4.3.
For any fixed time step such that , we have the following error estimates of (3.1):

There exists a constant , only depending on the exact solution, such that
when the time step size is sufficiently small. And there also holds

There exists a constant only depending on exact solution such that
when the time step size is sufficiently small.

There exists a constant only depending on exact solution such that
when the time step size is sufficiently small.
As the proof of the above theorem is similar to that of theorem 4.2, we omitted its details in this paper.
4.2. Proof of Theorem 4.2
The basic idea of our proof is to verify that satisfies the Gårding condition, and the truncation error is bounded by some norm of . Then the conclusion follows by Gronwall’s inequality.
First of all, we prove that satisfies the Gårding condition.
Lemma 4.1.
The sum of bilinear form satisfies the Gårding condition. That is, for any , there exists , such that
where and are positive constants that only depend on , , and the computation domain.
Proof.
By definition, we know that . Since , we choose . Noticing that , we get