[petsc-dev] [petsc-users] gamg failure with petsc-dev

[Stephan Kramer]

On 18/09/14 15:41, Mark Adams wrote:
> Thanks for updating me Stephan,
>
> This is a common problem.  I'm surprised we have not seen it more.
>
> Note, the problem with the solver dying when MatSetBlockSize is set
> (correctly).  Is some sort of bug.  You have Stokes and use FieldSplit,
> so GAMG is working on a 2D elasticity like problem. Is that correct?
>   And the equations are ordered node first:  [v0_x, v0_y, v1_x, v1_y,
> v2_x, ....]  correct?

Correct - although we're not using FieldSplit because of historic 
reasons (the code predates fieldsplit) - but yes, it's the velocity 
block K of a Stokes system (K G; G^T 0) with variable viscosity.

I finally figured out why the convergence was so much worse if we do set 
the block size: It turns out to do with the scaling of the number I put 
on the diagonal for the eliminated DOFs associated with strong bcs. This 
problem has free slip bcs and I eliminate the DOF associated with the 
normal component by zeroing out the row and column and putting an 
arbitrary number on the diagonal (I chose 1.0). Normally the number 
doesn't really matter as the row is completely decoupled from the rest 
of the system. However in this case, when using gamg with blocks, it 
affects the strong coupling criterion for *both* velocity components 
associated with the boundary node. This criterion is scaled with respect 
to the 1-norm of the diagonal block, thus if the diagonal for the normal 
component is chosen too big it also declares all connections to the dof 
associated with the tangential component weak. In this way all the 
boundary nodes (the pairs of velocity components at these nodes) became 
isolated clusters and were disregarded at the coarse level. After 
choosing a smaller value for the diagonals the solve now converges fine 
in roughly the same number of iterations as when not setting the block size.

The tricky thing here is the choice of the dummy diagonal number becomes 
very sensitive - with a typical gamg threshold of 0.01 it can't be more 
than an order out. Currently I use MatZeroRowsColumns() with a constant 
value for the diagonal. Because I have a variable viscosity it seems I 
can no longer do that. Now I can just go in and look at what the 
diagonal is for the other component at each boundary node and use that 
value instead (or at least a value that's smaller) - but it all seems a 
bit hackish to me. Perhaps we're following the wrong approach altogether 
and should be dealing with the bcs in an other way? Any thoughts would 
be more than welcome...

Cheers
Stephan