Hello Lawrence, Thanks for the tip on ipython. That is handy. I have tried using an LU solver and get a PETSc error, Traceback (most recent call last): File "linear_stommel_sw.py", line 57, in <module> sw_solver.solve() File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/variational_solver.py", line 196, in solve self.snes.solve(None, v) File "PETSc/SNES.pyx", line 537, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:170683) petsc4py.PETSc.Error: error code 56 [0] SNESSolve() line 4128 in /tmp/pip-SPrpSR-build/src/snes/interface/snes.c [0] SNESSolve_KSPONLY() line 40 in /tmp/pip-SPrpSR-build/src/snes/impls/ksponly/ksponly.c [0] KSPSolve() line 620 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest Again, the updated code can be found at the following in case someone wants to look at it. https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw.p... There is a fenics version of the code that does work that is located in the same directory. https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw_f... Any advice would be greatly appreciated. Best regards, 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, November 29, 2016 8:55 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] no normal flow Boundary Conditions On 29/11/16 13:49, Francis Poulin wrote:
Hello David,
Thanks a lot for the help. Some responses/follow ups below.
1) Sorry that I foolishly forgot to call the solver. That was a pretty silly mistake on my part.
2) Thanks for the advice on what solver to use. That makes a lot of sense.
3) The plot does show up. But it only lasts for a second before it vanishes. Any idea how to get it to last longer?
The plotting we have is designed for interactive use, you can run, for example with ipython: ipython -i my_script.py Then, you'll drop into a prompt when you're done, and you can interact then. See https://ipython.readthedocs.io/en/stable/interactive/plotting.html for how to interact with matplotlib within ipython.
4) I agree that the SW solver should also be modified to use LU. But at the moment I am using the generic solver and it doesn't converge. I suspect it's because the boundary conditions are not set up correctly. If I wanted to set no-normal flow BCs in my BDM velocity field I am doing the following, The fact that it doesn't converge makes me think that it is not set up correctly.
bc = DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3, 4))
If you have a linear problem that is not converging with an iterative method, switching to a direct solve will at least (if that then converges) give you confidence that your formulation is right, you just need to work harder on preconditioning. If the direct solver doesn't work, then you'll want to go back and figure out the problem with your formulation. Lawrence