Two-Level Finite Element Iterative Algorithm Based on Stabilized Method for the Stationary Incompressible Magnetohydrodynamics

In this paper, based on the stabilization technique, the Oseen iterative method and the two-level finite element algorithm are combined to numerically solve the stationary incompressible magnetohydrodynamic (MHD) equations. For the low regularity of the magnetic field, when dealing with the magnetic field sub-problem, the Lagrange multiplier technique is used. The stabilized method is applied to approximate the flow field sub-problem to circumvent the inf-sup condition restrictions. One- and two-level stabilized finite element algorithms are presented, and their stability and convergence analysis is given. The two-level method uses the Oseen iteration to solve the nonlinear MHD equations on a coarse grid of size H, and then employs the linearized correction on a fine grid with grid size h. The error analysis shows that when the grid sizes satisfy h=O(H2), the two-level stabilization method has the same convergence order as the one-level one. However, the former saves more computational cost than the latter one. Finally, through some numerical experiments, it has been verified that our proposed method is effective. The two-level stabilized method takes less than half the time of the one-level one when using the second class Nédélec element to approximate magnetic field, and even takes almost a third of the computing time of the one-level one when adopting the first class Nédélec element.


Introduction
Consider the following stationary incompressible MHD in Ω, div u = 0, div b = 0, in Ω, where Ω ∈ R d (d = 2, 3) is a bounded Lipschitz domain. R e and R m are the hydrodynamic and magnetic Reynolds numbers, respectively. S c is the coupling number, and f and g are source terms with ∇ · g = 0. n is the unit outward normal vector on ∂Ω.
Incompressible MHD describes the dynamics of a viscous, incompressible, electrically conducting fluid under an external magnetic field. The MHD (1) is a coupled multi-physical system of the classical Navier-Stokes equations and Maxwell's equations. MHD modelling has a number of applications in physics and engineering technology, such as radio wave propagation in ionosphere in geophysics, MHD engine, control of MHD boundary layer and liquid-metal MHD electricity generation (see [1]). Since MHD equations are strongly nonlinear and have many physical quantities, it is needed to find effective numerical methods to solve them.
For the MHD modelling (1) without the Lagrange multiplier r term, the early study of the exact penalty regularization finite element method on a convex domain is carried out in [2]. Based on this format, the nonconforming mixed finite element methods [3], the Stokes, Newton and Oseen finite element iterative methods [4,5], the penalty based finite element iterative methods [6], and the generalized Arrow-Hurwicz iterative methods [7] are investigated. In view of multi-physical coupling and nonlinearity of system (1), twolevel method and finite element iterative algorithms are combined by [8][9][10][11][12] to reduce the computing cost, and local and parallel finite element algorithms based on some iterations are proposed in [13][14][15][16]. On the other hand, a number of effective solvers based on the finite element methods are presented in [17][18][19]. To keep the physical property of the Gauss law of the magnetic field, the constrained transport divergence-free finite element method is designed in [20]. The coupled Stokes, Newton and Oseen-type iteration methods are studied and discussed for the (1) in [21] on a general Lipschitz domain. For the nonsmooth computational domain, the magnetic field belongs to a lower regularity space than H 1 (Ω), and the discrete finite element scheme with the Lagrange multiplier of (1) becomes a double-saddle points problem.
For the mixed finite element method, the component approximations must preserve the compatibility and satisfy the so-called inf-sup condition. It is well known that the lowest equal-order finite element pairs in engineering preferred do not satisfy the inf-sup condition. Numerical experiments show that the break of the inf-sup condition often leads to unphysical pressure oscillations. To avoid the instability problem and use the lowest equal-order elements, the popular stabilized methods based on local Gauss integrations are proposed and studied, for example, for the Stokes problem [22,23], the coupled Stokes-Darcy problem [24], the Stokes eigenvalue system [25], the Navier-Stokes equations [26][27][28][29] and the natural convection problem [30]. However, the stabilized finite element algorithm for MHD with respect to the Lagrange multiplier has not been reported.
In this paper, a two-level finite element iterative algorithm based on the stabilized method is proposed to numerically solve the stationary incompressible MHD equations. Compared to the existing literature, the stabilized scheme with the Lagrange multiplier proposed here have two main benefits. One is that the lowest equal-order finite element pairs can be used to approximate hydrodynamic subproblem, and the other is that our scheme preserve the physical property of Gauss law weakly for magnetic subproblem by adding the Lagrange multiplier. In the next section, the stabilized finite element discretization based on local Gauss integrations is designed and analyzed. To deal with the nonlinear term, the stabilized finite element method based on Oseen iteration is studied. The two-level stabilized finite element algorithm and its convergence are given in Section 3. In the last section, some numerical experiments are tested to support the theoretical analysis of our proposed method.

Stabilized Finite Element Discretization Based on Local Gauss Integrations
We will introduce some Sobolev spaces, and the norms of the product spaces: For all w, u, v ∈ X, b, c, d ∈ W, q ∈ Q, s ∈ S, let By the Lagrange multiplier technique, the variational form of system (1) is [31]: The compact form of (2) and (3) is read as The properties of the bilinear and trilinear forms from [32][33][34] are useful for our analysis. For all u, v ∈ X, b, c ∈ W, q ∈ Q, r ∈ S, there have C(w, d; u, b; u, b) = 0, whereN and λ 0 are positive constants that depend only on Ω. In the next content, we use C to represent a general positive constant independent of mesh sizes H and h. H and h(h H) are now two real positive parameters that tend to 0. T H is a uniformly regular partition of Ω into triangular(d = 2) or tetrahedral(d = 3) element K with diameters bounded by H, and T h is the fine mesh generated by a mesh refinement process to T H . Let T µ (µ = H, h) is a partition. P k (K) is the space of polynomials of degree k (positive integers) over K. P 1 element is utilized to approximate the velocity field, pressure and Lagrange multiplier, and two kinds of lowest order Nédélec elements are applied to approximate the magnetic field. The subspaces of X, W, Q, S are Here, Nédélec elements of the first family and the second one are as follows [35] N (1) where [P k (K)] d is the homogeneous polynomials of degree k.
W µ and S µ satisfy the discrete inf-sup condition [31] sup where the constantβ > 0 is independent of µ.
Denote P µ and R 0µ by the L 2 -orthogonal projectors Define the discrete Stokes operator by A 1µ = −P µ ∆ µ , in which ∆ µ is defined by (see [32,33] It is necessary to introduce some discrete estimates [33,34] The trilinear form C(·, ·, ·) has the properties [34]: for all w µ , It is apparent that the discrete inf-sup condition is not valid to the subspace X µ and Q µ . To meet the needs of this property, as in [22,26], a mixed stability term with the universal bilinear form is added: where for all p µ , q µ ∈ Q µ , K,i p µ q µ dξ means that makes use of an i-order(i = 1, 2) local Gauss integral to calculate it over the element K. Let Π µ : L 2 (Ω) → P 0 be a L 2 -projection with the properties as follows [22,36]: As a consequence, the local Gauss integral can be restated as: D µ (·, ·; ·, ·) satisfies the following important properties (see [22,26]): For all (u µ , The stabilized discrete scheme reads: where σ > 0 is an artificial viscosity parameter, (17) can be rewritten as: Then the bilinear form A µ satisfies the following coercive and continuous properties: where In order to derive error estimates, we introduce two projections. The Stokes projection is defined as follows [26,36] The Maxwell's projection is defined by [37]: By the property of Λ, it can be shown that Now, we will give the stability and error estimate for the problem (21).
The proof of Theorem 1 is shown in the section of Appendix A. In the following, the Oseen iteration is used to linearize the stabilized finite element discrete form (17). The stability and convergence is proven. The stabilized finite element algorithm based on Oseen iteration is stated as follows: Here, the initial value (u 0 Rewrite (31) and (32) in compact form as (33) and (34), we can access Using (20), we have (36) holds, it is sufficient to prove that it also holds for (33) and (34), using (8), we get and by applying (20) and (8), we derive that Thus, the proof of the first part of (36) has been finished. Applying (16) to (31), and using (10) in (32) with s = 0, we havê Combining with the stability of (u m µ , b m µ ) 1 , we can complete the proof of (36). Similar to (27), (37) can be derived. The proof ends.
Next, we will establish the upper bound of error ( Proof. Subtracting (35) from (21), there holds It follows from (20), (9) and Lemma 2 that Subtracting (31) from (16), we have Using the second equation of (17) minus equation (32) and choosing s = 0, we get Applying (16) and (10) to the above two equations, respectively, there holdŝ The result (39) can be obtained by (38) and Theorem 1. The proof ends.

Two-Level Stabilized Finite Element Algorithm
In this section, motivated by [38], two-level stabilized finite element algorithm for incompressible MHD equations is presented. The stability analysis and the optimal error estimation with respect to the mesh size H and h and the iterative step m are obtained.
To estimate the pressure, we rewrite (45) with s = 0 as . Applying (16), (10) and the standard technique to the above two equations, we can derive the second part of (44). We complete the proof.

Numerical Examples
In this section, some numerical experiments are shown to verify the correctness and effectiveness of the one-level stabilized finite element method and the two-level stabilized one. Here, the velocity, pressure and quasi-pressure are approximated by P 1 and the magnetic field by the first (or second) class Nédélec edge element. SFEM denotes by the stabilized finite element method (31)

Smooth solution in 2D:
Set Ω = [0, 1] 2 and R e = R m = S = 1, σ = 0.01. Given the source terms f, g such that the exact solution is

Tables 1 and 2 display the errors of SFEM and two-level SFEM for 2D MHD Equation (1).
It is shown that the corresponding errors are smaller and smaller along with the grid getting finer and finer, the convergence order is optimal. When h = O(H 2 ), the error accuracy of the two methods is almost the same. From CPU time, compared to SFEM, two-level SFEM save much computational cost.  The numerical results are listed in Tables 3 and 4 when the magnetic field is approximated by the second class Nédélec element. Clearly, the convergence order of || b − b h || 0 is one higher than that in Tables 1 and 2, which is consistent with the general theoretical analysis results of the Nédélec element.

Smooth solution in 3D:
Set Ω = [0, 1] 3 and R e = R m = S = 1, σ = 0.01. Given f, g such that the exact solution is: In Tables 5 and 6, the variable quantity (u, b, p, r) is approximated by P 1 , the first class Nédélec element, P 1 and P 1 , respectively. Tables 7 and 8 list the results when (u, b, p, r) is approximated by P 1 , the second class Nédélec element, P 1 and P 1 . It is observed that the numerical results agree well with the theoretical results of Theorems 1-3. On the other hand, from Tables 7 and 8, we find that SFEM (35) does not work when H = 1/16, however two-level SFEM (41) and (42) is valid in the current computing environment for our computer. On the other hand, the stability results of (43) are checked by Figure 1.
In the 2D case, there holds the regularity u ∈ H 1+λ (Ω), b ∈ H 2 3 (Ω) and p ∈ H λ (Ω). In Tables 9 and 10, (u, b, p, r) is approximated by P 1 , the second class Nédélec element, P 1 and P 1 . Because the regularity of velocity, magnetic field and pressure is low, the con-  Figure 2 we display the streamlines of the velocity field and magnetic field, and the contours of the pressure, which are consistent with the numerical results in the literature [40].

Conclusions
Based on the stabilization and the Lagrange multiplier techniques, the stabilized finite element algorithm is designed for the stationary incompressible MHD. The Lagrange multiplier technique idea helps us in dealing with the low regular magnetic field subproblem by H(curl; Ω)-conforming element. The stabilized one by using local Gauss integration allow us to adopt the lowest equal-order elements to approximate the flow field sub-problem. The stability and optimal convergence analysis are given. Furthermore, the two-level stabilized finite element algorithms are presented. In the first step we combine the stabilized finite element method with the Oseen iteration for the nonlinear MHD equations on a coarse grid. For the second step, we employ the linearized correction on a fine grid. We give the optimal error analysis, which shows that when the grid sizes satisfy h = O(H 2 ), the two-level stabilization method not only has the optimal convergence order, but also can save more computational cost than the one-level method. These theoretical analysis results have been verified by some numerical experiments.