Robust multigrid methods for isogeometric discretizations of the Stokes equations

TL;DR

This paper explores extending a robust multigrid solver, previously successful for Poisson problems, to efficiently handle the more complex Stokes equations in isogeometric analysis.

Contribution

It investigates the feasibility of adapting a proven multigrid method from Poisson to Stokes equations within isogeometric discretizations.

Findings

Multigrid method shows promise for Stokes equations

Convergence rates remain robust in grid size and polynomial degree

Extension from Poisson to Stokes is feasible

Abstract

In recent publications, the author and his coworkers have proposed a multigrid method for solving linear systems arizing from the discretization of partial differential equations in isogeometric analysis and have proven that the convergence rates are robust in both the grid size and the polynomial degree. So, far the method has only been discussed for the Poisson problem. In the present paper, we want to face the question if it is possible to extend the method to the Stokes equations.

Tables

Table 1. Table 1: Iteration counts for the unit square

		Taylor-Hood					Nédélec					Raviart-Thomas
	$\ell ╲p$		2	3	5	8		2	3	5	8		2	3	5	8
MINRES, preconditioned with SCMS-MG
	5		55	54	49	46		80	74	68	55		44	36	35	29
	6		54	58	53	51		76	76	70	63		44	37	36	32
	7		54	54	54	53		76	76	71	65		45	37	33	29
	8		50	51	55	OoM		71	71	67	65		41	37	33	29
MINRES, preconditioned with standard GS-MG
	5		64	167	>1k	>1k		84	213	>1k	>1k		124	219	>1k	>1k

Table 2. Table 2: Problem size and L 2 subscript 𝐿 2 L_{2} -errors for the unit square

Taylor-Hood

Nédélec

Raviart-Thomas

$p$

$\ell$

dof

$v$

$p$

dof

$v$

$p$

dof

$v$

$p$

2

4

2372

2e-5

1e-5

1637

2e-5

4e-5

869

3e-4

3e-2

5

9348

1e-6

6e-7

6341

1e-6

4e-5

3269

3e-5

2e-2

6

37124

7e-8

4e-6

24965

7e-7

9e-5

12677

7e-6

2e-2

7

147972

2e-8

7e-7

99077

8e-7

1e-4

49925

3e-6

2e-2

5

4

2891

2e-9

4e-8

2066

6e-8

2e-6

1202

9e-7

6e-3

5

10347

2e-9

1e-7

7154

8e-8

5e-6

3890

1e-6

3e-3

6

39083

3e-9

2e-7

26546

9e-7

2e-4

13876

6e-7

3e-3

7

151951

7e-9

4e-7

102194

2e-6

3e-4

52274

6e-7

4e-3

Table 3. Table 3: Iteration counts for the quarter annulus (direct composition)

		Taylor-Hood					Nédélec					Raviart-Thomas
	$\ell ╲p$		2	3	5	8		2	3	5	8		2	3	5	8
MINRES, preconditioned with SCMS-MG
	5		195	190	185	172		257	246	244	206		244	139	128	116
	6		208	217	213	199		295	296	280	241		192	170	142	129
	7		220	222	232	219		329	330	314	281		213	195	158	140
	8		231	239	244	OoM		333	342	333	306		223	200	168	149
MINRES, preconditioned with SCMS-MG-geo
	5		72	69	68	72		69	69	65	63		73	62	53	56
	6		77	75	73	79		76	74	64	70		71	69	59	63
	7		72	71	70	84		79	70	68	74		75	74	64	69
	8		74	73	72	OoM		73	73	71	78		71	70	68	74
MINRES, preconditioned with standard GS-MG
	5		70	173	>1k	>1k		110	225	>1k	>1k		182	220	>1k	>1k

Table 4. Table 4: Iteration counts for the quarter annulus (Piola transform)

		Taylor-Hood					Nédélec					Raviart-Thomas
	$\ell ╲p$		2	3	5	8		2	3	5	8		2	3	5	8
MINRES, preconditioned with SCMS-MG
	5		331	331	338	317		288	313	332	305		480	309	295	300
	6		407	400	402	371		361	387	405	374		368	344	323	299
	7		452	455	455	450		413	450	476	476		418	395	367	341
	8		487	485	500	OoM		458	494	556	568		441	438	411	361
MINRES, preconditioned with standard GS-MG
	5		70	165	>1k	>1k		69	164	>1k	>1k		206	199	>1k	>1k

Equations

$$-\Delta u+\nabla p=f\text{mbox}and\nabla \cdot u=0$$

$$\underset{a(u,v)}{\underbrace{(\nabla u,\nabla v)}}+\underset{b(v,p)}{\underbrace{(\nabla \cdot v,p)}}=(f,v)\forall v\in V,\underset{b(u,q)}{\underbrace{(\nabla \cdot u,q)}}=0\forall q\in Q.$$

$$\underset{q\in L_2(\Omega )}{\inf }\underset{v\in H^1(\Omega )}{\sup }\frac{(\nabla \cdot v,q)}{\Vert v{\Vert }_{H^1(\Omega )}\Vert q{\Vert }_{L_2(\Omega )}}\ge C,$$

$$\underset{q_h\in Q_h}{\inf }\underset{v_h\in V_h}{\sup }\frac{(\nabla \cdot v_h,q_h)}{\Vert v_h{\Vert }_{H^1(\Omega )}\Vert q_h{\Vert }_{L_2(\Omega )}}\ge C$$

$$A_{h}\,\underline{x}_{h}=\underline{f}_{h},\qquad\mbox{where}\qquad A_{h}=\left(\begin{array}[]{cc}K_{h}&D_{h}^{T}\\ D_{h}&0\end{array}\right)\qquad\mbox{and}\qquad\underline{x}_{h}=\left(\begin{array}[]{c}\underline{u}_{h}\\ \underline{p}_{h}\end{array}\right)$$

$$S_{p_1,p_2,h}^{q_1,q_2}:=S_{p_1,h}^{q_1}\otimes S_{p_2,h}^{q_2},$$

$$G:\hat{\Omega }=(0,1{)}^2\to \Omega ,$$

$$\begin{array}[]{lll}\hat{X}_{h}^{(\textnormal{TH})}:=\hat{V}_{h}^{(\textnormal{TH})}\times\hat{Q}_{h},&\qquad\hat{V}_{h}^{(\textnormal{TH})}:=S_{p+1,p+1}^{p-1,p-1}\times S_{p+1,p+1}^{p-1,p-1},&\\ \hat{X}_{h}^{(\textnormal{NE})}:=\hat{V}_{h}^{(\textnormal{NE})}\times\hat{Q}_{h},&\qquad\hat{V}_{h}^{(\textnormal{NE})}:=S_{p+1,p+1}^{p,p-1}\times S_{p+1,p+1}^{p-1,p},&\\ \hat{X}_{h}^{(\textnormal{RT})}:=\hat{V}_{h}^{(\textnormal{RT})}\times\hat{Q}_{h},&\qquad\hat{V}_{h}^{(\textnormal{RT})}:=S_{p+1,p}^{p,p-1}\times S_{p,p+1}^{p-1,p},&\qquad\hat{Q}_{h}:=S_{p,p}^{p-1,p-1},\\ \end{array}$$

$${\hat{V}}_h^{(\text{RT})}\subset {\hat{V}}_h^{(\text{NE})}\subset {\hat{V}}_h^{(\text{TH})}$$

$$V_h^{(X,\text{D})}:=\{v_h|v_h\circ G\in {\hat{V}}_h^{(X)}\},X\in \{\text{TH},\text{NE},\text{RT}\}.$$

$$V_h^{(X,\text{P})}:=\left\{v_h|\frac{1}{\text{mbox}det{J}_G}{J}_G{v}_h\circ G\in {\hat{V}}_h^{(X)}\right\},X\in \{\text{TH},\text{NE},\text{RT}\},$$

$$Q_h:=\{q_h|q_h\circ G\in {\hat{Q}}_h\}.$$

$$Q_{h}:=\left(\begin{array}[]{cc}K_{h}&\\ &\beta^{-1}M_{h}\end{array}\right),$$

$$u_{h}(x,y)=\left(\begin{array}[]{c}\cos(5x+5y)+\sin(5x-5y)\\ -1-\cos(5x+5y)+\sin(5x-5y)\end{array}\right),$$

Full text

1 Introduction

Isogeometric analysis (IgA) was introduced by Tom Hughes et al. in the 2005, cf. [11], aiming to improve the connection between computer aided design (CAD) and finite element (FEM) simulation. In IgA, as in CAD software, B-splines and non-uniform rational B-splines (NURBS) are used for representing both the geometrical objects of interest and the solution of the partial differential equation (PDE) to be solved.

In IgA, mostly B-splines or NURBS of maximum smoothness are used, i.e., having a spline degree of $p$, the functions are $p-1$ times continuously differentiable. Using such a function space, one obtains on the one hand the approximation power of high order functions, while on the other hand, unlike in standard high-order FEM, one does not suffer from a growth of the number of degrees of freedom.

From the computational point of view, the treatment of the linear systems arizing from the discretization with high spline degrees is still challenging as the condition number both of mass and stiffness matrices grows exponentially with the spline degree. In the early IgA literature, often finite element solvers have been transferred to IgA only with minimal adaptations. Numerical experiments indicate that such approaches result in methods that work well for small spline degrees, but their performance deteriorates as the degree is increased, often dramatically. In [10, 9], the author and his coworkers have proposed multigrid methods which are provable robust in the polynomial degree and the grid size. Numerical experiments indicated that the proposed approach of subspace corrected mass smoothers seems to pay off (compared to multigrid methods with a standard Gauss-Seidel smoother) for polynomial degrees of four or five.

In the present paper, we want to discuss the extension of the subspace corrected mass smoothers beyond the case of the Poisson problem to the Stokes flow problem. Unlike for the Poisson problem, for the Stokes problem already the setup of a stable isogeometric discretization is non standard. As there have already been results in the literature, we refer to paper [1], which serves as a basis of the present paper. Alternative approaches can be found in [7, 2, 5] and others. After introducing discretizations, we discuss the setup of the preconditioner.

For the Poisson problem, the multigrid solver has been applied directly and as a preconditioner for the conjugate gradient method. For the case of a non-trivial geometry transformation, in [9] a conjugate gradient method, preconditioned with the multigrid method for the parameter domain, has been used. It has been shown that in this case the resulting method is robust both the grid size and the polynomial degree, but not in the geometry transformation.

There are a few approaches how to carry this over to the Stokes equations. The first possibility is to apply the multigrid method directly to the problem of interest (all-at-once multigrid method), cf. [14] for a particularly popular method in standard FEM or [12] for a survey. As the results for the Poisson problem have indicated that a direct application of the multigrid method in the presence of a non-trivial geometry transformation is not optimal, we do not concentrate on that case.

So, we consider a Krylov space method with an appropriate preconditioner, living on the parameter domain. In principle, this could be the Stokes problem on the physical domain, but such a choice (an indefinite preconditioner for an indefinite problem) typically requires the use of a GMRES method, whose convergence cannot be easily proven by considering the spectrum of the preconditioned system, cf. [8]. So, we consider elliptic preconditioners, particularly block-diagonal preconditioners. As the Stokes equations are well-posed in the Sobolev space $H^{1}$ (velocity) and the Lebesgue space $L_{2}$ (pressure), we just setup preconditioners for those spaces (operator preconditioning).

As we observe that the subspace corrected mass smoothers suffer significantly from the geometry transformation, we propose a variant (by incorporating an approximation to the geometry transformation) which leads for several experiments to a significant speedup.

This paper is organized as follows. We will introduce the particular model problem in section 2 and discuss three kinds of discretizations for the mixed system in section 3. As a next step, in section 4, we propose a preconditioner. Finally, in section 5, we give the results of the numerical examples and draw some conclusions.

2 Model problem

Let $\Omega\subseteq\mathbb{R}^{2}$ be a simply connected domain with Lipschitz boundary $\partial\Omega$ and assume a force filed $f$ to be given on $\Omega$ and boundary data to be given on $\partial\Omega$. The Stokes flow model problem reads as follows. Find the velocity field $u$ and the pressure distribution $p$ such that

[TABLE]

hold on $\Omega$ and Dirichlet boundary conditions hold on $\partial\Omega$. After homogenization, we obtain a mixed variational form, which reads as follows. Find $u\in V:=H^{1}_{0}(\Omega)$ and $p\in Q:=L_{2}(\Omega)$ such that

[TABLE]

Here, and in what follows $L_{2}(\Omega)$, $H^{1}(\Omega)$ and $H^{1}_{0}(\Omega)$ are the standard Lebesgue and Sobolev spaces, and $(\cdot,\cdot)$ is the standard norm on $L_{2}(\Omega)$.

Existence and uniqueness of the solution and its dependence of the data follows from Brezzi’s theorem [3], which requires besides boundedness and $H^{1}$-coercivity of $a$ the inf-sup stability

[TABLE]

which is known to be satisfied for the Stokes problem, cf. [4].

3 Discretization

The discretization is done using a standard Galerkin approach, i.e., we replace the spaces $V$ and $Q$ by finite-dimensional subspaces $V_{h}$ and $Q_{h}$. As for the continuous problem, existence and uniqueness of the solution can be shown by Brezzi’s theorem. Boundedness and $H^{1}$-coercivity of $a$ follow directly from the continuous problem, but the inf-sup stability for the discrete problem does not. Therefore, we have to guarantee that the discrete inf-sup condition

[TABLE]

is satisfied, which is actually a condition on the discretization. In the subsection 3.2, we will discuss discretizations satisfying this condition.

Assuming a particular discretization and a basis for the chosen space, one ends up with a linear system to be solved: For given $\underline{f}_{h}$, find $\underline{x}_{h}$ such that

[TABLE]

and $K_{h}$ is a standard stiffness matrix and $D_{h}$ is a matrix representing the divergence.

3.1 Discretization in isogeometric analysis

Let $S_{p,h}^{q}$ be the space of all $q$ times continuously differentiable functions on $(0,1)$, which are piecewise polynomials of degree $p$ on a (uniform) grid of size $h=1/n$. As a basis for $S_{p,h}^{q}$ we choose the classical basis of B-splines, see, e.g., [6].

For computational domains $\Omega\subset\mathbb{R}^{2}$, we first define the spline spaces for the parameter domain $\hat{\Omega}=(0,1)^{2}$. On the parameter domain, we introduce the space of tensor-product splines, which reads as follows:

[TABLE]

where $A\otimes B$ denotes the linear span of all functions $(x,y)\mapsto u(x)v(y)$, where $u\in A$ and $v\in B$. Note that the restriction to two dimensions and to a uniform grid is only for ease of notation. The extension to three and more dimensions or to non-uniform grids is completely straight-forward. Assuming that physical domain $\Omega$ is the image of a B-spline or NURBS mapping

[TABLE]

we define the spline spaces on the physical domain typically using a classical pull-back principle. More complicated domains are represented patch-wise, where for each patch a separate geometry transformation $G$ exists. For simplicity, we do not discuss that in the present paper.

3.2 Stable discretizations for the Stokes problem

As mentioned above, it is required to set up the discretization such that the discrete inf-sup condition holds. We discuss this first for the parameter domain. Here, we follow the outline of the paper [1], where three spline space configurations have been proposed, which are variants of known stable spaces from standard finite elements: Taylor-Hood like splines $\hat{X}_{h}^{(\textnormal{TH})}$, Nédélec like splines $\hat{X}_{h}^{(\textnormal{NE})}$ and Raviart-Thomas like splines $\hat{X}_{h}^{(\textnormal{RT})}$. All of them utilize the same grid for both the velocity and the pressure, which makes the implementation significantly easier compared to approaches that are based on setting up two different grids (like IgA-variants of the macro elements, cf. [2]). All of these discretizations follow the spirit of IgA, allowing to freely choose the underlying polynomial degree $p$. For all of them, the smoothness is in the order of the polynomial degree, which preserves the feature that the number of degrees of freedom is basically not increased when the polynomial degree is increased.

For the case of two dimensions, the spaces are given by

[TABLE]

where $A\times B:=\{(a,b):a\in A,\;b\in B\}.$ Observe that these spline spaces are nested, i.e., we have

[TABLE]

and (for $n>>p$) a ratio of $9:5:3$ for the number of degrees of freedom. The extension of these definitions to three and more dimensions is straight-forward, cf. [1].

For all of these settings, the discrete inf-sup condition has been shown in [1]. For the Raviart-Thomas like splines, the discrete inf-sup condition cannot be proven if Dirichlet boundary conditions are present. As the method still seems to work well in practice, we include also the Raviart-Thomas discretization in our experiments.

The next step is to introduce the discretization for the physical domain. As outlined in the beginning of this section, the discretization, once introduced on the parameter domain, is typically defined on the physical domain just by direct composition:

[TABLE]

For the Stokes problem, as an alternative, the divergence preserving Piola transform has been proposed:

[TABLE]

where $J_{G}$ is the Jacobi matrix of $G$. The pressure distribution, which is a scalar quantity, is always mapped directly, i.e., in all cases we choose the direct composition

[TABLE]

In [1], the inf-sup stability has been shown if the Piola transform is used, for the Taylor-Hood like splines also if the direct composition is used. Again, we report also on the numerical results for the cases that are not covered by the convergence theory (direct composition for the Nédélec like and the Raviart-Thomas like splines).

4 Robust multigrid solvers

As outlined in the introduction, the multigrid preconditioner aims to represent the theoretical block-diagonal preconditioner

[TABLE]

where $K_{h}$ is the stiffness matrix, $M_{h}$ is the mass matrix and $\beta>0$ is an accordingly chosen scaling parameter. As mentioned above and as discussed in detail in [9], we use as preconditioner for the problem on the physical domain the corresponding preconditioner, say $\hat{Q}_{h}$, on the parameter domain. There, the matrices $M_{h}$ and $K_{h}$ are replaced by $\hat{M}_{h}$ and $\hat{K}_{h}$, their counterparts on the parameter domain. Note that the stiffness matrix acts on the velocity variable, a vector-valued quantity, and that this matrix is block-diagonal on the parameter domain and, iff the direct composition is used, on the physical domain. In all cases, $K_{h}$ and $\hat{K}_{h}$ are spectrally equivalent.

Instead of an exact inverse of the matrix $\hat{Q}_{h}$, we only need to realize an approximation to the application of $\hat{K}_{h}^{-1}$ and $\hat{M}_{h}^{-1}$ to any given vector. The approximation of $\hat{K}_{h}^{-1}$ is realized using one multigrid V-cycle with one pre- and one post-smoothing step of the subspace corrected mass smoother, as proposed in [9]. There, the algorithm was analyzed only for the case of splines of maximum smoothness, however it can be applied for any spline space and robustness in the polynomial degree can be guaranteed by a slight extension of the presented theory as long as the smoothness is in the order of the polynomial degree. As in the previous publications [10, 9], the grid hierarchy is set up for a fixed polynomial degree and a fixed smoothness by just uniformly refining the grid. Using this approach, one obtains nested spaces, so the setup of the coarse-grid correction is trivial.

One of the key observations which was leading to the results in [10, 9] was that the spectral equivalence of the mass matrix and its diagonal deteriorates if $p$ is increased. This has also to be taken into account when constructing the preconditioner for the pressure variable. Analogously to the smoother, we realize the application of $\hat{M}_{h}^{-1}$ exactly, based on the tensor-product structure of the mass matrix.

The preconditioner is symmetric and positive definite and can therefore be applied in the framework of a MINRES iteration.

5 Numerical results

The numerical experiments have been performed using the C++ library G+SMO, see [13], both for the unit square, i.e., for a problem without geometry transformation, and for a quarter annulus $\{(x,y)\in\mathbb{R}_{+}^{2}:1<x^{2}+y^{2}<4\}$. For both problems, the problem has been constructed (with inhomogeneous right-hand-side and inhomogeneous Dirichlet boundary conditions) such that the exact solution is

[TABLE]

and $p_{h}(x,y)=-(1+x)(1+y)+c$, where $c$ is chosen such that $\int_{\Omega}p_{h}\,\mbox{d}x=0$.

In Table 1, we see number of MINRES steps required for reducing the initial error (measured in the $\ell^{2}$-norm of the solution vector) by a factor of $10^{-6}$; cases where the memory was not enough are indicated with OoM. The discussion is done for all proposed discretization schemes. The need of the discussion of $p$-robust methods is easily observed when looking at the results for a standard preconditioner: We display the results if one multigrid V-cycle with Gauss-Seidel smoother is used for the velocity and one symmetric Gauss-Seidel sweep is used for the pressure (GS-MG). There, the number of iterations increases drastically if $p$ is increased. As the approach is perfectly robust in the grid size $h=2^{\ell}$, we omit the numbers for finer grids. Compared to that approach, the preconditioner proposed in Section 4 (SCMS-MG) shows results which are robust both in the grid size and the polynomial degree and which works well for all discretizations. Although the iteration numbers are smaller than for the GS-MG preconditioner, one has to consider that the costs of the SCMS-MG preconditioner are significantly higher than those of the GS-MG preconditioner, so the proposed method only pays of iff higher polynomial degrees (starting from $4$ or $5$) are considered. We have chosen $\beta=0.05$ and as damping parameter $\sigma$ of the underlying smoother, cf. [9], either $\sigma^{-1}=0.04\,\hat{h}^{2}$ (for Taylor-Hood and Nédélec) or $0.16\,\hat{h}^{2}$ (for Raviart-Thomas), where $\hat{h}$ is the grid size on the parameter domain. While some of the numbers might be improved by fine-tuning the parameters, the given tables for reasonable uniform choices show what one can expect for each of the methods.

In Table 2 we see how well the computed solution approximates the exact solution in terms of the $L^{2}$-norm. Here, we have used the abovementioned solver, where the the stopping criterion has been chosen to reach either a relative error of $10^{-10}$ or $100$ iterations. We present the error between the computed solution and the known exact solution (for the pressure after projecting into the space of functions with vanishing mean). We observe that, for the same choice of the polynomial order $p$ and the same grid size, the Taylor-Hood discretization yields to the smallest errors, for the cost of the largest number of degrees of freedom. For the Raviart-Thomas discretization (where the inf-sup condition cannot be shown for the chosen Dirichlet boundary conditions), we observe that the error for the velocity converges, while the error of the pressure stagnates at around $10^{-2}$. Observe moreover that for $p=5$, the approximation on the coarsest grid was fine enough such that the approximation error could not be improved by refinement.

For the case of the quarter annulus, we distinguish between the results obtained by the direct composition (Table 3) and for the Piola transform (Table 4). Again, we obtain first that GS-MG is robust in $h$, but the convergence deteriorates if the polynomial degree grows. As it leads to better results, we have set up the GS-MG on the physical domain. For the proposed SCMS-MG preconditioner, observe that the results behave similar to the results for the unit square, however the iteration counts are much larger, particularly if the Piola transform is used. For the direct composition, it is possible to improve the convergence significantly by replacing the mass and stiffness matrix on the parameter domain by a simple tensor-rank-one approximation of those matrices on the physical domain (SCMS-MG-geo). Note that the tensor-rank-one approximation does not lead to any additional computational costs after the assembling phase. The extension of such a rank-one geometry approximation to the Piola transform is not yet known. For the original SCMS-MG preconditioner, have chosen $\beta$ and $\sigma$ as for the first model problem, just for the Raviart-Thomas smoother for the case with Piola transformation, we have chosen $\beta=0.0025$. For the rank-one corrected version, we have chosen $\beta=0.01$; the damping has been chosen based on an approximation for constants of the inverse inequality on the physical domain.

As in the case of standard finite elements, there are several possibilities to discretize the mixed formulation of the Stokes equations. Our experiments indicate that it might pay off to use the (in terms of degrees of freedom) more expensive variant of Taylor Hood discretizations than the other variants, particularly because it is known that this discretization also works for direct composition. The $p$-robust smoothers which we have proposed for for the Poisson problem can be carried over also to the Stokes flow problem, however it seems that further investigation is necessary concerning its application in the framework of non-trivial geometry transformations.