exadg
Constraining the 0th DoF in case of continious Galerkin multigrid with singular operators
I was working on using global coarsening multigrid for the Taylor-Green Vortex application. But the multigrid implementation has the following Assert: https://github.com/exadg/exadg/blob/aec5721276aa012b3121b4daafc70b8c6505a0d5/include/exadg/solvers_and_preconditioners/multigrid/multigrid_preconditioner_base.cpp#L490
@peterrum suggested me to check how singular operators are addressed while using local smoothing multigrid. For the continuous case the following is done, whereas nothing is done for the DG case: https://github.com/exadg/exadg/blob/aec5721276aa012b3121b4daafc70b8c6505a0d5/include/exadg/solvers_and_preconditioners/multigrid/constraints.h#L145-L182 During our meeting today, @nfehn and @bergbauer raised the question if this is actually needed at all. This is why I am opening this issue. Unfortunately, I am not familiar with the reasoning or the implementation.
The reasoning is as follows: If we have a singular Laplace operator, the usual cause is that the kernel of A contains the vector z of all constants, i.e., A * z = 0 even though norm(z) != 0. When trying to solve a system of the kind
A * x = b
with singular A, we can have two cases. If b is in the image of A, i.e., there exists a vector y such that A * y = b, then every x = y + alpha * z is also a solution of A * x = b because A * z is zero and hence does not change the expression on the left hand side. This means that we have infinitely many solutions. In the second case, b is not in the image of A, so we cannot solve the given linear system.
If we want our numerical method to return a unique solution (that shows some appropriate measure of convergence), we must hence pose an additional condition and make sure b is compatible with A * x. In DG discretizations, our approach is to say that we want to find a solution with zero mean, which we achieve by manipulating the right hand side and solution at certain places in the algorithm. In theory, it would be enough to start the conjugate gradient method with an initial condition and right hand side that is compatible with the zero-mean solution, as any Krylov method will search the solution in the space spanned by the residual A * x - b, and all new information is created in the image of A (by successive multiplication of vectors by A, which is the only information added). However, in practice there are roundoff errors, especially for the coarse grid solver, so we might introduce some perturbations to this ideal case. Hence, we use the code https://github.com/exadg/exadg/blob/aec5721276aa012b3121b4daafc70b8c6505a0d5/include/exadg/solvers_and_preconditioners/multigrid/coarse_grid_solvers.h#L178-L179 which is the most sensitive place in multigrid regarding roundoff, besides the fix to the right hand side.
Continuous elements are more prone to roundoff effects that shift the solution around and create convergence issues. I don't recall if I understood it at some point in the past, but what I remember is that it has to do with the conservation properties of the method, which are stronger for DG than continuous FEM. In either case, we then manipulate a single entry in the matrix by a constraint to remove the singularity altogether. If it works, it is more robust than the DG setting, because we do no longer rely on the property of the Krylov method to preserve the invariants, but have a non-singular system to begin with. The downside is that we need to do it throughout all levels of the MG hierarchy in a consistent way. By consistent way I mean that these conditions should not contradict each other on different levels, because this choice means that we now have the constraint that the solution ought to be zero in a particular point, and if we request this on different points in physical space, the coarse grid corrections will get incompatible.
I think one could try to skip adding that single constraint and instead rely on the same mechanism as in the DG case, as long as the right hand sides are all fine. Does this help you understanding. There is literature about this Krylov projection property for singular systems you could google for. Can you try what happens?
Thanks for the detailed explanation. It helped a lot. I will also read into the chapter regarding this in Niklas's dissertation.
I will give it a try without the constraint and show the results here. But as far as I understand, removing the constraint might result in convergence problems caused by the roundoff effects. Wouldn't this be something problem dependent as well? Meaning, I could observe good convergence for the problem I choose, but this wouldn't rule out that this might cause convergence issues for other, more challenging problems.
By consistent way I mean that these conditions should not contradict each other on different levels, because this choice means that we now have the constraint that the solution ought to be zero in a particular point, and if we request this on different points in physical space, the coarse grid corrections will get incompatible.
This is the cause of the concern Peter had about applying the constraint on a random (0th) DoF using global coarsening multigrid, without extra care about the location of the point on different levels in the physical space. We might end up constraining different physical points on different levels, which would cause the coarse grid corrections to get incompatible.