Re: [firedrake] periodic boundary conditions?
Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future.
Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails.
If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain.
Is there another solver that might be a bit better at converging?
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: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk> Subject: [firedrake] periodic boundary conditions?
This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below).
Lawrence
Hello,
Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student.
I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at
https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph...
I have removed the jet and only have random initial conditions.
The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic.
I made two changes:
1) I asked it to be periodic in this line by switching "x" to "both".
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None )
2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0)
It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow?
bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )]
When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions?
I hope to make a nice simulation of this to see the structures merging and increasing in size.
What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time.
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916
Hello Colin, I agree that the solver should work and if you think the problem is something else I am happy to investigate other avenues. As for your 3 options: (1) The initial conditions are random. When I use the channel version it does converge fine. (2) Maybe. Below you will see how the mesh is defined in two cases. (3) I haven't touched the equations and thought they were set up correctly but maybe not. When I look at the diff between my two versions, channel vs doubly periodic, I get < mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None) ---
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="x", quadrilateral=True, reorder=None)
59c59 < psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) ---
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
It could be that the channel geometry is more forgiving to an error that is there but I am not sure why that would be or how to check it out. If anyone has any ideas of tests I can do please let me know. I can also share the code somehow if anyone wanted to take look at what it does. 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: Sunday, May 08, 2016 2:57 AM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello, Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future. Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails. If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain. Is there another solver that might be a bit better at converging? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>] Sent: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: [firedrake] periodic boundary conditions? This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below). Lawrence Hello, Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student. I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph... I have removed the jet and only have random initial conditions. The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic. I made two changes: 1) I asked it to be periodic in this line by switching "x" to "both". mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None ) 2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow? bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )] When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions? I hope to make a nice simulation of this to see the structures merging and increasing in size. What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg]
Hello again, Sorry to bother you again with this but I thought I might send a copy of my code that tries to enforce periodic BCs, in case someone wanted to look it over and see if I've done anything clearly wrong. I think it's done correctly but for some reason even when I reduce Dt by 10^5, it still does not converge. Any ideas of things I can try 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 ________________________________ From: Francis Poulin Sent: Sunday, May 08, 2016 8:52 AM To: firedrake@imperial.ac.uk Subject: RE: [firedrake] periodic boundary conditions? Hello Colin, I agree that the solver should work and if you think the problem is something else I am happy to investigate other avenues. As for your 3 options: (1) The initial conditions are random. When I use the channel version it does converge fine. (2) Maybe. Below you will see how the mesh is defined in two cases. (3) I haven't touched the equations and thought they were set up correctly but maybe not. When I look at the diff between my two versions, channel vs doubly periodic, I get < mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None) ---
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="x", quadrilateral=True, reorder=None)
59c59 < psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) ---
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
It could be that the channel geometry is more forgiving to an error that is there but I am not sure why that would be or how to check it out. If anyone has any ideas of tests I can do please let me know. I can also share the code somehow if anyone wanted to take look at what it does. 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: Sunday, May 08, 2016 2:57 AM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello, Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future. Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails. If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain. Is there another solver that might be a bit better at converging? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>] Sent: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: [firedrake] periodic boundary conditions? This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below). Lawrence Hello, Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student. I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph... I have removed the jet and only have random initial conditions. The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic. I made two changes: 1) I asked it to be periodic in this line by switching "x" to "both". mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None ) 2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow? bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )] When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions? I hope to make a nice simulation of this to see the structures merging and increasing in size. What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg]
Hi Francis, The first difference that jumped out at me is that you have set your Froude number to 0.0 - it is 1.0 in the original code. The code runs very nicely with F=1.0 again, although I'm not sure you should be setting Dirichlet BCs on a doubly periodic mesh. Hope that helps, Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 10 May 2016 17:14:10 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hello again, Sorry to bother you again with this but I thought I might send a copy of my code that tries to enforce periodic BCs, in case someone wanted to look it over and see if I've done anything clearly wrong. I think it's done correctly but for some reason even when I reduce Dt by 10^5, it still does not converge. Any ideas of things I can try 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 ________________________________ From: Francis Poulin Sent: Sunday, May 08, 2016 8:52 AM To: firedrake@imperial.ac.uk Subject: RE: [firedrake] periodic boundary conditions? Hello Colin, I agree that the solver should work and if you think the problem is something else I am happy to investigate other avenues. As for your 3 options: (1) The initial conditions are random. When I use the channel version it does converge fine. (2) Maybe. Below you will see how the mesh is defined in two cases. (3) I haven't touched the equations and thought they were set up correctly but maybe not. When I look at the diff between my two versions, channel vs doubly periodic, I get < mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None) ---
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="x", quadrilateral=True, reorder=None)
59c59 < psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) ---
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
It could be that the channel geometry is more forgiving to an error that is there but I am not sure why that would be or how to check it out. If anyone has any ideas of tests I can do please let me know. I can also share the code somehow if anyone wanted to take look at what it does. 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: Sunday, May 08, 2016 2:57 AM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello, Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future. Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails. If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain. Is there another solver that might be a bit better at converging? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>] Sent: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: [firedrake] periodic boundary conditions? This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below). Lawrence Hello, Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student. I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph... I have removed the jet and only have random initial conditions. The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic. I made two changes: 1) I asked it to be periodic in this line by switching "x" to "both". mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None ) 2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow? bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )] When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions? I hope to make a nice simulation of this to see the structures merging and increasing in size. What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg]
Btw, I realise F=1 might not be what you want to run, but that is where the change in the structure of your equations has occurred. ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 10 May 2016 17:22:26 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, The first difference that jumped out at me is that you have set your Froude number to 0.0 - it is 1.0 in the original code. The code runs very nicely with F=1.0 again, although I'm not sure you should be setting Dirichlet BCs on a doubly periodic mesh. Hope that helps, Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 10 May 2016 17:14:10 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hello again, Sorry to bother you again with this but I thought I might send a copy of my code that tries to enforce periodic BCs, in case someone wanted to look it over and see if I've done anything clearly wrong. I think it's done correctly but for some reason even when I reduce Dt by 10^5, it still does not converge. Any ideas of things I can try 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 ________________________________ From: Francis Poulin Sent: Sunday, May 08, 2016 8:52 AM To: firedrake@imperial.ac.uk Subject: RE: [firedrake] periodic boundary conditions? Hello Colin, I agree that the solver should work and if you think the problem is something else I am happy to investigate other avenues. As for your 3 options: (1) The initial conditions are random. When I use the channel version it does converge fine. (2) Maybe. Below you will see how the mesh is defined in two cases. (3) I haven't touched the equations and thought they were set up correctly but maybe not. When I look at the diff between my two versions, channel vs doubly periodic, I get < mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None) ---
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="x", quadrilateral=True, reorder=None)
59c59 < psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) ---
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
It could be that the channel geometry is more forgiving to an error that is there but I am not sure why that would be or how to check it out. If anyone has any ideas of tests I can do please let me know. I can also share the code somehow if anyone wanted to take look at what it does. 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: Sunday, May 08, 2016 2:57 AM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello, Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future. Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails. If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain. Is there another solver that might be a bit better at converging? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>] Sent: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: [firedrake] periodic boundary conditions? This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below). Lawrence Hello, Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student. I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph... I have removed the jet and only have random initial conditions. The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic. I made two changes: 1) I asked it to be periodic in this line by switching "x" to "both". mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None ) 2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow? bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )] When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions? I hope to make a nice simulation of this to see the structures merging and increasing in size. What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg]
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. 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 Shipton, Jemma [j.shipton@imperial.ac.uk] Sent: Tuesday, May 10, 2016 12:22 PM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, The first difference that jumped out at me is that you have set your Froude number to 0.0 - it is 1.0 in the original code. The code runs very nicely with F=1.0 again, although I'm not sure you should be setting Dirichlet BCs on a doubly periodic mesh. Hope that helps, Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 10 May 2016 17:14:10 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hello again, Sorry to bother you again with this but I thought I might send a copy of my code that tries to enforce periodic BCs, in case someone wanted to look it over and see if I've done anything clearly wrong. I think it's done correctly but for some reason even when I reduce Dt by 10^5, it still does not converge. Any ideas of things I can try 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 ________________________________ From: Francis Poulin Sent: Sunday, May 08, 2016 8:52 AM To: firedrake@imperial.ac.uk Subject: RE: [firedrake] periodic boundary conditions? Hello Colin, I agree that the solver should work and if you think the problem is something else I am happy to investigate other avenues. As for your 3 options: (1) The initial conditions are random. When I use the channel version it does converge fine. (2) Maybe. Below you will see how the mesh is defined in two cases. (3) I haven't touched the equations and thought they were set up correctly but maybe not. When I look at the diff between my two versions, channel vs doubly periodic, I get < mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None) ---
mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="x", quadrilateral=True, reorder=None)
59c59 < psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) ---
psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1)
It could be that the channel geometry is more forgiving to an error that is there but I am not sure why that would be or how to check it out. If anyone has any ideas of tests I can do please let me know. I can also share the code somehow if anyone wanted to take look at what it does. 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: Sunday, May 08, 2016 2:57 AM To: firedrake Subject: Re: [firedrake] periodic boundary conditions? Hi Francis, I'm very sceptical that it is anything to do with the solver. Either: (1) The input to the solver is already bad, or, (2) The mesh is messed up somehow, or (3) The equations are messed up somehow. Just checking, you didn't keep non-zero beta did you? That would not work with double periodic. all the best --cjc On 7 May 2016 at 22:33, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello, Sorry Lawrence for sending my email to the wrong address but thanks for forwarding it. I will be more careful in the future. Thanks for the suggestion Colin. I had tried reducing it by a factor of 10 and it still failed. Now I tried a factor of 10^5 and it still fails. If my implementation is correct, which is a big if, then I wonder if the Helmholtz problem that I'm solving is much more poorly conditioned on a doubly periodic domain than just a periodic domain. Is there another solver that might be a bit better at converging? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>] Sent: Saturday, May 07, 2016 11:41 AM To: <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: [firedrake] periodic boundary conditions? This message reached me as a list administrator from Francis Poulin (but I couldn't figure out how to release it onto the list, so instead, it's copied below). Lawrence Hello, Glad to see the posters are looking great. Sorry that I won't get to see them in person this year but glad I could contribute a little. Hopefully I'll get to attend next year with a grad student. I have a question on how to convert the QG code from a channel to periodic boundary conditions. You can find a demo version of the code at https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostroph... I have removed the jet and only have random initial conditions. The good news is that 2D turbulence does happen and structures merge but in a channel. To reproduce classical results I would like it to be doubly periodic. I made two changes: 1) I asked it to be periodic in this line by switching "x" to "both". mesh = PeriodicRectangleMesh(n0, n0, Lx, Ly, direction="both", quadrilateral=True, reorder=None ) 2) We should not impose Dirichlet boundary conditions. To do this I removed bcs=bc1 in the following psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0) It occurs to me that perhaps one needs to explicitly impose periodic boundary conditions. Is that the case? If yes then do I need to change the following somehow? bc1 = [DirichletBC(Vcg, 0., 1 ), DirichletBC(Vcg, 0., 2 )] When I run the code it does run but it doesn't converge (after 42 steps) and dies. This is for the same resolution and everything else that does converge for the channel case. Could it be that the periodic case is not as well conditioned if it is periodic in both directions? I hope to make a nice simulation of this to see the structures merging and increasing in size. What's really nice is that if I wanted to do 3D turbulence in the context of 3D the changes are pretty small. I need to define the grid to be 3D and impose Neumann BCs at the top and bottom. That I will ask about at a later time. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg]
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
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
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
PS In fact, it runs just fine if you delete the solver_parameters={blah} option. ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 11 May 2016 10:36:03 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? 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 _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 11 May 2016, at 10:39, Shipton, Jemma <j.shipton@imperial.ac.uk> wrote:
PS In fact, it runs just fine if you delete the solver_parameters={blah} option.
Beware! This is hiding a problem I think. The default solver parameters are to use GMRES, not CG. If you run with GMRES and 'ksp_monitor_true_residual': True, I get the following output. 0 KSP preconditioned resid norm 1.402486418742e-02 true resid norm 2.145922987616e-02 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02 So you can see that although the preconditioned norm is converging, the unpreconditioned norm is stagnating. If you provide a nullspace to the Krylov method and use GMRES, then the solution will be the minimum residual in the least squares sense in the space orthogonal to the nullspace. However, your solution may still have a large component in the nullspace (which is what I think is going on here). I'm not sure what the correct fix is though, unfortunately. Lawrence
I think you're problem might be that your right-hand side is inconsistent with your nullspace. The nullspace is not only a nullspace on the right-hand of the matrix, but also on the left-hand. If you consider the constant-1 test function, then the left-hand side vanishes (regardless of psi), but if q1 doesn't integrate to 0 over the domain, the right-hand side doesn't. At the level of the linear system, this is equivalent to left-multiplying with a constant vector on both sides: again the lhs disappears and the rhs doesn't. This means there is not solution to your system. If you think q1 should actually integrate to zero, but is not because of numerical error, you would have to project the constant part out before going into the solve. In petsc you could do that by setting a transpose nullspace, or explicitly calling MatNullSpaceRemove on the rhs vector - I'm not sure Firedrake has nice interface to do the same using the 'nullspace' you have Cheers Stephan On 11/05/16 10:55, Lawrence Mitchell wrote:
On 11 May 2016, at 10:39, Shipton, Jemma <j.shipton@imperial.ac.uk> wrote:
PS In fact, it runs just fine if you delete the solver_parameters={blah} option.
Beware! This is hiding a problem I think. The default solver parameters are to use GMRES, not CG.
If you run with GMRES and 'ksp_monitor_true_residual': True, I get the following output.
0 KSP preconditioned resid norm 1.402486418742e-02 true resid norm 2.145922987616e-02 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02
So you can see that although the preconditioned norm is converging, the unpreconditioned norm is stagnating. If you provide a nullspace to the Krylov method and use GMRES, then the solution will be the minimum residual in the least squares sense in the space orthogonal to the nullspace. However, your solution may still have a large component in the nullspace (which is what I think is going on here).
I'm not sure what the correct fix is though, unfortunately.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Stephan's email reminds me that when I was setting random initial conditions for turbulence in my PhD, we had to be a bit careful and ensure the mean PV over the domain was zero. I've just tried this in your code and it works: q0.dat.data[:] += 0.1*np.random.randn(q0.dof_dset.size) Area = Lx*Ly print assemble(q0*dx)/Area qmean = assemble(q0*dx)/Area q0 -= qmean print assemble(q0*dx)/Area Hope that helps! Jemma ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Stephan Kramer <s.kramer@imperial.ac.uk> Sent: 11 May 2016 11:21:05 To: firedrake Subject: Re: [firedrake] periodic boundary conditions? I think you're problem might be that your right-hand side is inconsistent with your nullspace. The nullspace is not only a nullspace on the right-hand of the matrix, but also on the left-hand. If you consider the constant-1 test function, then the left-hand side vanishes (regardless of psi), but if q1 doesn't integrate to 0 over the domain, the right-hand side doesn't. At the level of the linear system, this is equivalent to left-multiplying with a constant vector on both sides: again the lhs disappears and the rhs doesn't. This means there is not solution to your system. If you think q1 should actually integrate to zero, but is not because of numerical error, you would have to project the constant part out before going into the solve. In petsc you could do that by setting a transpose nullspace, or explicitly calling MatNullSpaceRemove on the rhs vector - I'm not sure Firedrake has nice interface to do the same using the 'nullspace' you have Cheers Stephan On 11/05/16 10:55, Lawrence Mitchell wrote:
On 11 May 2016, at 10:39, Shipton, Jemma <j.shipton@imperial.ac.uk> wrote:
PS In fact, it runs just fine if you delete the solver_parameters={blah} option.
Beware! This is hiding a problem I think. The default solver parameters are to use GMRES, not CG.
If you run with GMRES and 'ksp_monitor_true_residual': True, I get the following output.
0 KSP preconditioned resid norm 1.402486418742e-02 true resid norm 2.145922987616e-02 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02
So you can see that although the preconditioned norm is converging, the unpreconditioned norm is stagnating. If you provide a nullspace to the Krylov method and use GMRES, then the solution will be the minimum residual in the least squares sense in the space orthogonal to the nullspace. However, your solution may still have a large component in the nullspace (which is what I think is going on here).
I'm not sure what the correct fix is though, unfortunately.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Thank you very much Jemma, Lawrence and Stephan, Just to recap, the strong form of the problem that I am solving is \nabla_H^2 psi - F \psi = q Initially I am defining q to be a random variable so it mean is close to zero but not quite. -> In a periodic channel geometry the method will converge whether F is non-zero or not. -> In a doubly periodic domain it will only converge of F is non-zero. I did try 1e-5 and that still converged. When I set the mean to zero (thanks again Jemma) I agree that does work, without having to do anything with the nullspace. I guess the moral of the story is that if the domain is periodic it's best to have a PV with zero mean. I will hopefully not forget this anytime soon. I need to fix this up but in case you're interested, here is the starting to a move on 2D turbulence using firedrake. https://youtu.be/jhAI7TMN7_w 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 Stephan Kramer [s.kramer@imperial.ac.uk] Sent: Wednesday, May 11, 2016 6:21 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] periodic boundary conditions? I think you're problem might be that your right-hand side is inconsistent with your nullspace. The nullspace is not only a nullspace on the right-hand of the matrix, but also on the left-hand. If you consider the constant-1 test function, then the left-hand side vanishes (regardless of psi), but if q1 doesn't integrate to 0 over the domain, the right-hand side doesn't. At the level of the linear system, this is equivalent to left-multiplying with a constant vector on both sides: again the lhs disappears and the rhs doesn't. This means there is not solution to your system. If you think q1 should actually integrate to zero, but is not because of numerical error, you would have to project the constant part out before going into the solve. In petsc you could do that by setting a transpose nullspace, or explicitly calling MatNullSpaceRemove on the rhs vector - I'm not sure Firedrake has nice interface to do the same using the 'nullspace' you have Cheers Stephan On 11/05/16 10:55, Lawrence Mitchell wrote:
On 11 May 2016, at 10:39, Shipton, Jemma <j.shipton@imperial.ac.uk> wrote:
PS In fact, it runs just fine if you delete the solver_parameters={blah} option.
Beware! This is hiding a problem I think. The default solver parameters are to use GMRES, not CG.
If you run with GMRES and 'ksp_monitor_true_residual': True, I get the following output.
0 KSP preconditioned resid norm 1.402486418742e-02 true resid norm 2.145922987616e-02 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 1 KSP preconditioned resid norm 6.859114126588e-03 true resid norm 8.235483003565e-03 ||r(i)||/||b|| 3.837734648956e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 2 KSP preconditioned resid norm 4.372980631318e-03 true resid norm 4.523115353529e-03 ||r(i)||/||b|| 2.107771518191e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 3 KSP preconditioned resid norm 3.083878237363e-03 true resid norm 3.181122331413e-03 ||r(i)||/||b|| 1.482402840070e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 4 KSP preconditioned resid norm 2.328339850318e-03 true resid norm 2.482928092561e-03 ||r(i)||/||b|| 1.157044361280e-01 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 5 KSP preconditioned resid norm 1.850645396168e-03 true resid norm 2.076376448198e-03 ||r(i)||/||b|| 9.675913162685e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 6 KSP preconditioned resid norm 1.505245304176e-03 true resid norm 1.869875298565e-03 ||r(i)||/||b|| 8.713617913390e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 7 KSP preconditioned resid norm 1.228163738742e-03 true resid norm 1.716825645339e-03 ||r(i)||/||b|| 8.000406609403e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 8 KSP preconditioned resid norm 9.762580464613e-04 true resid norm 1.604598260859e-03 ||r(i)||/||b|| 7.477427056418e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 9 KSP preconditioned resid norm 7.187761512558e-04 true resid norm 1.509128698387e-03 ||r(i)||/||b|| 7.032538945228e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 10 KSP preconditioned resid norm 4.857940540733e-04 true resid norm 1.430687141460e-03 ||r(i)||/||b|| 6.667001330972e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 11 KSP preconditioned resid norm 3.400753386006e-04 true resid norm 1.370869815346e-03 ||r(i)||/||b|| 6.388252622562e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 12 KSP preconditioned resid norm 2.724827673558e-04 true resid norm 1.357871715518e-03 ||r(i)||/||b|| 6.327681484164e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 13 KSP preconditioned resid norm 2.023603200214e-04 true resid norm 1.348266511796e-03 ||r(i)||/||b|| 6.282921239844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 14 KSP preconditioned resid norm 1.405713028843e-04 true resid norm 1.342655679836e-03 ||r(i)||/||b|| 6.256774765844e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 15 KSP preconditioned resid norm 9.092448503001e-05 true resid norm 1.337633948725e-03 ||r(i)||/||b|| 6.233373501493e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 16 KSP preconditioned resid norm 5.750986093948e-05 true resid norm 1.335928348383e-03 ||r(i)||/||b|| 6.225425404790e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 17 KSP preconditioned resid norm 3.802056816287e-05 true resid norm 1.335095336743e-03 ||r(i)||/||b|| 6.221543570986e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 18 KSP preconditioned resid norm 2.223795106365e-05 true resid norm 1.334942349117e-03 ||r(i)||/||b|| 6.220830648727e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 19 KSP preconditioned resid norm 1.199219664791e-05 true resid norm 1.334849886266e-03 ||r(i)||/||b|| 6.220399771888e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 20 KSP preconditioned resid norm 6.045970214995e-06 true resid norm 1.334822540473e-03 ||r(i)||/||b|| 6.220272340509e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 21 KSP preconditioned resid norm 3.357136102289e-06 true resid norm 1.334814998460e-03 ||r(i)||/||b|| 6.220237194734e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 22 KSP preconditioned resid norm 2.109412994583e-06 true resid norm 1.334823835708e-03 ||r(i)||/||b|| 6.220278376303e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 23 KSP preconditioned resid norm 1.697698874597e-06 true resid norm 1.334823113503e-03 ||r(i)||/||b|| 6.220275010828e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 24 KSP preconditioned resid norm 1.512292251273e-06 true resid norm 1.334824567918e-03 ||r(i)||/||b|| 6.220281788403e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 25 KSP preconditioned resid norm 1.250309415572e-06 true resid norm 1.334825572707e-03 ||r(i)||/||b|| 6.220286470715e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 26 KSP preconditioned resid norm 7.983659954253e-07 true resid norm 1.334829784399e-03 ||r(i)||/||b|| 6.220306097201e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 27 KSP preconditioned resid norm 5.434753011846e-07 true resid norm 1.334830275105e-03 ||r(i)||/||b|| 6.220308383892e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 28 KSP preconditioned resid norm 3.583899393830e-07 true resid norm 1.334831849286e-03 ||r(i)||/||b|| 6.220315719574e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 29 KSP preconditioned resid norm 2.531881383802e-07 true resid norm 1.334830569064e-03 ||r(i)||/||b|| 6.220309753738e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 30 KSP preconditioned resid norm 1.983419151306e-07 true resid norm 1.334831092617e-03 ||r(i)||/||b|| 6.220312193496e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 31 KSP preconditioned resid norm 1.685855333188e-07 true resid norm 1.334830361112e-03 ||r(i)||/||b|| 6.220308784684e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 32 KSP preconditioned resid norm 1.289738938190e-07 true resid norm 1.334830636894e-03 ||r(i)||/||b|| 6.220310069824e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 33 KSP preconditioned resid norm 9.235117500873e-08 true resid norm 1.334830368242e-03 ||r(i)||/||b|| 6.220308817906e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 34 KSP preconditioned resid norm 6.503935947986e-08 true resid norm 1.334830235572e-03 ||r(i)||/||b|| 6.220308199666e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 35 KSP preconditioned resid norm 4.003439088177e-08 true resid norm 1.334830289372e-03 ||r(i)||/||b|| 6.220308450375e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 36 KSP preconditioned resid norm 2.933016192914e-08 true resid norm 1.334830350834e-03 ||r(i)||/||b|| 6.220308736789e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 37 KSP preconditioned resid norm 2.457785557100e-08 true resid norm 1.334830341474e-03 ||r(i)||/||b|| 6.220308693167e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 38 KSP preconditioned resid norm 2.003045563251e-08 true resid norm 1.334830452035e-03 ||r(i)||/||b|| 6.220309208385e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 39 KSP preconditioned resid norm 1.592346775337e-08 true resid norm 1.334830443617e-03 ||r(i)||/||b|| 6.220309169157e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 40 KSP preconditioned resid norm 1.285777369303e-08 true resid norm 1.334830523051e-03 ||r(i)||/||b|| 6.220309539319e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 41 KSP preconditioned resid norm 9.550269214853e-09 true resid norm 1.334830544472e-03 ||r(i)||/||b|| 6.220309639139e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 42 KSP preconditioned resid norm 6.311572416881e-09 true resid norm 1.334830559098e-03 ||r(i)||/||b|| 6.220309707295e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 43 KSP preconditioned resid norm 3.652454149992e-09 true resid norm 1.334830540583e-03 ||r(i)||/||b|| 6.220309621017e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 44 KSP preconditioned resid norm 2.037777932328e-09 true resid norm 1.334830540790e-03 ||r(i)||/||b|| 6.220309621981e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02 45 KSP preconditioned resid norm 1.128203860755e-09 true resid norm 1.334830532154e-03 ||r(i)||/||b|| 6.220309581738e-02
So you can see that although the preconditioned norm is converging, the unpreconditioned norm is stagnating. If you provide a nullspace to the Krylov method and use GMRES, then the solution will be the minimum residual in the least squares sense in the space orthogonal to the nullspace. However, your solution may still have a large component in the nullspace (which is what I think is going on here).
I'm not sure what the correct fix is though, unfortunately.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (5)
-
Colin Cotter
-
Francis Poulin
-
Lawrence Mitchell
-
Shipton, Jemma
-
Stephan Kramer