A robust solver for the finite element approximation of stationary incompressible MHD equations in 3D

# A robust solver for the finite element approximation of stationary incompressible MHD equations in 3D ††thanks:

Lingxiao Li Academy of Mathematics and System Sciences, Chinese Academy of Sciences; School of Mathematical Science, University of Chinese Academy of Sciences, Beijing, 100190, China.(lilingxiao@lsec.cc.ac.cn)    Weiying Zheng Academy of Mathematics and System Sciences, Chinese Academy of Sciences; School of Mathematical Science, University of Chinese Academy of Sciences, Beijing, 100190, China. This author is supported by China NSF under the grants 91430215 and by the National Magnetic Confinement Fusion Science Program 2015GB110003.(zwy@lsec.cc.ac.cn)
###### 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 Navier-Stokes equations and the quasi-static 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

 u⋅∇u+∇p−R−1eΔu−ScurlB×B=finΩ, (1a) curl(B×u+R−1mcurlB)=0% inΩ, (1b) divu=0,divB=0inΩ, (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

 u=g,B×n=Bs×nonΓ:=∂Ω. (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 well-posedness 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 two-dimensional 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 re-entrant angle, the magnetic field may not be in . It is preferable to use noncontinuous finite element functions to approximate , namely, the so-called 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 three-dimensional (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:

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

2. 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 Navier-Stokes equations are relatively well-studied in the literature (cf. e.g. [6, 7, 8, 29, 33, 34]). For moderate Reynolds number, the Picard iteration for stationary impressible Navier-Stokes equation is stable and efficient. At each iteration, one needs to solve the linearized problem, the Oseen equations

 w⋅∇u+∇p−R−1eΔu=f inΩ, (3a) divu=0 inΩ, (3b) u=g onΓ, (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 convection-diffusion (PCD) preconditioner [19], the least-squares 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 Picard-type 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 Picard-type 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 vector-valued 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:

 (u,v):=∫Ωu(x)v(x)dxand∥u∥L2(Ω):=(u,u)1/2.

Let the quotient space of be defined by

 L20(Ω):={v∈L2(Ω):∫Ωv(x)dx=0}=L2(Ω)/R.

Define where represents non-negative triple index. Let be the subspace of whose functions have zero traces on .

We define the spaces of functions having square integrable curl by

 H(curl,Ω) := {v∈L2(Ω):curlv∈L2(Ω)}, H0(curl,Ω) := {v∈H(curl,Ω):n×v=0onΓ},

which are equipped with the following inner product and norm

 (v,w)H(curl,Ω):=(v,w)+(curlv,curlw),∥v∥H(curl,Ω):=√(v,v)H(curl,Ω).

Here denotes the unit outer normal to .

With a Lagrange multiplier , we can rewrite (1) into an AL form

 u⋅∇u+∇p−γ∇divu−R−1eΔu−ScurlB×B=finΩ, (4a) Scurl(B×u+R−1mcurlB)+∇r=0inΩ, (4b) divu=0,divB=0inΩ, (4c) u=g,B×n=Bs×n,r=0onΓ. (4d)

where is the stabilization parameter or penalty parameter. Taking divergence on both sides of (4b) and using (4d) yields

 Δr=0inΩ,r=0inΓ.

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

 A((u,B),(v,φ))+O((u,B);(u,B),(v,φ))−B((p,r),(v,φ)) =(f,v), (5a) B((q,s),(u,B)) =0, (5b)

for all and , where the bilinear forms and trilinear form are defined respectively by

 A((u,B),(v,φ)) = R−1e(∇u,∇v)+γ(divu,divv)+SR−1m(∇×B,∇×φ), O((w,ψ);(u,B),(v,φ)) = (w⋅∇u,v)−S[(curlB,ψ×v)−(ψ×u,curlφ)], B((p,r),(v,φ)) = (p,divv)+(∇r,φ).

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 well-defined for high Reynolds number[30].

Now we introduce the finite element approximation to (5). Let be a quasi-uniform and shape-regular 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

 V(k,Th)={v∈H1(Ω):v|T∈Pk+1(T),∀T∈Th}.

First we choose the well-known Taylor-Hood - elements [3, Page 217-219] for the discretization of , namely,

 Vh:=V(2,Th)3∩H10(Ω),Qh:=V(1,Th).

From [3, Page 255-258], the discrete inf-sup condition holds

 sup0≠v∈Vh(q,divv)∥v∥H1(Ω)≥Cu∥q∥L2(Ω)∀q∈Qh, (6)

where is the inf-sup 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,

 ¯¯¯¯¯Ch={v∈H(curl,Ω):v|T∈P1(T),∀T∈Th},Ch=¯¯¯¯¯Ch∩H0(curl,Ω).

The finite element space for is defined by

 Sh=V(2,Th)∩H10(Ω).

Since , we easily get the inf-sup condition for the pair of finite element spaces

 sup0≠v∈Ch(∇s,v)∥v∥H(curl,Ω)≥|s|H1(Ω)≥Cb∥s∥H1(Ω)∀s∈Sh, (7)

where is the Poincáre constant depending only on .

The finite element approximation to (5) reads: Find and such that

 A((uh,Bh),(v,φ))+O((uh,Bh);(uh,Bh),(v,φ))−B((ph,rh),(v,φ)) =(f,v), (8a) B((q,s),(uh,Bh)) =0, (8b)

for all and . From (6) and (7) we know that the bilinear form satisfies the discrete inf-sup condition

 sup(vh,φh)∈Vh×ChB((qh,sh),(vh,φh))∥(vh,φh)∥Vh×Ch≥min(Cu,Cp)∥(qh,sh)∥Qh×Sh∀(qh,sh)∈Qh×Sh, (9)

where

 ∥(vh,φh)∥Vh×Ch:=√∥vh∥2H1(Ω)+∥φh∥2H(curl,Ω),∥(qh,sh)∥Qh×Sh:=√∥qh∥2L2(Ω)+∥sh∥2H1(Ω).

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 Picard-type 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 Navier-Stokes equations and the preconditioner for the Maxwell equations in mixed forms.

### 3.1 Picard-type 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

 SR−1m(curlδBk,curlφ)+(∇δrk,φ)+S(Bk×δuk,curlφ)=Rb(φ), (10a) (δBk,∇s)=Rr(s), (10b) −S(curlδBk,Bk×v)+F(uk;δuk,v)−(δpk,divv)=Ru(v), (10c) −(divδuk,q)=Rp(q), (10d)

where the trilinear form represents the convection-diffusion part of the fluid equation

 F(uk;δuk,v):=R−1e(∇δuk,∇v)+(uk⋅∇δuk,v)+γ(divδuk,divv),

and the residual functionals are defined by

 Rb(φ)=−SR−1m(curlBk,curlφ)−(∇rk,φ)−S(Bk×uk,curlφ), Rr(s)=−(Bk,∇s), Ru(v)=(f,v)−F(uk;uk,v)+S(curlBk,Bk×v)+(pk,divv), Rp(q)=(divuk,q).

After solving (10), the approximate solutions will be updated by

 Bk+1=Bk+θδBk,  rk+1=rk+θδrk,  uk+1=uk+θδuk,  pk+1=pk+θδpk (11)

with a relaxation factor .

To devise the preconditioner, we write problem (10) into an algebraic form

 Ax=b, (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

 x=⎛⎜ ⎜ ⎜⎝xbxrxuxp⎞⎟ ⎟ ⎟⎠,b=⎛⎜ ⎜ ⎜⎝RbRrRuRp⎞⎟ ⎟ ⎟⎠,A=⎛⎜ ⎜ ⎜ ⎜⎝CG⊤J⊤0G000−J0FB⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠. (13)

Let , , , be the bases of , , , and respectively. Then the entries of all block matrices are defined by

 Cij =SR−1m(curlφj,curlφi), Gij =S(φj,∇si), Jij =S(curlφj,Bk×vi) Fij =F(uk;vj,vi), Bij =−(divvj,qi).

Clearly the block matrices represents the differential operators appearing in the Navier-Stokes equations and the Maxwell equation on various finite element spaces

 C⇔SR−1mcurlcurl,G⇔−div onCh, G⊤⇔∇ onSh, F⇔(−R−1eΔ+uk⋅∇−γ∇div),B⇔−div onVh, B⊤⇔∇ onQh.

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

 J⊤⇔Scurl(Bk×w)onCh,J⇔Scurlψ×BkonVh. (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 post-multiply the second column of by and add it to the first column. This yields a matrix

 A1=⎛⎜ ⎜ ⎜ ⎜⎝C+σSrG⊤J⊤0G000−J0FB⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠,Sr:=G⊤L−1rG.

Next, pre-multiplying the first row of by and adding it to the second row, we get a matrix

 A2=⎛⎜ ⎜ ⎜ ⎜⎝C+σSrG⊤J⊤00−G(C+σSr)−1G⊤−G(C+σSr)−1J⊤0−J0FB⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠.

Note that represents the operator . Since

 div(SR−1mcurlcurl+σ∇Δ−1div)=σΔΔ−1div=σdiv,

formally we have

 div(SR−1mcurlcurl+σ∇Δ−1div)−1=σ−1div.

This means that

 G(C+σSr)−1≈σ−1G.

Since represents the coupling term , we have . Therefore,

 G(C+σSr)−1J⊤≈0.

Moreover, from [17], we know that is equivalent to in spectrum where is the mass matrix on . So one gets the approximation

 A2≈A3:=⎛⎜ ⎜ ⎜ ⎜⎝C+σMG⊤J⊤00−σ−1Lr00−J0FB⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠.

Next, pre-multiplying the first row of by and adding it to the third row, we get a matrix

 A4=⎛⎜ ⎜ ⎜ ⎜⎝C+σMG⊤J⊤00−σ−1Lr000J(C+σM)−1G⊤F+J(C+σM)−1J⊤B⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠.

Since , formally we have . Here is the identity operator. This means

 G(C+σM)−1J⊤≈0orJ(C+σM)−1G⊤≈0.

So one obtains an approximation of as follows

 A4≈A5:=⎛⎜ ⎜ ⎜ ⎜⎝C+σMG⊤J⊤00−σ−1Lr0000F+J(C+σM)−1J⊤B⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠. (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 field-fluid 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

 ^X=(C+σMJ⊤−JF). (17)

Remember from (14) that and represent, respectively, the two multiplication operators and for any given and . So is the algebraic representation of

 Scurl{(SR−1mcurlcurl+σI)−1Scurl(Bk×w)}×Bk (18)

Since commutates with , then (18) becomes

 S2{(SR−1mcurlcurl+σI)−1curlcurl(Bk×w)}×Bk

For sufficiently small, we adopt the following approximation

 curlcurl≈S−1Rm(SR−1mcurlcurl+σI).

This yields an approximation of (18)

 S2{(SR−1mcurlcurl+σI)−1curlcurl(Bk×w)}×Bk≈SRm(Bk×w)×Bk. (19)

Therefore, we get an approximation of the Schur complement

 F+J(C+σM)−1J⊤≈S,

where is the stiffness matrix associated with the bilinear form

 F(uk;u,v)+SRm(Bk×u,Bk×v).

This should yield a good preconditioner of , that is,

 (20)
###### Remark 3.1

In [31], Philips et al studied the preconditioner for two-dimensional 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

 βSRm(Bk×w)×Bk,

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

 −S(curlδB,B0×v)+F(u0;δu,v)=(f,v) ∀v∈Vh, (21a) SR−1m(curlδB,curlφ)+σ(δB,φ)+S(B0×δuk,curlφ)=0 ∀φ∈Ch, (21b)

where

 f=(1,sin(x),0),u0=(y,sin(x+z),1),B0=(sin(y)+cos(z),1−sin(x),1).

Clearly the stiffness matrix of (21) is which is given in (17). In the following, we test three cases of

 σ=1,10−2,10−4,

and three cases of physical parameters

 S=Rm=1,10,100.

The computational domain is the unit cube, namely, .

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

 Seu=ru,(C+σM)eb=rb−J⊤eu, (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 sub-problems in (22) are set by . From Table 13, 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 field-fluid 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

 (C+σMJ⊤−JF)∼(C+σMJ⊤0F). (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 Navier-Stokes equations

Combining (15) and (20), we get a preconditioner of , that is, the inverse of

 A6:=⎛⎜ ⎜ ⎜ ⎜⎝C+σMG⊤J⊤00−σ−1Lr0000SB⊤00B0⎞⎟ ⎟ ⎟ ⎟⎠. (24)

It is left to study the preconditioner for the lower right block of , namely,

 (SB⊤B0).

It amounts to solve the saddle point problem: Find such that

 F(uk;δu,v)+SRm(Bk×δu,Bk×v)−(δp,divv) =Ru(v)∀v∈Vh, (25a) −(divδu,q) =Rp(q)∀q∈Qh. (25b)

In [1, 2], Benzi et al studied the Osceen equation(namely ) and proposed to use the following preconditioner

 (FB⊤0−(R−1e+γ)−1Qp)−1, (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 time-dependent incompressible MHD, the work in [24] give useful insight using pressure mass matrix as a subblock. Inspired by them, we propose to precondition

 (SB⊤B0)by(SB⊤0−(R−1e+γ)−1Qp)−1. (27)

### 3.5 A robust preconditioner for the linearized MHD problem

Using (24) and (27), a preconditioner for is given by the inverse of

 ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝C+σMG⊤J⊤00−σ−1Lr0000SB⊤000−(R−1e+γ)−1Qp⎞⎟ ⎟ ⎟ ⎟ ⎟⎠, (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

 SR−1m[(curlu,curlv)+(u,v)].

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

 r(k)=b−Ax(k).

While do

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

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

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

4. Solve