- firedrake - imperial.ac.uk

Re: [firedrake] linear stability problem in firedrake
by Colin Cotter 03 May '16

03 May '16

Hi Francis, I'm hoping Patrick Farrell might be able to comment as he did something like this using dolfin-adjoint (I provided some test problems for the paper but wasn't involved in the implementation). all the best --cjc On 1 May 2016 at 17:01, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Dear Lawrence, David and Colin, > > Thanks a lot for the helpful information. I changed the subject since I > thought it might be helpful to have it match the topic of discussion. > > Now that I have an example on how to do this in firedrake I think that I > could probably adapt my QG stability code (which would go very well with > the nonlinear QG demo that I posted a while back). If there is an option > to use a direct eigenvalue method (lapack) then that would work for 1D > eigenvalue problems. Since the most unstable mode corresponds to the > eigenvalue with the largest imaginary part, and I don't suspect that > firedrake allows for imaginary targets, I don't imagine it would converge. > That's the problem I have in FEniCS. > > Colin: Thanks for the advice. For solving linear systems I know how to > decompose the problem into a real and imaginary part and that would be easy > enough to do. However, for an eigenvalue problem, I'm not sure how to work > around this issue. > > To illustrate the point, suppose that the linear system is something like > > du/dt = -i*A*u > > where A is a constant matrix. Formally, we can write the solution as > > u = u_0 exp(-i*A*t) > > and the eigenvalues that is going to grow the fastest is the one that has > the largest imaginary part. > > When we assume the ODE has a normal mode decomposition, exp(-i*lambda*t), > we get that the eigenvalue system becomes real, > > lambda*u= A*u > > but the eigenvalue that corresponds to the largest growth rates is the > imaginary part of the system. > > We could try a solution of the form exp(sigma*t), where the most unstable > mode is real (so we can use a real target) but then the system is > imaginary. A bit of a catch 22 as far as I can tell. > > There are two options that I think *should* work. > > 1) Set up the system in firedrake/fenics with real matrices, export the > matrices into PETSc format, and then learn how to compute the eigenvalues > directly in SLEPc. I imagine there is a steep learning curve but if I > build my only SLEPc with complex numbers then that *should* work. > > 2) Alternatively, I could propagate the identity matrix forward for some > time (how long is a bit of an issue) to find out what the > propagator/fundamental matrix looks like at some later time. Then, I can > compute the eigenvalues of that system, and the most unstable would be real > so it is something that I can do in firedrake/fenics. One problem is that > for a 2D problem, I need to integrate N^2 vectors so that would slow things > down, but at least it should work. > > Does anyone have any thoughts as to what might be easiest, or other clever > ideas? > > One goal that I have is to put together a demo to compute the linear > stability of a problem and then second part that actually solves the > nonlinear evolution. There are a variety of problems that I have in mind. > At the moment I am working on a one-layer QG case but then I'd like to do > the 3D problem that is almost identical. Would make for a very nice demo, > and also help one of my PhD students who is doing this as part of his > thesis right now. Any thoughts would be greatly appreciated. > > Cheers, Francis > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] Energy distribution
by Andrew McRae 03 May '16

03 May '16

Dear Will, Make a Function for outputting; e.g. energy_fn = Function(R) [if you want to output in the function space R] Then do energy_fn.interpolate(dot(u, u) + rho*rho) [for interpolation] or energy_fn.project(dot(u, u) + rho*rho) [for L^2 projection] followed by localenergyfile.write(energy_fn, time=t) Apologies for any typos in the above. Thanks, Andrew On 3 May 2016 at 14:00, William Booker <scwb(a)leeds.ac.uk> wrote: > Dear Firedrakers, > > > For the model, > > > https://github.com/booker213/compressible_acoustic/blob/master/Firedrake/1D… > > I'm interested in finding out how the energy in my system is distributed > at the end of my model. > > My naive approach, below, was to try and create a function, so I could > visualise it. > > > localenergyfile = File('./Results/local_energy.pvd') > E_local = Function(R, name="Local Energy") > E_local = 0.5*(inner(u,u) + rho**2 > localenergyfile.write(E_local , time =t) > > > But this does not work. How would I go about implementing this? > > Regards, > > Will >

1 0

Energy distribution
by William Booker 03 May '16

03 May '16

Dear Firedrakers, For the model, https://github.com/booker213/compressible_acoustic/blob/master/Firedrake/1D… I'm interested in finding out how the energy in my system is distributed at the end of my model. My naive approach, below, was to try and create a function, so I could visualise it. localenergyfile = File('./Results/local_energy.pvd') E_local = Function(R, name="Local Energy") E_local = 0.5*(inner(u,u) + rho**2 localenergyfile.write(E_local , time =t) But this does not work. How would I go about implementing this? Regards, Will

1 0

Converting FENICS code to Firedrake (iterations over mesh cells)
by George Ovchinnikov 03 May '16

03 May '16

Hi, I have the following code, which solves Poisson problem and prints solution in a cell and it neighbours also building cells adjacency matrix which is needed later on. The solving part maps to firedrake just fine, but I'm unsure about what to do with loops over cells and their neighbours. Here the code: u = TrialFunction(V) v = TestFunction(V) f = Constant(1) a = inner(nabla_grad(u), nabla_grad(v))*dx L = f*v*dx # Compute solution u = Function(V) solve(a == L, u, bcs) tdim = mesh.topology().dim() mesh.init(tdim - 1, tdim) cell_neighbors = {cell.index(): sum((filter(lambda ci: ci != cell.index(), facet.entities(tdim)) for facet in facets(cell)), []) for cell in cells(mesh)} print cell_neighbors print mesh.cells() for cell in cells(mesh): point = cell.midpoint() #print cell.index(), list(cells(mesh))[cell_neighbors[cell.index()][0]] nhs = [u(Cell(mesh, cx).midpoint()) for cx in cell_neighbors[cell.index()]] print u(point), nhs Sincerely, George

3 2

Re: [firedrake] Converting FENICS code to Firedrake (iterations over mesh cells)
by David Ham 03 May '16

03 May '16

Hi all, I believe the simplest way to find the cells which share a facet is: mesh.interior_facets.facet_cell Cheers, David On Sat, 30 Apr 2016 at 19:09 Homolya, Miklós <m.homolya14(a)imperial.ac.uk> wrote: > > Hi, > > To evaluate a field at the midpoint of all cells, you can just interpolate > onto a DG0 field. For example: > > midpoints = interpolate(u, FunctionSpace(mesh, 'DG', 0)) > > Details: http://firedrakeproject.org/interpolation.html#ufl-expressions > > Probably the simplest way to find cells that share a facet is to create a > RT1 (RTCF1 on non-simplices) function space on the mesh, and inspect its > cell_node_list attribute. It will be a #cells x #facets/cell matrix, if any > two row contains the same number, that means they share a facet. The > entries in midpoints.dat.data are ordered the same way as the rows of > cell_node_list. > > Regards, > Miklós > > ________________________________________ > From: firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of George Ovchinnikov <lives9(a)gmail.com> > Sent: 30 April 2016 15:54:10 > To: firedrake > Subject: [firedrake] Converting FENICS code to Firedrake (iterations over > mesh cells) > > Hi, > > I have the following code, which solves Poisson problem and > prints solution in a cell and it neighbours also building cells > adjacency matrix which is needed later on. > The solving part maps to firedrake just fine, but I'm unsure about what > to do with loops over cells and their neighbours. > > Here the code: > > u = TrialFunction(V) > v = TestFunction(V) > f = Constant(1) > a = inner(nabla_grad(u), nabla_grad(v))*dx > L = f*v*dx > # Compute solution > u = Function(V) > solve(a == L, u, bcs) > > tdim = mesh.topology().dim() > mesh.init(tdim - 1, tdim) > > cell_neighbors = {cell.index(): sum((filter(lambda ci: ci != cell.index(), > facet.entities(tdim)) > for facet in facets(cell)), []) > for cell in cells(mesh)} > print cell_neighbors > print mesh.cells() > > for cell in cells(mesh): > point = cell.midpoint() > #print cell.index(), > list(cells(mesh))[cell_neighbors[cell.index()][0]] > nhs = [u(Cell(mesh, cx).midpoint()) for cx in > cell_neighbors[cell.index()]] > > print u(point), nhs > > Sincerely, > George > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

1 0

Re: [firedrake] REXI time stepping in Firedrake
by David Ham 03 May '16

03 May '16

Hi Miklos, I suspect there is a little more to it than this. For example, dot products now require the complex conjugate to be taken of the first argument. Cheers, David On Sun, 1 May 2016 at 19:16 Homolya, Miklós <m.homolya14(a)imperial.ac.uk> wrote: > it probably wouldn't be a fatally hard amount of work to put this through > the form compiler. > > > In theory it is the same as PETSc, you just need to change a typedef. (not > tested though) > > > On 1 May 2016, at 12:52, David Ham <David.Ham(a)imperial.ac.uk> wrote: > > It's a little more subtle than this. PyOP2 does support complex and it > probably wouldn't be a fatally hard amount of work to put this through the > form compiler. The real issue is actually PETSc, because we *really* need > runtime switchable complex and that's an enormous amount of work in PETSc. > The PETSc guys know what to do, but the resource implications are > problematic. > > Regards, > > David > > On Sun, 1 May 2016 at 09:17 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > >> To rephrase that, the bigger pain is that Firedrake/FEniCS don't natively >> support complex data types. >> >> PETSc *already does* support complex data types, but with the annoyance >> that it has to be compiled separately for complex mode. >> >> On 1 May 2016 at 08:34, Colin Cotter <colin.cotter(a)imperial.ac.uk> wrote: >> >>> Hi Francis, >>> This was a specific thing about the formulation of rational >>> approximation of exponentials in the following paper: >>> http://arxiv.org/abs/1012.3196 >>> >>> Complex numbers are a pain because PETSc must be separately compiled >>> with complex data types. A possible workaround in current Firedrake/Fenics >>> is to have separate variables for real and imaginary parts and reinterpret >>> the problem as a real one of twice the size (although this limits the >>> solvers that can be used). >>> all the best >>> --cjc >>> >>> On 30 April 2016 at 21:14, Francis Poulin <fpoulin(a)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? >>>> >>>> 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(a)uwaterloo.ca >>>> Web: https://uwaterloo.ca/poulin-research-group/ >>>> Telephone: +1 519 888 4567 x32637 >>>> >>>> ------------------------------ >>>> *From:* firedrake-bounces(a)imperial.ac.uk [ >>>> firedrake-bounces(a)imperial.ac.uk] on behalf of Colin Cotter [ >>>> colin.cotter(a)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(a)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(a)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 >>> >>> >>> >> _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

1 0

Inconsistency between bilinear and residual solve
by Fryderyk Wilczynski 03 May '16

03 May '16

Hello all, I'm solving a system of pdes that are linear in (n+1) terms (in time disretised system). When I write the equations in Firedrake as a bilinear form and solve using: ... solve(a == L, out_u, bcs=[bc0, bc1, bc2, bc3], solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'}) I get a stable (energy dissipative) scheme, everything seems to be working fine. If I reformulate the equations so that I am solving the residual form: ... solve(F == 0, u, bcs=[bc0, bc1, bc2, bc3],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) all other things remainimg equal; the scheme breaks down after a number of timesteps, and throws an exception: Traceback (most recent call last): File "NSCH_residual_form.py", line 149, in <module> solve(F(u,u0,u,u0,W,n,myparameters) == 0, u, bcs=[bc0, bc1, bc2, bc3],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 164, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve solving_utils.check_snes_convergence(self.snes) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence %s""" % (snes.getIterationNumber(), msg)) RuntimeError: Nonlinear solve failed to converge after 17 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH Additionally, the results between the two methods are inconsistent (particularly the min, max ranges). What can I do to fix that? The scripts and pdf containing weak forms can be found on bitbucket: https://bitbucket.org/fryderyk216/navierstokes/src/9f0814f99fe23fa57ead8701… Best wishes, Fryderyk

2 2

Re: [firedrake] Inconsistency between bilinear and residual solve
by Andrew McRae 02 May '16

02 May '16

Dear Fryderyk, The usual advice here would be to build a "minimal failing example". 1) Can you reproduce the error on a 2x2 UnitSquareMesh? Then use that, not 100x100. 2) Run *one* timestep (no need for 100) and verify that the solutions are different (e.g. calculate some diagnostic). You should probably calculate that diagnostic on the initial conditions too, in case they weren't the same... 3) Forget about the full Navier--Stokes + Cahn--Hilliard problem. Throw away everything except one term on each side of the equation. Do 2). If the diagnostics were the same, add a term back in. Etc. Andrew On 2 May 2016 at 10:09, Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > Update: > > > In the residual form case I tried: > > > try: > solve(F(u,u0,u,u0,W,n,myparameters) == 0, u, bcs=[bc0, bc1, bc2, > bc3],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': > 'mumps'} ) # > except RuntimeError: > pass > > that works in the sense that the exception is no longer a problem. However > that slight inconsistency in results is still there. Is that to be expected > because of the different solvers? If that is the case, how can one minimise > that difference? > > > Thanks, > > Fryderyk > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Fryderyk Wilczynski <scfw(a)leeds.ac.uk> > *Sent:* 02 May 2016 09:44:16 > *To:* firedrake(a)imperial.ac.uk > *Subject:* [firedrake] Inconsistency between bilinear and residual solve > > > Hello all, > > > I'm solving a system of pdes that are linear in (n+1) terms (in time > disretised system). > > When I write the equations in Firedrake as a bilinear form and solve using: > > ... > > solve(a == L, out_u, bcs=[bc0, bc1, bc2, bc3], > solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': > 'mumps'}) > > > I get a stable (energy dissipative) scheme, everything seems to be working > fine. > > > If I reformulate the equations so that I am solving the residual form: > > ... > > solve(F == 0, u, bcs=[bc0, bc1, bc2, bc3],solver_parameters={'pc_type': > 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) > > > all other things remainimg equal; the scheme breaks down after a number of > timesteps, and throws an exception: > > Traceback (most recent call last): > File "NSCH_residual_form.py", line 149, in <module> > solve(F(u,u0,u,u0,W,n,myparameters) == 0, u, bcs=[bc0, bc1, bc2, > bc3],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': > 'mumps'} ) # > File > "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", > line 120, in solve > _solve_varproblem(*args, **kwargs) > File > "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", > line 164, in _solve_varproblem > solver.solve() > File "<decorator-gen-295>", line 2, in solve > File > "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", > line 203, in wrapper > return f(*args, **kwargs) > File > "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", > line 190, in solve > solving_utils.check_snes_convergence(self.snes) > File > "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", > line 62, in check_snes_convergence > %s""" % (snes.getIterationNumber(), msg)) > RuntimeError: Nonlinear solve failed to converge after 17 nonlinear > iterations. > Reason: > DIVERGED_LINE_SEARCH > > Additionally, the results between the two methods are inconsistent > (particularly the min, max ranges). > What can I do to fix that? > > The scripts and pdf containing weak forms can be found on bitbucket: > > > https://bitbucket.org/fryderyk216/navierstokes/src/9f0814f99fe23fa57ead8701… > > > Best wishes, > > Fryderyk > > >

1 0

Re: [firedrake] REXI time stepping in Firedrake
by David Ham 01 May '16

01 May '16

It's a little more subtle than this. PyOP2 does support complex and it probably wouldn't be a fatally hard amount of work to put this through the form compiler. The real issue is actually PETSc, because we *really* need runtime switchable complex and that's an enormous amount of work in PETSc. The PETSc guys know what to do, but the resource implications are problematic. Regards, David On Sun, 1 May 2016 at 09:17 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > To rephrase that, the bigger pain is that Firedrake/FEniCS don't natively > support complex data types. > > PETSc *already does* support complex data types, but with the annoyance > that it has to be compiled separately for complex mode. > > On 1 May 2016 at 08:34, Colin Cotter <colin.cotter(a)imperial.ac.uk> wrote: > >> Hi Francis, >> This was a specific thing about the formulation of rational >> approximation of exponentials in the following paper: >> http://arxiv.org/abs/1012.3196 >> >> Complex numbers are a pain because PETSc must be separately compiled with >> complex data types. A possible workaround in current Firedrake/Fenics is to >> have separate variables for real and imaginary parts and reinterpret the >> problem as a real one of twice the size (although this limits the solvers >> that can be used). >> all the best >> --cjc >> >> On 30 April 2016 at 21:14, Francis Poulin <fpoulin(a)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? >>> >>> 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(a)uwaterloo.ca >>> Web: https://uwaterloo.ca/poulin-research-group/ >>> Telephone: +1 519 888 4567 x32637 >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk [ >>> firedrake-bounces(a)imperial.ac.uk] on behalf of Colin Cotter [ >>> colin.cotter(a)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(a)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(a)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 >> >> >> >

2 1

linear stability problem in firedrake
by Francis Poulin 01 May '16

01 May '16

Dear Lawrence, David and Colin, Thanks a lot for the helpful information. I changed the subject since I thought it might be helpful to have it match the topic of discussion. Now that I have an example on how to do this in firedrake I think that I could probably adapt my QG stability code (which would go very well with the nonlinear QG demo that I posted a while back). If there is an option to use a direct eigenvalue method (lapack) then that would work for 1D eigenvalue problems. Since the most unstable mode corresponds to the eigenvalue with the largest imaginary part, and I don't suspect that firedrake allows for imaginary targets, I don't imagine it would converge. That's the problem I have in FEniCS. Colin: Thanks for the advice. For solving linear systems I know how to decompose the problem into a real and imaginary part and that would be easy enough to do. However, for an eigenvalue problem, I'm not sure how to work around this issue. To illustrate the point, suppose that the linear system is something like du/dt = -i*A*u where A is a constant matrix. Formally, we can write the solution as u = u_0 exp(-i*A*t) and the eigenvalues that is going to grow the fastest is the one that has the largest imaginary part. When we assume the ODE has a normal mode decomposition, exp(-i*lambda*t), we get that the eigenvalue system becomes real, lambda*u= A*u but the eigenvalue that corresponds to the largest growth rates is the imaginary part of the system. We could try a solution of the form exp(sigma*t), where the most unstable mode is real (so we can use a real target) but then the system is imaginary. A bit of a catch 22 as far as I can tell. There are two options that I think *should* work. 1) Set up the system in firedrake/fenics with real matrices, export the matrices into PETSc format, and then learn how to compute the eigenvalues directly in SLEPc. I imagine there is a steep learning curve but if I build my only SLEPc with complex numbers then that *should* work. 2) Alternatively, I could propagate the identity matrix forward for some time (how long is a bit of an issue) to find out what the propagator/fundamental matrix looks like at some later time. Then, I can compute the eigenvalues of that system, and the most unstable would be real so it is something that I can do in firedrake/fenics. One problem is that for a 2D problem, I need to integrate N^2 vectors so that would slow things down, but at least it should work. Does anyone have any thoughts as to what might be easiest, or other clever ideas? One goal that I have is to put together a demo to compute the linear stability of a problem and then second part that actually solves the nonlinear evolution. There are a variety of problems that I have in mind. At the moment I am working on a one-layer QG case but then I'd like to do the 3D problem that is almost identical. Would make for a very nice demo, and also help one of my PhD students who is doing this as part of his thesis right now. Any thoughts would be greatly appreciated. Cheers, Francis

1 0