For absent reaction (c ≡0) and decreasing diffusion (ν → 0) the problem (4.2) becomes convection dominated and looses its elliptic property, resulting in a nearly singular discrete problem, since the positive definiteness is solely grounded on ν. Problems with iterative solvers are the consequence, because their convergence rates directly depend on the condition of the matrix. Furthermore, the resulting discrete solutions then exhibit strong oscillations. Therefore, for practical applications, one needs to apply a stabilization method, that strengthens the ellipticity of the bilinear form.
Let in the following (4.3) - (4.4) be again discretized using a finite element spaceVh ⊂V, with an according triangulation Th, leading to a discrete problem as in (4.6).
The Streamline Diffusion/SUPG stabilization method (originally introduced in [HB79]), for exam-ple, adds artificial diffusion only in the direction of the convection b.
Definition 4.3.1 (SUPG method). The streamline upwind Petrov Galerkin (SUPG) or streamline diffusion method is defined as finding uh ∈Vh, such that
aSU P G(uh, vh) =fSU P G(vh), ∀ vh∈Vh, with
aSU P G(uh, vh) :=a(uh, vh) + XM
i=1
δTi(−ν∆uh+b· ∇uh+cuh, b· ∇vh)Ti, (4.11)
fSU P G(vh) :=f(vh) + XM
i=1
(f, b· ∇vh)Ti (4.12)
Here (·,·)Ti is the localL2-scalar product on the triangle/tetrahedron Ti: (u, v)Ti := (u, v)L2(Ti)=
Z
Ti
uv dx.
Especially for piece linear elements, where we have vanishing second order derivatives, and for c= 0 we see, that the SUPG method adds a small diffusive part to the variational formulation.
In the Galerkin least-squares method (see [HFH89]), the stabilization term is tested not only with the convection part but with the whole differential operator.
Definition 4.3.2 (GLS method). The Galerkin least squares (GLS) method is defined as finding uh∈Vh, such that
aGLS(uh, vh) =fGLS(vh), ∀ vh ∈Vh, with
aGLS(uh, vh) :=a(uh, vh) + XM
i=1
δTi(−ν∆uh+b· ∇uh+cuh,−ν∆vh+b· ∇vh+cvh)Ti, (4.13)
fGLS(vh) :=f(vh) + XM
i=1
(f,−ν∆vh+b· ∇vh+cvh)Ti (4.14)
Both of the above methods enhance the diffusive part of the bilinear form and therefore strengthen the positive definiteness of the stiffness matrix A. In PN S the stabilization parameter for both methods is chosen as
δTi := h2Ti
2ν(1 +P e2Ti)−12, (4.15)
wherehTi := diam(Ti) and
P eTi :=kbk∞,Ti
hTi
ν (4.16)
is the localPeclet number on elementTi. The Peclet number can be seen as a measure for the relation of convection and diffusion. A possible definition for the whole equation could be the following:
Definition 4.3.3. Problem (4.1) is called
purely convection dominated if min
i=1,...,MP eTi >1 (4.17)
purely diffusion dominated if max
i=1,...,MP eTi ≤1. (4.18)
Of course, mixtures of these types can occur, sinceP eTi is only a local measure, that can vary over the domain. For stability proofs and error estimations of the above stabilization methods, we refer to [RST96].
Algebraic multilevel methods
Since algebraic multigrid is heavily inspired by the concept of geometric multigrid and uses many terms thereof, it is hardly possible to describe the idea behind the former without first presenting the idea of the latter. But we will keep things as short as possible.
This chapter will give a brief overview first of geometric multigrid and then introduce the approach of algebraic multigrid.
5.1 Geometric multigrid
The foundations of (geometric) multilevel or multigrid methods were laid by Fedorenko1 [Fed61] and Bakhvalov [Bak66] in the 1960’s and further developed by Brandt (e.g. in [Bra73]) in the 1970’s.
It was only in the 1980’s when they became more popular since increasing CPU speed and larger memories then allowed to tackle bigger problems.
The motivation of geometric multigrid arose from the observation that simple solvers – splitting methods like Jacobi or Gauß-Seidel (see Section 8.2.1) – have a smoothing effect on the error of the approximate solution.
5.1.1 Smoothing
Let x =A−1b denote the exact solution and x(k) the approximation at step k. Remember that the according iteration reads
x(k):= (I−M−1A)xk+M−1b,
and thus, if e(k):=x(k)−x is the error in thek-th step we have for step k+ 1:
e(k+1) =x(k+1)−x=xk−x−M−1(Axk−b) = (I−M−1A)ek.
So the error components are damped with the iteration matrix I −M−1A. Looking, for example, at a simple boundary value problem like the 2-dimensional Poisson equation
−∆u=f, in Ω := (0,1)×(0,1) (5.1)
u= 0, on∂Ω, (5.2)
1Radii Petrovich Fedorenko, a Russian mathematician at the Keldysh Institute for Applied Mathematics, is always cited as the inventor of multi-level methods. Interestingly however in an article by Saad [Sv00], it is mentioned that the British engineer and mathematician Richard Vynne Southwell (1888–1970) in 1935 [Sou35] already had the the basic idea of a two-level method.
77
discretized using finite differences on a structured (N×N) grid produces a linear systemAx=bwith a symmetric positive definite stiffness matrix
A= 1
Solving this linear system with the (damped) Jacobi2 method using the relaxation parameter ω= 12 yields the iteration
In (5.5) the new component of the error vector is formed by the (weighted) average over the neighbour-ing points – this means that small oscillations are reduced quickly. The local information transport is very quick, whereas the global information transport (ifN is very large) is quite slow – which is why this method is a good smoother, but a poor solver.
Remark 5.1.1. Another approach is to look at the eigenvectors of the iteration matrix I −M−1A.
They are the same as those of A which in turn correspond to the discrete points of the eigenfunctions of the boundary value problem (5.1) – (5.2). Now the low and middle-frequent eigenfunctions (and thus the eigenvectors) belong to eigenvalues of I−M−1A which are close to 1 in contrast to the high frequent eigenvectors belonging to eigenvalues which are close to zero.
As a consequence of this, we only need to find a way to smooth the lower frequencies of the error as well, to have a good solver.
2Originally, the German mathematician Carl Gustav Jacob Javobi (1804–1851) in [Jac45] invented this iterative method, because the exact Gaussian elimination method was tedious and error-prone for large linear systems.
5.1.2 Coarse grid correction
Since partial differential equations are posed on a geometric domain Ω they are usually discretized using a grid of a certain mesh width, say h. Now it lies near to approximate the lower frequencies of the error with a dicretization on a coarser level, for instance on a grid of width 2h. On that level the same type of smoothing as above can be applied, and its low frequencies in turn can be smoothed on coarser grids recursively, which gives the classic multigrid method.
Precisely, we need ahierarchy of grids
Ωlmax ⊂Ωlmax−1 ⊂ · · ·Ω1⊂Ω0⊂Ω,
where Ωi ⊂ Ωj means that Ωj includes all points of Ωi. On each level l the problem is discretized using the finite element space Vl which consists of the according polynomial elements defined on grid Ωl, so we get a sequence of linear systems
Alxl=bl, with Al∈Rnl×nl, nl = dimVl (5.6) wherebl ∈Rnl is the dicrete version of the right hand sidef and xl∈Rnl the discrete approximation of u on level l. Especially we set A=A0, b=b0 and x=x0.
Since fine grids are often generated through the refinement of a simpler coarse grid, such hierarchies may be produced at no extra costs, if some rules are heeded.
Between the grids we need of course appropriate transfer operators for the interpolation.
Definition 5.1.2 (Restriction, Prolongation). The restrictionfrom a finer level l to the coarser level l+ 1 is denoted by
Rl :Vl −→Vl+1, while the reverse mapping from level l+ 1 to levell
Pl :Vl+1 −→Vl,
which has to be some suitable interpolation, is called the prolongation.
Usually we have Rl ∈Rnl+1×nl, Pl ∈Rnl×nl+1 and Rl =c(Pl)T with some constant c. Generally the prolongation is chosen such that a pointp∈Ωl is interpolated through its neighbouring points which lie in Ωl+1.
Remark 5.1.3. When the domain Ω has a simple geometry like a rectangle and is discretized using structured grids with regular refinement, all works very nicely. But if the domain is more complex and if we have unstructured grids generated by an adaptive refinement, then the construction of proper transfer operators is anything but trivial.
5.1.3 The algorithm
Let Sl be the smoother on level l and the sequence of linear equations as in (5.6) be given. Then, after kiterations withSl, the error is e(k)l =x(k)l −xl ⇐⇒xl=x(k)l −e(k)l and thus e(k)l would be the optimal correction for our current approximation x(k)l . Since now for the residual rl(k) the following relation holds
r(k)l =Alx(k)l −bl =Ale(k)l (5.7) we approximate e(k)l now by restrictingrl(k) on the coarser levell+ 1:
rl+1 :=Rlrl(k), (5.8)
where we solve the equation
Al+1el+1=rl+1. (5.9)
The resultel+1 is prolongated again to levell and then used to update our current solutionx(k)l : x(k,new)l :=x(k)l −Plel+1. (5.10) As this simple coarse-grid correction does not automatically guarantee convergence, additional smooth-ing is required after the correction (post-smoothsmooth-ing) or before (pre-smoothsmooth-ing).
Definition 5.1.4 (Two-grid-cycle, Multigrid-cycle). The process described above is called two-grid-cycleif only two levels are involved and the system on the coarsest level is solved exactly. If this process is again used recursively to solve (5.9), it is called a multigrid-cycle.
Finally we can summarize these ideas in the following general multigrid algorithm. The parameters ν1, ν2 control the number of pre- and post-smoothing steps, whileγldetermines the recursion pattern of the cycle:
γl≡1 : V-cycle, γl≡2 : W-cycle, γl=l : F-cycle.
Furthermore, let Smootherl(Al, xl, bl) be a function that applies the smoother Sl once to the System Alxl=bl.
Algorithm 1: MGCyclel(Al,xl,bl)
Input: MatrixAl, vector bl, smoother Sl, parameters ν1, ν2, γl, initial approximation xl
Output: An approximation xl to A−1l bl if l=lmaxthen
xl←A−1l bl else
for i= 1 toν1 do
xl← Smootherl(Al,xl,bl) end
rl+1 ←Rl+1(Alxl−bl) el+1 ←0
1
for i= 1 toγl do
el+1← MGCyclel+1(Al+1,el+1,rl+1) end
xl←xl−Plel+1 for i= 1 toν2 do
xl← Smootherl(Al,xl,bl) end
end
Remark 5.1.5. Algorithm 1 is the basic algorithm for all multigrid methods – even algebraic multigrid methods rely on this procedure since it is the core principle of every multigrid idea.
Remark 5.1.6. The runtime complexity of one multigrid-cycle is proportional to the total number of unknowns on the finest mesh, n0, if
lmax
X
l=0
nl≤Cnl0, (5.11)
with a constant C independent of n0. This is because the number of nonzeros in A is O(n0) when A comes from some discretization (FEM, FDM, FVM) of a PDE. The above equation is a reasonable condition, since the meshsize on a level l+ 1should be significantly smaller than on level l.
Remark 5.1.7. To get a sensible solving algorithm, we apply several multigrid-cycles, x(k):=MGCycle0(A, x(k−1), b)
untilx(k) satisfies some stopping criterion. In order to improve the initial valuex(o), we might use the so-called nested iterations: The system is solved directly on the coarsest level, i.e. xlmax :=A−1lmaxblmax, then prolongated to the next level (l=lmax−1), where some smoothing steps are applied, then again prolongated, and so on, until we get some x(o) on the finest level, where the normal multigrid solving algorithm can be started as above.
5.1.4 Convergence results
In order to show the convergence of a multigrid method, one needs to ensure that the smoother Sl really reduces the high-frequent error and that the error can be approximated on the coarser level as decribed through equations (5.8) – (5.10). The standard convergence theory for geometric multigrid is based on the following two properties.
Definition 5.1.8 (Smoothing property). The condition
kAlSlν1k ≤CSη(ν1)h−α, (5.12) with some function η withη(ν1)→0 for ν1→ ∞ and α >0 is called smoothing property.
Definition 5.1.9 (Approximation property). The condition
kA−1l −PlA−1l−1Rl+1k ≤CAh−α (5.13) is called approximation property.
For symmetric scalar problems there are already convergence results available as e.g. in [Hac85]
or [Man88]. Fewer results are available for nonsymmetric scalar problems like the convection-diffusion equation. For this type of problem Reusken was one of the first to show convergence for the the two-grid method and the W-cycle with damped Jacobi and Gauss-Seidel smoothing for finite differences [Reu00]. A similar result for a finite element discretization is given in [OR03].
Derived originally for elliptic problems, the (geometric) multigrid method delivers good results as well for other types of problems. For many cases, it only needs a constant number of multigrid-cycles to converge, independently of the meshsizeh, which results in an overall runtime complexity ofO(n0).
5.1.5 Disadvantages
Geometric multigrid are however not the optimal method in every case. For example, they show a lack of performance if the data of the underlying boundary value problem is not smooth enough. This is mainly due to the construction of the prolongation, which simply interpolates between points, and does not take into account the special properties of the (discretized) differential operator which forms the stiffness matrixA. As a first remedy, operator-dependend transfer operators have been developed by e.g. Alcouffe et al. [ABDP81] and de Zeeuw [dZ90].
In addition, as mentioned above, multigrid methods are not easy to implement independently of the problem. The interpolation is a crucial point, especially if we deal with complex geometries and unstructured grids. Adaptive mesh refinement with green refinement e.g. yields meshes that are not hierarchical and so are not applicable to the above approach.
Generally we need geometrical information about the underlying problem which can be missing
• because our software doesn’t allow access to the geometry information for some reason, maybe it’s because we have a legacy system whose mesh-generator is not capable of providing us with hierarchical grids or it is a third party product for which we don’t even have the source code
• or because the linear system does not stem from a discretized PDE at all, so we have no finite element spaces and no mesh information - e.g. optimization problems.
Speaking now in terms of computer science for geometric multigrid, the interface between the solver and the FEM-software doesn’t only consist of the linear system Ax=b, but also contains the meshes Ω0, . . . ,Ωlmax. Thus, a geometrical multigrid method isno black box solver.