Re: [firedrake] REXI time stepping in Firedrake
Hi Martin, I discussed this complex number issue with Terry Haut a while ago. He has an equivalent formulation that stays entirely in real numbers - this is because the complex values appear as conjugate pairs. So, it's probably easiest to move to this real-only formulation. By the way, the Coriolis term does make a big difference to how things are solved, I think that it will only be possible to efficiently solve this system using hybridisation techniques: this was the foundation of the proposal that we put in to NERC with Beth in January. To answer your other questions: we are converging towards a Firedrake dynamical core code, this is being developed on github: https://github.com/firedrakeproject/dcore but things are not really ready for integration with REXI yet, we still need to demonstrate that the basic discretisation works for 3D on the sphere. We have a plan to implement a multigrid solver for the the REXI scheme for FEM shallow water, and then develop an APINT scheme using Firedrake for shallow water on the sphere with Jemma, we should probably discuss to make sure that we don't tread on each other's toes. all the best --cjc On 30 April 2016 at 07:35, Martin SCHREIBER <M.Schreiber@exeter.ac.uk> wrote:
Dear Firedrakers,
I discussed yesterday with Andrew if Firedrake can be (ab)used to solve PDEs with a non-traditional time stepping method (we call it REXI) and he suggested me to drop an Email to this mailing list. (@David/Colin: It's related to our parallelization-in-time work).
Let L be the linear operator of the SWE on the sphere in advective form including the average height (by putting the perturbation into the non-linear part). L also contains the longitude-varying Coriolis term, but this shouldn't make a big difference here.
As one building block of the new time stepping method we like to solve the SoEs given by (\alpha - L) U = U0 with \alpha a complex shifted pole creating a non-singular SoE and U0 the DoFs to solve for.
I don't see any problem to write this in weak formulation as it is required in Firedrake.
Now the following questions arise:
* Is Firedrake able to assemble and solve these system of equations with the complex value \alpha? It's only important to keep the real values of "U". The imaginary values of U are not important if that's of any help.
* If not, is Firedrake still able to solve a reformulation which would have a 6x6 linear operator instead of a 3x3 and 6 DoFs (3 pseudo DoFs) per cell? (Basically reformulating the SoE with real values instead of complex ones).
* Is there any Firedrake code publicly available with the GungHo FEM dyn-core on the sphere with a C-grid?
* A question related to possible future work: Is Firedrake supporting Semi-Lagrangian methods efficiently?
Cheers,
Martin
-- Dr. rer.-nat. Martin Schreiber Lecturer, proleptic Scientific and High-Performance Computing University of Exeter
College of Engineering, Mathematics and Physical Sciences Harrison Building, Room 319, North Park Road, Exeter, EX4 4QF
_______________________________________________ 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
A good starting point to see how this should be done is the sw_linear_williamson2_triangle.py example in dcore/examples The solver there is not robust for large dt relative to the Coriolis timescale, but we are working on a hybridised solver for this that will be. If you just want to look into the REXI parallelism aspects without combining with spatial parallelism then you could use a direct solver on small problems for this. Let's talk more offline about it to avoid clogging up the list. all the best --cjc On 30 April 2016 at 20:40, Colin Cotter <colin.cotter@imperial.ac.uk> wrote:
Hi Martin, I discussed this complex number issue with Terry Haut a while ago. He has an equivalent formulation that stays entirely in real numbers - this is because the complex values appear as conjugate pairs. So, it's probably easiest to move to this real-only formulation.
By the way, the Coriolis term does make a big difference to how things are solved, I think that it will only be possible to efficiently solve this system using hybridisation techniques: this was the foundation of the proposal that we put in to NERC with Beth in January.
To answer your other questions: we are converging towards a Firedrake dynamical core code, this is being developed on github: https://github.com/firedrakeproject/dcore but things are not really ready for integration with REXI yet, we still need to demonstrate that the basic discretisation works for 3D on the sphere.
We have a plan to implement a multigrid solver for the the REXI scheme for FEM shallow water, and then develop an APINT scheme using Firedrake for shallow water on the sphere with Jemma, we should probably discuss to make sure that we don't tread on each other's toes.
all the best --cjc
On 30 April 2016 at 07:35, Martin SCHREIBER <M.Schreiber@exeter.ac.uk> wrote:
Dear Firedrakers,
I discussed yesterday with Andrew if Firedrake can be (ab)used to solve PDEs with a non-traditional time stepping method (we call it REXI) and he suggested me to drop an Email to this mailing list. (@David/Colin: It's related to our parallelization-in-time work).
Let L be the linear operator of the SWE on the sphere in advective form including the average height (by putting the perturbation into the non-linear part). L also contains the longitude-varying Coriolis term, but this shouldn't make a big difference here.
As one building block of the new time stepping method we like to solve the SoEs given by (\alpha - L) U = U0 with \alpha a complex shifted pole creating a non-singular SoE and U0 the DoFs to solve for.
I don't see any problem to write this in weak formulation as it is required in Firedrake.
Now the following questions arise:
* Is Firedrake able to assemble and solve these system of equations with the complex value \alpha? It's only important to keep the real values of "U". The imaginary values of U are not important if that's of any help.
* If not, is Firedrake still able to solve a reformulation which would have a 6x6 linear operator instead of a 3x3 and 6 DoFs (3 pseudo DoFs) per cell? (Basically reformulating the SoE with real values instead of complex ones).
* Is there any Firedrake code publicly available with the GungHo FEM dyn-core on the sphere with a C-grid?
* A question related to possible future work: Is Firedrake supporting Semi-Lagrangian methods efficiently?
Cheers,
Martin
-- Dr. rer.-nat. Martin Schreiber Lecturer, proleptic Scientific and High-Performance Computing University of Exeter
College of Engineering, Mathematics and Physical Sciences Harrison Building, Room 319, North Park Road, Exeter, EX4 4QF
_______________________________________________ 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
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916
Hello Colin, As a bit of an aside, is there somewhere I can read more about the "complex number issue" that you discussed with Terry Haut? I ask because I have been trying to develop code in FEniCS to compute the linear stability characteristics in different models, i.e. QG, SW and whatever else I can muster. I have had success using the direct solver. Unfortunately, the way the problem is formulated I want to find the eigenvalue with the largest imaginary part. as I'm sure you already know, the imaginary numbers in PETSc and SLEPc are not part of the default FEniCS package, or so I've been told. This seems to suggest that I need to use a direct solver, which is problematic when I want to solve 2D eigenvalue problems. I wonder whether the complex number issue you mentioned might help me to reformulate it in such a way that the eigenvalue I want is real, which FEniCS can solve? I would be very happy to do all of this in firedrake but I have not seen any eigenvalue examples and have heard this is on the wish list, but I am sure you guys have a long wish list. If there is anything I can do to help in testing this I would be happy to. Best regards, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Colin Cotter [colin.cotter@imperial.ac.uk] Sent: Saturday, April 30, 2016 3:40 PM To: firedrake Subject: Re: [firedrake] REXI time stepping in Firedrake Hi Martin, I discussed this complex number issue with Terry Haut a while ago. He has an equivalent formulation that stays entirely in real numbers - this is because the complex values appear as conjugate pairs. So, it's probably easiest to move to this real-only formulation. By the way, the Coriolis term does make a big difference to how things are solved, I think that it will only be possible to efficiently solve this system using hybridisation techniques: this was the foundation of the proposal that we put in to NERC with Beth in January. To answer your other questions: we are converging towards a Firedrake dynamical core code, this is being developed on github: https://github.com/firedrakeproject/dcore but things are not really ready for integration with REXI yet, we still need to demonstrate that the basic discretisation works for 3D on the sphere. We have a plan to implement a multigrid solver for the the REXI scheme for FEM shallow water, and then develop an APINT scheme using Firedrake for shallow water on the sphere with Jemma, we should probably discuss to make sure that we don't tread on each other's toes. all the best --cjc On 30 April 2016 at 07:35, Martin SCHREIBER <M.Schreiber@exeter.ac.uk<mailto:M.Schreiber@exeter.ac.uk>> wrote: Dear Firedrakers, I discussed yesterday with Andrew if Firedrake can be (ab)used to solve PDEs with a non-traditional time stepping method (we call it REXI) and he suggested me to drop an Email to this mailing list. (@David/Colin: It's related to our parallelization-in-time work). Let L be the linear operator of the SWE on the sphere in advective form including the average height (by putting the perturbation into the non-linear part). L also contains the longitude-varying Coriolis term, but this shouldn't make a big difference here. As one building block of the new time stepping method we like to solve the SoEs given by (\alpha - L) U = U0 with \alpha a complex shifted pole creating a non-singular SoE and U0 the DoFs to solve for. I don't see any problem to write this in weak formulation as it is required in Firedrake. Now the following questions arise: * Is Firedrake able to assemble and solve these system of equations with the complex value \alpha? It's only important to keep the real values of "U". The imaginary values of U are not important if that's of any help. * If not, is Firedrake still able to solve a reformulation which would have a 6x6 linear operator instead of a 3x3 and 6 DoFs (3 pseudo DoFs) per cell? (Basically reformulating the SoE with real values instead of complex ones). * Is there any Firedrake code publicly available with the GungHo FEM dyn-core on the sphere with a C-grid? * A question related to possible future work: Is Firedrake supporting Semi-Lagrangian methods efficiently? Cheers, Martin -- Dr. rer.-nat. Martin Schreiber Lecturer, proleptic Scientific and High-Performance Computing University of Exeter College of Engineering, Mathematics and Physical Sciences Harrison Building, Room 319, North Park Road, Exeter, EX4 4QF _______________________________________________ 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]
Dear Francis,
On 30 Apr 2016, at 21:14, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello Colin,
As a bit of an aside, is there somewhere I can read more about the "complex number issue" that you discussed with Terry Haut?
I ask because I have been trying to develop code in FEniCS to compute the linear stability characteristics in different models, i.e. QG, SW and whatever else I can muster. I have had success using the direct solver. Unfortunately, the way the problem is formulated I want to find the eigenvalue with the largest imaginary part. as I'm sure you already know, the imaginary numbers in PETSc and SLEPc are not part of the default FEniCS package, or so I've been told. This seems to suggest that I need to use a direct solver, which is problematic when I want to solve 2D eigenvalue problems. I wonder whether the complex number issue you mentioned might help me to reformulate it in such a way that the eigenvalue I want is real, which FEniCS can solve?
As long as your operator has real-valued entries, you can happily use SLEPc (via slepc4py) to compute eigenmodes with complex eigenvalues. If you want the largest imaginary part for your eigenvalues you just have to ask SLEPc to target the appropriate part of the spectrum. Although we have not hooked up a high-level interface for this, you can happily use SLEPc with the operators that Firedrake assembles. You can see one test doing this for the laplacian here: https://github.com/firedrakeproject/firedrake/blob/master/tests/regression/t... Modern versions of firedrake-install allow you to pass --slepc to install SLEPc compatibly with the rest of the distribution. This carries over to firedrake-update (so you can type firedrake-update --slepc to update your installation and build SLEPc as well).
I would be very happy to do all of this in firedrake but I have not seen any eigenvalue examples and have heard this is on the wish list, but I am sure you guys have a long wish list. If there is anything I can do to help in testing this I would be happy to.
Please let us know if the above approach works for you. Cheers, Lawrence
Hi Colin, On 30/04/16 20:40, Colin Cotter wrote:
Hi Martin, I discussed this complex number issue with Terry Haut a while ago. He has an equivalent formulation that stays entirely in real numbers - this is because the complex values appear as conjugate pairs. So, it's probably easiest to move to this real-only formulation.
I thought that the conjugate pairs (of \alpha) are important to reduce the sum of REXI terms to half of it but not for the reformulation to a real-only problem. This should be always possible - independent of the conjugate pairs (and also goes partly back to a document from Terry one year ago!).
By the way, the Coriolis term does make a big difference to how things are solved, I think that it will only be possible to efficiently solve this system using hybridisation techniques: this was the foundation of the proposal that we put in to NERC with Beth in January.
Back in 2015 I thought that it's a good idea to numerically investigate iterative solvers (Jacobi / CG / MG) for REXI in SWEET since I thought that we'd need iterative solvers for LFRic/GungHo to get REXI running. You can find a lynx document with some preliminary results about this here: https://github.com/schreiberx/sweet/tree/master/doc/notes_on_iterative_solve... It's just a sketch to start writing some documentation on the problem but it turned out to be by far more challenging. It was just a test implementation, but none of the implemented solvers seemed to converge fast enough or only with a very low (and useless) resolution. MG was implemented with a 2x2 <-> 1x1 restr./prolong. CG was only converging if implementing a non-standard residual computation. The "it doesn't make difference" referred to the weak formulation of the Coriolis term. It should be definitively challenging to write a GMG for this, but I don't see a major problem with an AMG solver. They should be able to cope with the varying Coriolis term and in particular with the \alpha poles generating large off-diagonal values. I'm more worried about numerical cancellation issues since we might get get several orders of magnitude different values in the matrix due to large imaginary \alpha values for large time steps.
To answer your other questions: we are converging towards a Firedrake dynamical core code, this is being developed on github: https://github.com/firedrakeproject/dcore but things are not really ready for integration with REXI yet, we still need to demonstrate that the basic discretisation works for 3D on the sphere.
Nice! Weather-drake ;-)
We have a plan to implement a multigrid solver for the the REXI scheme for FEM shallow water, and then develop an APINT scheme using Firedrake for shallow water on the sphere with Jemma, we should probably discuss to make sure that we don't tread on each other's toes.
Thanks for letting me know about the NERC project. Apart from the APinT stuff it sounds very much like what I am/was partly working on (also involving other researchers...). Maybe next time all the people who are *really* working on this should hand in a project... Let's discuss the rest off-list. Cheers, Martin
all the best --cjc
On 30 April 2016 at 07:35, Martin SCHREIBER <M.Schreiber@exeter.ac.uk <mailto:M.Schreiber@exeter.ac.uk>> wrote:
Dear Firedrakers,
I discussed yesterday with Andrew if Firedrake can be (ab)used to solve PDEs with a non-traditional time stepping method (we call it REXI) and he suggested me to drop an Email to this mailing list. (@David/Colin: It's related to our parallelization-in-time work).
Let L be the linear operator of the SWE on the sphere in advective form including the average height (by putting the perturbation into the non-linear part). L also contains the longitude-varying Coriolis term, but this shouldn't make a big difference here.
As one building block of the new time stepping method we like to solve the SoEs given by (\alpha - L) U = U0 with \alpha a complex shifted pole creating a non-singular SoE and U0 the DoFs to solve for.
I don't see any problem to write this in weak formulation as it is required in Firedrake.
Now the following questions arise:
* Is Firedrake able to assemble and solve these system of equations with the complex value \alpha? It's only important to keep the real values of "U". The imaginary values of U are not important if that's of any help.
* If not, is Firedrake still able to solve a reformulation which would have a 6x6 linear operator instead of a 3x3 and 6 DoFs (3 pseudo DoFs) per cell? (Basically reformulating the SoE with real values instead of complex ones).
* Is there any Firedrake code publicly available with the GungHo FEM dyn-core on the sphere with a C-grid?
* A question related to possible future work: Is Firedrake supporting Semi-Lagrangian methods efficiently?
Cheers,
Martin
-- Dr. rer.-nat. Martin Schreiber Lecturer, proleptic Scientific and High-Performance Computing University of Exeter
College of Engineering, Mathematics and Physical Sciences Harrison Building, Room 319, North Park Road, Exeter, EX4 4QF
_______________________________________________ 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>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr. rer.-nat. Martin Schreiber Lecturer, proleptic Scientific and High-Performance Computing University of Exeter College of Engineering, Mathematics and Physical Sciences Harrison Building, Room 319, North Park Road, Exeter, EX4 4QF
On 1 May 2016, at 09:51, Martin SCHREIBER <M.Schreiber@exeter.ac.uk> wrote:
It was just a test implementation, but none of the implemented solvers seemed to converge fast enough or only with a very low (and useless) resolution. MG was implemented with a 2x2 <-> 1x1 restr./prolong. CG was only converging if implementing a non-standard residual computation.
I don't know quite what you mean by this, but note that CG minimises the A-norm of the error irrespective of what residual calculation you choose to terminate the iteration. So be careful! Lawrence
participants (4)
-
Colin Cotter
-
Francis Poulin
-
Lawrence Mitchell
-
Martin SCHREIBER