Local Fourier analysis for mixed finite-element methods for the Stokes equations

Local Fourier analysis for mixed finite-element methods for the Stokes equations

Yunhui He yunhui.he@mun.ca Scott P. MacLachlan smaclachlan@mun.ca Department of Mathematics and Statistics, Memorial University of Newfoundland, St. John’s, NL A1C 5S7, Canada

In this paper, we develop a local Fourier analysis of multigrid methods based on block-structured relaxation schemes for stable and stabilized mixed finite-element discretizations of the Stokes equations, to analyze their convergence behavior. Three relaxation schemes are considered: distributive, Braess-Sarazin, and Uzawa relaxation. From this analysis, parameters that minimize the local Fourier analysis smoothing factor are proposed for the stabilized methods with distributive and Braess-Sarazin relaxation. Considering the failure of the local Fourier analysis smoothing factor in predicting the true two-grid convergence factor for the stable discretization, we numerically optimize the two-grid convergence predicted by local Fourier analysis in this case. We also compare the efficiency of the presented algorithms with variants using inexact solvers. Finally, some numerical experiments are presented to validate the two-grid and multigrid convergence factors.

Monolithic multigrid, Block-structured relaxation, local Fourier analysis, mixed finite-element methods, Stokes Equations
65N55, 65F10, 65F08, 76M10

1 Introduction

In recent years, substantial research has been devoted to efficient numerical solution of the Stokes and Navier-Stokes equations, due both to their utility as models of (viscous) fluids and their commonalities with many other physical problems that lead to saddle-point systems (see, for example elman2006finite (), and many of the other references cited here). In the linear (or linearized) case, solution of the resulting matrix equations is seen to be difficult, due to indefiniteness and the usual ill-conditioning of discretized PDEs. In the literature, block preconditioners (cf. elman2006finite () and the references therein) are widely used, due to their easy construction from standard multigrid algorithms for scalar elliptic PDEs, such as algebraic multigrid ruge1987algebraic (). However, monolithic multigrid approaches adler2016monolithic (); MR3439774 (); MR558216 (); MR2833487 (); MR2001083 () have been shown to outperform these preconditioners when algorithmic parameters are properly chosen MR3639325 (); farrell2018augmented (). The focus of this work is on the analysis of such monolithic multigrid methods in the case of stable and stabilized finite-element discretizations of the Stokes equations.

Local Fourier analysis (LFA) MR1807961 (); wienands2004practical () has been widely used to predict the convergence behavior of multigrid methods, to help design relaxation schemes and choose algorithmic parameters. In general, the LFA smoothing factor provides a sharp prediction of actual multigrid convergence, see MR1807961 (), under the assumption of an “ideal” coarse-grid correction scheme (CGC) that annihilates low-frequency error components and leaves high-frequency components unchanged. In practice, the LFA smoothing and two-grid convergence factors often exactly match the true convergence factor of multigrid applied to a problem with periodic boundary conditions brandt1994rigorous (); stevenson1990validity (); MR1807961 (). Recently, the validity of LFA has been further analysed rodrigo2017validity (), extending this exact prediction to a wider class of problems. However, the LFA smoothing factor is also known to lose its predictivity of the true convergence in some cases HM2018LFALaplace (); MR2840198 (); friedhoff2013local (). In particular, the smoothing factor of LFA overestimates the two-grid convergence factor for the Taylor-Hood () discretization of the Stokes equations with Vanka relaxation MR2840198 (). Even for the scalar Laplace operator, the LFA smoothing factor fails to predict the observed multigrid convergence factor for higher-order finite-element methods HM2018LFALaplace ().

Two main questions interest us here. First, we look to extend the study of MR2840198 () to consider LFA of block-structured relaxation schemes for finite-element discretizations of the Stokes equations. Secondly, we consider if the LFA smoothing factor can predict the convergence factors for these relaxation schemes. Recently, LFA for multigrid based on block-structured relaxation schemes applied to the marker-and-cell (MAC) finite-difference discretization of the Stokes equations was shown to give a good prediction of convergence NLA2147 (), in contrast to the results of MR2840198 (). Thus, a natural question to investigate is whether the contrasting results between NLA2147 () and MR2840198 () is due to the differences in discretization or those in the relaxation schemes considered. Here, we apply the relaxation schemes of NLA2147 () to the discretization from MR2840198 (), as well as an “intermediate” discretization using stabilized approaches.

In recent decades, many block relaxation schemes have been studied and applied to many problems, including Braess-Sarazin-type relaxation schemes adler2016monolithic (); braess1997efficient (); MR1685458 (); MR3439774 (); MR1810326 (), Vanka-type relaxation schemes adler2016monolithic (); MR848451 (); MR3439774 (); MR2840198 (); MR2263039 (); MR3488076 (); MR2001083 (), Uzawa-type relaxation schemes john2016analysis (); MR3217219 (); MR3580776 (); MR2833487 (); MR833993 (), distributive relaxation schemes bacuta2011new (); MR558216 (); wittum1989multi (); MR3071182 (); oosterlee2006multigrid () and other types of methods MR3427891 (); MR3432850 (). Even though LFA has been applied to distributive relaxation MR1049395 (); wienands2004practical (), Vanka relaxation MR2840198 (); MR3488076 (); sivaloganathan1991use (); MR1131567 (), and Uzawa-type schemes MR3217219 () for the Stokes equations, most of the existing LFA has been for relaxation schemes using (symmetric) Gauss-Seidel (GS) approaches, and for simple finite-difference and finite-element discretizations. Considering modern multicore and accelerated parallel architectures, we focus on schemes based on weighted Jacobi relaxation with distributive, Braess-Sarazin, and Uzawa relaxation for common finite-element discretizations of the Stokes equations.

Some key conclusions of this analysis are as follows. First, while the LFA smoothing factor gives a good prediction of the true convergence factor for the stabilized discretizations with distributive weighted Jacobi and Braess-Sarazin relaxation, it does not for the Uzawa relaxation (in contrast to what is seen for the MAC discretization NLA2147 (); MR1049395 ()). For no cases does the LFA smoothing factor offer a good prediction of the true convergence behaviour for the (stable) discretization, suggesting that the discretization is responsible for the lack of predictivity, consistent with the results in HM2018LFALaplace (); MR2840198 (). For both stable and stabilized discretizations, we see that standard distributive weighted Jacobi relaxation loses some of its high efficiency, in contrast to what is seen for the MAC scheme NLA2147 (); MR1049395 () but that robustness can be restored with an additional relaxation sweep. Exact Braess-Sarazin relaxation is also highly effective, with LFA-predicted convergence factors of in the stabilized cases and in the stable case. To realize these rates with inexact cycles, however, requires nested W-cycles to solve the approximate Schur complement equation accurately enough in the stabilized case, although simple weighted Jacobi on the approximate Schur complement is observed to be sufficient in the stable case. For Uzawa-type relaxation, we see a notable gap between predicted convergence with exact inversion of the resulting Schur complement, versus inexact inversion, although some improvement is seen when replacing the approximate Schur complement with a mass matrix approximation, as is commonly used in block-diagonal preconditioners wathen2009chebyshev (); wathen1993fast (); silvester1994fast (). Overall, however, we see that distributive weighted Jacobi (DWJ) (with 2 sweeps of Jacobi relaxation on the pressure equation) outperforms both Braess-Sarazin relaxation (BSR) and Uzawa relaxation, for the stabilized discretizations, while DWJ and inexact BSR offer comparable performance for the stable discretization.

We organize this paper as follows. In Section 2, we introduce two stabilized and the stable mixed finite-element discretizations of the Stokes equations in two dimensions (2D). In Section 3, we first review the LFA approach, then discuss the Fourier representation for these discretizations. In Section 4, LFA is developed for DWJ, BSR, and Uzawa-type relaxation, and optimal LFA smoothing factors are derived for the two stabilized methods with DWJ and BSR. Multigrid performance is presented to validate the theoretical results. Section 5 exhibits optimized LFA two-grid convergence factors and measured multigrid convergence factors for the discretization. Furthermore, a comparison of the cost and effectiveness of the relaxation schemes is given. Conclusions are presented in Section 6.

2 Discretizations

In this paper, we consider the Stokes equations,


where is the velocity vector, is the (scalar) pressure of a viscous fluid, and represents a (known) forcing term, together with suitable boundary conditions. Because of the nature of LFA, we validate our predictions against the problem with periodic boundary conditions on both and . Discretizations of (1) typically lead to a linear system of the following form:


where corresponds to the discretized vector Laplacian, and is the negative of the discrete divergence operator. If the discretization is naturally unstable, then is the stabilization matrix, otherwise . In this paper, we discuss two stabilized and the stable finite-element discretizations.

The natural finite-element approximation of Problem (1) is: Find and such that



and , are finite-element spaces. Here, satisfies homogeneous Dirichlet boundary conditions in place of any non-homogenous essential boundary conditions on . Problem (3) has a unique solution only when and satisfy an inf-sup condition (see elman2006finite (); MR2373954 (); brezzi1988stabilized (); dohrmann2004stabilized ()).

2.1 Stabilized discretizations

The standard equal-order approximation of (3) is well-known to be unstable brezzi1988stabilized (); elman2006finite (). To circumvent this, a scaled pressure Laplacian term can be added to (3); for a uniform mesh with square elements of size , we subtract

for . With this, the resulting linear system is given by

where is the Laplacian operator for the pressure. Denote , and , where . From elman2006finite (), the red-black unstable mode , can be moved from a zero eigenvalue to a unit eigenvalue ( giving stability without loss of accuracy) by choosing so that


where is the mass matrix. Substituting the bilinear stiffness and mass matrices into (4), we find . We refer to this method as the Poisson-stabilized discretization (PoSD).

An projection to stabilize the discretization, proposed in dohrmann2004stabilized (), stabilizes with


where is the projection from into the space of piecewise constant functions on the mesh. We refer to this method as the projection stabilized discretization (PrSD). The element matrix of (5) is given by

where is the element mass matrix for the bilinear discretization and . In the projection stabilized method, we can write , where is given by the 9-point stencil

Applying (4) to , we find that is the optimal choice.

2.2 Stable discretizations

In order to guarantee the well-posedness of the discrete system (2) with , the discretization of the velocity and pressure unknowns should satisfy an inf-sup condition,

where is a constant. Taylor-Hood () elements are well known to be stable MR2373954 (); elman2006finite (), where the basis functions associated with these elements are biquadratic for each component of the velocity field and bilinear for the pressure.

3 LFA preliminaries

3.1 Definitions and notations

In many cases, the LFA smoothing factor offers a good prediction of multigrid performance. Thus, we will explore the LFA smoothing factor and true (measured) multigrid convergence for the three types of relaxations considered here. We first introduce some terminology of LFA, following MR1807961 (); wienands2004practical (). We consider the following two-dimensional infinite uniform grids,


The coarse grids, , are defined similarly.

Figure 1: At left, the mesh used for discretization. At right, the mesh used for discretization. Points marked by correspond to , those marked by correspond to , those marked by correspond to and those marked by 9  correspond to .

Let be a scalar Toeplitz operator defined by its stencil acting on as follows:


with constant coefficients , where is a function in . Here, is a finite index set. Because is formally diagonalized by the Fourier modes , where and , we use as a Fourier basis with . High and low frequencies for standard coarsening (as considered here) are given by

Definition 3.1.

If, for all functions ,

we call the symbol of .

In what follows, we consider linear systems of operators, which read


where depends on which discretization we use.

For the stabilized approximations, the degrees of freedom for both velocity and pressure are only located on as pictured at left of Figure 1. In this setting, the in (7) are scalar Toeplitz operators. Denote as the symbol of . Each entry in is computed as the (scalar) symbol of the corresponding block of , following Definition 3.1. Thus, is a matrix. All blocks in are diagonalized by the same transformation on a collocated mesh.

However, for the discretization, the degrees of freedom for velocity are located on , containing four types of meshpoints as shown at right of Figure 1. The Laplace operator in (7) is defined by extending (6), with taken to be a finite index set of values, with , , , and . With this, the (scalar) Laplace operator is naturally treated as a block operator, and the Fourier representation of each block can be calculated based on Definition 3.1, with the Fourier bases adapted to account for the staggering of the mesh points. Thus, the symbols of and are matrices. For more details of LFA for the Laplace operator using higher-order finite-element methods, refer to HM2018LFALaplace (). Similarly to the Laplace operator, both terms in the gradient, and , can be treated as ()-block operators. Then, the symbols of and are matrices, calculated based on Definition 3.1 adapted for the mesh staggering. The symbols of and are the conjugate transposes of those of and , respectively. Finally, . Accordingly, is a matrix for the discretization.

Definition 3.2.

The error-propagation symbol, , for a block smoother on the infinite grid satisfies

for all , and the corresponding smoothing factor for is given by

where is an eigenvalue of .

In Definition 3.2, for the stabilized case (and is a matrix) and for the stable case (where is a matrix).

The error-propagation symbol for a relaxation scheme, represented by matrix , applied to either the stabilized or stable scheme is written as

where represents parameters within , the block approximation to , is an overall weighting factor, and and are the symbols for and , respectively. Note that is a function of some parameters in Definition 3.2. In this paper, we focus on minimizing with respect to these parameters, to obtain the optimal LFA smoothing factor.

Definition 3.3.

Let be the set of allowable parameters and define the optimal smoothing factor over as

If the standard LFA assumption of an “ideal” CGC holds, then the two-grid convergence factor can be estimated by the smoothing factor, which is easy to compute. However, as expected, we will see that this idealized CGC does not lead to a good prediction for some cases we consider below. When the LFA smoothing factor fails to predict the true two-grid convergence factor, the LFA two-grid convergence factor can still be used. Thus, we give a brief introduction to the LFA two-grid convergence factor in the following.


We use the ordering of for the four harmonics. To apply LFA to the two-grid operator,


we require the representation of the CGC operator,

where is the multigrid interpolation operator and is the restriction operator. The coarse-grid operator, , can be either the Galerkin or rediscretization operator.

Inserting the representations of into (8), we obtain the Fourier representation of two-grid error-propagation operator as


in which stands for the block diagonal matrix with diagonal blocks, , and .

Here, we use the standard finite-element interpolation operators and their transposes for restriction. For , the symbol is well-known MR1807961 () while, for the nodal basis for , the symbol is given in HM2018LFALaplace ().

Definition 3.4.

The asymptotic two-grid convergence factor, , is defined as

In what follows, we consider a discrete form of , denoted by , resulting from sampling over only a finite set of frequencies. For simplicity, we drop the subscript throughout the rest of this paper, unless necessary for clarity.

3.2 Fourier representation of discretization operators

3.2.1 Fourier representation of the stabilized discretization

By standard calculation, the symbols of the stiffness and mass stencils are

respectively. The stencils of the partial derivative operators and are

respectively, and the corresponding symbols are

where denotes the conjugate transpose. Thus, the symbols of the stabilized finite-element discretizations of the Stokes equations are given by

For the Poisson-stabilized discretization, the symbol of is . For the projection stabilized method, following (5), the symbol of is


For convenience, we write for the last block of Equation (2), and its symbol as in the rest of this paper.

3.2.2 Fourier representation of stable discretizations

The symbols of the stiffness and mass stencils for the discretization using nodal basis functions in 1D are

respectively HM2018LFALaplace (). Here, we note that the entries correspond to the symbols associated with basis functions at the nodes of the mesh, while the entries correspond to the symbols associated with cell-centre (bubble) basis functions. The off-diagonal entries express the interaction between the two types of basis functions. Then, the Fourier representation of in 2D can be written as a tensor product,

The tensor product preserves block structuring; that is, is a matrix, ordered as mesh nodes, -edge midpoints, -edge midpoints, and cell centres. Each row of reflects the connections between one of the four types of degrees of freedom with each of these four types. Similarly, there are four types of stencils for and .

The stencils and the symbols of for the nodal, -edge, -edge, and cell-centre degrees of freedom are

respectively. Denote .

Similarly to , the symbol of the stencil of can be written as

Thus, the Fourier representation of the finite-element discretization of the Stokes equations can be written as


Note that the Fourier symbol for the discretization is a matrix, and that the LFA smoothing factor for the approximation generally fails to predict the true two-grid convergence factor HM2018LFALaplace (); MR2840198 (). The same behavior is seen for the relaxation schemes considered here. Therefore, we do not present smoothing factor analysis for this case and only optimize two-grid LFA predictions numerically.

4 Relaxation for discretizations

4.1 DWJ relaxation

Distributive GS (DGS) relaxation MR558216 (); oosterlee2006multigrid () is well known to be highly efficient for the MAC finite-difference discretization MR1807961 (), and other discretizations MR3427891 (); chen2015multigrid (). Its sequential nature is often seen as a significant drawback. However, Distributive weighted Jacobi (DWJ) relaxation was recently shown to achieve good performance for the MAC discretization NLA2147 (). Thus, we consider DWJ relaxation for the finite-element discretizations considered here. The discretized distribution operator can be represented by the preconditioner

Then, we apply blockwise weighted-Jacobi relaxation to the distributed operator


where we note that the operators and are formed by taking products of the discrete derivative operators and, thus, do not satisfy the identity .

The discrete matrix form of is

where is the Laplacian operator discretized at the pressure points. For standard distributive weighted-Jacobi relaxation (with weights ), we need to solve a system of the form


then distribute the updates as . We use in the block of (12), because the diagonal entries of the block will be of the form of a constant times (up to boundary conditions), for both stabilization terms. The error propagation operator for the scheme is, then, .

The symbol of the blockwise weighted-Jacobi operator, , is

By standard calculation, the eigenvalues of the error-propagation symbol, , are


where and

Noting that is very simple, we first consider a lower bound on the optimal LFA smoothing factor corresponding to .

Lemma 4.1.

and this value is achieved if and only if .


It is easy to check that for . The minimum of is with or and the maximum is with or . Thus, under the condition . ∎

Remark 4.1.

The optimal smoothing factor for damped Jacobi relaxation for the finite-element discretization of the Laplacian is with . Thus, this offers an intuitive lower bound on the possible performance of block relaxation schemes that include this as a piece of the overall relaxation.

From (13), we see that the only difference between the eigenvalues of DWJ relaxation for the Poisson-stabilized and projection stabilized methods is in the third eigenvalue, which depends on and, consequently, on the stabilization term.

4.1.1 Poisson-stabilized discretization with DWJ relaxation

For the Poisson-stabilized case, with and . By standard calculation, , with , and with or .

Theorem 4.1.

The optimal smoothing factor for the Poisson-stabilized discretization with DWJ relaxation is , that is,

and is achieved if and only if


with the condition that . Because , we need to require for all to achieve this factor. It follows that . ∎

4.1.2 Projection stabilized discretization with DWJ relaxation

For the projection stabilized discretization, depends on given in (9), and standard calculation gives with and with or .

Theorem 4.2.

The optimal smoothing factor for the projection stabilized discretization with DWJ relaxation is , that is,

and is achieved if and only if