A robust solver for the finite element approximation of stationary incompressible MHD equations in 3D ^{†}^{†}thanks:
Abstract
In this paper, we propose a robust solver for the finite element discrete problem of the stationary incompressible magnetohydrodynamic (MHD) equations in three dimensions. By the mixed finite element method, both the velocity and the pressure are approximated by conforming finite elements, while the magnetic field is approximated by conforming edge elements. An efficient preconditioner is proposed to accelerate the convergence of the GMRES method for solving the linearized MHD problem. We use three numerical experiments to demonstrate the effectiveness of the finite element method and the robustness of the discrete solver. The preconditioner contains the least undetermined parameters and is optimal with respect to the number of degrees of freedom. We also show the scalability of the solver for moderate physical parameters.
Key words. Incompressible magnetohydrodynamic equations, mixed finite element method, preconditioner, parallel computing.
1 Introduction
Magnetohydrodynamics (MHD) has broad applications in our real world. It describes the interaction between electrically conducting fluids and magnetic fields. It is used in industry to heat, pump, stir, and levitate liquid metals. Incompressible MHD model also governs the terrestrial magnetic filed maintained by fluid motion in the earth core and the solar magnetic field which generates sunspots and solar flares[5]. The incompressible MHD model consists of the incompressible NavierStokes equations and the quasistatic Maxwell equations. The magnetic field influences the momentum of the fluid through Lorentz force, and conversely, the motion of fluid influences the magnetic field through Faraday’s law. In this paper, we are studying the efficient iterative solver for the stationary MHD equations
(1a)  
(1b)  
(1c) 
where is the velocity of the fluid, is the hydrodynamic pressure, is the magnetic flux density or the magnetic field provided with constant permeability, is the fluid Reynolds number, is the magnetic Reynolds number, is the coupling constant concerning the Lorentz force, and stands for the external force. We assume that is a bounded Lipschitz domain. The system of equations are complemented with Dirichlet boundary conditions
(2) 
There are extensive papers in the literature to study numerical solutions of incompressible MHD equations (cf. e.g. [14, 16, 18, 13, 21, 24, 26, 27, 28, 31] and the references therein). In [18], Gunzburger et al studied wellposedness and the finite element method for the stationary incompressible MHD equations. The magnetic field is discretized by the conforming finite element method. Strauss et al studied the adaptive finite element method for twodimensional MHD equations [21]. Very recently, based on the nodal finite element approximation to the magnetic field, Philips et al proposed a block preconditioner based on an exact penalty formulation of stationary MHD equations [31]. We also refer to [14] for a systematic analysis on finite element methods for incompressible MHD equations. When the domain has reentrant angle, the magnetic field may not be in . It is preferable to use noncontinuous finite element functions to approximate , namely, the socalled edge element method[10, 25]. In 2004, Schötzau proposed a mixed finite element method to solve the stationary incompressible MHD equations where edge elements are used to solve the magnetic field. To our knowledge, efficient solvers for threedimensional (3D) MHD equations are still rare in the literature, particularly, for large Reynolds number and large coupling number . An efficient solver should possess two merits:

the convergence rate is independent of the mesh or the number of degrees of freedom (DOFs);

the algorithm is robust with respect to the physical parameters.
The objective of this paper is to propose a preconditioned GMRES method to solve the linearized discrete problem of (1)–(2). We shall adopt the mixed finite element method proposed in [32] and study efficient preconditioners for the linearized problem.
Over the past three decades, fast solvers for incompressible NavierStokes equations are relatively wellstudied in the literature (cf. e.g. [6, 7, 8, 29, 33, 34]). For moderate Reynolds number, the Picard iteration for stationary impressible NavierStokes equation is stable and efficient. At each iteration, one needs to solve the linearized problem, the Oseen equations
(3a)  
(3b)  
(3c) 
where is the approximate solution at the previous step. Iterative methods for discrete Oseen equations mainly consist of Krylov subspace methods, multigrid methods, or their combinations. In terms of parallel computing and practical implementation, it is preferable to use Krylov subspace method combined with an effective preconditioner. Among them, the pressure convectiondiffusion (PCD) preconditioner [19], the leastsquares commutator (LSC) preconditioner [6, 7, 8], and the augmented Lagrangian(AL) preconditioner [1, 2, 30] prove robust and efficient for relatively large Reynolds number. In this paper, we shall study the AL finite element method for the stationary MHD equations. Based on a Picardtype linearization of the disrecte problem, we develop an efficient preconditioner for solving the linear problem. The preconditioner proves to be robust when the Reynolds number and the coupling number are relatively large and to be optimal with respect to the number of DOFs.
The paper is organized as follows: In section 2, we introduce some notations for Sobolev spaces. A mixed finite element method is proposed to solve the AL formulation of the stationary MHD equations. In section 3, we introduce a Picardtype linearization for the discrete MHD problem and devise an efficient preconditioner for solving the linear discrete problem. A preconditioned GMRES algorithm is also presented for the implementation of the discrete solver. In Section 4, we present three numerical experiments to verify the optimal convergence rate of the mixed finite element method, to demonstrate the optimality and the robustness of the MHD solver, and to demonstrate the scalability for parallel computing. Throughout the paper we denote vectorvalued quantities by boldface notation, such as .
2 Mixed finite element method for the MHD equations
First we introduce some Hilbert spaces and Sobolev norms used in this paper. Let be the usual Hilbert space of square integrable functions equipped with the following inner product and norm:
Let the quotient space of be defined by
Define where represents nonnegative triple index. Let be the subspace of whose functions have zero traces on .
We define the spaces of functions having square integrable curl by
which are equipped with the following inner product and norm
Here denotes the unit outer normal to .
With a Lagrange multiplier , we can rewrite (1) into an AL form
(4a)  
(4b)  
(4c)  
(4d) 
where is the stabilization parameter or penalty parameter. Taking divergence on both sides of (4b) and using (4d) yields
This means in actually. Note that , thus (4) is equivalent to (1)–(2). In the rest of this paper, we are going to study the augmented problem (4) instead of the original problem.
A weak formulation of (4) reads: Find and such that and on and
(5a)  
(5b) 
for all and , where the bilinear forms and trilinear form are defined respectively by
Assuming small data, Schötzau proved the existence and uniqueness of the solution to (5) without the penalized term . The purpose of this paper is to propose a robust solver for the discrete problem. This extra term in the new formula makes the discrete problem more welldefined for high Reynolds number[30].
Now we introduce the finite element approximation to (5). Let be a quasiuniform and shaperegular tetrahedral mesh of . Let denote the maximal diameter of all tetrahedra on the mesh. For any , let be the space of polynomials of degree on and be the corresponding space of vector polynomials. Define the Lagrange finite element space of the th order by
First we choose the wellknown TaylorHood  elements [3, Page 217219] for the discretization of , namely,
From [3, Page 255258], the discrete infsup condition holds
(6) 
where is the infsup constant independent of the mesh size. We shall also use .
The finite element space for is chosen as Nédélec’s edge element space of the first order in the second family[25], namely,
The finite element space for is defined by
Since , we easily get the infsup condition for the pair of finite element spaces
(7) 
where is the Poincáre constant depending only on .
The finite element approximation to (5) reads: Find and such that
(8a)  
(8b) 
for all and . From (6) and (7) we know that the bilinear form satisfies the discrete infsup condition
(9) 
where
Based on the assumption of small data, we can prove that the discrete problem (8) has a unique solution. Again we do not elaborate on the details and pay our attention to fast solvers of the discrete solution (8).
3 A preconditioner for the linearized finite element problem
In this section, we are going to study the solution of the nonlinear discrete problem (8). First we propose a Picardtype iterative method for solving (8). At each nonlinear iteration, the linearized problem consists of an AL Oseen equation with Lorentz force and a Maxwell equation coupled with the fluid. The preconditioner for the linearized MHD equation depends crucially on the preconditioner for the penalized NavierStokes equations and the preconditioner for the Maxwell equations in mixed forms.
3.1 Picardtype method for the discrete MHD equations
In this subsection, we consider the Picard linearization of (8). For convenience, we rearrange the order of variables as in the linearized problem. Let be the approximate solutions of (8) from the previous iteration. The error equation for these approximate solutions reads: Find such that
(10a)  
(10b)  
(10c)  
(10d) 
where the trilinear form represents the convectiondiffusion part of the fluid equation
and the residual functionals are defined by
After solving (10), the approximate solutions will be updated by
(11) 
with a relaxation factor .
To devise the preconditioner, we write problem (10) into an algebraic form
(12) 
where the solution vector consists of the degrees of freedom for respectively, is the residual vector, and is the stiffness matrix. In block forms, they can be written as
(13) 
Let , , , be the bases of , , , and respectively. Then the entries of all block matrices are defined by
Clearly the block matrices represents the differential operators appearing in the NavierStokes equations and the Maxwell equation on various finite element spaces
Here is understood as the dual operator of or . Moreover, are algebraic representations of the two multiplication operators which couple the magnetic field and the conducting fluid. For any given and , we have
(14) 
The relationships between these operators play an important role in devising a robust preconditioner for the linearized problem.
3.2 Preconditioning for the linearized MHD equations
Let be the stiffness matrix of on and let be a constant. First we postmultiply the second column of by and add it to the first column. This yields a matrix
Next, premultiplying the first row of by and adding it to the second row, we get a matrix
Note that represents the operator . Since
formally we have
This means that
Since represents the coupling term , we have . Therefore,
Moreover, from [17], we know that is equivalent to in spectrum where is the mass matrix on . So one gets the approximation
Next, premultiplying the first row of by and adding it to the third row, we get a matrix
Since , formally we have . Here is the identity operator. This means
So one obtains an approximation of as follows
(15) 
This means that is actually a natural preconditioner for , except for the difficulties in computing the inverse of the block . In the next subsection, we are going to derive a good approximation of so that its approximation inverse is easy to compute iteratively.
3.3 An efficient preconditioner for the magnetic fieldfluid coupling block
From (15), the key step to compute is how to precondition the block
(16) 
We note that is the precise Schur complement of the following matrix which accounts for the coupling between and
(17) 
Remember from (14) that and represent, respectively, the two multiplication operators and for any given and . So is the algebraic representation of
(18) 
Since commutates with , then (18) becomes
For sufficiently small, we adopt the following approximation
This yields an approximation of (18)
(19) 
Therefore, we get an approximation of the Schur complement
where is the stiffness matrix associated with the bilinear form
This should yield a good preconditioner of , that is,
(20) 
Remark 3.1
In [31], Philips et al studied the preconditioner for twodimensional MHD equations where the magnetic field is discretized with conforming finite elements. Based on an exact penalty formulation, they propose to approximate the coupling effect between magnetic field and velocity by
where is a parameter depending on the mesh size and the magnitude of .
Compared with [31], for the 3D MHD equations and the conforming approximation of , our approximation to the coupling effect in the preconditioner level does not need the extra parameter and has more advantages in practical computations such as adaptive computing, though they are similar.
Now we demonstrate numerically the robustness of the preconditioner in (20) with respect to the parameter and the mesh size . Since it is the approximation that is concerned here, we fix and and test the efficiency of the preconditioner for different values of and .
We consider the system of linear equations: Find and such that
(21a)  
(21b) 
where
Clearly the stiffness matrix of (21) is which is given in (17). In the following, we test three cases of
and three cases of physical parameters
The computational domain is the unit cube, namely, .
0.216506  4  14  91 

0.108253  4  14  74 
0.054127  4  14  64 
0.027063  4  14  61 
0.216506  4  14  91 

0.108253  4  14  73 
0.054127  4  14  64 
0.027063  4  14  61 
0.216506  4  14  91 

0.108253  4  14  73 
0.054127  4  14  64 
0.027063  4  14  61 
0.216506  3  14  82 

0.108253  4  14  92 
0.054127  4  15  97 
0.027063  4  15  98 
We use preconditioned GMRES method to solve (21) and the preconditioner is set by (20). This means that we need solve the residual equation at each GMRES iteration
(22) 
where stand for the residual vectors and stand for the error vectors. The tolerance for the relative residual of the GMRES method is set by . The tolerances for solving the two subproblems in (22) are set by . From Table 1–3, we find that the convergence of the preconditioned GMRES is uniform with respect to both and . An interesting observation is that, for large , the number of GMRES iterations even decreases when . In this case, the magnetic fieldfluid coupling becomes strong.
Table 4 shows the number of preconditioned GMRES iterations where the approximate Schur complement in (20) is replaced with the matrix . It amounts to devise a preconditioner of by dropping its left lower block
(23) 
This is the classical Riesz map preconditioning in [20] or the operator preconditioning in [11]. Comparing Table 3 with Table 4, we find that, for large , the convergence of GMRES method with this preconditioner becomes slower and deteriorates when . This becomes even more apparent when solving the whole MHD system (see Table 7 for the computation of driven cavity flow).
3.4 A preconditioner for the augmented NavierStokes equations
Combining (15) and (20), we get a preconditioner of , that is, the inverse of
(24) 
It is left to study the preconditioner for the lower right block of , namely,
It amounts to solve the saddle point problem: Find such that
(25a)  
(25b) 
In [1, 2], Benzi et al studied the Osceen equation(namely ) and proposed to use the following preconditioner
(26) 
where is the mass matrix on . It is proved that the above preconditioner is efficient for relatively large Reynolds number. With , namely the Stokes equations, we refer to [22, 23] for similar arguments. And with for the timedependent incompressible MHD, the work in [24] give useful insight using pressure mass matrix as a subblock. Inspired by them, we propose to precondition
(27) 
3.5 A robust preconditioner for the linearized MHD problem
Using (24) and (27), a preconditioner for is given by the inverse of
(28) 
where the parameter can be used to tune the efficiency of the preconditioner. According to our experience, any value works well for high Reynolds and moderate and . Moreover, to fix the parameter , we set so that is associated with the bilinear form
This yields our final preconditioner for the stiffness matrix , defined by
(29) 
Now we are in the position to present the preconditioned GMRES algorithm for solving the linear system (12) of the MHD problem. The idea is to use an approximation of to precondition . For convenience in notation, given a vector which has the same size as one column vector of , we let be the vectors which consist of entries of and correspond to respectively.
Algorithm 3.2 (Preconditioned GMRES Algorithm)
Given the tolerances and , the maximal number of GMRES iterations , and the initial guess for the solution of (12). Set and compute the residual vector
While do

Solve by the CG method with the diagonal preconditioning and tolerance for the relative residual.

Solve by preconditioned GMRES method with tolerance . The preconditioner is the one level additive Schwarz method with overlap = 2 [4].

Solve by preconditioned CG method with tolerance . The preconditioner is the algebraic multigrid method (AMG) solver [9].

Solve