May 2015 - firedrake - imperial.ac.uk

Re: [firedrake] Mapping nodal values from parent to extruded mesh
by David Ham 04 Jun '15

04 Jun '15

Hi Tuomas, Concerning the sets of spaces which this will work for, I think that the current implementation probably works for any horizontal space but only for vertical CG. One could just try relaxing the test in that case. The reason for the restriction to vertical CG is that fs.bt_masks uses the topological association of nodes with mesh entities in order to work out which nodes are on the top or bottom of the cell. In order to allow DG in the vertical, it would be necessary to support the geometric definition (ie which basis functions do not vanish on the top/bottom). If you look at functionspace.py:78 you can see where the bottom and top masks are generated. This uses entity_closure_dofs() from FIAT. In order to use the geometric definition of dofs one would need to support using facet_support_dofs() (which is how BC maps are set up at functionspace.py:375). Currently FIAT TensorFiniteElement objects do support entity_closure_dofs() but nobody has done the legwork to get them to support facet_support_dofs. Concerning documentation, COFFEE is deplorably underdocumented (although the author is on this list so maybe this will change ;). The actual available AST nodes are only "documented" by reading the source: https://github.com/coneoproject/COFFEE/blob/master/coffee/base.py However, you are fundamentally writing a PyOP2 kernel, and the C api for that *is* documented at: http://op2.github.io/PyOP2/kernels.html so that hopefully helps somewhat. Otherwise, keep asking here or in IRC (IRC tends to get faster responses, at least between about 10am and 8pm UK time). BTW, any chance of seeing you at FEniCS '15? Cheers, David On 10 April 2015 at 02:44, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi all, > > A couple of questions regarding those extruded mesh -> parent mesh copy > operations, > > A while back I got the reverse 2d->3d copy working with a hard-coded pyop2 > kernel (like in mesh extrusion). I'm not familiar with COFFEE syntax, is > there documentation/examples somewhere? > > The 3d->2d copy routine in extrusion_extraction branch checks that the > function space is CG. Does this method easily generalize to other spaces? > DGxCG prisms seem to work OK. I'm also interested in using DGxDG and > RTxCG/DG spaces. > > > Cheers, > > Tuomas > > > On 03/05/2015 10:24 AM, Tuomas Karna wrote: > > Thanks David, > > This is great, I'll try doing the reverse operation. > > - Tuomas > > On 03/05/2015 08:20 AM, David Ham wrote: > > Hi Tuomas, > > This wasn't there this morning but I've implemented one of the cases > (pulling out the top and bottom maps). The result is in the > extrusion_extraction branches of both PyOP2 and Firedrake. The 2d->3d > operation would be rather similar except that you'd have to interrogate the > fiat_element on the extruded Function in order to determine which extruded > nodes are "on top of" the 2d nodes you have. Unfortunately I won't have > time to do that one soon (huge amounts to do before SIAM CSE next week) but > feel free to have a try and complain when it doesn't work! > > Cheers, > > David > > On 4 March 2015 at 22:23, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> Hi all, >> >> I'd need to copy nodal values between fields on parent and extruded >> meshes. For example, copy 2d->3d (constant over vertical) or 3d->2d >> (extract surface/bottom level). The horizontal function space is the >> same. Is there an easy way to do this? Is there a map of extruded nodes >> somewhere? >> >> Thanks, >> >> Tuomas >> >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > > > > -- > Dr David Ham > Departments of Mathematics and Computing > Imperial College London > > http://www.imperial.ac.uk/people/david.ham > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

4 6

Projection on DG0 space fails on large core counts
by Eike Mueller 03 Jun '15

03 Jun '15

Dear firedrakers, I now re-ran my code on up to 1536 cores on ARCHER, but I get a problem when I try to project an expression onto a DG0 function space on an extruded grid. The full (very large) log is here https://gist.github.com/eikehmueller/83a5fc139e1fedb5306c but as far as I can tell the following crashes: r_p.project(expression,solver_parameters={'ksp_type':'cg','pc_type':'jacobi’}) and here is the relevant part of the trace that I attempted to reconstruct: File "/work/n02/n02/eike/git_workspace/firedrake/firedrake/function.py", line 157, in project return projection.project(b, self, *args, **kwargs) File "/work/n02/n02/eike/git_workspace/firedrake/firedrake/projection.py", line 94, in project […] solving_utils.check_snes_convergence(self.snes) File "/work/n02/n02/eike/git_workspace/firedrake/firedrake/variational_solver.py", line 163, in solve File "/work/n02/n02/eike/git_workspace/PyOP2/pyop2/profiling.py", line 199, in wrapper %s""" % (snes.getIterationNumber(), msg)) File "/work/n02/n02/eike/git_workspace/firedrake/firedrake/solving_utils.py", line 62, in check_snes_convergence return f(*args, **kwargs) RuntimeError: Nonlinear solve failed to converge after 1 nonlinear iterations. It does work fine on smaller processor numbers. Maybe the PETSc integers overflow again, the number of cells is 5242880 x 64 = 335544320 ~ 2^{28}, which is not too far from 2^{32}, but I thought I check in case you’ve seen something similar before. I thought I had managed to run problems of this size in the past (i.e. earlier this year). Thanks, Eike

2 1

EUCLID BoomerAMG smoother on ARCHER
by Eike Mueller 29 May '15

29 May '15

Hi, I just tried to run the PETSc fieldsplit solver with the EUCLID smoother for the BoomerAMG preconditioner on ARCHER, and I get this error message: ============= error stack trace ==================== [18] ERROR: global index 112 not found in map Numbering_dhGlobalToLocal file= Numbering_dh.c line= 175 [18] called from: Factor_dhSolveSetup file= Factor_dh.c line= 555 [18] called from: Euclid_dhSetup file= Euclid_dh.c line= 282 [18] called from: HYPRE_EuclidSetup file= HYPRE_parcsr_Euclid.c line= 278 any idea what might cause this? I can run this without problems on my laptop. Thanks, Eike

2 12

FEniCS '15 registration. Imperial people: yes, you do have to register.
by David Ham 28 May '15

28 May '15

Dear Imperial FEniCS '15 participants (and apologies for spamming the other people on these lists), In response to a couple of enquiries this week, please be aware that Imperial staff and students who are intending to attend FEniCS '15 *do* have to register. The only exception is if you only plan to participate in the poster session and/or the public lecture. If you intend to attend the rest of the workshop then you need to register. Registration is comparatively cheap (given that you get food + beverages for 2.5 days) but it does get more expensive after this weekend! Register online at: http://bit.ly/1J5jIkJ Cheers, David

1 0

Coordinates of nodes
by Colin Cotter 27 May '15

27 May '15

Dear all, What's the easiest way to get global coordinates evaluated at node points for normal evaluation nodes of H(div) spaces? Cheers Cjc

1 0

FEniCS '15 Early registration closes on 31 May
by David Ham 27 May '15

27 May '15

Dear all, This is a gentle reminder that FEniCS '15 early registration closes this coming Sunday. From Monday, registration will cost £200 instead of £150. Note that this is not merely an attempt to fleece you: I need to make catering and other orders so I really need to know fairly firm numbers. If you haven't registered yet, please head to the link below and do so! http://bit.ly/1J5jIkJ Cheers, David

1 0

Additional poster spaces, and poster board size.
by David Ham 22 May '15

22 May '15

Dear all, A few people have emailed me asking if they can submit a late poster abstract. I have now sourced the poster boards and there isn't a problem fitting in a few more posters. I have therefore emailed those people to tell them how to submit an abstract. If there is anyone else out there who would like to put in a poster and who has not yet done so, please email me and I'll make it possible for you to submit an abstract. For those of you who will be bringing posters to the workshop, the poster boards will accommodate an A0 portrait format poster. The boards are perhaps 2m tall and don't have legs, so there is plenty of vertical space, but A0 landscape posters won't fit. Regards, David

1 0

Re: [firedrake] 3d momentum equation with RT spaces
by Andrew McRae 20 May '15

20 May '15

Hi Tuomas, Are you sure you want to use RT1xP1, as opposed to RT1xP0[DG]? This would be the usual element for 'mimetic' discretisations. The code for RT1xP0 would be [possibly with typos...] --------------------------- from firedrake import * mesh2d = UnitSquareMesh(10, 10) layers = 6 mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) Uh_elt = FiniteElement('RT', triangle, 1) Uv_elt = FiniteElement('DG', interval, 0) U_elt = HDiv(OuterProductElement(Uh_elt, Uv_elt)) U = FunctionSpace(mesh, U_elt) P1 = FunctionSpace(mesh, 'CG', 1) # automatically gives the appropriate P1 x P1 element f = Function(P1) f.interpolate(Expression('-x[0]')) sol = Function(U) test = TestFunction(U) # sol + grad(f) = 0 F = inner(sol, test)*dx + dot(grad(f), test)*dx solve(F == 0, sol) print sol.dat.data.min(), sol.dat.data.max() --------------------------- FIAT will throw an error if you try to do this with RT1xP1, however (though it will let you build a curl-conforming 'edge' element by using HCurl). If you genuinely need RT1xP1 then I can probably turn off some FIAT validation so that it can be used. I think it's more likely that you want RT1xP0, however, and a corresponding vertical velocity element of P0xP1, HDiv'd. On 20 May 2015 at 02:22, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi all, > > I'm implementing 3D momentum equation for horizontal velocity components > (u,v) that live in RT1xP1 space. So the space is 3D but the field only has > 2 components. I'm getting some errors from FFC related to missing > quadrature rules: > > File > "/home/tuomas/.local/lib/python2.7/site-packages/FFC-1.6.0.dev0-py2.7-linux-x86_64.egg/ffc/quadrature/quadraturetransformerbase.py", > line 1031, in _create_mapping_basis > name, non_zeros, zeros, ones = self.name_map[name] > KeyError: 'FE1_C2' > > If I interpret that correctly, FE1_C2 is a quad rule for the third > component of the vector space, which doesn't exist in this case. If I > replace the RT space with DG (dim=2) everything works. Is this a FFC issue > or am I doing this wrong? > > - Tuomas > > --- > > from firedrake import * > > mesh2d = UnitSquareMesh(10, 10) > layers = 6 > mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) > > #U = VectorFunctionSpace(mesh, 'DG', 1, vfamily='CG', vdegree=1, dim=2) > U = FunctionSpace(mesh, 'RT', 1, vfamily='CG', vdegree=1) > P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) > > f = Function(P1) > f.interpolate(Expression('-x[0]')) > sol = Function(U) > > test = TestFunction(U) > > # sol + grad(f) = 0 > F = inner(sol, test)*dx + (grad(f)[0]*test[0] + grad(f)[1]*test[1])*dx > > solve(F == 0, sol) > print sol.dat.data.min(), sol.dat.data.max() > >

2 1

3d momentum equation with RT spaces
by Tuomas Karna 20 May '15

20 May '15

Hi all, I'm implementing 3D momentum equation for horizontal velocity components (u,v) that live in RT1xP1 space. So the space is 3D but the field only has 2 components. I'm getting some errors from FFC related to missing quadrature rules: Konsole output File "/home/tuomas/.local/lib/python2.7/site-packages/FFC-1.6.0.dev0-py2.7-linux-x86_64.egg/ffc/quadrature/quadraturetransformerbase.py", line 1031, in _create_mapping_basis name, non_zeros, zeros, ones = self.name_map[name] KeyError: 'FE1_C2' If I interpret that correctly, FE1_C2 is a quad rule for the third component of the vector space, which doesn't exist in this case. If I replace the RT space with DG (dim=2) everything works. Is this a FFC issue or am I doing this wrong? - Tuomas --- from firedrake import * mesh2d = UnitSquareMesh(10, 10) layers = 6 mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) #U = VectorFunctionSpace(mesh, 'DG', 1, vfamily='CG', vdegree=1, dim=2) U = FunctionSpace(mesh, 'RT', 1, vfamily='CG', vdegree=1) P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) f = Function(P1) f.interpolate(Expression('-x[0]')) sol = Function(U) test = TestFunction(U) # sol + grad(f) = 0 F = inner(sol, test)*dx + (grad(f)[0]*test[0] + grad(f)[1]*test[1])*dx solve(F == 0, sol) print sol.dat.data.min(), sol.dat.data.max()

2 2

Pushing branch to mapdes/petsc in preparation for pull request
by Eike Mueller 08 May '15

08 May '15

Dear firedrakers, I finished making the changes to PETSc which allow setting more complex smoothers (in particular the Euclid ILU(k) smoother) for hypre’s BoomerAMG. With Euclid’s ILU(k) smoother I now also get a converging iteration for higher orders. With my changes it becomes possible to control the options for Euclid (mainly drop tolerance and fill-in level k) via the PETSc options database. For this I cloned mapdes/petsc and created a branch extend_hypre_boomer_amg_interface. When I’m now trying to push this to the remote I get an authentication failure: eikemueller@Eikes-MBP $ git push -u origin extend_hypre_boomeramg_interface -vv Pushing to https://bitbucket.org/mapdes/petsc.git Username for 'https://bitbucket.org': em459 Password for 'https://em459@bitbucket.org': fatal: Authentication failed for 'https://bitbucket.org/mapdes/petsc.git/' Is this because I’m not allowed to push to this repository (I appear to have access level ‘Read'), or simply an issue with my connection to bitbucket? I was going to push to mapdes and then create a pull request from there to get this into the PETSc upstream master, this seems easier than me forking off PETSc, making the changes in my fork and then create a pull request from this. Once this has gone into the main PETSc repository, what should I do to get the changes in mapdes’ master branch? Ask one of you to sync mapdes/master with the upstream master? Thanks, Eike

2 5