April 2020 - firedrake - imperial.ac.uk

very slow code when using a lot of functions belonging to Real spaces
by Roberto Federico Ausas 04 May '20

04 May '20

Hello, I have recently posted some questions related to the use of Real function spaces for some FSI problems I'm trying to solve. The problem I'm facing now is that the solving phase is consuming absurd amount of time. It does not seem to be the PETSc resolution itself, but the assembly phase, since, once the KSP iterations begin they proceed reasonably fast. The piece of code I believe is relevant is: #------------------------------------------------------------------------ . . . EU = VectorElement("Lagrange", 'triangle', orderU) EP = FiniteElement("Lagrange", 'triangle', orderP) EL = FiniteElement("Real", 'triangle', 0) ES = [] ndofs = 20 for i in range(ndofs): ES.append(FiniteElement("Real", 'triangle', 0)) EM = MixedElement([EU,EP,EL,ES]) W_f = FunctionSpace(mesh, EM) # mesh has a few hundred nodes sol = Function(W_f) auxs = split(sol) u = auxs[0] p = auxs[1] l = auxs[2] s = [] for i in range(ndofs): s.append(auxs[i+3]) auxT = TestFunctions(W_f) v = auxT[0] q = auxT[1] r = auxT[2] d = [] for i in range(ndofs): d.append(auxT[i+3]) uT = u for i in range(ndofs): uT = uT + s[i]*phifluid[i] # ( phifluid[:] are functions belonging to W_f that are zero everywhere except on some boundary nodes ) print(' |> - Fluid part') Res = momentum_residual(uT, p, l, s, v, mesh) # Stokes operator Res += mass_residual(uT, p, l, s, q, r, mesh) # Incompressibility restriction print(' |> - Beam part') aream1 = 1/assemble(Constant(1) * dx(mesh)) for i in range(ndofs): Res += momentum_residual(uT, p, l, s, d[i]*phifluid[i], mesh) force_b = Constant(force_beam(vphi, dphibeam[i])) Res -= d[i]*force_b*aream1*dx # Until this point the code runs reasonably fast # Here starts the bottleneck idsplit = '2' for i in range(ndofs): idsplit = idsplit+','+str(i+3) # idsplit ends up being 2,3,4,5,... sp = { "mat_type": "matfree", "snes_type": "newtonls", "snes_monitor": None, "snes_converged_reason": None, "snes_linesearch_type": "basic", "snes_max_it": 8, "snes_rtol": 1e-10, "ksp_type": "fgmres", "ksp_monitor_true_residual": None, "ksp_max_it": 10, "pc_type": "fieldsplit", "pc_fieldsplit_type": "schur", "pc_fieldsplit_schur_fact_type": "full", "pc_fieldsplit_0_fields": "0,1", "pc_fieldsplit_1_fields": idsplit, "fieldsplit_0_ksp_type": "preonly", "fieldsplit_0_pc_type": "python", "fieldsplit_0_ksp_monitor_true_residual": None, "fieldsplit_0_pc_python_type": "firedrake.AssembledPC", "fieldsplit_0_assembled_pc_type": "lu", "fieldsplit_0_assembled_pc_factor_mat_solver_type": "mumps", "fieldsplit_0_assembled_mat_mumps_icntl_14": 200, "mat_mumps_icntl_14": 200, "fieldsplit_1_ksp_type": "gmres", "fieldsplit_1_ksp_monitor_true_residual": None, "fieldsplit_1_ksp_max_it": ndofs+1, "fieldsplit_1_ksp_convergence_test": "skip", "fieldsplit_1_pc_type": "none" } solve(Res == 0, sol, bcs=bcs, J=Jac, solver_parameters=sp) #------------------------------------------------------------------------ Any tip on what could be going on would be welcome. Regards -- Roberto F. Ausas Instituto de Ciências Matemáticas e de Computação - ICMC/USP São Carlos, SP - Brasil *http://www.lmacc.icmc.usp.br/~ausas* <http://www.lmacc.icmc.usp.br/~ausas>

4 7

Nelder fellow lecture course on Multigrid Methods by Luke Olson
by Cotter, Colin J 29 Apr '20

29 Apr '20

Dear all, Prof Luke Olson from University of Illinois at Urbana-Champaign is a leading expert in algebraic multigrid methods, visiting the department for six months on a Nelder Fellowship. Unfortunately this was interrupted by the covid19 pandemic and he and his family have returned to Illinois. However, Luke is still able to give his scheduled course on multigrid methods in a virtual format. This will comprise 4x2 hours sessions (Luke anticipates that he won't use the whole 2 hours, probably something just over an hour per session with discussion afterwards). These will be: Wednesday 6th May 4pm (BST) Friday 8th May 4pm (BST) Monday 11th May 4pm (BST) Friday 15th May 4pm (BST) There is a sign-up sheet for the course here: https://docs.google.com/spreadsheets/d/13Z5VSijFgn3SGng9YJMfR8nRF8ATs4wpy7s… and course materials will be uploaded to this git repository: https://github.com/lukeolson/imperial-multigrid The course will be a mixture of analysis and practical aspects. Luke is a great communicator and teacher so I hope that you will be able to join this excellent course. See you there! all the best --Colin Professor Colin Cotter (he/him) Department of Mathematics 755, Huxley Building Imperial College London South Kensington Campus United Kingdom of Great Britain and Northern Ireland +44 2075943468

1 0

warning about "Creating new TransferManager to transfer data to coarse grids"
by Ed Bueler 20 Apr '20

20 Apr '20

Dear Firedrake -- When I run the demo code https://www.firedrakeproject.org/demos/geometric_multigrid.py.html I get output with the warning firedrake:WARNING Creating new TransferManager to transfer data to coarse grids firedrake:WARNING This might be slow (you probably want to save it on an appctx) The same warning appears for my own codes, which often use -pc_type mg. This is new behavior for me, i.e. as of a firedrake update about 2 months ago. How should I update this demo to avoid the warning message? Is there an example where it is added to an appctx? I'm curious about the cause of the warning, because I would not want mesh-specific info in a PETSc appctx, generally, I think ... Thanks! Ed -- Ed Bueler Dept of Mathematics and Statistics University of Alaska Fairbanks Fairbanks, AK 99775-6660 306C Chapman

4 5

Re: [firedrake] Issue solving problem involving function in R space
by Roberto Federico Ausas 17 Apr '20

17 Apr '20

Hello, Thanks for your answer. The solution you provided worked. Effectively, I was using Fenics, but then I needed to switch to firedrake because of the possibility of using Hermite elements to solve the structure part (a Kirchhoff beam). In the FSI interaction scheme I'm thinking, I will probably need a few dozens of variables belonging to R spaces (these essentially are the dofs associated with the beam problem). I hope the strategy you provided will work as well in such a case. Another difficulty I will face is that in order to construct the residual of my problem, I will need to integrate over the fluid and the beam meshes, the coupling being done through these Real space variables. A fully monolithic approach is probably out of question. Regards On Fri, Apr 17, 2020 at 8:00 AM <firedrake-request(a)imperial.ac.uk> wrote: > Send firedrake mailing list submissions to > firedrake(a)imperial.ac.uk > > To subscribe or unsubscribe via the World Wide Web, visit > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > or, via email, send a message with subject or body 'help' to > firedrake-request(a)imperial.ac.uk > > You can reach the person managing the list at > firedrake-owner(a)imperial.ac.uk > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of firedrake digest..." > > > Today's Topics: > > 1. Issue solving problem involving function in R space > (Roberto Federico Ausas) > 2. Re: Issue solving problem involving function in R space > (Patrick Farrell) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Fri, 17 Apr 2020 00:22:02 -0300 > From: Roberto Federico Ausas <rfausas(a)icmc.usp.br> > To: firedrake(a)imperial.ac.uk > Subject: [firedrake] Issue solving problem involving function in R > space > Message-ID: > < > CAH6Ec4msgbi57NxOuLHB-W3LZm1yzhpa86DORwdobSWr2KmEUg(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > Hello, > > I need to solve a problem involving functions belonging to an R space. > > For robustness, I need to use direct solvers for the linear part, but I'm > getting the following: > > ValueError: Monolithic matrix assembly is not supported for systems with > R-space blocks. > > Is there any way of doing this? > Any help would be appreciated > Regards > > -- > Roberto F. Ausas > Instituto de Ci?ncias Matem?ticas e de Computa??o - ICMC/USP > S?o Carlos, SP - Brasil > *http://www.lmacc.icmc.usp.br/~ausas* <http://www.lmacc.icmc.usp.br/~ausas > > >

2 1

Issue solving problem involving function in R space
by Roberto Federico Ausas 17 Apr '20

17 Apr '20

Hello, I need to solve a problem involving functions belonging to an R space. For robustness, I need to use direct solvers for the linear part, but I'm getting the following: ValueError: Monolithic matrix assembly is not supported for systems with R-space blocks. Is there any way of doing this? Any help would be appreciated Regards -- Roberto F. Ausas Instituto de Ciências Matemáticas e de Computação - ICMC/USP São Carlos, SP - Brasil *http://www.lmacc.icmc.usp.br/~ausas* <http://www.lmacc.icmc.usp.br/~ausas>

2 1

Re: [firedrake] SymPy to UFL
by Ham, David A 16 Apr '20

16 Apr '20

Hi Andrew, I think you can just use the `.subs()` method on a sympy expression and provide mappings from Sympy symbols to Firedrake objects such as Functions, TestFunctions, or TrialFunctions. Alternatively, you can write your manufactured solution in UFL directly and just use UFL’s differentiation operations to compute the required source term, bypassing sympy altogether. Cheers, David From: <firedrake-bounces(a)imperial.ac.uk> on behalf of Andrew Hicks <ahick17(a)lsu.edu> Date: Thursday, 16 April 2020 at 01:35 To: firedrake <firedrake(a)imperial.ac.uk> Subject: [firedrake] SymPy to UFL Hi all, I am currently trying to work on projects such as manufacturing solutions for PDE’s. I will mainly be using SymPy for this. My question is, is there a way to take a SymPy expression in my code and apply some function to it to convert it to UFL format, so that I can then use it in my firedrake functions? Thanks, Andrew Hicks

1 0

SymPy to UFL
by Andrew Hicks 16 Apr '20

16 Apr '20

Hi all, I am currently trying to work on projects such as manufacturing solutions for PDE’s. I will mainly be using SymPy for this. My question is, is there a way to take a SymPy expression in my code and apply some function to it to convert it to UFL format, so that I can then use it in my firedrake functions? Thanks, Andrew Hicks

1 0

Changing mesh coordinates restrictions
by Thomas Ranner 15 Apr '20

15 Apr '20

Hi, I am interested in using firedrake to investigate using different "higher order" meshes. By this I mean replacing mesh coordinates with piecewise polynomial function from a particular function space. A typical example would be a higher order Lagrange space but I am also interested in investigating other function spaces too (e.g. Hermite). According to the documentation I can find online this should be possible (https://www.firedrakeproject.org/mesh-coordinates.html) I've written what I think is about the simplest test case - solving a Poisson equation on a unit square (taken from https://www.firedrakeproject.org/demos/helmholtz.py.html) and used L2 projection to define the new coordinates of the identity map (I understand that interpolation isn't always defined for different spaces). Test case simple.py attached. I can run successfully with a Lagrange space of different orders: ``` (firedrake) Singularity> python simple.py Lagrange 3 Lagrange (3) mesh with 200 elements created successfully. problem on mesh defined successfully problem solved successfully error 0.053347460759624356 ``` But the equivalent code for a Hermite basis doesn't work (and gives the long errors attached). ``` (firedrake) Singularity> python simple.py Hermite &> hermite.out ``` Can anyone help with this? Is there any reason why this code shouldn't run? Thanks, Tom

3 3

The UK's Leading Supplier of Face Masks, Latex Gloves, Hand Sanitizer
by Adam Martin 15 Apr '20

15 Apr '20

[DEFAULT_HEADER][ We still have stock of Face Masks! While stocks last! Protect yourself and your employees now! CLICK HERE! ]( http://r.indigosmtp.com/mk/cl/f/7Vs0AnXPbVMvlb0RH3Blmnx-yswF-7LAsxTakZCltQG… ) Call 01268 768 768 NOW! PROTECTU have been supplying PPE for over 25 years We have a full range of Face Masks, Hand Sanitizer, Latex and Nitrile Gloves and more in stock now Which mask is best for Coronavirus? We recommend the following Masks FFP1 GOOD (filters 80% particles) FFP2 (N95) EXCELLENT (filters 94% of particles) FFP3 (N99) THE BEST (filters 99% of particles) [ CLICK HERE TO BUY NOW ]( http://r.indigosmtp.com/mk/cl/f/LBYi9i3E_zkjqsTdgM0Ltbe5MsYjxCEdIsC6WtNLKuQ… ) © 2020 PROTECT U - UK'S LEADING SUPPLIER OF PPE We hope you dont mind receiving the odd email from us This email was sent to firedrake(a)imperial.ac.uk You received this email because you are registered with the Indigo Group [ Unsubscribe here ]( http://r.indigosmtp.com/mk/un/g4A1HlAo2vIzXA_eISZdqoEWgw4AVcFA8IT0MKIW1v-xr… )

1 0

Re: [firedrake] demo update queries
by Ham, David A 15 Apr '20

15 Apr '20

Hi Onno, See below for answers. From: <firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove <O.Bokhove(a)leeds.ac.uk> Date: Tuesday, 14 April 2020 at 18:45 To: firedrake <firedrake(a)imperial.ac.uk> Subject: [firedrake] demo update queries Dear Firedrakers, Some queries: (a1) I have been updating one of the tutorials by adapting the python script of the BL-demo; how can I test its page, and where do I submit it as candidate update? The demos are tested using the following command: py.test tests/demos/test_demos_run.py In order to contribute a new or modified demo (thanks!), you need to make a branch with your changes, push it to GitHub and make a pull request. You need push access to the Firedrake repository on GitHub, which I will give you if you email me your github username. (a2) In essence I updated the demo so that it has reduced output and reproduces the Figure 3 test of an associated publication, which figure I included as target for reproduction by any user, because then the user would know whether the outcome is right; also adding a small python script to visualise the energy (lower panel of Figure 3 attached). Since Figure 3a and 3b contain line plots at multiple times: rather than using an animation in Paraview: how can one with "Plot over line" keep multiple line plots at different time levels in one plot (since this would literally reproduce Figure 3a and 3b)? If you pass any of the plotting commands the keyword argument “axes=a” where “a” is a matplotlib axes object, then the plot will happen on those axes. This can be used to overplot. Also the old demo had some old expired links which I fixed. (b) Finally: with h_n0 = Function(V) on a 2D mesh, how can I find the maximum value on the mesh? This is in preparation for a potential flow test with a nonlinear, 8.4-fold amplification. I now use/added: ... from firedrake.petsc import PETS ... with h_n0.dat.vec_ro as hh: L_inf_err = hh.max()[1] ... but that seems a bit overkill? You do have to do this or the PyOP2 equivalent (using .dat.data_ro) because there isn’t in general a straightforward way of finding the maximum value of a Function from its DoFs. There is for DG0 and CG1, but in other cases the function maximum may not occur at a node. Cheers, David MPI next!?! Thank you and kind regards, Onno

1 0