A GMG solver for a simple Poisson test problem on a hexahedron fails when assembled as a linear problem, but succeeds when assembled as a non-linear problem.
This is due to the special treatment linear problems get in the P1-continuity constraint: While non-linear problems assemble zero defect entries and unity lines in the Jacobian (which effectively takes constrained DoFs out of the calculation), linear problems assemble zero right-hand side and an interpolating line in the matrix (thus keeping the constrained DoFs as real DoFs).
This does not seem to go down well in the multigrid hierarchy. But even if the restriction matrix (and prolongation) is adjusted so that the constrained DoFs' defects stay where they are, this does not seem to be a good idea, because the constrained defects are scaled independently from the rest and could - on small-scale geometries - dominate the total defect.
I would much rather treat the linear case just like the non-linear one by setting unity rows in the matrix.
The main problem with this approach is that this will make it necessary for the solution in the hanging nodes to be enforced separately AFTER solving the linear system, which is not a big deal per se, but will require an update in all applications currently using hanging nodes in linear problems.
A GMG solver for a simple Poisson test problem on a hexahedron fails when assembled as a linear problem, but succeeds when assembled as a non-linear problem.
This is due to the special treatment linear problems get in the P1-continuity constraint: While non-linear problems assemble zero defect entries and unity lines in the Jacobian (which effectively takes constrained DoFs out of the calculation), linear problems assemble zero right-hand side and an interpolating line in the matrix (thus keeping the constrained DoFs as real DoFs). This does not seem to go down well in the multigrid hierarchy. But even if the restriction matrix (and prolongation) is adjusted so that the constrained DoFs' defects stay where they are, this does not seem to be a good idea, because the constrained defects are scaled independently from the rest and could - on small-scale geometries - dominate the total defect.
I would much rather treat the linear case just like the non-linear one by setting unity rows in the matrix. The main problem with this approach is that this will make it necessary for the solution in the hanging nodes to be enforced separately AFTER solving the linear system, which is not a big deal per se, but will require an update in all applications currently using hanging nodes in linear problems.
testScriptAndGeom.zip