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

TL;DR

This paper introduces a robust finite element solver for 3D stationary incompressible MHD equations, featuring an efficient preconditioner that enhances convergence and scalability for practical computations.

Contribution

It develops a new preconditioner for the finite element discretization of 3D stationary incompressible MHD equations, improving solver robustness and efficiency.

Findings

Preconditioner accelerates GMRES convergence.

Solver demonstrates robustness across numerical experiments.

Scalability is confirmed for moderate parameters.

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 H1-conforming finite elements, while the magnetic field is approximated by H(curl)-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.

Tables

Table 1. Table 1: Number of preconditioned GMRES iterations for σ = 1 𝜎 1 \sigma=1 .

$h$	$S=R_m=1$	$S=R_m=10$	$S=R_m=100$
0.216506	4	14	91
0.108253	4	14	74
0.054127	4	14	64
0.027063	4	14	61

Table 2. Table 2: Number of preconditioned GMRES iterations for σ = 10 − 2 𝜎 superscript 10 2 \sigma=10^{-2} .

$h$	$S=R_m=1$	$S=R_m=10$	$S=R_m=100$
0.216506	4	14	91
0.108253	4	14	73
0.054127	4	14	64
0.027063	4	14	61

Table 3. Table 3: Number of preconditioned GMRES iterations for σ = 10 − 4 𝜎 superscript 10 4 \sigma=10^{-4} .

$h$	$S=R_m=1$	$S=R_m=10$	$S=R_m=100$
0.216506	4	14	91
0.108253	4	14	73
0.054127	4	14	64
0.027063	4	14	61

Table 4. Table 4: Number of preconditioned GMRES iterations for σ = 10 − 4 𝜎 superscript 10 4 \sigma=10^{-4} with 𝕊 = 𝔽 𝕊 𝔽 \mathbb{S}=\mathbb{F} .

$h$	$S=R_m=1$	$S=R_m=10$	$S=R_m=100$
0.216506	3	14	82
0.108253	4	14	92
0.054127	4	15	97
0.027063	4	15	98

Table 5. Table 5: Convergence rate of the finite element discrete problem.

$h$	${\Vert 𝒖-{𝒖}_h\Vert }_{H^1}$	order	${\Vert p-p_h\Vert }_{L^2}$	order	${\Vert 𝑩-{𝑩}_h\Vert }_{𝑯\left(\mathrm{𝐜𝐮𝐫𝐥}\right)}$	order
0.4330	2.893e-03	—	1.848e-03	—	4.811e-02	—
0.2165	7.071e-04	2.033	3.908e-04	2.242	2.378e-02	1.017
0.1083	1.745e-04	2.019	9.197e-05	2.087	1.182e-02	1.009
0.0541	4.335e-05	2.009	2.214e-05	2.055	5.893e-03	1.004

Table 6. Table 6: The mesh sizes and the numbers of DOFs.

Mesh	$h$	DOFs for $(𝑩,r)$	DOFs for $(𝒖,p)$
${𝒯}_1$	0.2165	13281	15468
${𝒯}_2$	0.1083	97985	112724
${𝒯}_3$	0.0541	752001	859812
${𝒯}_4$	0.0271	5890817	6714692

Table 7. Table 7: Average GMRES iteration number: R e = 100.0 , S = 100.0 , R m = 1.0 formulae-sequence subscript 𝑅 𝑒 100.0 formulae-sequence 𝑆 100.0 subscript 𝑅 𝑚 1.0 R_{e}=100.0,S=100.0,R_{m}=1.0 .

Mesh	$N_{\mathrm{picard}}\times N_{\mathrm{gmres}}$ with $\mathrm{𝐁𝐮𝐁𝐯}$	$N_{\mathrm{picard}}\times N_{\mathrm{gmres}}$ without $\mathrm{𝐁𝐮𝐁𝐯}$
${𝒯}_1$	$6\times 51.5$	$7\times 108.1$
${𝒯}_2$	$6\times 43.5$	$7\times 102.3$
${𝒯}_3$	$6\times 36.8$	$7\times 192.6$
${𝒯}_4$	$6\times 32.5$	$>7\times 200.0$

Table 8. Table 8: Scalability for the discrete solver(Average Iter and Total time).

Mesh	Total DOFs	Cores	$N_{\mathrm{gmres}}$	Time (s)	Efficiency
${𝒯}_1$	210709	1	16.3	106.6	—
${𝒯}_2$	403221	2	14.5	128.2	83.2%
${𝒯}_3$	850965	4	15.8	235.7	45.2%
${𝒯}_4$	1611813	8	15.8	300.6	35.5%
${𝒯}_5$	3151909	16	14.3	299.7	35.6%

Equations

$$\mathbf{u}\cdot \nabla \mathbf{u}+\nabla p-R_e^{-1}\Delta \mathbf{u}-S\mathrm{curl}\mathbf{B}\times \mathbf{B}=\mathbf{f}in\Omega ,$$

$$\mathrm{curl}\left(\mathbf{B}\times \mathbf{u}+R_m^{-1}\mathrm{curl}\mathbf{B}\right)=0in\Omega ,$$

$$\mathrm{div}\mathbf{u}=0,\mathrm{div}\mathbf{B}=0in\Omega ,$$

$$\mathbf{u}=\mathbf{g},\mathbf{B}\times \mathbf{n}={\mathbf{B}}_s\times \mathbf{n}on\Gamma :=\partial \Omega .$$

$$\mathbf{w}\cdot \nabla \mathbf{u}+\nabla p-R_e^{-1}\Delta \mathbf{u}=\mathbf{f}$$

$$\mathrm{div}\mathbf{u}=0$$

$$\mathbf{u}=\mathbf{g}$$

$$(u,v):={\int }_{\Omega }u(\mathbf{x})v(\mathbf{x})\mathrm{d}\mathbf{x}and{\Vert u\Vert }_{L^2(\Omega )}:=(u,u{)}^{1/2}.$$

$$L_0^2(\Omega ):=\left\{v\in L^2(\Omega ):{\int }_{\Omega }v(\mathbf{x})\mathrm{d}\mathbf{x}=0\right\}=L^2(\Omega )/\mathbb{R}.$$

$$\mathbf{H}(\mathrm{curl},\Omega )$$

$${\mathbf{H}}_0(\mathrm{curl},\Omega )$$

$$(\mathbf{v},\mathbf{w}{)}_{\mathbf{H}(\mathrm{curl},\Omega )}:=(\mathbf{v},\mathbf{w})+(\mathrm{curl}\mathbf{v},\mathrm{curl}\mathbf{w}),{\Vert \mathbf{v}\Vert }_{\mathbf{H}(\mathrm{curl},\Omega )}:=\sqrt{(\mathbf{v},\mathbf{v}{)}_{\mathbf{H}(\mathrm{curl},\Omega )}}.$$

$$\mathbf{u}\cdot \nabla \mathbf{u}+\nabla p-\gamma \nabla \mathrm{div}\mathbf{u}-R_e^{-1}\Delta \mathbf{u}-S\mathrm{curl}\mathbf{B}\times \mathbf{B}=\mathbf{f}in\Omega ,$$

$$S\mathrm{curl}\left(\mathbf{B}\times \mathbf{u}+R_m^{-1}\mathrm{curl}\mathbf{B}\right)+\nabla r=0in\Omega ,$$

$$\mathrm{div}\mathbf{u}=0,\mathrm{div}\mathbf{B}=0in\Omega ,$$

$$\mathbf{u}=\mathbf{g},\mathbf{B}\times \mathbf{n}={\mathbf{B}}_s\times \mathbf{n},r=0on\Gamma .$$

$$\Delta r=0in\Omega ,r=0in\Gamma .$$

$$\mathcal{A}((\mathbf{u},\mathbf{B}),(\mathbf{v},\phi ))+\mathcal{O}((\mathbf{u},\mathbf{B});(\mathbf{u},\mathbf{B}),(\mathbf{v},\mathbf{\phi }))-\mathcal{B}((p,r),(\mathbf{v},\mathbf{\phi }))$$

$$\mathcal{B}((q,s),(\mathbf{u},\mathbf{B}))$$

$$\mathcal{A}((\mathbf{u},\mathbf{B}),(\mathbf{v},\phi ))$$

$$\mathcal{O}((\mathbf{w},\mathbf{\psi });(\mathbf{u},\mathbf{B}),(\mathbf{v},\phi ))$$

$$\mathcal{B}((p,r),(\mathbf{v},\mathbf{\phi }))$$

$$V(k,{\mathcal{T}}_h)=\left\{v\in H^1(\Omega ):v{|}_T\in P_{k+1}(T),\forall T\in {\mathcal{T}}_h\right\}.$$

$${\mathbf{V}}_h:=V(2,{\mathcal{T}}_h{)}^3\cap {\mathbf{H}}_0^1(\Omega ),Q_h:=V(1,{\mathcal{T}}_h).$$

$$\underset{0\ne \mathbf{v}\in {\mathbf{V}}_h}{\sup }\frac{(q,\mathrm{div}\mathbf{v})}{{\Vert \mathbf{v}\Vert }_{H^1(\Omega )}}\ge C_u{\Vert q\Vert }_{L^2(\Omega )}\forall q\in Q_h,$$

$${\overline{\mathbf{C}}}_h=\left\{\mathbf{v}\in \mathbf{H}(\mathrm{curl},\Omega ):\mathbf{v}{|}_T\in {\mathbf{P}}_1(T),\forall T\in {\mathcal{T}}_h\right\},{\mathbf{C}}_h={\overline{\mathbf{C}}}_h\cap {\mathbf{H}}_0(\mathrm{curl},\Omega ).$$

$$S_h=V(2,{\mathcal{T}}_h)\cap H_0^1(\Omega ).$$

$$\underset{0\ne \mathbf{v}\in {\mathbf{C}}_h}{\sup }\frac{(\nabla s,\mathbf{v})}{{\Vert \mathbf{v}\Vert }_{\mathbf{H}(\mathrm{curl},\Omega )}}\ge {\left|s\right|}_{H^1(\Omega )}\ge C_b{\Vert s\Vert }_{H^1(\Omega )}\forall s\in S_h,$$

$$\mathcal{A}(({\mathbf{u}}_h,{\mathbf{B}}_h),(\mathbf{v},\mathbf{\phi }))+\mathcal{O}(({\mathbf{u}}_h,{\mathbf{B}}_h);({\mathbf{u}}_h,{\mathbf{B}}_h),(\mathbf{v},\mathbf{\phi }))-\mathcal{B}((p_h,r_h),(\mathbf{v},\mathbf{\phi }))$$

$$\mathcal{B}((q,s),({\mathbf{u}}_h,{\mathbf{B}}_h))$$

$$\underset{({\mathbf{v}}_h,{\mathbf{\phi }}_h)\in {\mathbf{V}}_h\times {\mathbf{C}}_h}{\sup }\frac{\mathcal{B}((q_h,s_h),({\mathbf{v}}_h,{\mathbf{\phi }}_h))}{{\Vert ({\mathbf{v}}_h,{\mathbf{\phi }}_h)\Vert }_{{\mathbf{V}}_h\times {\mathbf{C}}_h}}\ge \min (C_u,C_p){\Vert (q_h,s_h)\Vert }_{Q_h\times S_h}\forall (q_h,s_h)\in Q_h\times S_h,$$

$${\Vert ({\mathbf{v}}_h,{\mathbf{\phi }}_h)\Vert }_{{\mathbf{V}}_h\times {\mathbf{C}}_h}:=\sqrt{{\Vert {\mathbf{v}}_h\Vert }_{{\mathbf{H}}^1(\Omega )}^2+{\Vert {\mathbf{\phi }}_h\Vert }_{\mathbf{H}(\mathrm{curl},\Omega )}^2},{\Vert (q_h,s_h)\Vert }_{Q_h\times S_h}:=\sqrt{{\Vert q_h\Vert }_{L^2(\Omega )}^2+{\Vert s_h\Vert }_{H^1(\Omega )}^2}.$$

$$S{R}_m^{-1}(\mathrm{curl}\delta {\mathbf{B}}_k,\mathrm{curl}\mathbf{\phi })+(\nabla \delta r_k,\mathbf{\phi })+S({\mathbf{B}}_k\times \delta {\mathbf{u}}_k,\mathrm{curl}\mathbf{\phi })=R_b(\mathbf{\phi }),$$

$$(\delta {\mathbf{B}}_k,\nabla s)=R_r(s),$$

$$-S(\mathrm{curl}\delta {\mathbf{B}}_k,{\mathbf{B}}_k\times \mathbf{v})+\mathcal{F}({\mathbf{u}}_k;\delta {\mathbf{u}}_k,\mathbf{v})-(\delta p_k,\mathrm{div}\mathbf{v})=R_u(\mathbf{v}),$$

$$-(\mathrm{div}\delta {\mathbf{u}}_k,q)=R_p(q),$$

$$\mathcal{F}({\mathbf{u}}_k;\delta {\mathbf{u}}_k,\mathbf{v}):=R_e^{-1}(\nabla \delta {\mathbf{u}}_k,\nabla \mathbf{v})+({\mathbf{u}}_k\cdot \nabla \delta {\mathbf{u}}_k,\mathbf{v})+\gamma (\mathrm{div}\delta {\mathbf{u}}_k,\mathrm{div}\mathbf{v}),$$

$$R_b(\mathbf{\phi })=-S{R}_m^{-1}(\mathrm{curl}{\mathbf{B}}_k,\mathrm{curl}\mathbf{\phi })-(\nabla r_k,\mathbf{\phi })-S({\mathbf{B}}_k\times {\mathbf{u}}_k,\mathrm{curl}\mathbf{\phi }),$$

Full text

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

[TABLE]

where ${\boldsymbol{u}}$ is the velocity of the fluid, $p$ is the hydrodynamic pressure, ${\boldsymbol{B}}$ is the magnetic flux density or the magnetic field provided with constant permeability, $R_{e}$ is the fluid Reynolds number, $R_{m}$ is the magnetic Reynolds number, $S$ is the coupling constant concerning the Lorentz force, and ${\boldsymbol{f}}\in{\boldsymbol{L}}^{2}({\Omega})$ stands for the external force. We assume that $\Omega$ is a bounded Lipschitz domain. The system of equations are complemented with Dirichlet boundary conditions

[TABLE]

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 ${\boldsymbol{H}}^{1}({\Omega})$-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 ${\boldsymbol{H}}^{1}({\Omega})$. It is preferable to use noncontinuous finite element functions to approximate ${\boldsymbol{B}}$, 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 $R_{e}$ and large coupling number $S$. An efficient solver should possess two merits:

the convergence rate is independent of the mesh or the number of degrees of freedom (DOFs); 2. 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

[TABLE]

where ${\boldsymbol{w}}$ 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 ${\boldsymbol{L}}^{2}({\Omega}):=(L^{2}({\Omega}))^{3}$.

2 Mixed finite element method for the MHD equations

First we introduce some Hilbert spaces and Sobolev norms used in this paper. Let $L^{2}(\Omega)$ be the usual Hilbert space of square integrable functions equipped with the following inner product and norm:

[TABLE]

Let the quotient space of $L^{2}({\Omega})$ be defined by

[TABLE]

Define $H^{m}(\Omega):=\{v\in L^{2}(\Omega):D^{\xi}v\in L^{2}(\Omega),|\xi|\leq m\}$ where $\xi$ represents non-negative triple index. Let $H^{1}_{0}(\Omega)$ be the subspace of $H^{1}(\Omega)$ whose functions have zero traces on $\Gamma$.

We define the spaces of functions having square integrable curl by

[TABLE]

which are equipped with the following inner product and norm

[TABLE]

Here ${\boldsymbol{n}}$ denotes the unit outer normal to $\Gamma$.

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

[TABLE]

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

[TABLE]

This means $r=0$ in $\Omega$ actually. Note that $\operatorname{div}{\boldsymbol{u}}=0$, 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 $({\boldsymbol{u}},{\boldsymbol{B}})\in{\boldsymbol{H}}^{1}({\Omega})\times\boldsymbol{H}(\operatorname{{\bf curl}},{\Omega})$ and $(p,r)\in L^{2}_{0}(\Omega)\times H_{0}^{1}({\Omega})$ such that ${\boldsymbol{u}}={\boldsymbol{g}}$ and ${\boldsymbol{B}}\times{\boldsymbol{n}}={\boldsymbol{B}}_{s}\times{\boldsymbol{n}}$ on $\Gamma$ and

[TABLE]

for all $({\boldsymbol{v}},\varphi)\in{\boldsymbol{H}}_{0}^{1}({\Omega})\times\boldsymbol{H}_{0}(\operatorname{{\bf curl}},{\Omega})$ and $(q,s)\in L^{2}_{0}(\Omega)\times H_{0}^{1}({\Omega})$, where the bilinear forms and trilinear form are defined respectively by

[TABLE]

Assuming small data, Schötzau proved the existence and uniqueness of the solution to (5) without the penalized term $\gamma(\operatorname{div}{\boldsymbol{u}},\operatorname{div}{\boldsymbol{v}})$. 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 $\mathcal{T}_{h}$ be a quasi-uniform and shape-regular tetrahedral mesh of $\Omega$. Let $h$ denote the maximal diameter of all tetrahedra on the mesh. For any $T\in\mathcal{T}_{h}$, let $P_{k}(T)$ be the space of polynomials of degree $k\geq 0$ on $K$ and ${\boldsymbol{P}}_{k}(T)=\left(P_{k}(T)\right)^{3}$ be the corresponding space of vector polynomials. Define the Lagrange finite element space of the $k$-th order by

[TABLE]

First we choose the well-known Taylor-Hood $P_{2}$-$P_{1}$ elements [3, Page 217-219] for the discretization of $({\boldsymbol{u}},p)$, namely,

[TABLE]

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

[TABLE]

where $C_{u}$ is the inf-sup constant independent of the mesh size. We shall also use $\overline{\boldsymbol{V}}_{h}=V(k,\mathcal{T}_{h})^{3}$.

The finite element space for ${\boldsymbol{B}}$ is chosen as Nédélec’s edge element space of the first order in the second family[25], namely,

[TABLE]

The finite element space for $r$ is defined by

[TABLE]

Since $\nabla S_{h}\subset{\boldsymbol{C}}_{h}$, we easily get the inf-sup condition for the pair of finite element spaces ${\boldsymbol{C}}_{h}\times S_{h}$

[TABLE]

where $C_{b}>0$ is the Poincáre constant depending only on $\Omega$.

The finite element approximation to (5) reads: Find $({\boldsymbol{u}}_{h},{\boldsymbol{B}}_{h})\in\overline{\boldsymbol{V}}_{h}\times\overline{\boldsymbol{C}}_{h}$ and $(p_{h},r_{h})\in Q_{h}\times S_{h}$ such that

[TABLE]

for all $({\boldsymbol{v}},\boldsymbol{\varphi})\in{\boldsymbol{V}}_{h}\times{\boldsymbol{C}}_{h}$ and $(q,s)\in Q_{h}\times S_{h}$. From (6) and (7) we know that the bilinear form $\mathcal{B}(\cdot,\cdot)$ satisfies the discrete inf-sup condition

[TABLE]

where

[TABLE]

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 $({\boldsymbol{B}}_{h},r_{h},{\boldsymbol{u}}_{h},p_{h})$ in the linearized problem. Let $({\boldsymbol{B}}_{k},r_{k},{\boldsymbol{u}}_{k},p_{k})\in\overline{\boldsymbol{C}}_{h}\times S_{h}\times\overline{\boldsymbol{V}}_{h}\times Q_{h}$ be the approximate solutions of (8) from the previous iteration. The error equation for these approximate solutions reads: Find $(\delta{\boldsymbol{B}}_{k},\delta r_{k},\delta{\boldsymbol{u}}_{k},\delta p_{k})\in{\boldsymbol{V}}_{h}\times Q_{h}\times{\boldsymbol{C}}_{h}\times S_{h}$ such that

[TABLE]

where the trilinear form $\mathcal{F}$ represents the convection-diffusion part of the fluid equation

[TABLE]

and the residual functionals are defined by

[TABLE]

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

[TABLE]

with a relaxation factor $0<\theta\leq 1$.

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

[TABLE]

where the solution vector $\mathbf{x}$ consists of the degrees of freedom for $(\delta{\boldsymbol{B}}_{k},\delta r_{k},\delta{\boldsymbol{u}}_{k},\delta p_{k})$ respectively, $\mathbf{b}$ is the residual vector, and $\mathbb{A}$ is the stiffness matrix. In block forms, they can be written as

[TABLE]

Let $\{{\boldsymbol{v}}_{i}:\;1\leq i\leq N_{V}\}$, $\{\boldsymbol{\varphi}_{i}:\;1\leq i\leq N_{C}\}$, $\{q_{i}:\;1\leq i\leq N_{Q}\}$, $\{s_{i}:\;1\leq i\leq N_{S}\}$ be the bases of ${\boldsymbol{V}}_{h}$, ${\boldsymbol{C}}_{h}$, $Q_{h}$, and $S_{h}$ respectively. Then the entries of all block matrices are defined by

[TABLE]

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

[TABLE]

Here $-\operatorname{div}$ is understood as the dual operator of $\nabla|_{S_{h}}$ or $\nabla|_{Q_{h}}$. Moreover, $\mathbb{J},\mathbb{J}^{\top}$ are algebraic representations of the two multiplication operators which couple the magnetic field and the conducting fluid. For any given ${\boldsymbol{w}}\in{\boldsymbol{V}}_{h}$ and $\boldsymbol{\psi}\in{\boldsymbol{C}}_{h}$, we have

[TABLE]

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 $\mathbb{L}_{r}$ be the stiffness matrix of $-\Delta$ on $S_{h}$ and let $\sigma>0$ be a constant. First we post-multiply the second column of $\mathbb{A}$ by $\sigma\mathbb{L}_{r}^{-1}\mathbb{G}$ and add it to the first column. This yields a matrix

[TABLE]

Next, pre-multiplying the first row of $\mathbb{A}_{1}$ by $-\mathbb{G}(\mathbb{C}+\sigma\mathbb{S}_{r})^{-1}$ and adding it to the second row, we get a matrix

[TABLE]

Note that $\mathbb{G}(\mathbb{C}+\sigma\mathbb{S}_{r})$ represents the operator $\operatorname{div}\left(SR_{m}^{-1}\operatorname{{\bf curl}}\operatorname{{\bf curl}}+\sigma\nabla\Delta^{-1}\operatorname{div}\right)$. Since

[TABLE]

formally we have

[TABLE]

This means that

[TABLE]

Since $\mathbb{J}^{\top}$ represents the coupling term $\operatorname{{\bf curl}}(S{\boldsymbol{B}}_{k}\times{\boldsymbol{v}}_{h})$, we have $\mathbb{G}\mathbb{J}^{\top}\approx 0$. Therefore,

[TABLE]

Moreover, from [17], we know that $\mathbb{C}+\sigma\mathbb{M}$ is equivalent to $\mathbb{C}+\sigma\mathbb{G}^{\top}\mathbb{L}_{r}^{-1}\mathbb{G}$ in spectrum where $\mathbb{M}$ is the mass matrix on ${\boldsymbol{C}}_{h}$. So one gets the approximation

[TABLE]

Next, pre-multiplying the first row of $\mathbb{A}_{3}$ by $\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}$ and adding it to the third row, we get a matrix

[TABLE]

Since $\operatorname{div}\left(SR_{m}^{-1}\operatorname{{\bf curl}}\operatorname{{\bf curl}}+\sigma{\mathbf{I}}\right)=\sigma\operatorname{div}$, formally we have $\operatorname{div}\left(SR_{m}^{-1}\operatorname{{\bf curl}}\operatorname{{\bf curl}}+\sigma{\mathbf{I}}\right)^{-1}=\sigma^{-1}\operatorname{div}$. Here ${\mathbf{I}}$ is the identity operator. This means

[TABLE]

So one obtains an approximation of $\mathbb{A}_{4}$ as follows

[TABLE]

This means that $\mathbb{A}_{5}^{-1}$ is actually a natural preconditioner for $\mathbb{A}$, except for the difficulties in computing the inverse of the block $\mathbb{F}+\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}\mathbb{J}^{\top}$. In the next subsection, we are going to derive a good approximation of $\mathbb{F}+\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}\mathbb{J}^{\top}$ so that its approximation inverse is easy to compute iteratively.

3.3 An efficient preconditioner for the magnetic field-fluid coupling block

From (35), the key step to compute $\mathbb{A}_{5}^{-1}$ is how to precondition the $2\times 2$ block

[TABLE]

We note that $\mathbb{F}+\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}\mathbb{J}^{\top}$ is the precise Schur complement of the following matrix which accounts for the coupling between ${\boldsymbol{B}}_{h}$ and ${\boldsymbol{u}}_{h}$

[TABLE]

Remember from (26) that $\mathbb{J}^{\top}$ and $\mathbb{J}$ represent, respectively, the two multiplication operators $S\operatorname{{\bf curl}}({\boldsymbol{B}}_{k}\times{\boldsymbol{w}})$ and $S\operatorname{{\bf curl}}\boldsymbol{\psi}\times{\boldsymbol{B}}_{k}$ for any given ${\boldsymbol{w}}$ and $\boldsymbol{\psi}$. So $\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}\mathbb{J}^{\top}{\boldsymbol{w}}$ is the algebraic representation of

[TABLE]

Since $(SR_{m}^{-1}\operatorname{{\bf curl}}\operatorname{{\bf curl}}+\sigma\mathbf{I})^{-1}$ commutates with $\operatorname{{\bf curl}}$, then (42) becomes

[TABLE]

For $\sigma>0$ sufficiently small, we adopt the following approximation

[TABLE]

This yields an approximation of (42)

[TABLE]

Therefore, we get an approximation of the Schur complement

[TABLE]

where $\mathbb{S}$ is the stiffness matrix associated with the bilinear form

[TABLE]

This should yield a good preconditioner of $\mathbb{X}$, that is,

[TABLE]

Remark 3.1

In [31], Philips et al studied the preconditioner for two-dimensional MHD equations where the magnetic field is discretized with $H^{1}({\Omega})$-conforming finite elements. Based on an exact penalty formulation, they propose to approximate the coupling effect between magnetic field and velocity by

[TABLE]

where $\beta>0$ is a parameter depending on the mesh size and the magnitude of ${\boldsymbol{B}}_{k}$.

Compared with [31], for the 3D MHD equations and the $\boldsymbol{H}(\operatorname{{\bf curl}},{\Omega})$-conforming approximation of ${\boldsymbol{B}}$, 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 (44) with respect to the parameter $\sigma$ and the mesh size $h$. Since it is the approximation $\mathbb{J}(\mathbb{C}+\sigma\mathbb{M})^{-1}\mathbb{J}^{\top}\approx\mathbb{S}$ that is concerned here, we fix $R_{e}=1.0$ and $\gamma=1.2$ and test the efficiency of the preconditioner for different values of $S$ and $R_{m}$.

We consider the system of linear equations: Find $\delta{\boldsymbol{u}}\in{\boldsymbol{V}}_{h}$ and $\delta{\boldsymbol{B}}\in{\boldsymbol{C}}_{h}$ such that

[TABLE]

where

[TABLE]

Clearly the stiffness matrix of (45) is $\hat{\mathbb{X}}$ which is given in (41). In the following, we test three cases of $\sigma$

[TABLE]

and three cases of physical parameters

[TABLE]

The computational domain is the unit cube, namely, $\Omega=(0,1)^{3}$.

We use preconditioned GMRES method to solve (45) and the preconditioner is set by (44). This means that we need solve the residual equation at each GMRES iteration

[TABLE]

where ${\mathbf{r}}_{b},{\mathbf{r}}_{u}$ stand for the residual vectors and ${\mathbf{e}}_{b},{\mathbf{e}}_{u}$ stand for the error vectors. The tolerance for the relative residual of the GMRES method is set by $10^{-6}$. The tolerances for solving the two sub-problems in (46) are set by $10^{-3}$. From Table 1–3, we find that the convergence of the preconditioned GMRES is uniform with respect to both $\sigma$ and $h$. An interesting observation is that, for large $S=R_{m}$, the number of GMRES iterations even decreases when $h\to 0$. In this case, the magnetic field-fluid coupling becomes strong.

Table 4 shows the number of preconditioned GMRES iterations where the approximate Schur complement $\mathbb{S}$ in (44) is replaced with the matrix $\mathbb{F}$. It amounts to devise a preconditioner of $\mathbb{X}$ by dropping its left lower block

[TABLE]

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 $S=R_{m}$, the convergence of GMRES method with this preconditioner becomes slower and deteriorates when $h\to 0$. 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 (35) and (44), we get a preconditioner of $\mathbb{A}$, that is, the inverse of

[TABLE]

It is left to study the preconditioner for the lower right $2\times 2$ block of $\mathbb{A}_{6}$, namely,

[TABLE]

It amounts to solve the saddle point problem: Find $(\delta{\boldsymbol{u}},\delta p)\in{\boldsymbol{V}}_{h}\times Q_{h}$ such that

[TABLE]

In [1, 2], Benzi et al studied the Osceen equation(namely ${\boldsymbol{B}}_{k}=\mathbf{0}$) and proposed to use the following preconditioner

[TABLE]

where $\mathbb{Q}_{p}$ is the mass matrix on $Q_{h}$. It is proved that the above preconditioner is efficient for relatively large Reynolds number. With ${\boldsymbol{u}}_{k}={\boldsymbol{B}}_{k}=\mathbf{0}$, namely the Stokes equations, we refer to [22, 23] for similar arguments. And with ${\boldsymbol{u}}_{k}=\mathbf{0}$ 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

[TABLE]

3.5 A robust preconditioner for the linearized MHD problem

Using (56) and (63), a preconditioner for $\mathbb{A}$ is given by the inverse of

[TABLE]

where the parameter $\gamma$ can be used to tune the efficiency of the preconditioner. According to our experience, any value $\gamma\sim O(1)$ works well for high Reynolds and moderate $S$ and $R_{m}$. Moreover, to fix the parameter $\sigma$, we set $\sigma=SR^{-1}_{m}$ so that $\mathbb{C}+\sigma\mathbb{M}$ is associated with the bilinear form

[TABLE]

This yields our final preconditioner for the stiffness matrix $\mathbb{A}$, defined by

[TABLE]

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 $\mathbb{P}$ to precondition $\mathbb{A}$. For convenience in notation, given a vector ${\mathbf{x}}$ which has the same size as one column vector of $\mathbb{A}$, we let $({\mathbf{x}}_{b},{\mathbf{x}}_{r},{\mathbf{x}}_{u},{\mathbf{x}}_{p})$ be the vectors which consist of entries of ${\mathbf{x}}$ and correspond to $({\boldsymbol{B}}_{h},r_{h},{\boldsymbol{u}}_{h},p_{h})$ respectively.

Algorithm 3.2** (Preconditioned GMRES Algorithm)**

Given the tolerances $\varepsilon\in(0,1)$ and $\varepsilon_{0}\in(\varepsilon,1)$, the maximal number of GMRES iterations $N>0$, and the initial guess $\mathbf{x}^{(0)}$ for the solution of (12). Set $k=0$ and compute the residual vector

[TABLE]

While $\left(k<N\;\;\&\;\;\left\|{\mathbf{r}^{(k)}}\right\|_{2}>\varepsilon\left\|{\mathbf{r}^{(0)}}\right\|_{2}\right)$ do

Solve $\mathbb{Q}_{p}{\mathbf{e}}_{p}=-\left(R_{e}^{-1}+\gamma\right){\mathbf{r}}_{p}^{(k)}$ by the CG method with the diagonal preconditioning and tolerance $\varepsilon_{0}$ for the relative residual. 2. 2.

Solve $\mathbb{S}{\mathbf{e}}_{u}={\mathbf{r}}^{(k)}_{u}-\mathbb{B}^{\top}{\mathbf{e}}_{p}$ by preconditioned GMRES method with tolerance $10^{-3}$. The preconditioner is the one level additive Schwarz method with overlap = 2 [4]. 3. 3.

Solve $\mathbb{L}_{r}{\mathbf{e}}_{r}=-SR^{-1}_{m}{\mathbf{r}}^{(k)}_{r}$ by preconditioned CG method with tolerance $\varepsilon_{0}$. The preconditioner is the algebraic multigrid method (AMG) solver [9]. 4. 4.

Solve $(\mathbb{C}+SR^{-1}_{m}\mathbb{M}){\mathbf{e}}_{b}={\mathbf{r}}^{(k)}_{b}-\mathbb{J}^{\top}{\mathbf{e}}_{u}-\mathbb{G}^{\top}{\mathbf{e}}_{r}$ by preconditioned CG method with tolerance $\varepsilon_{0}$ and the Hiptmair-Xu preconditioner [12]. 5. 5.

Update the solution: $\mathbf{x}^{(k+1)}:=\mathbf{x}^{(k)}+\mathbf{e}^{(k)}$ . 6. 6.

Set $k:=k+1$ and compute the residual vector $\mathbf{r}^{(k)}=\mathbf{b}-\mathbb{A}\mathbf{x}^{(k)}$.

End while.**

4 Numerical experiments

In this section, we present three numerical experiments to verify the convergence rate of finite element approximation to the augmented Lagrangian(AL) formulation of the MHD, to demonstrate the robustness of the preconditioner, and to demonstrate the scalability of the parallel solver. The parallel code is developed based on the finite element package—Parallel Hierarchical Grids (PHG) [35, 36].

example 4.1

This example is to verify the convergence rate of the finite element discrete problem (8). The analytic solutions are chosen as

[TABLE]

The parameters are set by $R_{e}=S=R_{m}=1$ and $\gamma=1$.

From Table 5, we find that the convergence rates for ${\boldsymbol{u}}_{h}$, $p_{h}$, ${\boldsymbol{B}}_{h}$ are given by

[TABLE]

Remember that we are using the second-order Lagrangian finite elements for discretizing ${\boldsymbol{u}}$, the first-order Lagrangian finite elements for discretizing $p$, and the first-order Nédélec’s edge elements in the second family for discretizing ${\boldsymbol{B}}$. This means that the optimal convergence rates are obtained for all variables.

example 4.2** (Driven Cavity Flow)**

The example is the benchmark problem of a driven cavity flow. The righthand side is given by ${\boldsymbol{f}}=\mathbf{0}$ and the boundary conditions are given by ${\boldsymbol{g}}=(g_{1},0,0)^{\top}$ and ${\boldsymbol{B}}_{s}=(1,0,0)^{\top}$ where $g_{1}=g_{1}(z)$ is a continuous function and satisfies

[TABLE]

where $h$ is the mesh size. The parameters are set by

[TABLE]

The purpose of this example is to verify the effectiveness of the mixed finite element method for engineering benchmark problem and demonstrate the robustness of the discrete solver with respect to the mesh size $h$. The computational domain is set by $\Omega=(0,1)^{3}$. We set the tolerances by $\varepsilon=10^{-6}$ and $\varepsilon_{0}=10^{-3}$ in Algorithm 3.2. The tolerance is $10^{-4}$ for the relative residual of the nonlinear iterations. Table 6 shows the mesh sizes and the numbers of DOFs on the meshes which we are using.

Let $N_{\rm picard}$ denote the number of nonlinear iterations to reduce the relative residual by a factor $10^{-4}$. Let $N_{\rm gmres}$ denote the number of preconditioned GMRES iterations for solving the linearized problem (12). Therefore, $N_{\rm picard}\times N_{\rm gmres}$ represents the total computational quantity for solving the nonlinear problem (8). Remember that the approximate Schur complement $\mathbb{S}$ in (65) is defined by the bilinear form

[TABLE]

As mentioned in the last paragraph of Subsection 3.3, dropping the second term gives $\mathbb{S}=\mathbb{F}$. In Table 7, we show the effectiveness of the preconditioner $\mathbb{P}$ with and without the term $\mathbf{BuBv}:=SR_{m}({\boldsymbol{B}}_{k}\times{\boldsymbol{u}},{\boldsymbol{B}}_{k}\times{\boldsymbol{v}})$ in $\mathbb{S}$. An interesting observation is that, with $\mathbf{BuBv}$, the number of GMRES decays when the mesh is refined successively. However, without this term, the number of GMRES iterations increases considerably.

Now we give the visualization of the simulation results. Firstly, Fig. 1 give the grayscale figure of the magnitude of ${\boldsymbol{u}}_{h}$. Three 2D projection velocity streamlines of ${\boldsymbol{u}}_{h}$ are show in Fig. 2. Fig. 3 shows the contour of the pressure $p_{h}$.

Finally, Fig. 4 shows the distribution of of $\left|{{\boldsymbol{B}}_{h}}\right|$ on three cross-sections of $\Omega$ at $x=0.5$, $y=0.5$, and $z=0.5$ respectively.

example 4.3** (Scalability)**

This example also computes the driven cavity flow and is used to test the scalability of the solver for moderate parameters. The physical parameters $R_{e}=S=R_{m}=1$, and $\gamma=0.1$.

The error tolerance of the nonlinear iteration is $10^{-4}$ relative to the initial residual. In step 2 of Algorithm 3.2, we replace the one level additive Schwarz preconditioner with the BoomerAMG preconditioner [9]. This yields better parallel efficiency of the solver for moderate parameters. We carried out the computations on 5 successively refined meshes. Table 8 shows the scalability of the discrete solver by parallel computing on five successively refined meshes. The parallel efficiencies are above $35\%$ and good for such a complex problem.

For large Reynolds number, the main challenge for the scalability lies in step 2 of Algorithm 3.2, that is, the solution of the system of algebraic equations

[TABLE]

We should admit that neither the one level additive Schwarz method nor the classical AMG method can provide ideal scalability solely with our code. Multilevel-based preconditioning and stabilizations for the convection term should be resorted to. This will be our future work and will not be discussed here.