On 10/05/16 17:29, Francis Poulin wrote:
Hello Jemma,
Thanks for the suggestion, that helps a lot. You're right that if F = 1.0 it does work without a problem. I tried F = 0.0001 and that also seems to work very well. There is something about the purely barotropic (rigid-lid) problem that the solver does not like. In the case of a channel the F=0.0 case works without a problem.
In pseudo-spectral codes I know that inverting an ellitpic operator has a problem with F = 0 because of the zero wavenumber. Not sure if there is an analogue in FE.
With this I will happily consider F = 1.0 or really, really small but not zero. I will try to make a nice animation of 2D turbulence. When I do I will send it along, in case people are curious to see the results.
If F is zero then your elliptic operator has a nullspace of the constants. If your right hand side is consistent (i.e. orthogonal to the constant vector) then all you need to do is tell the Krylov solver to project the constants out of the solution. See http://firedrakeproject.org/solving-interface.html#solving-singular-systems for details on how to specify this correctly. Cheers, Lawrence