4.1 Velocity Systems
As mentioned in Section 2.4.2, the linear systems for the velocity in the segregated approach are highly diagonally dominant. Therefore we do not apply any specialized AMG strategies for those systems, but stick with a simple BiCGStab2 scheme, that has been proven to be sufficiently fast and stable, see Section 4.5.1. This section will also give some examples regarding the diagonal dominance of those systems and the computational effort to solve them by using BiCGStab2 versus using a sophisticated AMG method. It turns out that in most cases the overhead that comes with using an AMG method does not outweigh the lower iteration count.
4.2.2 Coarsening Strategies
Regarding the coarsening strategy, we recall from Section 2.4.1 that compared to standard discretizations of Poisson’s equation the FPM matrices A are fairly dense.
While with a standard Finite Difference method in 3D the Laplace stencil would include 7 entries, in the FPM we have 40 entries per row. When constructing the coarse level Galerkin operator
AH =RAP, (4.3)
this increased density becomes even more exaggerated, cf. Section 5.4.1. This means that
• the coarse level operatorAH consumes more memory than in comparable, mesh-based systems and
• since AH has a lot of entries, the smoothing steps on the coarse level are more expensive, again compared to mesh-based systems.
To counter this effect, we coarsen more aggressively compared to what would be the standard AMG approach [138]: Where in the basic AMG algorithm, for a variable that we choose to be a coarse level variable (C-variable), all its strongly connected, di-rect neighbors become fine level variable (F-variable), in the so-called (l, k)-aggressive coarsening strategy this definition is extended to all variables that are either strongly connected direct neighbors or to which at leastlpaths of lengthkalong strong connec-tions exist. In our case, we used a (1,2)-aggressive coarsening strategy, meaning that for every C-point iC, all its strongly connected neighbors jk as well as the points that are strongly connected to thosejkbecome F-points. This leads to a lot less C-points on the coarse level and therefore to an overall smaller number of non-zero entries. When using an aggressive coarsening strategy, we need to employ the so-called multi-pass in-terpolation [138] which proceeds in several passes, using the direct inin-terpolation [138]
where possible and then use interpolation formulas at neighboring F-variables for those F-variables that are not connected to any C-variable. For the aggressive coarsening strategy we are using here, this process will need at most four passes before every F-variable has been assigned an interpolation formula [138]. We use this aggressive coarsening strategy to create the second level in our hierarchy. The other levels – if needed – are created in the standard Ruge-Stüben fashion. Comparisons between the standard and the aggressive AMG coarsening can be found in Section 5.4.1.
Positive couplings
AMG methods have originally been developed with symmetric, positive definite M-matrices in mind. Stüben [138] shows how for these M-matrices the error of a given approximation, after applying some iterations of a relaxation scheme, varies smoothly along large negative couplings. In an M-matrix, all off-diagonal entries, i.e. couplings, are negative. This does not hold true for the GFDM matrices. But because of the measures being taken to make the matrix weakly diagonally dominant where possible, cf. Section 2.4.1, most of the couplings are in fact negative, cf. Section 4.5.4. In
addition to that, since there are 40 neighbors in almost every neighborhood1, we expect to have at least one large negative coupling in every spatial direction, avoiding possible one-sided interpolations. Therefore, in our AMG algorithm we regard positive couplings as weak couplings, i.e. we neither coarsen nor interpolate along them.
4.2.3 Interpolation and Restriction
We employ a standard interpolation as described in [138] and let the restriction be the transpose of that interpolation. The idea to use the transpose of the interpolation as the restriction originates from the case of symmetric matrices but has shown to work sufficiently well in our case, too. We did conduct experiments with the modified lAIR restriction operator described by Southworth [135] but did not notice any sig-nificant improvements, possibly because Southworth’s work focuses on naturally non-symmetric problems like convection dominated phenomena. Another reason is that with our preprocessing of the matrix and the AMG techniques we use the convergence rates of our new methods are already fairly competitive (as the numerical experiments show) and could not be improved much, if at all, by the modified restriction operator.
4.2.4 Smoothing Schemes
Although the pressure matrices arising in the FPM do not have some of the properties that most FD matrices have, they are still weakly diagonally dominant for the most part, cf. Section 4.5.4. Also, the pressure systems do not have any special block structure with scaling issues like the saddle point problems in the coupled approach have. Hence we can apply a standard smoother. In our case this is the Gauss-Seidel relaxation scheme applied in forward order in the pre-smoothing step and in reverse order in the smoothing step. We perform one iteration in the pre- and post-smoothing step each.
4.2.5 Accelerator
As explained in Chapter 3, Multigrid methods are usually used as a preconditioning method for Krylov subspace methods. Then, in AMG terminology, the Krylov sub-space method being used is called an accelerator for the AMG method.
A standard accelerator choice for a Poisson-like problem discretized using Finite Dif-ference Methods would be the CG method. Since we are dealing with non-symmetric matrices in the case of a GFDM discretization, we use the BiCGStab method in-stead, cf. Section 3.2.2.
4.2.6 Coarse Level Solvers
In most situations, we employ a direct solver on the coarsest level ΩH of our multigrid hierarchy. The only exception to that is if the coarsest level is diagonally dominant
1Except for at the boundary.
with
maxi∈ΩHδi = maxi
P
i6=j|aij|
|aii| ≤0.9, (4.4)
in which case we use a simple one-level relaxation method, the default being Gauss-Seidel. Regarding the direct solver, we have two main options: Intel’s PARDISO and a standard LU-decomposition2. In the former case, the linear system on the coarsest level is collected on one single process, solved there and the solution is distributed back to all other processes. In the latter case, all processes involved in the AMG algorithm solve the coarse level system as well. Which of these methods gives a faster overall algorithm depends on the problem at hand, see Section 5.4.2. Both are direct methods though, so they do not influence the convergence of the overall AMG method. For the relaxation scheme in case of diagonally dominant coarse levels, we iterate until a relative residual reduction of 0.01 or a maximum number of 200 iterations is reached.
The latter case would be considered as an error though but does not occur when the coarse level matrix is diagonally dominant as it is the case here.
4.2.7 Setup Re-use Within a Time Step
Multigrid methods always consist of two separate phases: The setup phase and the solution phase. In the setup phase, the coarse levels are created, along with the corresponding interpolation and restriction operators. Since the coarse level system is also known after these steps, any decomposition or other form of setup for the coarse level solver can also be counted towards the setup phase. Then, in the solution phase, the actual iterations are carried out. This phase consists of the smoothing, applying the transfer operators, solving on the coarsest level and applying the accelerator, see Section 3.3.3. The important observation here is that if we have to solve two linear systems
Au =f, (4.5)
Au =g (4.6)
in sequence, then we only need to compute the setup phase once and can re-use all the computed operators to also solve the second system.
Remark 4.1. Even if we had to solve
Au=f, (4.7)
Bu=g (4.8)
with B ≈A,A, B ∈Rn×n, we still have a number of options in order to avoid having to compute the full setup forB:
• We can keep the full setup forAwhen solvingB, including all transfer and coarse level operators. Note that when applying the smoother and the Krylov iteration on the finest level, we still useB. But on all the coarse levels, the operators that were originally constructed for solving Au =f are used. In this method, we do
2BLAS implementation [79] [34] [18].
not need to compute anything in the setup phase, but the convergence of the overall method may not be very good, ifA and B are too different.
• Another option is to use the C/F-Splitting (cf. Section 3.3.3) from A but re-compute all the transfer and coarse level operators using B. This way we save the cost of finding a C/F-Splitting for B, but we do need to re-compute all the operators. This gives an algorithm that is much better suited for solvingBu=g but the setup phase is not as cheap as in the first option.
Unfortunately it is very hard to tell for any given A and B if one of the two options above or a completely new setup gives the fastest overall solution method. In some cases it can be a good option to perform a certain number of test cycles when solving Bu=g and then either continue with the iterations if the convergence rate is similar to the one that was observed when solving Au=f or to stop iterating and to create a new setup. Note that this approach may introduce some overhead and its success in terms of trying to save run-time is highly dependent on the application and also the specific model in some cases.
In the case of the FPM at least we know that the hydrostatic pressure system (2.51) and the correction pressure system (2.58) both yield the same matrixA and the only difference is in the corresponding right-hand sides. Therefore we only need to compute one setup in order to solve both of these systems.
For the dynamic pressure system (2.52) the same holds in the default case. However, it is possible to impose different boundary conditions on the dynamic pressure systems for reasons of stability of the Chorin projection approach, see Section 2.3.3. In this case, the matrixB corresponding to the dynamic pressure system changes in those rows that correspond to boundary conditions that are being changed. Not only can this lead to a decreased convergence rate when re-using the setup from the hydrostatic system but these changes may also introduce or remove independent subsystems within the linear system which is important for our AMG method in the parallel case, see Section 5.2. In this case, the setup cannot be re-used because the additional couplings introduced or removed by the change in the boundary conditions that are responsible for merging or creating independent subsystems also have an impact on the communication pattern of our AMG method. Therefore, the setup cannot be re-used without substantial computational effort in the update of the communication patterns which is why in those cases we opt to compute a completely new setup. For these reasons, we only re-use the setup from the hydrostatic pressure system for the dynamic pressure system if the boundary conditions did not change.