mantle_convection test inverts a (near) singular matrix

[James Maddison]

The mantle_convection test fails to apply a reference pressure node or supply null space information to the linear solvers. As a result the test is inverting a singular matrix. Rounding errors mean that the assembled matrix is non-singular, and invertible with a direct solver, but very ill-conditioned.

[Marie Rognes (meg-simula)]

Thanks for pointing at this.

The Stokes system(s) in that demo are indeed rather singular. However, the null space is not the usual constant pressure space (because the velocity is not specified on the entire boundary), but is rather spanned by three velocity eigenmodes.

I'm not quite sure what the best fix is. Is this one of the tests included in the dolfin-adjoint paper? If so, I think we should just leave it as is for the sake of reproducibility. If not, we could replace the direct solves with appropriately orthogonalized Krylov solves.

[Patrick Farrell (pefarrell)]

We can fix it, I think -- I need to replace those figures anyway for other reasons. The addition of a sentence on specifying the nullspace of the operator should be possible, I hope.

[Marie Rognes (meg-simula)]

Ok, I'll replace the LUs by appropriately nullspaced Krylov solves in mantle_convection.py then.

[Marie Rognes (meg-simula)]

Sorry, please disregard my original comment -- I forgot to include the boundary conditions when playing around with the testing and forgot to think as well ;-) Yes, there is a pressure mode and I will fix it.

[James Maddison (jamesmadd)]

> Yes, there is a pressure mode and I will fix it.

See timestepping/tests/long/mantle_convection for a simple fix.

[Marie Rognes (meg-simula)]

Thanks. I'm not a big fan of pinning pressures at points, but I couldn't get the adjoint Krylov solve to converge with the constant pressure null space (perhaps the transpose nullspace is something else?) at the moment , so I have pushed a fix using this approach.

Could someone else run mantle_convection/paper.py to check that those results still make sense? (I didn't have time to now.) The quick test (mantle_convection/mantle_convection.py) results seemed to.

[Patrick Farrell (pefarrell)]

Did you manage to take a look at the SVD of the assembled system? The transpose nullspace should just be the left singular vector where the right singular vector is the constant pressure mode. (I imagine you won't have time for this; if not, I'll do it next week.)