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?
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))
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
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 David Ham [David.Ham@imperial.ac.uk]
Sent: Tuesday, November 29, 2016 4:51 AM
To: firedrake
Subject: Re: [firedrake] no normal flow Boundary Conditions
Sent: Tuesday, November 29, 2016 4:51 AM
To: firedrake
Subject: Re: [firedrake] no normal flow Boundary Conditions
Hi Francis,
On Tue, 29 Nov 2016 at 04:29 Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
Some progress since last I wrote:
1) I was able to get the QG linear solver working. I do have a couple of questions though.
https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_qg.py
First, the solver that Jemma and I used for the nonlinear QG code seems to yield the trivial solution. Any idea why that is?
I don't think it's finding a trivial solution. I think that the issue is that you never solve the system. Having set up your linear solver, you are missing:
psi_solver.solve()
However, the linear solver in your code has another problem.
The equation isn't symmetric (the beta*phi*psi.dx(0)*dx term) so CG is unlikely to work. For 2D problems, using a direct solver is usually the right start:
solver_parameters={
'ksp_type':'preonly',
'pc_type':'lu'
}
Second, should we be able to use matplotlib to plot the solutions? I was able to get a demo to work but I can't do it for this script.
There are two problems here. The first is that you are calling the plot command from Matplotlib. Unsurprisingly, this does not know about Firedrake Functions. You need to call the Firedrake plot command.
p = plot(psi_soln)
p.show()
However, more importantly, the matplotlib functionality is mostly useful for debugging small systems. Your 50 x 50 mesh is already going to be very slow via this route.
2) The SW problem is not working much better I'm afraid. It does not converge to a solution, which is I suppose better than converging to an incorrect solution. I have a couple of concerns.
https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw.py
First, i am trying to impose no-normal flow BCs. I am trying to do what I understood Colin suggested but since it's not converging, maybe I misunderstood.
Second, I presume there is a better solver that I should be using?
Same problem as above. You are attempting to use CG for a non-symmetric system. In this case it's a saddle point system. Try using a direct solver instead.
Regards,
David
Thanks in advance for your help.
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 Francis Poulin [fpoulin@uwaterloo.ca]
Sent: Saturday, November 26, 2016 6:54 PMHello,
Minor progress.
I found out that the problem with my QG code was how I imposed the boundary conditions. Now I am happy to say that they both run. Unfortunately, they both give me a zero solution even though the winds are not zero. I plotting the winds and so I know that they're not zero. Maybe I"ll have more luck tomorrow.
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 Francis Poulin [fpoulin@uwaterloo.ca]
Sent: Saturday, November 26, 2016 12:58 PM
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] no normal flow Boundary Conditions
Hello again,
Minor update on my wind-driven gyre problem.
I just put together a SW version of the code that solve the Linear Stommel problem. You can find the code in the same github repo,
https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw.py
To my surprise the solver does not complain, and it is exactly the same thing as in the QG version so maybe that part of it is ok. Unfortunately, the solution seems to be zero. So this code is at least giving me an answer, but a wrong one.
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 Francis Poulin [fpoulin@uwaterloo.ca]
Sent: Saturday, November 26, 2016 7:36 AM
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] no normal flow Boundary Conditions
Hello,
Thank you Colin and everyone else for your responses. That was very helpful.
We are working on a bunch of problems right now, both QG and SW. We have many different codes working in FEniCS but I am going to start translating them into Firedrake. Maybe make a demo if people seem interested?
The first code is to solve the Linear Stommel problem in one-layer QG. I was able to translate most of it ok but I seem to be having problems with the solver for some reason.
If anyone had any advice on what I am doing wrong that would be greatly appreciated. You can find a link on github
https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_qg.py
After I get this working I will like to solve a variety of other problems. One is the Munk problem, then look at non-linear gyre solutions, then find these solutions in SW.
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 Colin Cotter [colin.cotter@imperial.ac.uk]
Sent: Friday, November 25, 2016 4:15 AM
To: firedrake
Subject: Re: [firedrake] no normal flow Boundary Conditions
Hi Francis,That depends on the formulation of the equations. If you are using a streamfunction-vorticity formulation as for QG, then you need to set a zero boundary condition for the streamfunction, as a scalar. If you are using a velocity formulation based on H(div) spaces (BDM, RT etc.) then there are only u.n DOFs on the boundary, so setting:bcs = [DirichletBC(Z.sub(0), Constant((0, 0)), (blah,)), ...actually only sets the normal component. That's one of the reasons that I like these spaces for GFD.cheers--cjc
On 25 November 2016 at 01:17, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
Sorry for the bother but I have a simple question to ask.
I am currently developing some code to solve for wind-driven gyres and basin modes using QG and SW. My question is how is it best to impose no-normal flow BCs,
$$\vec u \cdot \hat n = 0$$
for any boundary. I found this in the Navier-Stokes example which says that we should impose zero on particular boundaries. I guess that would work but if one had a more complicated geometry I imagine that would have it's problems.
bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))]
Any suggestions would be greatly appreciated.
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
--
