The node (or point) based ordering is a variant that is only possible with equal-order elements, like P1-P1-stab, which is used in PN S. Here, those variables that belong to the same degree of freedom (but different physical quantities) are grouped together. The new stiffness matrixFπ can be derived fromF by applying an according permutation matrix Π:
Fπx= ΠTFΠΠT u would have the form:
as in (7.34) or (7.37). This type of ordering is also referred to as Vanka-type of ordering. However, the original method in [Van86] consisted of an element-wise ordering, in contrast to the point-wise ordering used here.
Basically one has the possibility to solve the arising systems with iterative solvers such as Krylov methods, since the point based ordering obviously inhibits Schur complement methods. However, this type of ordering also has some advantages. For example, the data structure for the sparse matrixFπ doesn’t need to store the indices of all single scalar entries but only the indices1 of each 3×3-block
1see Section 2.5.2 in Chapter 2 about sparse storage formats
(4×4-block in 3D).
Numerically, a Krylov method would lead to more or less the same result for (8.2) as well as for (8.5) – except for the permutation.
However, since the size of the blocks is known in advance, the matrix-vector multiplication can be better optimized and vectorized by the compiler (cf. Section 3.8). This has direct consequences for most Krylov subspace methods, which are strongly based on this operation. Thus, one and the same Krylov method for 8.2 and 8.5 would exhibit a greater execution speed for the latter (see also [DLH00]).
Furthermore, we are able to use specialized block versions of numerical algorithms, as we show in the following.
8.2.1 Simple iteration and preconditioning methods for the point based ordering In order to increase the convergence speed of iterative methods, one is interested in a preconditioner M ∈RK×KforFπ (orF), such thatM−1Fπ ≈IorFπM−1 ≈I. Generally,M is some approximation to Fπ. Traditional preconditioners arise from splitting methods, which split the matrixFπ additively intoFπ =M−N whereM is easy to invert. Usual choices are
M =D (Jacobi),
M =ω−1D (damped Jacobi), M =L+D (Gauß-Seidel), M =ω−1(ωL+D) (SOR),
where L is the strict lower triagonal part L, D is the diagonal part and ω ∈ (0,2) some relaxation parameter. Because of their slow convergence rate, they are of course outdated as a pure solver when compared to e.g. Krylov methods, but still are important as a smoother for multigrid methods (cf.
chapter 5).
For saddle point problems however, these traditional iteration schemes above show comparatively poor convergence rates, regardless of the ordering of the unknowns. This changes when exploiting the block structure of the matrix.
Generalized SOR
The relaxed Gauß-Seidel or SOR (Successive Over-Relaxation) algorithm can be formulated as a generalized version for block matrices of the following form. Let A∈RN×N be partitioned inton·n submatrices such that
A= (Aij)ni,j=1, Aij ∈Rni×nj, n1, . . . , nn∈N, Xn
i=1
ni =N. (8.6)
Then we can formulate the generalized SOR method for block matrices like this:
Following the generic programming paradigm, Algorithm 5 is parameterized with invert, which shall depict an appropriate solution method, that determines or approximates A−1ii yi. It might be chosen different types of solvers for different diagonal entries, for example exact solvers like an LU -decomposition if Aii is relatively small, or some iterative solver again, if Aii is relatively large and sparse, leading to a hybrid method. For Aii ∈R it is the simple inversion of a real number, leading to the standard SOR2.
For the point based ordering withP1 elements for each unknown, we haven1 =. . .=nn= 3 for a 2D problem (4 in 3D), in any case a small constant number, such that a direct solver can be employed.
2Using SOR again for the inversion of the diagonal blocks would result in method that is mathematically equivalent to the standard SOR that simply is applied to the whole matrixA.
Algorithm 5: Generalized SOR
Input: MatrixA∈Rn×n, vector b∈Rn, initial approximation x0 ∈Rn Output: An approximation xk to A−1b
repeat
for i= 1 ton do
1
yi ←bi−Pi−1
j=1Aijxkj −Pi−1
j=1Aijxk−1j
2
zi ←invert(Aii,yi) xki ←ωizi(Ini−ωi)xk−1i end
until convergence
In BLANC and in our library MiLTON, we directly compute the inverses of allAii in advance, since it obviously can be reused during the iteration. We will refer to this method also asblock SOR.
The relaxation weights ωi∈Rni×ni are now matrices instead of scalars, however in BLANC, they are restricted to diagonal matrices. Numerical experiments show that it is important to be able to vary the relaxation parameters for different physical quantities, see [M¨ul01]. Experimental studies support that the choice ofωvel ≥1 (overrelaxation of the velocity part) andωp≪1 (underrelaxation of the pressure part), with the relaxation weight
ω :=
ωvel
ωvel ωp
(8.7)
is beneficial for the 2D Oseen equation.
Backward SOR, SSOR
In Algorithm 5, the components ofxk are updated fromxk1 to xkn. Therefore, the inner loop starting at line 5 it is also called a forward SOR step. If the directions of the outer loop in line 5 and of the summation in line 5 are reversed, such that the components of xk are updated from xkn to xk1 it is called abackward SOR step.
Combing these two alternating steps together leads to the SSOR method, which exhibits greater stability for convection dominated problems in practical computations. It is used as a smoother for our algebraic multigrid methods throughout this thesis.
Block ILU
Another class of methods that are widely used as preconditioners (and sometimes as multigrid smoothers) are the incomplete factorizations, the most famous being the ILU (Incomplete LU) decomposition.
For matrices like (8.6) ILU and its variants can be formulated as block algorithms (often referred to as BILU) ifn1 =. . .=nn. Since the inverses of the diagonal entries have to be computed explicitly, the dimension of the submatrices cannot be arbitrarily high.
However, simple ILU(0) preconditioners often fail, especially when applied to indefinite systems as described in [Saa92]. Our own experience with ILU supports this observation, even when an ILU decomposition is feasible, it often shows a much worse convergence behaviour than e.g. an SSOR preconditioner. Variants of ILU, that allow fill-in (ILU(p), or that employ pivoting (ILUP), drop-tolerance (ILUT), or that add entries to the diagonal instead of dropping them (MILU), however seem to promise better convergence, see e.g. [CS97]. A short but comprehensive overview of these methods is given in [CvdV94].
Algebraic multigrid for systems of equations
The methods introduced in Chapter 5 are well suited for scalar elliptic problems like the Poisson equation or the general (stabilized) advection diffusion reaction equation. However they become rather ineffective when applied to a matrix arising from a system of equations like a coupled Poisson, an Oseen or Navier-Stokes problem in two or three dimensions. Here, the classical AMG algorithm could find connectivity between physical quantities where there is none. Our own numerical computations support this – the standard AMG method from Section 5.2 applied to an Oseen problem may still converge (in simple cases) but the convergence rate heavily slows down1. A smoother may still work as a solver, but then the smooth error components (consisting of eigenvectors of Sl that belong to large eigenvalues) are not well enough approximated on the next coarser level.
For saddle point problems like Oseen or Stokes, there is another difficulty: the incompressibility constraint (7.2) causes the indefiniteness of the system matrix for most of the standard discretizations.
Only with according pressure stabilization schemes, a positive definite matrix can be gained.
So there is obviously a demand for adapted algebraic multigrid multigrid methods for saddle point problems. One of the first adaptions of algebraic multigrid ideas for the Navier-Stokes equations (for a finite volume discretization) was presented in [Raw95], where a Additive Correction Multigrid was used to construct the prologation and restriction operators. In [Web01] the smoothed and unsmoothed aggregation approaches where applied.
Furthermore there also exist hybrid methods, based on Schur complement iteration schemes, which use AMG only for inverting the Schur complement or the velocity part (see e.g. [EHST02] or [GNR98]).
In the following two sections, we present customized, fully coupled, algebraic multigrid methods for mixed problems. These methods basically differ in the ordering of the unknowns, that were introduced in the previous chapter. All of these methods however, are mainly based on the classical AMG approach presented by Ruge and St¨uben in [RS87].
9.1 Algebraic multigrid for the unknown based ordering
An algebraic multigrid method for systems of equations that are assembled into a matrix using un-known based ordering as in (8.2) was already suggested in [RS87]. At least, it is recommended to prolongate and restrict the physical quantities separately, using a block-diagonal prolongation operator like
Pvar =
P1 0
. ..
0 Pn
, (9.1)
1This deterioration of convergence was e.g. also reported in [Oel01] for coupled Laplace equations.
139
withn= 3 for the two-dimensional Oseen problem (n= 4 in 3D). The prolongation matricesP1, . . . , P3
are constructed like in the scalar case out of Axx, Ayy, Cpp.
Problems then arise in the detail. Since equal-order elements are not LBB-stable, mixed problems are often discretized usingPk−Pk−1 (Taylor-Hood) elements. Different polynomial degrees however lead to different dimensions for Cpp and the rest of the matrices on the diagonal. This leads to problems when no pressure stabilization matrix is assembled (how to construct Pn in this case?).
Even if Cpp 6= 0, and Pn could be constructed out of it, it is not ensured, that Pn coarsens away (geometrical) points of the underlying mesh, that P1, . . . , Pn−1 don’t. This could lead to a situation, where the velocity variables on the coarse levels would live on totally different points than the pressure variables.
Another issue is the correct choice of the smoother. Standard SOR or SSOR smoothing has little effect on these equations. One suggestion is therefore to employ one step of a Schur complement method, such as the Braess-Sarazin smoother [BS97].
Various (successful) suggestions to overcome these problems are explored in [Wab03] for the Navier-Stokes, respectively the Oseen and Stokes equation. The method forP2-P1 elements introduced there however uses additional geometrical information to construct the coarser levels, therefore it is not purely algebraic.