Consistency enforcement for the iterative solution of weak Galerkin finite element approximation of Stokes flow

We consider the lowest-order weak Galerkin (WG) finite element approximation of Stokes flow problems in the form

where $\Omega\subset\mathbb{R}^{d}$ $(d=2,3)$ is a bounded polygonal/polyhedral domain, $\mu>0$ is the fluid kinematic viscosity, $\mathbf{u}$ is the fluid velocity, $p$ is the fluid pressure, $\mathbf{f}$ is a body force, $\mathbf{g}$ is a non-zero boundary datum of the velocity satisfying the compatibility condition $\int_{\partial\Omega}\mathbf{g}\cdot\mathbf{n}=0$, and $\mathbf{n}$ is the unit outward normal to the boundary of the domain. The resulting linear algebraic system is a singular saddle point system where the singularity reflects the nonuniqueness of the pressure solution. While it is well known (e.g., see [4, Remark 6.12] and [16, Section 10.2]) that Krylov subspace methods, such as the minimal residual method (MINRES) and the generalized minimal residual method (GMRES), work well for consistent singular systems, the underlying system is nonconsistent in general when $g$ is not identically zero (cf. Section 2). This inconsistency can cause MINRES and GMRES to fail to converge.

The objective of this work is to study a simple consistency enforcement strategy by modifying the right-hand side of the linear system that results from the WG discretization of (1). We shall prove that the optimal-order convergence is not affected by the modification. Moreover, we consider inexact block diagonal and triangular Schur complement preconditioners for the efficient iterative solution of the singular but consistent saddle point system resulting from the modified scheme. Bounds for the eigenvalues and the residual of MINRES and GMRES for the corresponding preconditioned systems are established. These bounds show that the convergence of MINRES and GMRES with the inexact block diagonal and triangular Schur complement preconditioning is independent of the fluid kinematic viscosity $\mu$ and mesh size. It is worth mentioning that the preconditioned system associated a block diagonal Schur complement preconditioner is diagonalizable and this MINRES can be used for system solving and spectral analysis can be used to analyze the convergence of MINRES. On the other hand, the preconditioned system with a block triangular Schur complement preconditioner is not diagonalizable and we need to use GMRES for the iterative solution of the system. Moreover, the spectral analysis cannot be used directly to analyze the convergence of GMRES. Instead, it needs to be combined with Lemmas A.1 and A.2 of [7] that provide a bound for the residual of GMRES in terms of the norm of the off-diagonal blocks in the preconditioned system and the performance of GMRES to the preconditioned Schur complement.

Numerical solution of Stokes flow problems has continuously gained attention from researchers. Particularly, a variety of finite element methods have been studied for those problems; e.g., see [1] (mixed finite element methods), [3, 17] (virtual element methods), [8, 15] (hybrid discontinuous Galerkin methods), and [11, 14, 18] (weak Galerkin (WG) finite element methods). We use here the lowest-order WG method for the discretization of Stokes flow problems. It is known (cf. Lemma 2.1 or [19]) that the lowest-order WG method, without using stabilization terms, satisfies the inf-sup condition (for stability) and has the optimal-order convergence. Moreover, the error in the velocity is independent of the error in the pressure (pressure-robustness) and the error in the pressure is independent of the viscosity $\mu$ ($\mu$-semi-robustness). On the other hand, little work has been done so far for the efficient iterative solution of the saddle point system arising from the WG approximation of Stokes problems. Recently, the authors considered in [7] the iterative solution of the lowest-order WG approximation of Stokes problems through regularization and provided a convergence analysis for MINRES and GMRES with inexact block diagonal and triangular Schur complement preconditioners, respectively.

The rest of this paper is organized as follows. In Section 2, the weak formulation for Stokes flow and its discretization by the lowest-order WG method are described. Consistency enforcement and related error analysis are discussed in Section 3. In Section 4, the block diagonal and block triangular Schur complement preconditioning and convergence of MINRES and GMRES for the modified system are studied. Section 5 presents two- and three dimensional numerical experiments to verify the theoretical findings and showcase the effectiveness of the preconditioners. Finally, the conclusions are drawn in Section 6.

In this section we describe the lowest-order weak Galerkin approximation of the Stokes flow problem (1).

The weak formulation of (1) is to find $\mathbf{u}\in H^{1}(\Omega)^{d}$ and $p\in L^{2}(\Omega)$ such that $\mathbf{u}|_{\partial\Omega}=\mathbf{g}$ (in the weak sense) and

Consider a connected quasi-uniform simplicial mesh $\mathcal{T}_{h}=\{K\}$ on $\Omega$. A mesh is called to be connected if any pair of its elements is linked by a chain of elements that share an interior facet with each other. We introduce the following discrete weak function spaces on $\Omega$:

where $P_{0}(K)$ and $P_{0}(e)$ denote the set of constant polynomials defined on element $K$ and facet $e$, respectively. Note that $\mathbf{u}_{h}\in\mathbf{V}_{h}$ is approximated on both interiors and facets of mesh elements while $p_{h}\in W_{h}$ is approximated only on element interiors. Then, we can define the lowest-order WG approximation of the Stokes problem (1) as: finding $\mathbf{u}_{h}\in\mathbf{V}_{h}$ and $p_{h}\in W_{h}$ such that $\mathbf{u}_{h}^{\partial}|_{\partial\Omega}=Q_{h}^{\partial}\mathbf{g}$ and

where $Q_{h}^{\partial}$ is a $L^{2}$-projection operator onto $\mathbf{V}_{h}$ restricted on each facet and the lifting operator $\mathbf{\Lambda}_{h}:\mathbf{V}_{h}\to RT_{0}(\mathcal{T}_{h})$ is defined [10, 19] as

Here, $RT_{0}(K)$ is the lowest-order Raviart-Thomas space, viz.,

Notice that $\mathbf{\Lambda}_{h}\mathbf{v}$ depends on $\mathbf{v}^{\partial}$ but not on $\mathbf{v}^{\circ}$.

The discrete weak gradient and divergence operators in (5) are defined as follows. For a scalar function or a component of a vector-valued function, $u_{h}=(u_{h}^{\circ},u_{h}^{\partial})$, the discrete weak gradient operator $\nabla_{w}:W_{h}\rightarrow RT_{0}(\mathcal{T}_{h})$ is defined as follows

where $\mathbf{n}$ is the unit outward normal to $\partial K$ and $(\cdot,\cdot)_{K}$ and $(\cdot,\cdot)_{\partial K}$ are the $L^{2}$ inner product on $K$ and $\partial K$, respectively. For a vector-valued function $\mathbf{u}_{h}$, $\nabla_{w}\mathbf{u}_{h}$ is viewed as a matrix with each row representing the weak gradient of a component. By choosing $\mathbf{w}$ properly in (7) and using the fact that $\nabla_{w}u_{h}\in RT_{0}(K)$, we can obtain (e.g., see [6])

where $\varphi_{K}^{\circ}$ and $\varphi_{K,i}^{\partial}$ denote the basis functions of $P_{0}(K)$ and $P_{0}(e_{K,i})$, respectively, $e_{K,i}$ denotes the $i$-th facet of $K$, $\mathbf{n}_{K,i}$ is the unit outward normal to $e_{K,i}$,

and $\mathbf{x}_{K,i}$, $i=1,...,d+1$ denote the vertices of $K$.

The discrete weak divergence operator $\nabla_{w}\cdot:\mathbf{V}_{h}\to\mathcal{P}_{0}(\mathcal{T}_{h})$ is defined independently through

Note that $\nabla_{w}\cdot\mathbf{u}|_{K}\in P_{0}(K)$. Moreover, by taking $w=1$ we have

where $\langle\mathbf{u}\rangle_{e_{K,i}}$ denotes the average of $\mathbf{u}$ on facet $e_{K,i}$ and $|e_{K,i}|$ is the $(d-1)$-dimensional measure of $e_{K,i}$.

The scheme (5) achieves the optimal-order convergence as shown in the following lemma.

To cast (5) into a matrix-vector form, hereafter we use $\mathbf{v}_{h}$ interchangeably for the WG approximation of $\mathbf{v}$ in $\mathbf{V}_{h}$ and the vector formed by its values $(\mathbf{v}_{h,K}^{\circ},\mathbf{v}_{h,K,i}^{\partial})$ for $i=1,...,d+1$ and $K\in\mathcal{T}_{h}$, excluding those on $\partial\Omega$. Similarly, $\mathbf{q}_{h}$ is used to denote both the WG approximation of $q$ and the vector formed by $q_{h,K}$ for all $K\in\mathcal{T}_{h}$. Then, we can write (5) into

where the matrices $A$ and $B^{\circ}$ and vectors $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ are defined as

Notice that (16) is a singular saddle point system where the pressure solution is unique up to a constant. Moreover, the system is inconsistent in general when $\mathbf{g}$ is not identically zero. To explain this, from (20) we have

This gives rise to

where we have used the fact that $\partial\Omega$ consists of element facets because $\Omega$ is assumed to be a polygon or a polyhedron. If we choose $Q_{h}^{\partial}\mathbf{g}|_{e}=\langle\mathbf{g}\rangle_{e}$, $\forall e\in\partial\Omega$, from the compatibility condition we have

Unfortunately, in general the averages $\langle\mathbf{g}\rangle_{e}$ need to be approximated numerically. The approximation error can cause a non-zero $\sum_{K\in\mathcal{T}_{h}}\mathbf{b}_{2}(K)$ and thus the inconsistency of the system. This can also happen if $Q_{h}^{\partial}$ is defined differently, for instance, $Q_{h}^{\partial}\mathbf{g}|_{e}$ is defined as the value of $\mathbf{g}$ at the barycenter of $e$. Krylov subspace methods, like MINRES and GMRES, may fail to converge when applied to this inconsistent singular system although they are known to work well for consistent singular systems, at least when the initial guess is taken as zero; e.g., see [4, Remark 6.12] and [16, Section 10.2].

In this section we consider an approach of enforcing the consistency/compatibility condition by modifying the right-hand term $\mathbf{b}_{2}$ and study the effects of this modification on the optimal-order convergence of the numerical solution.

We propose to modify $\mathbf{b}_{2}$ to enforce the consistency. Specifically, we define

where $N$ is the number of elements in $\mathcal{T}_{h}$ and

By definition, we have $\sum_{K}\tilde{\mathbf{b}}_{2}(K)=0$. Moreover, from the compatibility condition $\int_{\partial\Omega}\mathbf{g}\cdot\mathbf{n}dS=0$, we can rewrite $\alpha_{h}$ as

The modified scheme reads as

where $A$, $B^{\circ}$, and $\mathbf{b}_{1}$ are given in (17), (18), and (19), respectively. Note that (25) is still singular but consistent.

An immediate question about the modified scheme is how much the accuracy of the numerical solution is affected by the modification. To answer this, we provide an error analysis for the modified scheme (25) in the following.

Let $\mathbf{u}\in H^{2}(\Omega)^{d}$ and $p\in H^{1}(\Omega)$ be the solutions of the Stokes problem (1) and $\mathbf{u}_{h}\in\mathbf{V}_{h}$ and $p_{h}\in W_{h}$ be the numerical solutions of (25). We split the error into the projection and discrete error as

where the $L^{2}$ projection operator is defined as $Q_{h}=(Q_{h}^{\circ},Q_{h}^{\partial})$. The operator $Q_{h}$ is considered to be componentwise when applied to vector-valued functions. The projection error $\mathbf{u}-Q_{h}\mathbf{u}$ and $p-Q_{h}p$ are determined by the approximation capacity of the finite element spaces. Our primary focus is on $\mathbf{e}_{h}$ and $e_{h}^{p}$.