# A closed-form multigrid smoothing factor for an additive Vanka-type smoother applied to the Poisson equation

We consider an additive Vanka-type smoother for the Poisson equation discretized by the standard finite difference centered scheme. Using local Fourier analysis, we derive analytical formulas for the optimal smoothing factors for two types of smoothers, called vertex-wise and element-wise Vanka smoothers, and present the corresponding stencils. Interestingly, in one dimension the element-wise Vanka smoother is equivalent to the scaled mass operator obtained from the linear finite element method, and in two dimensions the element-wise Vanka smoother is equivalent to the scaled mass operator discretized by bilinear finite element method plus a scaled identity operator. Based on these discoveries, the mass matrix obtained from finite element method can be used as an approximation to the inverse of the Laplacian, and the resulting mass-based relaxation scheme features small smoothing factors in one, two, and three dimensions. Advantages of the mass operator are that the operator is sparse and well conditioned, and the computational cost of the relaxation scheme is only one matrix-vector product; there is no need to compute the inverse of a matrix. These findings may help better understand the efficiency of additive Vanka smoothers and develop fast solvers for numerical solutions of partial differential equations.

## Authors

• 4 publications
• 16 publications
11/09/2021

### Novel mass-based multigrid relaxation schemes for the Stokes equations

In this work, we propose three novel block-structured multigrid relaxati...
11/04/2021

### A note on using the mass matrix as a preconditioner for the Poisson equation

We show that the mass matrix derived from finite elements can be effecti...
11/30/2021

### A novel multigrid method for elliptic distributed control problems

Large linear systems of saddle-point type have arisen in a wide variety ...
12/02/2020

### Explicit geometric construction of sparse inverse mass matrices for arbitrary tetrahedral grids

The geometric reinterpretation of the Finite Element Method (FEM) shows ...
06/03/2020

### Easy and Efficient preconditioning of the Isogeometric Mass Matrix

This paper deals with the fast solution of linear systems associated wit...
01/16/2017

### Poisson Vector Graphics (PVG) and Its Closed-Form Solver

This paper presents Poisson vector graphics, an extension of the popular...
07/14/2016

### Finite Element Integration with Quadrature on the GPU

We present a novel, quadrature-based finite element integration method f...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Consider the Poisson equation in one dimension (1D) or two dimensions (2D):

 −Δu=f, (1)

where and in 1D or and in 2D. The function is assumed to be sufficiently smooth so that finite-difference discretizations of provide effective approximations of the solution. The differential operator stands for the Laplacian: in 1D and in 2D. We assume a uniform mesh discretization with meshsize .

We apply the standard three- and five-point finite difference discretizations for the Laplacian; the corresponding stencils are given by

 Ah=1h2[−12−1] (2)

and

 Ah=1h2⎡⎢⎣−1−14−1−1⎤⎥⎦, (3)

respectively.

Let us denote the corresponding linear system by

 Ahuh=bh. (4)

The numerical solution of (4) is one of the most extensively explored topics in the numerical linear algebra literature. The discrete Laplacian

is a symmetric positive definite M-matrix, its eigenvalues are explicitly known, and it is used as a primary benchmark problem for the development of fast solvers. When the problem is large, iterative solvers that allow a high level of parallelism are often preferred.

One of the most efficient methods for solving (4) is multigrid [stuben1982multigrid, MR1156079]. The choice of multigrid components, such as a relaxation scheme as a smoother, plays an important role in the design of fast methods.

The choice of additive Vanka as relaxation scheme, which is suitable for parallel computing, has recently drawn considerable attention. Vanka-type smoothers have been applied to the Navier–Stokes equations [john2002higher, SManservisi_2006a, SPVanka_1986a], the Poisson equation using continuous and discontinuous finite elements methods [he2021local], the Stokes equations with finite element methods [farrell2021local, voronin2021low], poroelasticity equations [adler2021monolithic] in monolithic multigrid, and other problems. A restricted additive Vanka for Stokes using discretizations is presented in [saberi2021restricted], which shows its competitiveness with the multiplicative Vanka smoother. Nonoverlapping block smoothing using different patches has been discussed in [claus2021nonoverlapping] for the Stokes equations discretized by the marker-and-cell scheme. Muliplicative Vanka smoothers in combination with multigrid methods are discussed in [MR4024766, de2021two, franco2018multigrid]. Vanka-type relaxation has been used in many contexts, and for more details, we refer the reader to [JAdler_etal_2015b, john2002higher, VJohn_LTobiska_2000a].

In the literature, Vanka-type relaxation schemes demonstrate their high efficiency in a multigrid setting, but there seems no theoretical analysis for the convergence speed even for the simple Poisson equation. In this work we take steps towards closing this gap by considering the additive Vanka relaxation for the Poisson equation, and exploring stencils for the Vanka patches. We derive analytical optimal smoothing factors for two types of additive Vanka patches used in [he2021local] for hybridized and embedded discontinuous Galerkin methods. Moreover, we also find the corresponding stencils for the Vanka operator, and show that they are closely related to the scaled mass matrix obtained from the finite element method. Based on this discovery, we propose the mass-based relaxation scheme, which yields rapid convergence. This mass-based relaxation is very simple: the computation cost is only matrix-vector product and there is no need to solve the subproblems needed in an additive Vanka setting. Another advantage of the mass matrices obtained from (bi)linear elements is sparsity.

Solvers for the Poisson equation often form the first step for designing fast solvers for more complex problems, such as the Stokes equations, and Navier-Stokes equations. We therefore believe that the findings in this work can give some hints in future for designing fast numerical methods for these complex problems.

The remainder of this paper is organized as follows. In Section 2 we introduce the two types of additive Vanka smoothers for the Poisson equation. In Section 3, we present our theoretical analysis of optimal smoothing factors in one and two dimensions. Based on our analysis we also propose a mass-based smoother for the three-dimensional problem, where the mass matrix is obtained from the trilinear finite element method. In Section 4, we numerically validate our analytical observations and present an LFA two-grid convergence factor. Finally, in Section 5 we discuss our findings and draw some conclusions.

## 2 Vanka-type smoother

We are interested in exploring the structure of an additive Vanka-type smoother for solving the linear system (4) using multigrid. In general, this type can be thought of as related to the family of block Jacobi smoothers, which are suitable for parallel computation and are typically highly efficient within the context of multigrid smoothing.

Let the degrees of freedom (DoFs) of

be the set such that . is a restriction operator that maps the vector onto the vector in . Define

 Ai=ViAhVTi.

Then, we update current approximation by a single Vanka relaxation given by:
For ,

 Aiδui=Vi(bh−Ahuj)

and

 uj+1=uj+N∑i=1VTiWiδui.

A single iteration of Vanka can be represented as

 M=N∑i=1VTiWiA−1iVi, (5)

where the weighting matrix is given by the natural weights of the overlapping block decomposition. Each diagonal entry is equal to the reciprocal of the number of patches that the corresponding degree of freedom appears in. We refer to the Vanka operator.

For a single additive Vanka relaxation process, the relaxation error operator is given by

 S=I−ωMAh. (6)

A key factor is the choice of the patch, that is, the . Here, following [he2021local], we consider two patches, shown in Figure 1. We refer to the left patch in Figure 1 as an element-wise patch and the right one as a vertex-wise patch. We denote the corresponding relaxation error operators in (6) as and , respectively. The collection of circles indicate the number of DoFs in one patch . This means that the resulting subproblem is associated with a small matrix whose size is or . In the remainder of this work, for simplicity and clarity we use subscripts and to distinguish between the corresponding operators for element-wise Vanka and vertex-wise Vanka, respectively.

## 3 Local Fourier analysis

Local Fourier analysis (LFA) [MR1807961, wienands2004practical] is a useful tool for predicting and analyzing the convergence behavior of multigrid and other numerical algorithms. In this section, we employ LFA to examine the spectrum or spectral radius of the underlying operators to better understand the proposed Vanka-type smoother.

When applying LFA to multigrid methods, high and low frequencies [MR1807961] for standard coarsening () are introduced:

 θ∈Tlow=[−π2,π2)d,θ∈Thigh=[−π2,3π2)d\[−π2,π2)d,

where is the dimension of the underlying problem.

To quantitatively assess the performance of multigrid, LFA-predicted two-grid convergence and smoothing factors are used. In many cases the LFA smoothing factor is used due to its simplicity. It works under the assumption that the smoothing process reduces the high frequencies and leaves the low frequencies unchanged. The LFA smoothing factor typically offers a rather sharp bound on the actual two-grid performance.

###### Definition 3.1.

Let be the relaxation error operator. Then, the corresponding LFA smoothing factor for is given by

 μ(ω)=maxθ∈Thigh{ρ(˜S(θ,ω))}, (7)

where is the symbol of , is algorithmic parameter, and denotes the spectral radius of the matrix .

Note that the LFA smoothing factor is a function of . Often, one can minimize (7) with respect to to obtain fast convergence speed. We define the optimal smoothing factor as

 μopt=minωμ. (8)

In this work, the symbol for the Laplacian considered is a scalar, so the spectral radius is reduced to the maximum of a scalar function. In the following, we use LFA to identify the optimal smoothing factor for the additive Vanka-type relaxation schemes and explore the structure of the Vanka operator defined in (5). Before providing our detailed analysis of the smoothing factor for different relaxation schemes, we summarize our results in Table 1. The table provides a review of quantitative results of additive Vanka smoothers, and for comparison we include results for the standard point-wise damped Jacobi smoother.

Note that a general form of the symbol of additive Vanka operator for the Stokes equations has been discussed in [farrell2021local], which gives , where is called the relative Fourier matrix. Here, we can directly apply the formula of to our additive Vanka operator; see [farrell2021local].

In the analysis that follows, we will use to denote the imaginary scalar satisfying

 ι2=−1.

### 3.1 Symbols of Vanka patches in 1D

In this subsection, we first consider the analytical symbol of the element-wise patch, then the vertex-wise patch for the Laplacian in 1D. We discuss the optimal smoothing factor for each case and derive the corresponding stencil for the Vanka operator.

#### 3.1.1 Element-wise Vanka patch in 1D

It can easily be shown that the symbol of , see (2), is given by

 ˜Ah=1h22(1−cosθ). (9)

Moreover, for the element-wise patch the subproblem matrix is

 Ai=1h2(2−1−12).

Following [farrell2021local], the relative Fourier matrix is

 Φ=(100eιθ).

Then, the symbol of is given by , where

 ˜V=(11),˜W=12(1001),A−1i=h23(2112).

Based on the above formulas, we obtain

 ˜Me=h26(4+eιθ+e−ιθ). (10)

Formula (10) indicates that the element-wise Vanka patch corresponds to the stencil

 Me=h26[141]. (11)

Recall that the mass stencil in 1D using linear finite elements is given by

 Mfe=h6[141]. (12)

This means that the element-wise Vanka operator is equivalent to a scaled mass matrix obtained from the linear finite element method, .

Next, we give the optimal smoothing factor for the element-wise Vanka relaxation scheme.

###### Theorem 3.1.

The optimal smoothing factor of for the vertex-wise Vanka in 1D is

 μe,opt=minωmaxθ∈Thigh|1−ω˜Me˜Ah|=117≈0.059, (13)

where the minimum is uniquely achieved at .

###### Proof.

When ,

 ˜Me˜Ah=23(2−cosθ−cos2θ)∈[43,32].

Thus,

 μ(ω)=max|1−ω˜Me˜Ah|=max{∣∣∣1−ω43∣∣∣,∣∣∣1−ω32∣∣∣}.

To minimize , we require

 ∣∣∣1−ω43∣∣∣=∣∣∣1−ω32∣∣∣,

which gives . Then, . ∎

It is well known that the optimal smoothing factor for damped Jacobi relaxation for the Laplacian in 1D (with ) is . This suggests that using the additive Vanka smoother for multigrid achieves much faster convergence.

#### 3.1.2 Vertex-wise Vanka patches in 1D

We now consider the vertex-wise patch. The subproblem matrix is

 Ai=1h2⎛⎜⎝2−10−12−10−12⎞⎟⎠.

Again, following [farrell2021local], the relative Fourier matrix is

 Φ=⎛⎜⎝e−ιθ0001000eιθ⎞⎟⎠.

Then, the symbol of is given by , where

 ˜V=⎛⎜⎝111⎞⎟⎠,˜W=13⎛⎜⎝100010001⎞⎟⎠,A−1i=h24⎛⎜⎝321242123⎞⎟⎠.

Using the above formulas, we have

 ˜M=˜VT˜WΦTA−1iΦ˜V=h212(10+4eιθ+4e−ιθ+eι2θ+4e−ι2θ). (14)

Based on (14), the stencil of is

 Mv=h212[141041]. (15)

Compared with (11), the vertex-wise Vanka uses a wider stencil.

###### Theorem 3.2.

The optimal smoothing factor of for the vertex-wise Vanka in 1D is

 μv,opt=minωmaxθ∈Thigh|1−ω˜Mv˜Ah|=126≈0.039, (16)

where the minimum is uniquely achieved at .

###### Proof.

From (9) and (14), we have

 ˜Mv˜Ah =13(1−cosθ)(5+4cosθ+cos(2θ)), =23(1−cosθ)(2+2cosθ+cos2θ), =23(1−cosθ)((cosθ+1)2+1).

Note that when . Let , where . To identify the range of , we first compute its derivative:

 g′(x)=−23x(3x+2)<0,forx∈[−1,0].

It follows that

 (109)2=g(−2/3)≤g(x)≤g(0)=g(1)=43.

That is, for . To minimize , we require . Then, it follows that

Again, the optimal smoothing factor for vertex-wise Vanka is significantly smaller (and hence better) than that of the damped Jacobi relaxation scheme.

### 3.2 Symbols of Vanka patches in 2D

Similarly to previous subsection, we first consider the analytical symbol for the element-wise patch, then for the vertex-wise patch for the Laplacian in 2D.

#### 3.2.1 Element-wise Vanka patch in 2D

The symbol of Laplace operator discretized by five-point stencil, see (3), is

 ˜Ah=1h2(4−eιθ1−eιθ2−e−ιθ1−e−ιθ2)=1h2(4−2cosθ1−2cosθ2). (17)

For the vertex-wise patch, following [farrell2021local], the relative Fourier matrix is

 Φ=⎛⎜ ⎜ ⎜ ⎜⎝10000eιθ10000eιθ20000eι(θ1+θ2)⎞⎟ ⎟ ⎟ ⎟⎠.

Then, the symbol of is given by , where

 ˜W=14⎛⎜ ⎜ ⎜⎝1000010000100001⎞⎟ ⎟ ⎟⎠,˜V=⎛⎜ ⎜ ⎜⎝1111⎞⎟ ⎟ ⎟⎠,Ai=1h2⎛⎜ ⎜ ⎜⎝4−1−10−140−1−104−10−1−14⎞⎟ ⎟ ⎟⎠.

We have

 A−1i=h2⎛⎜ ⎜ ⎜ ⎜⎝7/241/121/121/241/127/241/241/121/121/247/241/121/241/121/127/24⎞⎟ ⎟ ⎟ ⎟⎠,

and it follows that

 ˜Me =˜VT˜WΦTA−1iΦ˜V, =h296(28+4(eιθ1+e−ιθ1+eιθ2+e−ιθ2)+eι(θ1+θ2)+e−ι(θ1+θ2)+(eιθ1e−ιθ2+e−ιθ1eιθ2)), =h224(7+2(cosθ1+cosθ2)+cosθ1cosθ2). (18)

Next, we give the optimal smoothing factor for the element-wise Vanka relaxation in 2D.

###### Theorem 3.3.

The optimal smoothing factor of for the element-wise Vanka relaxation in 2D is given by

 μe,opt=minωmaxθ∈Thigh|1−ω˜Me˜Ah|=725≈0.280, (19)

where the minimum is uniquely achieved at .

###### Proof.

From (17) and (18), we have

 ˜Me˜Ah =124(7+2(cosθ1+cosθ2)+cosθ1cosθ2)(4−2cosθ1−2cosθ2), =112(7+2(cosθ1+cosθ2)+cosθ1cosθ2)(2−(cosθ1+cosθ2)),

where . Let

 g(η1,η2)=112(7+2(η1+η2)+η1η2)(2−η1−η2).

We have and is continuous in . By the Extreme Value Theorem, achieves its extremal values at the boundary of or its derivatives are zeros. We first consider the derivatives,

 gη1 =112(−3−4η1−2η1η2−η22−2η2), gη2 =112(−3−4η2−2η1η2−η22−2η1).

Solving gives . Thus, is a possible global extreme value.

Next, we compute the extremal values of at the boundary . Due to the symmetric of , we only need to consider the following two cases.

• Case 1: and . We have

 g(−1,η2)=112(5+η2)(3−η2).

Thus,

 1=g(−1,1)≤g(−1,η2)≤g(−1,−1)=43.
• Case 2: and . We have

 g(1,η2)=14(3+η2)(1−η2).

Thus,

 0=g(1,1)≤g(1,η2)≤g(1,−1)=1.

If we restrict , we find that the maximum and minimum of are given by

 g(−1,−1)=43,g(1,0)=34, (20)

respectively. It follows that and . ∎

It is well known that the optimal smoothing factor for damped Jacobi relaxation for the Laplacian in 2D is with [MR1807961]. This suggests that using the additive Vanka smoother for multigrid method, convergence is faster compared to the damped Jacobi relaxation scheme.

Based on the symbol of , we obtain the stencil of ,

 Me=h296⎡⎢⎣1414284141⎤⎥⎦. (21)

Recall that the mass matrix stencil using bilinear finite elements is

 Mfe=h236⎡⎢⎣1414164141⎤⎥⎦. (22)

Now, we can make a connection between and :

 Me=38Mfe+h28I, (23)

where

 I=⎡⎢⎣000010000⎤⎥⎦. (24)

The relationship (23) is interesting, and suggests that the mass matrix obtained from bilinear elements might be a good approximation to the inverse of . Let us, then, move to consider the mass matrix (22) as an approximation to the inverse of .

###### Theorem 3.4.

Given the mass stencil in (22) and the relaxation scheme , the corresponding optimal smoothing factor is

 μfe,opt=minωmaxθ∈Thigh|1−ω˜Mfe˜Ah|=13≈0.333, (25)

where the minimum is uniquely achieved at .

###### Proof.

It can easily be shown that for . Thus, the optimal is . Then, . ∎

From Theorem 3.4, we see that the optimal smoothing factor of 0.333 for the mass-based relaxation is close to the optimal smoothing factor 0.280 for the element-wise Vanka patch, and it is better than 0.391 obtained from the vertex-wise Vanka, see (30), discussed in the next subsection. Thus, mass matrix could be used as a good approximation to the inverse of the Laplacian considered here.

#### 3.2.2 Vertex-wise Vanka patches in 2D

Now, we analyse the smoothing factor for the vertex-wise patch. Following [farrell2021local], the relative Fourier matrix is

 Φ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝e−ιθ200000e−ιθ100000100000eιθ100000eιθ2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (26)

The symbol of is given by with

 ˜W=15I,˜V=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝11111⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,Ai=1h2⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝40−10004−100−1−14−1−100−14000−104⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (27)

where

stands for identity matrix with size

.

It can be shown that

 A−1i=h2⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝13/481/481/121/481/481/4813/481/121/481/481/121/121/31/121/121/481/481/1213/481/481/481/481/121/4813/48⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (28)

From (26), (27) and (28), we have

 ˜Mv =˜VT˜WΦTA−1iΦ˜V, =h2240(8(e−ιθ2+eιθ2+e−ιθ1+eιθ1)+68+(e−ιθ2+e−ιθ1)2+(eιθ2+eιθ1)2+2(e−ιθ2eιθ1+eιθ2e−ιθ1)), =h2240(16(cosθ1+cosθ2)+68+2((cosθ1+cosθ2)2−(sinθ1+sinθ2)2)+4cos(θ1−θ2)) =h2120((cosθ1+cosθ2+4)2+(cosθ1+cosθ2)2+16).

Based on the symbol of , we can write down the corresponding stencil of as follows:

 Mv=h2240⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣00100028201868810282000100⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (29)

Now, we are able to give the optimal smoothing factor for the vertex-wise Vanka relaxation scheme.

###### Theorem 3.5.

The optimal smoothing factor of for the vertex-wise Vanka relaxation in 2D is

 μopt=minωmax|1−ω˜Mv˜Ah|=923≈0.391, (30)

where the minimum is uniquely achieved at .

###### Proof.

We first compute

 ˜Mv˜Ah =1120((cosθ1+cosθ2+4)2+(cosθ1+cosθ2)2+16)(4−2cosθ1−2cosθ2) =160((4+η)2+η2+16)(2−η),

where with .

Let . We find that

This means that is a decreasing function. Thus, for , we have

 0=g(η=2)≤g(η)≤g(η=−2)=85.

This means that .

If we restrict , then . In this situation, we have

 710=g(η=1)≤g(η)≤g(η=−2)=85.

Since for , we have and . ∎

###### Remark 3.1.

Note that in 2D the element-wise Vanka, see (21), uses fewer points than that of vertex-wise Vanka, see (29). However, the corresponding optimal smoothing factor of element-wise Vanka, see (19), is smaller than that of vertex-wise Vanka, see (30), which is different than the case in 1D.

### 3.3 Extension to the 3D case

While we do not include a smoothing analysis of Vanka-type solvers for the 3D case, we can still make a few interesting observations. In particular, motivated by our findings on the potential role of the mass matrix for relaxation, we further explore the scaled mass matrix in 3D as an approximation to the inverse of the Laplacian. Let

 M3=h−4Me⊗Me⊗Me,

where is defined in (11). The symbol of

can be obtained by tensor product given by

 ˜M3=h227(2+cosθ1)(2+cosθ2)(2+cosθ3).

The symbol of in 3D is

 ˜Ah=1h22(3−cosθ1−cosθ2−cosθ3).

Let , which is the Jacobi matrix.

###### Theorem 3.6.

If we consider the point-wise damped Jacobi as a preconditioner for the Laplacian, then the corresponding optimal smoothing factor of is

 μJ,opt=minωmaxθ∈Thigh|1−ω˜MJ˜Ah|=57≈0.714, (31)

where the minimum is uniquely achieved at .

###### Proof.

Note that . For , . Thus, to minimize , we require . It follows . ∎

Next, we consider mass-based relaxation scheme for the Laplacian in 3D.

###### Theorem 3.7.

Given the scaled mass stencil in (3.3) and mass-based relaxation scheme in 3D, the corresponding optimal smoothing factor is

 μm,opt=minωmaxθ∈Thigh|1−ω˜M3˜Ah|=131212≈0.618, (32)

where the minimum is uniquely achieved at .

###### Proof.

Let

 T=˜M3˜Ah=227(2+cosθ1)(2+cosθ2)(2+cosθ2)(3−cosθ1−cosθ2−cosθ3).

Define with . To find the extremal values of , we will consider its derivatives and the function values at the boundary of underlying domain. We compute the derivatives of with respect to and , given by

 gx =227(2+y)(2+z)(1−2x−y−z), gy =227(2+x)(2+z)(1−2y−x−z), gz =227(2+x)(2+y)(1−2z−x−y).

Solving with gives . However, such that does not belong to .

Let us define , , and . Note that corresponds to