6.2 Convection diffusion problems
6.2.2 A 3D example
We would like to test the AMG performance also in a three dimensional setting. We examine a 3D rotation flow, which is described by the following problem:
Find u(x, y, z) :R3→R, such that the following mixed boundary conditions:
~n· ∇u= 0 for x= 0 (Neumann)
This problem is discretized on a structured 3D grid (tetrahedrons) using linear elements und GLS stabilization.
For all problems with ν ∈ {1,10−2,10−4,10−6}, the stabilization parameter δTi was chosen as in formula (4.15) in Chapter 4.
The vector fieldb(see Figure 6.14) induces a circular convection around thez-axis. This convection drags the solution profile of the inflow boundary (Figure 6.13) at y = 0 through the whole domain and transports it to the outflow boundary at x= 0. Slice plots of the three-dimensional solution are shown in the next four subsections.
0
Figure 6.13: Solution profile at the inflow border
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 6.14: Convection field viewed from above (z≡const)
h 1/4 1/8 1/16 1/32 1/64 1/129 inner points (= matrix dimension) 36 392 3600 30752 254016 2064512
nonzero entries 290 4174 42662 382294 3230390 26547574 Table 6.8: mesh widths and matrix properties for the 3D problem
In Table 6.8 we see the properties of the stiffness matrices resulting from aP1 discretization on a structured mesh inPN S.
0
Figure 6.15: Solution profile at the outflow bor-der, ν = 1
Figure 6.17: Solution profile at the outflow bor-der, ν = 10−2
The profile at the outflow boundary is determined by the amount of diffusion controlled by the vis-cosity parameterν. For ν= 1 the diffusion dominates over the convection, the information stemming from the inflow boundary is transported slowly in every space dimension (nearly) uniformly (Figure
6.16). The convection is comparably weak, such that the solution at the outflow boundary is nearly zero (Figure 6.15). At ν = 10−2, more of the input profile is transported into the direction of the outflow boundary, however blurred by the diffusion (Figures 6.17 and 6.18).
0
Figure 6.19: Solution profile at the outflow bor-der, ν= 10−4
Figure 6.21: Solution profile at the outflow bor-der, ν= 10−6
With diminishing diffusion, at ν = 10−4 and ν = 10−6, we see that the sharp input profile is transported through the whole domain (Figures 6.20 and 6.22) to the output boundary (Figures 6.19 and 6.21) nearly without any losses.
The effect of the strong coupling parameter θ on the coarsening process
Again we would like to look on the way the strong coupling determinates the generated levels. The coarsening parameter θ is again chosen out of the set Θ from (6.4). The according setup times are shown in Tables 6.9 to 6.12. Using the legend from Figure 6.5 we display the data of the generated levels in the Figures 6.23 to 6.26 for the four values ofν, exemplarily for mesh widthh = 1/65. The upper curves again show the nonzeros of the matrices, the lower curves indicate their dimension.
As we can see in Tables 6.9 to 6.12, the time factor from a mesh sizeh to h2 is at best about 9 to 10. This is somewhat bigger than the factor 8 that the problem sizes differ (cf. Table 6.8). However, since the number of generated levels is in all cases restricted to five, the SuperLU inverter has to treat a relatively large matrix on the coarsest level.
θ h= 18 h = 161 h= 321 h= 641
0.2 0.02 0.31 3.31 34.22
0.3 0.02 0.21 2.09 18.82
0.4 0.02 0.21 2.08 18.89
0.5 0.02 0.23 2.18 19.36
0.6 0.02 0.24 2.29 20.44
0.7 0.02 0.22 2.36 23.4
0.8 0.02 0.24 2.26 21.86
Table 6.9: AMG setup times [s],ν = 1
θ h= 18 h= 161 h= 321 h= 641
0.2 0.02 0.17 1.68 17.1
0.3 0.02 0.17 1.7 15.53
0.4 0.01 0.18 1.55 13.43
0.5 0.02 0.16 1.34 13.92
0.6 0.01 0.14 1.63 17.5
0.7 0.02 0.15 1.65 19.58
0.8 0.01 0.14 1.5 21.49
Table 6.10: AMG setup times [s], ν = 10−2 θ h= 18 h = 161 h= 321 h= 641
0.2 0.02 0.23 2.53 23.39
0.3 0.02 0.21 2.14 20.86
0.4 0.01 0.2 2.02 19.74
0.5 0.02 0.18 1.75 16.76
0.6 0.01 0.17 1.62 16.09
0.7 0.02 0.14 1.33 13.59
0.8 0.01 0.13 1.24 12.69
Table 6.11: AMG setup times [s], ν= 10−4
θ h= 18 h= 161 h= 321 h= 641
0.2 0.02 0.23 2.71 25.11
0.3 0.02 0.21 2.1 20.32
0.4 0.02 0.21 2.1 20.14
0.5 0.02 0.17 1.74 16.71
0.6 0.02 0.17 1.6 16.12
0.7 0.01 0.14 1.39 13.32
0.8 0.02 0.12 1.22 11.74
Table 6.12: AMG setup times [s], ν = 10−6 For the diffusion dominated problems, we see that values between 0.3 and 0.5 for θ are optimal because they not only reduce the C variables, but also the nonzero entries. This is the same effect as mentioned in the 2D context, at first the levels doesn’t differ much, lower values for θ introduce more nonzeros in the coarse level matrix, but this increased neighbourhood is likely to be coarsened drastically in the next step.
For convection dominated problems, values between 0.5 and 0.8 seem to be more appropriate, since they result in a stronger reduction of the nonzero entries. Again, as in the two-dimensional setting, bigger values for θ decrease the elements in the set of interpolation variables Pi, therefore yielding prolongation and coarse level matrices with a smaller number of nonzero entries. This is also supported by the setup times listed in the Tables below. However, this effect is not so evident as for the 2D problem.
In the three dimensional setting, as for the two dimension case, we observe, that the number of variables from level to level is about halved in each step, independent of the space dimension. Since it takes into account also the direction of the stream, this type of coarsening is also referred to as
levels
dimensions,nonzeros
0 1 2 3 4 5
102 103 104 105 106
Figure 6.23: Level hierarchy forν= 1
levels
dimensions,nonzeros
0 1 2 3 4 5
102 103 104 105 106
Figure 6.24: Level hierarchy for ν = 10−2
levels
dimensions,nonzeros
0 1 2 3 4 5
102 103 104 105 106
Figure 6.25: Level hierarchy forν = 10−4
levels
dimensions,nonzeros
0 1 2 3 4 5
102 103 104 105 106
Figure 6.26: Level hierarchy for ν = 10−6 directional coarsening or semi-coarsening. In the AMG context, this is simply a byproduct of the strong neighbourhood relationship and the C/F-splitting.
For geometric multigrid, the factor from a fine level to a coarse level dimension is normally 1/4 in 2D and 1/8 in 3D, which results in a lesser computational complexity. However, for anisotropic grids or strong convection, the semi-coarsening approach often has superior convergence properties (cf. [Mul89]).
Nevertheless, in 3D, semi-coarsening has the drawback, that the overall complexity is not reduced strong enough, the AMG setup phase takes significantly longer than in the 2D case. A remedy can be to apply agressive coarsening as introduced in [St¨u99], where not only direct connections, but also indirect connections (i.e. neighbours of neighbours) are considered.
The effect of relaxation on the solution process
For the 3D problem, we seek the optimal relaxation parameter(s). Again, ω is taken from (6.6). As we see in the Tables 6.13 and 6.14, the algorithm converges over a wide range of values for diffusion dominated problem.
However, the smaller ν is, the more important becomes the underrelaxation of the smoother. In fact, many AMG level hierarchies (for differentθ) didn’t yield convergence in the convection dominated case for ω≥1. This sustains the fact that underrelaxation is inevitable for these types of problem.
ν = 1 ν = 10−2
ν = 10−4 ν= 10−6
Convergence speed of the AMG method
Again, we compare the AMG method with a GMRES(m) Krylov solver, using a restart length of m= 20, in combination with a SSOR preconditioner (with the optimal relaxation parameterω= 0.7).
time [s]
relativeresidualnorm
0 5 10 15
10-10 10-8 10-6 10-4 10-2 100
GMRES
AMG
Figure 6.27: Convergence of GMRES vs. AMG forh= 641
time [s]
relativeresidualnorm
0 25 50 75 100 125 150 175
10-10 10-8 10-6 10-4 10-2 100
GMRES
AMG
Figure 6.28: Convergence of GMRES vs. AMG forh= 1281
In Figures 6.27 and 6.28 we exemplarily see the convergence of the residual in the L2-norm (6.7) for the strongly convection-dominated problem with ν = 10−6 at h = 641 and h = 1281 , the two finest mesh-widths, that were possible in 3D on a single workstation with 4GB of memory.
In order to shorten the setup time for the AMG solver, we increase the number of levels generated to 7 forh= 641 and 8 forh= 1281 . Since we have hardly any setup cost for the GMRES solver and the SSOR preconditioner, the Krylov solver has a headstart (red curves). The AMG setup time is much higher, roughly 5 sec. forh =641 and 38 sec. at h= 1281 .
variables
variablespersecond
102 103 104 105
1000 3000 5000 7000 9000 11000 13000
nu = 1 nu = 0.01 nu = 0.0001 nu = 0.000001
Figure 6.29: Calculation speed for the 3D problem
We can see, that the GMRES solver is clearly superior on the coarser mesh-width. At this problem size, the stiffness matrix’ condition number is small enough to allow the Krylov solver to converge quickly. Also, the convection field and the boundary conditions are comparatively simple, and thus doesn’t cause the problems with boundary layers as in the 2D example.
Only at the finer mesh-width, the AMG solver shows a more typical multigrid behaviour, it can compete with (and even outperform) the Krylov solver. Although the convergence curves doesn’t differ much, we see that the scaling of AMG between the problem sizes is very much better (150 sec.
to 16 sec. versus 160 sec to 4 sec.).
The diagram in Figure 6.29 finally shows the calculation speeds on the different problem sizes.
Here the optimal parameters from the previous empirical studiess were used. In contrast to the 2D case, however, we observe a degradation of the speed for small h. The computational effort due to the increased densitity of the matrices slows down the setup (cf. Proposition 5.3.9) as well as the multigrid iteration.
AMG for mixed problems
123
Mixed problems in computational fluid dynamics
The physical discipline of fluid dynamics is the study of fluids (that is liquids and gases) in motion.
Its fundamental mathematical model, that describes continuous flows, is given by the Navier-Stokes equations1. Since a general solution can’t be derived in closed form, one is dependent upon numerical approximations that the scientific area of computational fluid dynamics (CFD) is trying to give.
In this Chapter, we shortly describe the (incompressible) Navier-Stokes problem, its linearization, the finite element discretization and the Oseen equation which emerges in this connection. Finally, we introduce the stabilization methods that needs to be applied for a successful computation.
7.1 The Navier-Stokes equations
We start with the definition of the main CFD problem, which is the background and reason for our research on numerical methods.
Definition 7.1.1 (Navier-Stokes problem). One of the most important and difficult problems in com-putational fluid dynamics consists of the evolutionary incompressible Navier-Stokes equations, which describes the flow and the pressure of a fluid in motion. More exactly, we are interested in the distri-bution of the pressurep(t, x) and the velocity u(t, x) = (u1, . . . , ud)T(t, x) in the time-space cylinder QT := (0, T)×Ω with respect to a given source-term f(t, x) = (f1, . . . , fd)T(t, x):
∂u
∂t −ν∆u+ (u· ∇)u+∇p=f in QT, (7.1)
∇ ·u= 0 in QT, (7.2)
with inital and boundary conditions
u(0, x) =u0(x), (7.3)
u= 0 on (0, T)×∂Ω (7.4)
where 0 < ν ≤ 1, (t, x) ∈ QT, and Ω ⊂ Rd, d ∈ {1,2,3}. The scalar ν > 0 is the viscosity parameter while equation (7.2) is the incompressibility constraint.
A main obstacle in solving these equations is the nonlinearity, which prevents us from giving an analytical solution in the classical sense for arbitrary space dimension d. Leray ([Ler34]) in 1934 showed, that weak solutions exist under certain conditions. However, even if a weak solution exists,
1Named after Claude-Louis Navier (1785–1836), French engineer and physicist, and George Gabriel Stokes (1819–
1903), Irish mathematician and physicist.
125
its uniqueness can not be guaranteed, especially in the three-dimensional case and for arbitrary large time-intervals (0, T). Therefore, one is interested in computing a discrete numerical approximation in order to be able to simulate practical industrial applications.
This is usually achieved by the discretization in time and space, leading to systems of linear equa-tions with a relatively high dimension number. Solving these linear systems, most Navier-Stokes solvers spend the biggest amount of time during the whole computing process. Here, efficient algorithms and a thorough implementation have to be applied.
Consequently, each component of the CFD software has to be chosen with care. An important topic, beside the mesh generation, is the discretization in time and space.
• The spatial (semi-) discretization on the domain Ω is usually done by using finite element (or finite difference or finite volume) methods.
• For the discretization of the time interval (0, T) we normally use a time-stepping method like a single-step θ-scheme.
We will apply the time discretization before the space discretization, resulting in an outer time loop, which seems to be a quite natural approach to an instationary problem, however other methods are also common. In each time step a stationary Navier-Stokes type of equation in the form of a (discrete) nonlinear and indefinite saddle point problem has to be solved.
A closer look at the corresponding discretization methods will be taken in the Subsections 7.1.1 and 7.1.2.
Now to solve the arising nonlinear problems, we have basically two alternatives:
• First treat the nonlinearity by a linearization technique for the whole problem like Newton or a fixed point iteration which leads to linear subproblems. These can be solved by using a fully coupled or a decoupled approach.
• Another possibility is to first split the problem into equations for the velocity and the pres-sure components, which then in turn can be treated with a linearization method for nonlinear problems.
We will focus on the first approach, using a fully coupled AMG solver for the linear problems (see Chapter 9). Decoupling methods, such as Schur complement methods, will shortly be discussed in Chapter 8.
In any case, our CFD-software is forced to calculate a solution to a linear system with a large-dimensional sparse matrix in the end. Mainly two variants of solver algorithms are eligible for this task:
• Krylov subspace methods like GMRES, BiCGSTAB, etc. are usually reliable methods. Though they need an appropriate preconditioning.
• Multigrid methods, especially when adapted to the concrete problem, promise even better per-formance.
Geometric multigrid however, is not trivial to implement, especially the interpolation for abitrary unstructured, anisotropic meshes. Also, for decreasingν these methods often show degraded conver-gence rates. Instead, the focus of this work will lie on the implementation and application of algebraic multigrid methods for saddle point problems, to which Chapter 9 is dedicated.