Optimal Convergence Analysis of Two-Level Nonconforming Finite Element Iterative Methods for 2D/3D MHD Equations

Several two-level iterative methods based on nonconforming finite element methods are applied for solving numerically the 2D/3D stationary incompressible MHD equations under different uniqueness conditions. These two-level algorithms are motivated by applying the m iterations on a coarse grid and correction once on a fine grid. A one-level Oseen iterative method on a fine mesh is further studied under a weak uniqueness condition. Moreover, the stability and error estimate are rigorously carried out, which prove that the proposed methods are stable and effective. Finally, some numerical examples corroborate the effectiveness of our theoretical analysis and the proposed methods.


Introduction
Magnetohydrodynamics (MHD) describes the interaction of electrically conductive fluids with electromagnetic fields. The governing equations describing the MHD system is a strong coupling nonlinear system coupled with the Navier-Stokes equations and Maxwell equations. Magnetohydrodynamics has become widespread in such areas of astrophysics, controlled thermonuclear reactions and industry. For the study of the MHD problem, a large amount of research and analysis have been carried out in recent decades. The well-posedness of weak form solution of MHD equations can be guaranteed in [1,2]. With regard to theoretical analysis, the regularity, long-time behaviors of solution of MHD problems and the error estimation of FEM are studied in [3,4]. We can refer to [5][6][7][8][9] and their references for many Galerkin finite element methods (FEM) analysis and studied on the MHD system. For the MHD problem, a series of one-level iteration methods and their error estimation are studied in [10], and some coupled type iteration methods are designed and discussed in [11] on a general Lipschitz domain.
The present paper mainly focuses on the study of nonconforming finite element. Low-order nonconforming FEM has the advantages over the conforming FEM in terms of simplicity and small support sets of basis functions. For the Stokes and Navier-Stokes equation, the nonconforming FEM are studied in [12,13], in which the discretization of velocity space uses the nonconforming element, and the discretization of pressure space uses the piecewise constant element. In addition, a nonconforming FEM are also proposed in [14][15][16], which differs from the discrete pressure space using a piecewise linear element. Due to the limitation of inf-sup condition, it has the advantage of simple structure and has been well applied in solving various problems. For instance, an Oseen iterative algorithm for the conduction-convection equations with nonconforming FEM is carefully studied in [17]. In addition, the low-order nonconforming FEM is used to solve 3D MHD system, then they make deep and systematic analysis and research in [18,19].
In order to raise the efficiency of computation, much work has sought to answer two-level methods for solving the nonlinear problems with conforming FEM. For the Navier-Stokes equations, some two-grid schemes and their error estimation are presented in [20][21][22][23][24]. Regarding the two-level methods approach to discussing the MHD problems, we can refer to the literature [25][26][27]. Furthermore, the MHD problem is studied well enough by a two-level Newton method at a small magnetic Reynolds number and a twolevel method under the hypotheses of a small data in [28,29]. In addition, some two-level iterative methods are used to solve the MHD problem, then rigorous systematic analysis and numerical tests are carried out in [30].
In our previous work, three linearized and one level nonconforming discretization for the 2D/3D MHD equations were proposed in [31]. In order to solve it more efficiently and to make our work more completely, the algorithms we study are looking to link two-level iterative methods with the nonconforming FEM to solve the MHD problem. The study is based on approximating velocity space by the nonconforming element, and the conforming linear elements are used in the discretization of magnetic field and pressure space. We analyze the error estimates of three classical one-level iterative methods under different uniqueness conditions. Then, more comprehensive and diverse two-level iterative methods based on the iterative solution firstly calculated by Stokes, Newton and Oseen iterative methods on a coarse grid and then the correction solution calculated by Stokes, Newton and Oseen corrections on a fine grid are proposed. Moreover, we perform a systematic and in-depth analysis of our proposed methods under the different strong uniqueness conditions. Finally, several numerical texts are executed and the accuracy of the proposed methods are proved.
The following describes the components of this article. The mathematical setting of the MHD equations is introduced and the nonconforming FEM is proposed in Sections 2 and 3. In Section 4, three classical one-level iterative methods and error analysis results are given. In Section 5, we put forward some two-level algorithms and deduce more comprehensively theoretical analysis. In Section 6, several numerical simulations are tested to verify the accuracy of the previous results. In the last part, the summary and prospects of the paper are given.
In this next paper, C represents a real constant, which represents different values in different cases and is independent of the coefficients of the system equations, the grid sizes h and H. Notations without special interpretation are used for their usual meaning.

Preliminaries
Let Ω be a bounded convex region in R d (d = 2, 3) with boundary ∂Ω, to study the 2D/3D stationary incompressible MHD problem, which is modeled as listed below: find the velocity u, magnetic field b and pressure p such that where f ∈ H −1 (Ω) d , g ∈ L 2 (Ω) d denote prescribed body terms, and n represents outward unit normal vector. The nonlinear system includes three different coefficients: hydrodynamic Reynolds number R e , magnetic Reynolds number R m , and coupling coefficient S c . It should still be noted out that the actual physically meaningful induction equation corresponds only to the special case g = 0.
In addition, we have the following properties from [10,11,18]: for all u, κ, e ∈ H 1 0 (Ω), and b, D, E ∈ H 1 n (Ω), where λ 1 is the constant from Finally, to analyze the error estimates in the following sections, we can obtain the important theorem as follows [2,11,18]: is well-posed and the unique solution (u, b) satisfies: To obtain the H 2 stability of the solution to (2), we give the following assumption where the domain Ω satisfies regular properties as follows. Assumption 1. First, the steady Stokes problem is introduced as follows: −∆κ + ∇r = f , ∇ · κ = 0, in Ω, κ| ∂Ω = 0, and the unique solution (κ, r) satisfies Then, we introduce Maxwell's equations Similarly, the unique solution D satisfies

Remark 1.
If Ω is a convex polygon or polyhedron, or if ∂Ω is of C 2 , the conclusion is that the assumption is tenable [32,33].

Nonconforming Discretization
Here, we consider the regular triangulation T µ = {K} that partitions the domain Ω into triangles or quadrangles. Here, we defined the positive parameter µ = max K∈T µ {µ K : µ K = diam(K)}, the boundary edge Γ j = ∂K j ∩ ∂Ω and the interior boundary Γ jk = Γ kj = ∂K j ∩ ∂K k . Denote the centers of Γ j and Γ jk by υ j and υ jk , respectively. This is followed by the discrete spaces as shown below, in which the discretization of velocity space is approximated by a nonconforming element, and the discretization of pressure and magnetic field space is approximated by a conforming linear element. X 1µ = e µ ∈ L 2 (Ω) d : e µ K ∈ P 1 (K) d , ∀K ∈ T µ ; e µ (υ jk ) = e µ (υ kj ), e µ (υ j ) = 0, ∀j, k , where P 1 (K) stands for the continuous piecewise polynomial space. In fact, the finite element pair(X 1µ , M µ ) we studied in this paper satisfies the inf-sup condition, which has been rigorously proven in [34]. The finite element spaces X 1µ , X 2µ and M µ satisfy the interpolation theory: for any [13,19], To give the discrete variational form of (2), we introduce the following compatibility conditions [14,15]: where [e µ ] = e µ | Γ jk − e µ | Γ kj , representing e µ through an interface Γ jk .
Combining the above spaces and compatibility conditions, the discrete variational form of (2) is recast: For all κ µ ∈ X 1µ , D µ ∈ X 2µ , we define

Lemma 2.
Here are two important properties of the trilinear forms [27]: The important lemma followed there, which has been proved in [11,14,37].

Lemma 3.
There exists a constant β > 0 such that Combining with the above conclusions, we can obtain the existence and uniqueness results of (6) from [18].

Theorem 3. Suppose that
is well-posed and its unique solution satisfies For the convenience of the interval, it is expected to be easier to estimate the error of the solutions (u µ , b µ ), p µ ∈ X 1µ × X 2µ × M µ . Then, the projection operators ((S µ , Q µ ), P µ ) are given as follows: As the finite element spaces X 1µ , X 2µ and M µ satisfy property (5), the above projection operators satisfy the following conclusion [16,17]: Next, we draw two important conclusions from this section. One is the error estimate of H 1 -norm of finite element solution without proof, which has been given in [19]. Another important result is the error estimate of L 2 -norm of finite element solution obtained by duality theory [38], which has been proved and can be found in [19].

Iterative Methods
Recently, for the MHD problem, three iteration methods under different uniqueness conditions are presented in [10] and three iteration methods based on nonconforming FEM are designed in [31] on a Lipschitz domain. Therefore, we further give the important conclusion of this section, three iterative methods on convex region and their error estimates. Three iterative methods appear as follows: Method 2: the Newton iteration method.
The initial value ((u 0 Next, we further find the stability of these iterative methods, which can be proved by using similar methods in [10,31]. and (u µ − u n µ , b µ − b n µ ) and p µ − p n µ have the following bounds: for all n ≥ 0.
and (u µ − u n µ , b µ − b n µ ) and p µ − p n µ have the following bounds: for all n ≥ 0.
According to Theorems 4-7, we draw the important conclusions of this section as follows.
for all n ≥ 0.
For further clarification about all the methods, we describe them as follows.

Remark 2.
From the above discussion, we can see clearly that method 1 is the simplest. In the condition 0 < σ < 5 11 , method 2 is stable and has exponential convergence rate for n, which is faster than methods 1 and 3. In addition, if 0 < σ < 1, method 3 is stable and convergent and has the widest application scope.

Two-Level Iterative Methods
Some two-level methods for solving the MHD problem with different Reynolds numbers are proposed in [30]. In this section, based on nonconforming FEM, to further study more effective two-level methods, several two-level methods are proposed in different conditions as follows.
If 0 < σ < 2 5 , we know that three iterative methods are all stable and convergent. We firstly propose nine two-level methods to derive the iterative solution ((u m H , b m H ), p m H ) using methods 1-3 on coarse mesh T H , and then to solve the correction solution ((u mh , b mh ), p mh ) using corrections 1-3 on fine mesh T h . Here, H and h are set as positive numbers tending to zero (0 < h ≤ H). All of the proposed methods can be shown below: Step 1. Solve the MHD equations on the coarse grid, Step 2. One step correction on the fine grid: find ((u mh , b mh ), p mh ) ∈ X h × M h satisfy: Correction 1: Stokes correction.
In what follows, we set out to give the following results for the above two-level iterative methods by introducing projection operators. where Then, we can derive the error equation by subtract (33) from (36), Subtracting (37) from (12), we obtain (38), and applying the property (13), we deduce Then, using (13), (39) and the triangle inequality, we can hit bottom Next, using the property (7) and (9) in (38) yields Finally, we can obtain the conclusion by applying (13), (40) and the triangle inequality, Thus, the proof is done.
for all m ≥ 0.
(42) (42), and making the use of (13), we obtain Then, using (13), (43) and the triangle inequality, we can come to the conclusion that Next, we can deduce the conclusion by using properties (7) and (9) in (42), Finally, we can obtain the conclusion by applying (13), (44) and the triangle inequality, Thus, the proof is done.  ((u mh , b mh ), p mh ) calculated by correction 3. Then, (u − u mh , b − b mh ) and p − p mh satisfy: for all m ≥ 0.
(46) (46), and making the use of (13), we can derive Then, using (13), (47) and the triangle inequality, we can come to the conclusion that Next, by using (7) and (9) in (46), we deduce that Finally, we can obtain the conclusion by applying (13), (48) and the triangle inequality, Thus, the proof is done.

Two-Level Iterative Method with
If 2 5 < σ < 5 11 , we know that iterative methods 2 and 3 are stable and convergent, while method 1 is not convergent. We propose six two-level iterative methods to solve the iterative solution ((u m H , b m H ), p m H ) using methods 2 and 3 on coarse mesh T H , and then to solve the correction solution ((u mh , b mh ), p mh ) using corrections 1-3 on fine mesh T h . All of the two-level methods are shown below: Step 1. Solve the MHD equations on the coarse grid, ((u m H , b m H ), p m H ) ∈ X H × M H provided by methods 2 and 3, respectively.
Step 2. One step correction on the fine grid, ((u mh , b mh ), p mh ) ∈ X h × M h provided by corrections 1-3, respectively.
Then, we further give the following important results about the estimation of the above two-level methods.
and ((u mh , b mh ), p mh ) calculated by correction 2, (u − u mh , b − b mh ) and p − p mh satisfy: and ((u mh , b mh ), p mh ) calculated by correction 3, (u − u mh , b − b mh ) and p − p mh satisfy: for all m ≥ 0.

Two-Level Iterative Method with
we can see that method 3 is the only stable and convergent method. We propose three two-level iterative methods to solve the iterative solution ((u m H , b m H ), p m H ) using method 3 on coarse mesh T H , and then to solve the correction solution ((u mh , b mh ), p mh ) using corrections 1-3 on fine mesh T h . Three two-level methods are shown below: Step 1. Solve the MHD equations on the coarse grid, ((u m H , b m H ), p m H ) ∈ X H × M H provided by method 3.
Step 2. One step correction on the fine grid, ((u mh , b mh ), p mh ) ∈ X h × M h provided by corrections 1-3, respectively.
Then, we further obtain the following theoretical results about the estimation of the above two-level methods.
Theorem 13. Under the conditions of Theorem 4 and the condition 5 and ((u mh , b mh ), p mh ) calculated by correction 2, (u − u mh , b − b mh ) and p − p mh satisfy: and ((u mh , b mh ), p mh ) calculated by correction 3, (u − u mh , b − b mh ) and p − p mh satisfy: for all m ≥ 0.

One-Level Oseen Iterative
Method with 1 − ( F * F 0 ) 1 2 < σ < 1 It can be analyzed from Theorem 8 that method 3 is the only stable and convergence method with the weaker condition 1 − ( F * F 0 ) 1 2 < σ < 1. We present a one-level method to solve the iterative solution ((u mh , b mh ), p mh ) on coarse mesh T h , as shown below.
Step 1. Solve the MHD equations on the fine grid, ((u mh , b mh ), p mh ) ∈ X h × M h obtained by method 3.
and (u h − u mh , b h − b mh ) and p h − p mh satisfy: for all n ≥ 0.
Combining Theorems 4 and 14, we further derive the final conclusion of this section for all n ≥ 0.

Remark 4.
In terms of the convergence rate, we know that two-level iterative methods that combine method i (i = 1, 2, 3) with correction j (j = 1, 3) are linear convergent. Moreover, two-level methods that combine method i (i = 1, 2, 3) with correction 2 have an exponential convergence rate. Therefore, a two-level method that combined method 2 with correction 2 has a faster convergence speed under the unique condition 0 < σ < 5 11 . In case of 5 11 , a two-level iterative method that combined method 3 with correction 2 converges the fastest.

Numerical Experiments
In this part, three numerical tests are rendered to substantiate the good performance of our proposed methods for the MHD equations. Taking a fluid problem with smooth true solution and Hartmann flow as examples, the optimal convergence rate and computational cost of the proposed scheme are tested. The last one of the driven cavity flow shows good simulated fluid motion results. Moreover, we use low order nonconforming finite element pair P 1 nc-P 1 b-P 1 for the velocity, magnetic field and pressure. Throughout this section, we denote M1, M2, M3, C1, C2, and C3 as the abbreviations of methods 1-3 and corrections 1-3.
where α > 0 satisfies the uniqueness conditions of two-level methods.
Firstly, set R e = R m = S c = 1. The error values and order of convergence data calculated by our presented two-level methods are shown in Table 1. The CPU times are showed in Figure 1. From this table, we can see that the orders of u − u n h 1 , b − b n h 1 and p − p n h 0 are almost all equal to one. It is clear that the two-level algorithms can guarantee the stability of the stabilization method for fluid problems with smooth true solutions. We can see from Table 1 and Figure 1 that Mi + C2 (i = 1, 2, 3) has the best accuracy, while Mi + C1 (i = 1, 2, 3) spends the least computational time with the smaller mesh division. Since the trilinear terms of Mi + C1 (i = 1, 2, 3) are the easiest, the trilinear terms of Mi + C2 (i = 1, 2, 3) are most complicated in our presented two-level algorithms.  Then, numerical results for calculating velocity using the nonconforming element P 1 nc are given in Table 1, and numerical results for calculating velocity using the piecewise linear element P 1 are shown in Table 2. By comparing and analyzing the data in Tables 1 and 2, we can find that the errors calculated by our proposed two-level methods are smaller.

The Hartmann Flow
In this case, we will test a classical MHD problem, the 2D Hartmann flow. It is affected by a steady unidirectional flow in the channel Ω = [0, 10] × [−1, 1]. Set Ha = √ R e R m S c in 2D Hartmann flow. Meanwhile, the transverse magnetic field b D = (0, 1) is inflicted on the boundaries of this system. The analytic solutions are given by: Furthermore, the boundary conditions are given indicated below: where p d (x, y) = p(x, y), p 0 is zero, and I is the identity matrix.
Taking G = 0.1 and considering two different schemes: (1) Ha = 1 (R e = 1, R m = 1, S c = 1), (2) Ha = 10 (R e = 10, R m = 1, S c = 10). The first components of the analytic solutions u(y), B(y) and the numerical ones u(y k ), B(y k )(y k = −1 + 0.1k, k = 0, · · · , 20) provided by Mi + Cj (i, j = 1, 2, 3) are showed in Figures 2-4. Since some two-level methods do not converge at high Hartmann numbers, we present only convergent images.    Then, to better illustrate the precision of our proposed algorithms, the error values and order of convergence data calculated by our presented two-level algorithms are shown in Table 3 when scheme (1) is selected. Here, we only show the discretization errors of velocity, magnetic field and pressure of one method, since the numerical results of other methods are roughly the same in numerical values. We can seen from the table that the convergence orders of velocity and magnetic field in H 1 -norm reach the first order, and the pressure in L 2 -norm achieves a better approximate convergence rate higher than the first order. The results show that the two-level algorithms can guarantee the accuracy of numerical experiment for the Hartmann flow problem. Next, we test the convergence of Mi + Cj (i, j = 1, 2, 3). Setting G = 0.1 and choosing three different schemes: (1) Ha = 1 (R e = 1, R m = 1, S c = 1), (2) Ha = 10 (R e = 10, R m = 1, S c = 10), (3) Ha = 10 √ 15 (R e = 150, R m = 1, S c = 10). All the nine proposed methods are available for scheme (1) in Figure 5, which means that all the proposed methods are convergent for 0 < σ < 2 5 . However, Figure 6 shows M1 + Cj (j = 1, 2, 3) are not convergent for 2 5 < σ < 5 11 because the parameter selection is not in 2 5 < σ < 5 11 . Figure  7 shows If Ha = 10 √ 15 (R e = 150, R m = 1, S c = 10), only M3 + Cj (j = 1, 2, 3) is convergent while Mi + Cj (i = 2, 3, j = 1, 2, 3) is not convergent, since the uniqueness condition is not met.
We can also see from those figures that M1 + Cj (j = 1, 2, 3) only applies to the case of a small Hartmann number. In some cases with large Hartmann numbers, M3 + Cj (j = 1, 2, 3) can be chosen. Figure 5 shows that M2 + Cj (j = 1, 2, 3) converges faster than Mi+Cj (i = 1, 2, j = 1, 2, 3) since M2 has exponential convergence. All in all, these numerical experimental results demonstrate the effectiveness of our theoretical analysis and the proposed methods.

Driven Cavity Flow
In the last case, the numerical simulation of a classical fluid problem with driven cavity flow is showed. We consider the flow in the 2D domain Ω = [−1, 1] × [−1, 1] with Γ D = ∂Ω, and there is no analytical solution. Let the external force terms f and g be zero; then, their boundary conditions are defined as follows: From the above two examples and theoretical analysis, M3 + Cj (j = 1, 2, 3) has a wider application. Therefore, we mainly simulate the effect of the driven cavity flow by M3 + C2 in this part. In Figure 8, the velocity streamline for three different parameters R e = 1, 100 and 1000 with R m = 1, S c = 1 are showed. As R e enhances, the number of vortices produced by velocity streamlines increases to three. We use the same change pattern in Figure 9, and the streamlines of velocity for S c = 50, 500, 5000 with R e = 1, R m = 1, from which we can notice that the vortices produced by velocity streamlines also divide into three vortices and move upward.  Then, Figures 10-12 show the variation trend of numerical streamlines of the magnetic by M3 + C2 as parameters change. We vary R e = 10, 100, 1000 with R m = 5, S c = 1 in Figure 10. Analogously, we vary S c from 1 in (a), to 100 in (b), then to 5000 in (c) with R e = 1, R m = 10 in Figure 12, from which we can all observe that the streamlines have a tendency to change straight.  Figure 10. Numerical streamlines of the magnetic field drawn using data obtained from M3 + C2, wherein R e is set as 10, 100 and 1000 respectively.  Figure 11. Numerical streamlines of the magnetic field drawn using data obtained from M3 + C2, wherein R m is set as 1, 5 and 15 respectively. Conclusively, Figure 11 shows the trends of magnetic field streamlines for R m = 1, 5, 15 with R e = 10, S c = 1. As R m enhances, the shape of streamline changes from a straight line to a curve. The above phenomenon indicates an increase in curvature.  (c) R e = 1, R m = 10, S c = 5000 Figure 12. Numerical streamlines of the magnetic field drawn using data obtained from M3 + C2, wherein S c is set as 1, 100 and 5000 respectively.

Conclusions
Based on nonconforming FEM, several two-level methods for solving the stationary incompressible MHD equations have been presented under different unique conditions in this paper. Combining theoretical analysis with numerical experiments, the two-level method that combined method 2 with correction 2 has faster convergence speed and better calculation accuracy under the unique condition 0 < σ < 5 11 . In case of 5 11 < σ < 1 − ( F * F 0 ) 1 2 , the two-level method that combined method 3 with correction 2 has a good advantage. If 1 − ( F * F 0 ) 1 2 < σ < 1, one-level Oseen iterative method is a unique scheme.

Conflicts of Interest:
The authors declare no conflict of interest.