An AMG saddle point preconditioner with application to mixed Poisson problems on adaptive quad/cube meshes

TL;DR

This paper introduces a specialized algebraic multigrid (AMG) preconditioner for saddle point problems in mixed Poisson discretizations, demonstrating near mesh-independent convergence on adaptive meshes with heterogeneous coefficients.

Contribution

The paper develops a novel AMG saddle point preconditioner with a stabilized prolongation operator tailored for mixed Poisson problems on adaptive meshes.

Findings

SPAMG preconditioner achieves mesh-independent iteration counts.

Effective on both 2D and 3D adaptive meshes.

Handles heterogeneous coefficients efficiently.

Abstract

We investigate various block preconditioners for a low-order Raviart-Thomas discretization of the mixed Poisson problem on adaptive quadrilateral meshes. In addition to standard diagonal and Schur complement preconditioners, we present a dedicated AMG solver for saddle point problems (SPAMG). A key element is a stabilized prolongation operator that couples the flux and scalar components. Our numerical experiments in 2D and 3D show that the SPAMG preconditioner displays nearly mesh-independent iteration counts for adaptive meshes and heterogeneous coefficients.

Tables

Table 1. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	459	299	22	9	8	8
$5$	$>$1000	773	22	10	8	8
$6$	-	$>$1000	24	12	9	10
$7$	-	-	24	14	10	10
$8$	-	-	24	14	10	11
$9$	-	-	24	14	10	11

Table 2. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	459	299	22	9	8	8
$5$	$>$1000	773	22	10	8	8
$6$	-	$>$1000	24	12	9	10
$7$	-	-	24	14	10	10
$8$	-	-	24	14	10	11
$9$	-	-	24	14	10	11

Table 3. (b)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4-7$	$>$1000	$>$1000	87	10	7	7
$5-8$	-	-	87	10	7	8
$6-9$	-	-	99	13	10	10
$7-10$	-	-	98	13	10	10
$8-11$	-	-	112	15	11	11
$9-12$	-	-	149	15	11	11

Table 4. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	20	9	7	7
$4$	22	11	8	8
$5$	22	13	9	9
$6$	24	14	9	10
$7$	24	15	11	12

Table 5. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	20	9	7	7
$4$	22	11	8	8
$5$	22	13	9	9
$6$	24	14	9	10
$7$	24	15	11	12

Table 6. (b)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3-6$	238	11	8	8
$4-7$	322	11	8	8
$5-8$	338	13	9	9
$6-9$	383	14	10	10
$7-10$	390	17	12	12

Table 7. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	$>$1000	$>$1000	85	18	12	11
$5$	-	-	102	20	13	13
$6$	-	-	113	23	14	14
$7$	-	-	130	24	16	15
$8$	-	-	142	25	16	17
$9$	-	-	169	25	17	18

Table 8. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	$>$1000	$>$1000	85	18	12	11
$5$	-	-	102	20	13	13
$6$	-	-	113	23	14	14
$7$	-	-	130	24	16	15
$8$	-	-	142	25	16	17
$9$	-	-	169	25	17	18

Table 9. (b)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4-7$	$>$1000	$>$1000	250	18	11	10
$5-8$	-	-	255	19	13	12
$6-9$	-	-	415	22	14	13
$7-10$	-	-	338	24	15	15
$8-11$	-	-	394	24	16	16
$9-12$	-	-	539	24	16	17

Table 10. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	65	19	13	13
$4$	84	21	14	14
$5$	103	23	15	16
$6$	113	24	17	17
$7$	128	26	18	19

Table 11. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	65	19	13	13
$4$	84	21	14	14
$5$	103	23	15	16
$6$	113	24	17	17
$7$	128	26	18	19

Table 12. (b)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3-6$	330	22	14	14
$4-7$	456	23	15	14
$5-8$	486	24	16	15
$6-9$	471	25	17	18
$7-10$	493	27	18	20

Table 13. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	$>$1000	$>$1000	36	12	8	9
$5$	-	-	56	14	9	8
$6$	-	-	55	14	9	10
$7$	-	-	47	15	10	11
$8$	-	-	43	15	10	11
$9$	-	-	41	15	11	11

Table 14. (a)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4$	$>$1000	$>$1000	36	12	8	9
$5$	-	-	56	14	9	8
$6$	-	-	55	14	9	10
$7$	-	-	47	15	10	11
$8$	-	-	43	15	10	11
$9$	-	-	41	15	11	11

Table 15. (b)

Level	# Iterations
	NoPC	Diag	Schur	SPAMG
				Uzawa	Vanka one	Vanka scale
$4-7$	$>$1000	$>$1000	379	12	8	8
$5-8$	-	-	547	13	8	8
$6-9$	-	-	930	13	9	9
$7-10$	-	-	$>$1000	15	10	10
$8-11$	-	-	-	15	10	10
$9-12$	-	-	-	16	10	11

Table 16. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	40	13	9	9
$4$	45	13	8	9
$5$	56	15	9	9
$6$	53	15	10	10
$7$	47	15	11	11

Table 17. (a)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3$	40	13	9	9
$4$	45	13	8	9
$5$	56	15	9	9
$6$	53	15	10	10
$7$	47	15	11	11

Table 18. (b)

Level	# Iterations
	Schur	SPAMG
		Uzawa	Vanka one	Vanka scale
$3-6$	637	15	11	11
$4-7$	855	14	10	11
$5-8$	$>$1000	15	11	11
$6-9$	-	16	12	12
$7-10$	-	20	14	15

Equations

$$\Vert v{\Vert }_m^2:=\sum _{|\alpha |\le m}{\int }_{\Omega }|D^{\alpha }v{|}^2\mathrm{d}x,|v{|}_m^2:=\sum _{|\alpha |=m}{\int }_{\Omega }|D^{\alpha }v{|}^2\mathrm{d}x.$$

$$\mathbf{H}(\mathrm{div},\Omega ):=\{\mathbf{v}\in (L^2(\Omega ){)}^d:\mathrm{div}\mathbf{v}\in L^2(\Omega )\},$$

$$\Vert \mathbf{v}{\Vert }_{\mathrm{div}}^2=\Vert \mathbf{v}{\Vert }_0^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_0^2.$$

$$H_{0,D}^1(\Omega ):=\{\varphi \in H^1(\Omega ):\varphi {|}_{{\Gamma }_D}=0\}.$$

$$-\mathrm{div}\mathcal{K}(\mathbf{x})\nabla p=f$$

$$p=p_0$$

$$\mathcal{K}(\mathbf{x})\nabla p\cdot \mathbf{n}=g$$

$$f\in L^2(\Omega ),p_0\in H^{1/2}({\Gamma }_D)\text{and}g\in L^2({\Gamma }_N).$$

$$\mathbf{u}=\mathcal{K}\nabla p$$

$$-\mathrm{div}\mathbf{u}=f$$

$$p=p_0$$

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

$${\mathbf{H}}_{0,N}(\Omega ):=\{\mathbf{\tau }\in \mathbf{H}(\mathrm{div},\Omega ):⟨\mathbf{\tau }\cdot \mathbf{n},\varphi ⟩=0\text{ for all }\varphi \in H_{0,D}^1(\Omega )\}.$$

$${\int }_{\Omega }{\mathcal{K}}^{-1}\mathbf{u}\cdot \mathbf{v}\mathrm{d}x+{\int }_{\Omega }p\mathrm{div}\mathbf{v}\mathrm{d}x$$

$${\int }_{\Omega }q\mathrm{div}\mathbf{u}\mathrm{d}x$$

$$a(\mathbf{u},\mathbf{v}):={\int }_{\Omega }{\mathcal{K}}^{-1}\mathbf{u}\cdot \mathbf{v}\mathrm{d}x,b(\mathbf{v},p):={\int }_{\Omega }p\mathrm{div}\mathbf{v},$$

$$\mathbb{A}(\mathbf{v}):={\int }_{{\Gamma }_D}{p}_0(\mathbf{v}\cdot \mathbf{n})\mathrm{d}s,\mathbb{B}(q):=-{\int }_{\Omega }fq\mathrm{d}x.$$

$$a({\mathbf{u}}_0,\mathbf{v})+b(\mathbf{v},p)$$

$$b(\mathbf{u},q)$$

$$a({\mathbf{u}}_h,{\mathbf{v}}_h)+b({\mathbf{v}}_h,p_h)$$

$$b({\mathbf{u}}_u,q_h)$$

$$\mathbf{v}{|}_{{\Omega }_e}=\left\{\begin{array}{cc}(a_0+b_{0}x,a_1+b_{1}y{)}^T& \text{if }d=2\\ (a_0+b_{0}x,a_1+b_{1}y,a_2+b_{2}z{)}^T& \text{if }d=3\end{array}\right\},p{|}_{{\Omega }_e}=c_0,$$

$$\Vert \mathbf{u}-{\mathbf{u}}_h{\Vert }_0$$

$$\Vert p-p_h{\Vert }_0$$

$${\mathcal{A}}_h:=\left[\begin{array}{cc}\mathbf{A}& {\mathbf{B}}^{\top }\\ \mathbf{B}& 0\end{array}\right]\sim \left[\begin{array}{cc}{\mathrm{Id}}_h& -{\mathrm{grad}}_h\\ {\mathrm{div}}_h& 0\end{array}\right],$$

$${\mathcal{A}}_h=\left[\begin{array}{cc}\mathbf{I}& 0\\ \mathbf{B}{\mathbf{A}}^{-1}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& -\mathbf{S}\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& {\mathbf{A}}^{-1}{\mathbf{B}}^{\top }\\ 0& \mathbf{I}\end{array}\right],\mathbf{S}=\mathbf{B}{\mathbf{A}}^{-1}{\mathbf{B}}^{\top },$$

$$\mathcal{B}=\left[\begin{array}{cc}\mathbf{M}& 0\\ 0& -\mathbf{N}\end{array}\right],$$

$$\mathbf{M}\mathbf{A}\approx \mathbf{I},$$

$$\mathbf{N}\mathbf{S}\approx \mathbf{I}.$$

$${\mathbf{S}}_{BFBt}^{-1}:=({\mathbf{BC}}^{-1}{\mathbf{B}}^{\top }{)}^{-1}({\mathbf{BC}}^{-1}{\mathbf{AD}}^{-1}{\mathbf{B}}^{\top })({\mathbf{BD}}^{-1}{\mathbf{B}}^{\top }{)}^{-1},$$

$${\mathbf{P}}_l\in {\mathbb{R}}^{n_l\times n_{l+1}},{\mathbf{P}}_l:=\left[\begin{array}{c}{\mathbf{P}}_{l,F}\\ {\mathbf{I}}_{l,C}\end{array}\right].$$

$${\mathbf{\Lambda }}_{l+1}:={\mathbf{P}}_l^{\top }{\mathbf{\Lambda }}_l{\mathbf{P}}_l.$$

$$e_i\approx -\frac{1}{{\lambda }_{ii}}\sum _{j\ne i}{\lambda }_{ij}{e}_j.$$

$$\left[\begin{array}{cc}\mathbf{A}& {\mathbf{B}}^{\top }\\ \mathbf{B}& -\mathbf{C}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ p\end{array}\right]=\left[\begin{array}{c}\mathbf{v}\\ q\end{array}\right].$$

$${\mathcal{A}}_l=\left[\begin{array}{cc}{\mathbf{A}}_l& {\mathbf{B}}_l^{\top }\\ {\mathbf{B}}_l& -{\mathbf{C}}_l\end{array}\right],$$

Full text

1 Introduction

In many applications one is more interested in the gradient of the solution of the Poisson equation than in the solution itself. One practical example is solving the Richards equation for subsurface flow, where the scalar variable represents the water pressure [34, 6]. One approach is to use standard methods such as finite differences (FD), finite elements (FE) or finite volumes (FV), to obtain an approximate solution and then to compute itself gradient via numerical differentiation, which may lead to a loss in accuracy [3]. A mixed finite element (MFE) discretization then appears to be a more natural choice since the gradient of the solution is part of the unknown variables of the formulation. Generally, in MFE methods the vector valued quantity is approximated at least with the order of accuracy of the scalar unknown [13]. In addition, the computed solution satisfies mass conservation at the element level, a property highly relevant for simulating fluid flow problems where the governing partial differential equations (PDE) are derived from mass balance laws [9].

MFE discretizations lead to symmetric indefinite systems of algebraic equations that may be solved by iterative methods. In order to obtain a solution with a reasonable investment of computational resources, the use of optimal preconditioners becomes mandatory. Although the existence of an optimal preconditioner for general (anisotropic) coefficients remains an open question, there has been a significant amount of work in this direction. [2] proposes two preconditioners for the operator $(I-\operatorname{\mathrm{grad}}\operatorname{\mathrm{div}})$ in two dimensions, one based on domain decomposition and another on multigrid. The latter has been generalized to operators of the form $(\rho I-\mu\operatorname{\mathrm{grad}}\operatorname{\mathrm{div}})$ and three-dimensional problems [4]. From the observation that a MFE discretization of a generalized diffusion equation is well posed in a product of two different discrete function spaces, [32] and [31] propose two block-diagonal preconditioners. One of them is what we will refer to as Schur complement preconditioner and consists of a lumped diagonal approximation of the $(1,1)$ block and an algebraic multigrid V-cycle to approximate a Schur complement on the $(2,2)$ block. The preconditioners introduced in these papers are shown to be optimal with respect to the mesh size and various classes of coefficients (such as a conductivity tensor). Recent work [26] introduces a new preconditioner whose key ingredient is the approximation of the $(1,1)$ block with an auxiliary space multigrid method [27] that offers optimality with respect to a larger class of coefficients. The authors include numerical evidence for two dimensional problems posed on uniform rectangular meshes.

For some problems we may require a very fine mesh in order to correctly resolve the phenomena we are trying to model. A uniform mesh might be undesirable or even impractical given the computational resources it requires. One solution to this problem is to use locally refined meshes, which use the correct resolution only in the portion of the domain that really requires it [5, 8, 15]. There are essentially two possibilities: either we allow the so called hanging nodes or not. In the context of MFE, implementations making use of adaptively refined meshes with hanging nodes are not the standard, nevertheless the theoretical framework has been established in the early 90s [19, 20]. Unfortunately, none of the above mentioned methods for preconditioning address the case of adaptively refined meshes. Hence, the goal of this paper is to present a multigrid preconditioner that retains its robustness in such a case for both two and three dimensional problems. In addition, we aim to allow for variable matrix-valued coefficients. In our approach, we acknowledge the different scales in the saddle point structure and introduce a prolongation operator that suitably couples the vector and scalar unknowns. We demonstrate superiority of this approach to the Schur complement method using a variety of numerical examples.

2 Problem formulation

In this section we briefly review the mixed formulation of a Poisson-type partial differential equation. The material covered here is standard and found in many textbooks; see e.g. [9, 11], augmented with a few recent results.

Let $\Omega\in\mathbb{R}^{d}$ for $d\in\{2,3\}$ be a bounded Lipschitz domain with boundary $\Gamma=\partial\Omega$. $H^{m}(\Omega)$ will denote the standard Hilbert space of functions in $L^{2}(\Omega)$ whose weak derivatives up to order $m\geq 0$ are also square integrable. $H^{m}(\Omega)$ is endowed with the usual norm and seminorm,

[TABLE]

We define the velocity space

[TABLE]

which is a Hilbert space with the norm

[TABLE]

We assume that $\Gamma=\Gamma_{D}\cup\Gamma_{N}$, with $\Gamma_{D}\cap\Gamma_{N}=\emptyset$, and that $\Gamma_{D}$ has nonzero $(d-1)$-dimensional Lebesgue measure. Lastly, define (in the weak sense)

[TABLE]

In the following, we use bold mathematical symbols to denote vectors and matrices over $\mathbb{R}$ and the normal font for scalar quantities.

We consider the equation

[TABLE]

where $\boldsymbol{n}$ is the outer normal vector to the boundary $\Gamma$. The conductivity tensor $\boldsymbol{\mathcal{K}}(\boldsymbol{x})$ is a $d\times d$ symmetric positive definite matrix whose smallest eigenvalue is bounded uniformly away from zero. Furthermore, the data is required to satisfy

[TABLE]

Introducing the variable $\boldsymbol{u}=\boldsymbol{\mathcal{K}}\nabla p$ leads to the mixed first-order system

[TABLE]

2.1 Weak formulation

To derive the mixed weak formulation of (7) we introduce the following space,

[TABLE]

The dual paring $\langle\,\cdot\,,\,\cdot\,\rangle$ is defined via Green’s formula; see [9, pp. 50]. Multiplying the first equation in (7) by $\boldsymbol{\mathcal{K}^{-1}}$ and then by a test function $\boldsymbol{v}\boldsymbol{\in}\boldsymbol{H}_{0,N}(\Omega)$, the second by some $q\in L^{2}(\Omega)$, integrating over $\Omega$ and using Green’s formula on the gradient term yields the following weak formulation: Find $(\boldsymbol{u},p)\in\boldsymbol{H}_{0,N}(\Omega)\times L^{2}(\Omega)$ such that

[TABLE]

Let us define the bilinear forms $a:\boldsymbol{H}(\operatorname{\mathrm{div}},\Omega)\times\boldsymbol{H}(\operatorname{\mathrm{div}},\Omega)\to\mathbb{R}$ and $b:\boldsymbol{H}(\operatorname{\mathrm{div}},\Omega)\times L^{2}(\Omega)\to\mathbb{R}$,

[TABLE]

and linear functionals $\mathbb{A}:\boldsymbol{H}(\operatorname{\mathrm{div}},\Omega)\to\mathbb{R}$, $\mathbb{B}:L^{2}(\Omega)\to\mathbb{R}$,

[TABLE]

To enforce the Neumann boundary condition (5c), (7d), let $\tilde{p}$ a classical solution of (5) with $f=0$ and $p_{0}=0$. Then, taking $\boldsymbol{\tilde{u}}=\boldsymbol{\mathcal{K}}\nabla\tilde{p}$ leads to the following mixed weak formulation of (7): Find $\boldsymbol{u}=\boldsymbol{\tilde{u}}+\boldsymbol{u}_{0}$ with $\boldsymbol{u}_{0}\in\boldsymbol{H}_{0,N}(\Omega)$ and $p\in L^{2}(\Omega)$ such that

[TABLE]

Existence and uniqueness of a solution for the problem (12) follows from standard arguments, namely establishing coercivity of the bilinear form $a(\,\cdot\,,\,\cdot\,)$, and the Ladyženskaja-Babuška-Brezzi (LBB) condition; see e.g. [13].

2.2 Discretization

We pick finite dimensional subspaces $\boldsymbol{X}_{0}^{h}\subset\boldsymbol{H}_{0,N}(\Omega)$ and $V^{h}\subset L^{2}(\Omega)$ and define the following problem: Find $(\boldsymbol{u}_{h},p_{h})\in\boldsymbol{X}_{0}^{h}\times V^{h}$ such that

[TABLE]

To ensure that (13) is well posed, the pair of spaces $(\boldsymbol{X}_{0}^{h},V^{h})$ must be chosen such that the LBB condition is fulfilled for the discrete problem. Many spaces with this property have been developed since the early seventies, such as the Raviart-Thomas [33] and Brezzi-Douglas-Marini [12] spaces. In this paper we choose the lowest order Raviart-Thomas discretization $\mathcal{RT}_{0}$ defined on rectangles/hexahedra. Hence, in an element $\Omega_{e}$ the velocity and pressure test functions take the form

[TABLE]

respectively, where $a_{i},b_{i}$ for $i=0,\ldots,d-1$ and $c_{0}$ are constants. The degrees of freedom are shown in Figure 1.

The following approximation properties are well known for $\mathcal{RT}_{0}$ in the context of affine elements defined on uniform meshes [9],

[TABLE]

where we assume that the pair $(\boldsymbol{u},p)$ fulfils the regularity requirements required by the right hand side of (15).

2.3 Adaptivity

For well behaved (smooth) problems posed on convex domains, the use of a uniform mesh usually offers satisfactory results when computing a numerical solution, that is, there is an optimal trade-off between numerical effort (computational resources) invested and effective reduction of the error. Nevertheless, there are situations in which the mesh resolution required to accurately reproduce the physical behavior of the underlying PDE becomes computationally impractical if imposed on the whole domain.

Adaptive mesh refinement (AMR) provides a valuable tool in order to reduce the computational complexity in such situations by increasing the mesh resolution only locally where is required. As stated in section Section 2.2, we will work with meshes composed of rectangles/hexahedra. The refinement schemes used in this paper include the case of a mesh with hanging nodes, that is, we allow (a nonempty) intersection of two elements to be a complete side of a neighboring element. Additionally, we choose to use 2:1 balanced meshes: The length ratio between a coarse and a fine element is at most of factor two; see Figure 2. Non-balanced meshes would be possible as well but require more technical work regarding the definition of MFE spaces and parallelization.

Continuity of fluxes across interfaces for this kind of meshes can be enforced in several ways. One is to eliminate te degrees of freedom at hanging nodes [19]. Another approach is to use Mortar finite elements [1]. We will follow the former. Due to our assumption of 2:1 balance, for a $\mathcal{RT}_{0}$ discretization given a hanging node with flux value $\tilde{u}$, there is only one non-hanging node $u$ lying in the same edge/face; see Figure 2.

Estimates for locally refined meshes using hanging nodes have been studied in [19, 20]. Essentially, it is shown that the $\mathcal{RT}_{0}$ spaces still respect the LBB condition after introducing locally refined grids, and hence the estimates (15) remain valid.

2.4 Preconditioning

The discretization of (7) via stable MFE leads to a saddle point problem defined by the following block matrix,

[TABLE]

where $\boldsymbol{A}$ is the vector mass matrix and $\boldsymbol{B}$ is the discrete divergence operator. This is a well studied problem and there are several solution methods available; see e.g. [7] for a comprehensive review. They range from Uzawa algorithms and its variants [10], projection methods [23], to block factorization methods [16, 29]. In contrast to most methods that treat flux and pressure individually, we introduce a monolithic multigrid method originally developed for the Stokes equations [30] that is to our knowledge yet unpublished.

To prepare the following exposition, let us briefly discuss notation and some background. We consider the case that the matrix (16) is symmetric and indefinite. The factorization

[TABLE]

implies that $\mathcal{A}_{h}$ is congruent to a block diagonal matrix [21]. This fact motivates the use of a preconditioner of the form

[TABLE]

where $\boldsymbol{M}$ and $\boldsymbol{N}$ satisfy

[TABLE]

A simple choice is to take $\boldsymbol{M}$ as the inverse of the lumped mass matrix $\boldsymbol{A}$. The Schur complement $\boldsymbol{S}$ represents the operator $-\Delta_{h}$. Hence, (19b) suggests $\boldsymbol{N}\approx\Delta_{h}^{-1}$, and the second block of the preconditioner should approximate the inverse of a discrete Laplacian in pressure space. Then, we can use a solver for elliptic operators such as multigrid to apply $\boldsymbol{N}$. Nevertheless, it is a concern that our pressure belongs to $L^{2}(\Omega)$, and strictly speaking we do not have the required regularity to apply $\Delta^{-1}_{h}$. An option to deal with this problem is to use the auxiliary space technique [51], in which the idea is to use a multigrid preconditioner for continuous pressure elements and then project it in to the desired space of discontinuous pressure elements in combination with a suitable smoothing operator. We note that some approaches only use a one-sided factorization of (17), which leads to a block-triangle form of $\mathcal{B}$; e.g. [18, 28]. For problems in which the (1,1) block of $\mathcal{A}_{h}$ has a non-symmetric term, this variant offers a faster convergence of the iterative solver at the price of evaluating an additional sparse matrix vector product compared to the block diagonal preconditioner; see [22].

In the context of the (Navier-)Stokes equations, dedicated approximations to the inverse of the Schur complement have been proposed [18, 17]. These include the pressure Schur complement methods [44] and the least squares commutator preconditioner a.k.a. BFBt and its extensions [16]. Following the presentation from [36], the BFBt approximation of the inverse of the Schur complement can be written

[TABLE]

where $\boldsymbol{C}$ and $\boldsymbol{D}$ are diagonal and symmetric positive definite matrices. The original choice sets $\boldsymbol{C}$ and $\boldsymbol{D}$ to the lumped velocity mass matrix of the system. New modifications have been introduced in an effort to improve the effectiveness of the preconditioners in cases where the equations present high variability in the (scalar) coefficients [35, 36]. For Poisson instead of Stokes, $\boldsymbol{A}$ is the velocity mass matrix instead of a discrete Laplacian. Hence, with $\boldsymbol{C}=\boldsymbol{D}\approx\boldsymbol{A}$ the method reduces to $\boldsymbol{S}_{BFBt}\approx\boldsymbol{S}$, the usual Schur complement.

3 Multigrid for Saddle Point Problems

Multigrid (MG) methods provide efficient preconditioners for important classes sparse of linear systems, in particular if the matrix system arises from the discretization of an elliptic PDE. Their main advantage is that they scale linearly in the number of unknowns $N$, i.e., that they require only $O(N)$ computational work and memory. Multigrid has been primarily developed for symmetric positive M-matrices as they typically arise from FD/FV/FE discretizations of (scalar) second order elliptic PDEs.

As already indicated in the previous section, multigrid algorithms can be used inside Schur complement preconditioners to invert one (or both, depending on the application) of the diagonal blocks. While this approach allows for an easy reuse of existing and well-established techniques, it does no longer guarantee linear convergence for the block system as a whole. The multigrid cycles are applied only to the sub-problem(s), while couplings between the unknowns are handled by the outer iteration (for example GMRES, BiCGStab or (inexact) Uzawa). In consequence, the outer iteration essentially determines the overall convergence speed, even if the inner sub-blocks can be solved quickly.

The question is whether it is possible to build a multigrid hierarchy for the coupled system to take into account the cross-couplings on all levels. Several developments have been made in this direction. Geometric multigrid methods for (Navier-)Stokes have been proposed in e.g. [45, 39]. In [46, 47, 48], this has been extended to a “semi-algebraic” AMGe, where the coarse levels are still determined geometrically. In this section, we will present the fully algebraic method from the original thesis [30]. As in classical AMG for M-matrices, this method only requires the matrix for building a robust multigrid hierarchy and hence can adapt itself to difficulties such as anisotropic or jumping coefficients as well as non-uniform meshes that cannot be coarsened geometrically.

3.1 Algebraic multigrid (AMG)

The whole hierarchy of grids $\left\{\Omega_{l}\right\}_{l=1}^{L}$, differential operators $\left\{{\boldsymbol{\Lambda}}_{l}\right\}_{l=1}^{L}$, and transfer operators $\left\{{\boldsymbol{P}}_{l}\right\}_{l=1}^{L-1}$, $\left\{{\boldsymbol{R}}_{l}\right\}_{l=1}^{L-1}$ needed in the multigrid cycle (the solution phase) is computed from the system matrix ${\boldsymbol{\Lambda}}$ during a setup phase. In this document we will understand the term “grid” as a set of indices, since AMG methods require no geometric mesh information. The strength of AMG is its ability to deal with difficulties in the operator, such as heterogeneous, strongly varying coefficients as well as unstructured meshes, for which a hierarchy cannot be (easily) identified. To set the stage for our proposed method, let us briefly recapitulate the classical Ruge-Stüben AMG setup [37, 42]. We start on the finest level $l=1$ using the fine system matrix ${\boldsymbol{\Lambda}}_{1}:={\boldsymbol{\Lambda}}$ and the finest set of degrees of freedom $\Omega_{1}=\{1,\ldots,n_{1}\}$.

We decompose the grid $\Omega_{l}$ into the set of fine grid points $F_{l}$ and coarse grid points $C_{l}$. The latter form the next coarser grid $\Omega_{l+1}:=C_{l}$ of size $n_{l+1}$. 2. 2.

We compute the interpolation matrix

[TABLE]

The submatrix $\boldsymbol{P}_{l,F}$ contains the interpolation weights for all fine grid points $i\in F_{l}$. For the coarse grid points $i\in C_{l}$, interpolation is just the identity. 3. 3.

We compute the next coarser matrix by the Galerkin product

[TABLE]

The algorithm is then recursively applied to the input matrix $\boldsymbol{\Lambda}_{l+1}$. If $n_{l}$ is small enough such that $\boldsymbol{\Lambda}_{l}$ can be efficiently factored by a direct solver, the recursion is terminated. In our experiments we use a redundant serial solver if $n_{l}\leq 1000$.

We recall that for a fast multigrid algorithm, the error components not efficiently reduced by smoothing (the smooth vectors) must be well represented within the range of the interpolation operator $\boldsymbol{P}_{l}$. In AMG, we take a simple smoothing scheme like Jacobi or Gauss-Seidel relaxation. For symmetric positive definite M-matrices, an investigation of the smooth error components $\boldsymbol{e}=(e_{1},\ldots,e_{n})^{\top}$ reveals that these satisfy $\boldsymbol{D}^{-1}\boldsymbol{\Lambda}\boldsymbol{e}\approx 0$, where $\boldsymbol{D}:=\mathrm{diag}(\boldsymbol{\Lambda})$. If we denote the coefficients of $\boldsymbol{\Lambda}$ by $\lambda_{ij}$, this means

[TABLE]

Equation (23) already delivers us a template for the interpolation formula: The value at the fine grid point $i$ should be approximated by the values at those points $j$ for which $-\lambda_{ij}$ is relatively large, while those (scaled) entries $\lambda_{ij}$ also provide the interpolation weights. In consequence, a substantial amount of those large negative connections should lead (directly or indirectly) to coarse grid points $j\in\boldsymbol{C}_{l}$. This imposes certain conditions on the selection of the coarse grid points $\boldsymbol{C}_{l}\subset\Omega_{l}$. Many algorithms to the coarse grid selection and construction of the interpolation have been proposed. We do not describe them in detail here, but refer to [37, 42, 24, 41, 40]. In our experiments we use Falgout coarsening to select the set of coarse grid points $\mathfrak{C}_{l}$. We compute the interpolation weights $\boldsymbol{P}_{l,\mathfrak{F}}$ using the modified classical interpolation scheme and dropped interpolation weights smaller than 1/20 of the largest absolute weight per row. The restriction matrix is taken as the transpose of the interpolation.

The Galerkin ansatz for the coarse grid matrix (22) has two benefits: First, for symmetric positive definite $\boldsymbol{\Lambda}_{l}$, $\boldsymbol{\Lambda}_{l+1}$ is also symmetric positive definite for any full (column) rank interpolation operator $\boldsymbol{P}_{l}$. Second, the resulting two-grid correction operator $\boldsymbol{I}-\boldsymbol{P}_{l}\boldsymbol{\Lambda}_{l+1}^{-1}\boldsymbol{P}_{l}\boldsymbol{\Lambda}_{l}$ is an orthogonal projector, which (extended recursively over all levels and combined with smoothing) ensures that the multigrid cycle converges [42].

We now present an algebraic multigrid approach for a monolithic solution of saddle point problems of the form

[TABLE]

We will construct a multigrid hierarchy of saddle point matrices indexed by $l$,

[TABLE]

as well as block interpolation matrices $\mathcal{P}_{l}$. Note that we include a lower right block $\boldsymbol{C}_{l}$, which will be zero on the finest level, $\boldsymbol{C}_{1}=0$, and non-zero for $l\geq 2$.

3.2 Smoothers for saddle point problems

Classical relaxation schemes like the Jacobi or Gauss-Seidel iteration are not suitable for saddle point systems, since the smoothing properties for these schemes rely on the positive definiteness of the matrix. We present two dedicated smoothing schemes for saddle point systems [39]: First, a predictor-corrector algorithm, which combines segregated sweeps over the physical components, and second, a box relaxation scheme, where small saddle point subsystems are solved within a global Schwarz method.

Uzawa relaxation.

Our first option is a symmetric inexact Uzawa relaxation scheme. Within each iteration, a predictor $\boldsymbol{u}^{\ast}$ for the flux is computed first, which is used to update $\boldsymbol{u}$. Finally, employing the updated values of $p_{l}^{it+1}$, the next iterate $\boldsymbol{u}_{l}^{it+1}$ is computed,

[TABLE]

The matrices $\boldsymbol{\hat{A}}_{l}$ and $\boldsymbol{\hat{S}}_{l}$ are chosen such that $\boldsymbol{\hat{A}}_{l}-{\boldsymbol{A}}_{l}$ and $\boldsymbol{\hat{S}}_{l}-{\boldsymbol{B}}_{l}\boldsymbol{\hat{A}}_{l}^{-1}{\boldsymbol{B}}_{l}^{\top}-{\boldsymbol{C}}_{l}$ are symmetric positive definite. Furthermore, they should be easily invertible. For example, we can choose to use the scaled diagonals of $\boldsymbol{A}_{l}$ and $\boldsymbol{B}_{l}\boldsymbol{\hat{A}}_{l}^{-1}\boldsymbol{B}_{l}^{\top}+\boldsymbol{C}_{l}$, respectively. The magnitude of the scaling can be obtained with help of the power iteration on the latter matrices. For convergence and further properties of this smoother we refer to [39].

Vanka smoothing.

The second alternative requires us to first decompose the computational domain into small overlapping patches. To this end, we employ the non-zero structure of $\boldsymbol{B}_{l}$ and construct a patch $\Omega_{l,j}$ for each row of $\boldsymbol{B}_{l}$: For the $j$-th row, $\Omega_{l,j}$ consists of the index $j$ as well as all indices $i$ such that there exists an entry $b_{ji}\neq 0$,

[TABLE]

The transfer between the global domain $\Omega_{l}$ and the subdomains $\Omega_{l,j}$ is accomplished by (optionally scaled) injection operators $\boldsymbol{V}_{l,j}$ (flux) and $\boldsymbol{W}_{l,j}$ (pressure),

[TABLE]

where $\boldsymbol{J}_{j}$ and $\boldsymbol{\tilde{J}}_{j}$ are binary matrices that map the subdomain into the global domain for velocity and pressure, respectively; i.e., each of their columns contains exactly one unit entry.

First, the residuals are restricted to the subdomain,

[TABLE]

Then, on each subdomain, a small saddle point problem of the form

[TABLE]

is solved using a sparse direct solver. Finally, the updates are prolongated to the global domain,

[TABLE]

The small matrices $\boldsymbol{\hat{A}}_{l,j}$, $\boldsymbol{B}_{l,j}$ and $\boldsymbol{\hat{S}}_{l,j}$, where the latter is just a scalar, are defined by

[TABLE]

Again, $\boldsymbol{\hat{A}}_{l}$ is a scaled version of the diagonal of ${\boldsymbol{A}}_{l}$ such that $\boldsymbol{\hat{A}}_{l}-{\boldsymbol{A}}_{l}$ is positive definite. With the aid of a power iteration, $\beta>0$ is computed such that

[TABLE]

is symmetric positive definite.

The iteration can be performed either additively (i.e., the residuals are computed once, then all subdomain solves are performed independently) or multiplicatively (after each subdomain solve the residuals are updated). In the latter case, it is beneficial to perform a symmetric sweep, i.e., after one complete sweep the subdomain solves are performed again in reverse order. If we choose

[TABLE]

the additive smoother coincides with the Uzawa method described in the previous section [39]. On the other hand, in the multiplicative case the simple choice

[TABLE]

may result in faster convergence.

3.3 AMG setup for saddle point systems

The starting point for our saddle point AMG method (SPAMG) is the block diagonal matrix

[TABLE]

This operator is symmetric positive definite. Let us assume that we can apply the classical AMG setup algorithm as described in Section 3.1 to each of the blocks $\boldsymbol{A}_{l}$ and $\boldsymbol{B}_{l}\boldsymbol{\hat{A}}_{l}^{-1}\boldsymbol{B}_{l}^{\top}+\boldsymbol{C}_{l}$. We construct coarse grids and interpolation operators $\boldsymbol{P}_{l,u}$ and $\boldsymbol{P}_{l,p}$ for each of these blocks and obtain a block interpolation operator

[TABLE]

Now, one idea would be to construct the coarse grid operator in the usual Galerkin way, $\mathcal{A}_{l+1}:=\mathcal{P}^{\top}_{l}\mathcal{A}_{l}\mathcal{P}_{l}$. Unlike in the symmetric positive definite case, however, we cannot be sure that $\mathcal{A}_{l+1}$ is invertible, thus we might attain an unusable multigrid hierarchy. To prevent this, we modify the interpolation operator $\mathcal{P}_{l}$ by the multiplication of a stabilization term. To this end, we re-write the velocity prolongation block such that the rows corresponding to the fine grid points come first,

[TABLE]

where $\boldsymbol{I}_{l,u,C}$ is the identity injecting the values from level $l+1$ to level $l$ and $\boldsymbol{P}_{l,u,CF}$ contains the interpolation weights from coarse to fine computed from the matrix $\boldsymbol{A}_{l}$. Analogously, we write

[TABLE]

The stabilized interpolation operator is computed by

[TABLE]

Now, we use the modified Galerkin product

[TABLE]

to compute the coarse grid operator. For this matrix, we can show an inf-sup-condition if (a) the fine grid matrix $\mathcal{A}_{l}$ fulfils such a condition and (b) the interpolation operator $\boldsymbol{P}_{l,u}$ for the velocity satisfies certain approximation properties, which most usual AMG interpolation schemes do [30, Lemma 4.6]. The invertibility of the coarse grid matrix (41) is ensured by its lower right block

[TABLE]

This can be understood as a partial Schur complement that ensures the stability of the coarse grid matrix.

In reference to recent work, the results published in [49] and [50] indicate that it is neither always necessary to stabilize the interpolation operator on every level, nor that interpolation and restriction need to be stabilized simultaneously. A further investigation of these techniques, however, is beyond the scope of this paper: It may be quite involved to determine algorithmically and separately for each level whether the stabilization is required.

4 Numerical Results

In this section, we evaluate the effectiveness of SPAMG compared to the diagonal and Schur preconditioners for uniform and adaptive meshes. We choose the examples based on manufactured solutions, that is, we prescribe a target pressure field and compute the velocity field, boundary conditions, and right hand sides in order to satisfy (7). We also choose synthetic adaptive refinement tailored to the reference solutions. Examples 1 to 3 are typical benchmarks with smooth coefficients. In example 4 we challenge the numerical solver and the preconditioner by using a conductivity tensor with strong coefficient variation. We employ three different flavors of smoothers inside SPAMG: The inexact Uzawa scheme (26) and two variants of the multiplicative Vanka smoother, either using the scaling (34) (“Vanka Scale”) or the unscaled injection (35) (“Vanka One”).

We delegate the parallel mesh management to the octree-based AMR software library p4est [14, 25]. This library provides a collection of algorithms that implement scalable parallel AMR operations. In particular, we employ p4est to create and modify a hexahedral triangulation of the unit square/cube and to introduce a numbering of the degrees of freedom suitable for a $\mathcal{RT}_{0}$ discretization. Regarding solvers and preconditioners, we use the GMRES [38] implementation provided by the software library hypre [43]. In order to approximate the product $\boldsymbol{N}\boldsymbol{r}$ for a given residual $\boldsymbol{r}$, we employ one SPAMG V-cycle. We use the parallel multigrid implementation BoomerAMG from hypre as the frame to which our saddle-point AMG is added.

4.1 Homogeneous Dirichlet boundary conditions

We solve (7) on the unit cube using an identity conductivity tensor and homogeneous Dirichlet boundary conditions. We compute the right hand side based on the manufactured solution

[TABLE]

in 2D and

[TABLE]

in 3D, respectively.

With this example we verify the correctness of our implementation of the $\mathcal{RT}_{0}$ discretization. The discretization error converges the predicted rates as confirmed in Figure 3 and Figure 4. The iteration counts displayed in Table 1 and Table 2 confirm the (well known) robustness of the Schur preconditioner for uniform meshes. For adaptive meshes, the Schur complement preconditioner incurs high iteration counts, particularly for the 3d case. The three variants of the SPAMG preconditioner retain mesh independent iteration counts for both uniform and adaptive meshes.

4.2 Inhomogeneous Dirichlet/Neumann boundary conditions

We solve (7) with an identity tensor. We impose homogeneous mixed homogeneous Dirichlet / Neumann boundary conditions and compute the right hand side based on the exact solution

[TABLE]

in 2D and

[TABLE]

in 3D, respectively. The Neumann boundary is set at $y=0$ and $y=1$ in both cases.

As in the previous section, the theoretical converge rates agree with the bounds (15). Due to the smoothness solution and the equation constant coefficient, no additional benefit is expected from local adaptation of the mesh; see Figure 5 and Figure 6. The iteration counts (not shown) are similar to the previous example.

4.3 A non-trivial conductivity tensor

In this example we approximate the solution of (7) for the case of a non-diagonal conductivity tensor $\boldsymbol{\mathcal{K}}$. We impose non-zero Dirichlet boundary conditions. The manufactured solution is

[TABLE]

in 2D and

[TABLE]

in 3D, respectively. The conductivity tensor is given by

[TABLE]

in 2D and

[TABLE]

in 3D.

The three variants of the SPAMG preconditioner offer mesh independent iteration counts for uniform and adaptive meshes. The Schur complement again produces growing iteration counts, in particular for the three dimensional case and even for uniform meshes. See Table 3 and Table 4.

4.4 High conductivity contrast

To further examine the robustness of the SPAMG preconditioner, we solve (7) with a conductivity tensor that exhibits strong coefficient variation. We enforce inhomogeneous mixed Dirichlet boundary conditions and compute the right hand side terms based on the manufactured solution

[TABLE]

in 2D and

[TABLE]

in 3D. The conductivity Tensor $\boldsymbol{\mathcal{K}}$ is given by the identity matrix scaled pointwise with a continuously differentiable function $m(\boldsymbol{x};\boldsymbol{x}_{0},a,b,c)$ constructed to fulfil the following properties:

• •

$m(\boldsymbol{x};\boldsymbol{x}_{0},a,b,c)=1\>$ if $\lVert\boldsymbol{x}-\boldsymbol{x}_{0}\rVert\geq b$,

• •

$m(\boldsymbol{x};\boldsymbol{x}_{0},a,b,c)=1-c\>$ if $\lVert\boldsymbol{x}-\boldsymbol{x}_{0}\rVert\leq a$ and

• •

$m(\boldsymbol{x};\boldsymbol{x}_{0},a,b,c)\in(1-c,1)\>$ if $a<\lVert\boldsymbol{x}-\boldsymbol{x}_{0}\rVert<b$.

Such a function can be constructed by defining

[TABLE]

and

[TABLE]

see Figure 7.

Hence, the first three arguments define the location and radius of the support of $m$.

The parameter $c$ allows us to tune the coefficients to vary across the domain. From an application point of view, if $c$ is close to one the function $m$ models a medium in which almost no flow is allowed within a circle (sphere) centered at $\boldsymbol{x}_{0}$ with radius $a$. Given this information, it is clear that the velocity field is likely to have a strong gradient within the ring-shaped region $a<\lVert\boldsymbol{x}-\boldsymbol{x_{0}}\rVert<b$. We provide Figure 8 and Figure 9 to illustrate this, refining a given element on the mesh whenever it overlaps this region.

In fact, adaptive refinement is necessary in this example to obtain optimal convergence, which we show in Figure 10 (2D) and Figure 11 (3D). The plot of accuracy versus degrees of freedom in Figure 12 supports this observation.

We display the iteration counts for various preconditioners in Table 5 (2D) and Table 6 (3D). While the Schur preconditioner still functions in 2D, it fails completely in three space dimensions. All three variants of SPAMG, on the other hand, lead to almost mesh independent iteration counts. The “Vanka One” variant seems best with an iteration count of just 10 to reduce the relative error of the linear system of equations by 6 orders of magnitude. Thus, SPAMG is robust against coefficient functions whose magnitude varies by at least a factor of 1000 in a narrow region.

5 Conclusion

The purpose of this paper is to describe the SPAMG multigrid preconditioner for saddle point systems. It is a monolithic AMG method applicable to the block system arising from mixed Poisson and Stokes discretizations. One key element of the construction is a coupled prolongation operator that stabilizes the Galerkin product from one level to the next coarser one.

Our numerical examples are obtained from a Raviart-Thomas discretization of the mixed Poisson system. We show the correctness of the solver and preconditioner by reproducing optimal convergence rates and, more importantly, demonstrate mesh-independent iteration counts of a preconditioned GMRES solve. We see that the SPAMG preconditioner offers robustness against highly graded adaptive meshes in both two and three space dimensions and tolerates strongly varying coefficients of the PDE as well as a non-trivial, matrix-valued conductivity. In addition, we show that such is not achieved by a standard Schur complement block preconditioner.

Essentially, mixed adaptive discretizations are handled well by SPAMG. Some open questions remain, such as studying the robustness of the preconditioner for even higher conductivity contrast on the order of $10^{6}$ and extending it to largely disparate eigenvalues in the conductivity tensor (i.e., problems involving anisotropy). Investigating the behavior of SPAMG for higher-order RT discretizations is another possible extension.

Acknowledgments

Authors B. and F. gratefully acknowledge financial support by the SFB/TR 32 “Patterns in Soil-Vegetation-Atmosphere Systems: Monitoring, Modeling, and Data Assimilation” funded by the Deutsche Forschungsgemeinschaft (DFG), as well as travel support from the DFG-funded Hausdorff Center for Mathematics, Bonn, Germany.

Author M. acknowledges the support of the former SFB611 “Singular Phenomena and Scaling in Mathematical Models” funded by the Deutsche Forschungsgemeinschaft (DFG).