Hi Francis, I think you have the syntax correct. When I tried this it at least changed the error! I now get DIVERGED_INDEFINITE_PC instead of DIVERGED_DTOL - is that what you see too? Jemma ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 11 May 2016 02:44:42 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Thank you Lawrence! That helps a great deal in understanding the problem. I looked at the example you suggested and they point out that you need to add these lines nullspace = VectorSpaceBasis(constant=True) solve(a == L, u, nullspace=nullspace) I am guessing that first line I presume should be added as is. The second line, which defines the solver, I would think would have to be changed appropriately. # Set up Elliptic inverter psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) psi_solver = LinearVariationalSolver(psi_problem, solver_parameters={ 'ksp_type':'cg', 'pc_type':'sor' }) I tried adding nullspace=nullspace after the } but that didn't seem to work. What should the syntax be? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk] Sent: Tuesday, May 10, 2016 2:44 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] periodic boundary conditions? 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 _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake