- firedrake - imperial.ac.uk

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

Pointwise calculations
by Colin Cotter 07 May '15

07 May '15

Dear Firedrakers, One of our favourite testcases, the Galewsky jet, is annoying in that the initial conditions for one of the variables is expressed as a 1-dimensional integral which must be evaluated numerically. In our Fortran code, we used the Fluidity Python interface to evaluate this function pointwise at the node points of the coordinate field, and were thus able to use Numpy quadrature routines to evaluate the integral. What is the best way to do this in Firedrake? all the best --Colin

2 1

Update on using hypre's BoomerAMG
by Eike Mueller 03 May '15

03 May '15

Dear all, having done some runs on ARCHER now, hypre’s BoomerAMG is doing very well at lowest order. On a full node (=24 cores) I get for a problem with 81920 x 64 = 5E6 cells (all times including setup times, which are ~2/3 or the total runtime): matrixfree, geometric multigrid: 11.15 PETSc, GAMG (old results): 13.80 PETSc, BoomerAMG: 8.41 This is for the BoomerAMG standard parameters, I haven’t attempted any tuning, since looking back through my notes from two years ago this will only give maybe 10% more performance (it could, however, make a difference for the scalbility). Note also that I use the default smoother (which is symmetric-SOR/Jacobi), as I don’t think I can use block ILU with BoomerAMG. In the tests we did in the past with hypre we did the same, i.e. use a pointwise smoother in BoomerAMG. Since the predominant coupling is in the vertical direction, AMG will first coarsen in this direction until the problem is isotropic, and then a pointwise smoother works fine. This is very different to what’s done in the tensor-product geometric multigrid, which coarses in the horizontal only and so needs the columnwise smoother. This might also explain why the GAMG (which did use ILU) gives worse performance than BoomerAMG. Unfortunately it doesn’t look good at higher order order. While the matrixfree solver converges as before, once I use BoomerAMG the convergence rate is close to 1.0 I will try tweaking the parameters (there are lots!) to see if that helps. Eike

1 0

unexpected keyword argument 'search'
by Colin Cotter 30 Apr '15

30 Apr '15

Dear Firedrake list, Upon executing the following: from firedrake import * n = 10 mesh = UnitSquareMesh(n,n) V1 = FunctionSpace(mesh,"BDM",2) V2 = FunctionSpace(mesh,"DG",1) W = MixedFunctionSpace((V1,V2)) F = Function(V1).project(Expression(("sin(x[0]*2*pi)*cos(x[1]*2*pi)", "-cos(x[0]*2*pi)*sin(x[1]*2*pi)"))) I get: File "/home/cjc1/firedrake/firedrake/firedrake/function.py", line 157, in project return projection.project(b, self, *args, **kwargs) File "/home/cjc1/firedrake/firedrake/firedrake/projection.py", line 69, in project f.interpolate(v) File "/home/cjc1/firedrake/firedrake/firedrake/function.py", line 221, in interpolate subset) File "/home/cjc1/firedrake/firedrake/firedrake/function.py", line 275, in _interpolate op2.par_loop(*args) File "<string>", line 2, in par_loop File "/home/cjc1/firedrake/PyOP2/pyop2/versioning.py", line 150, in modifies_arguments retval = func(*args, **kwargs) File "/home/cjc1/firedrake/PyOP2/pyop2/op2.py", line 276, in par_loop return backends._BackendSelector._backend.par_loop(kernel, iterset, *args, **kwargs) File "/home/cjc1/firedrake/PyOP2/pyop2/base.py", line 4252, in par_loop return _make_object('ParLoop', kernel, it_space, *args, **kwargs).enqueue() File "/home/cjc1/firedrake/PyOP2/pyop2/backends.py", line 72, in _make_object return _BackendSelector._backend.__dict__[obj](*args, **kwargs) File "<string>", line 2, in __init__ File "/home/cjc1/firedrake/PyOP2/pyop2/utils.py", line 122, in wrapper return f(*args, **kwargs) File "/home/cjc1/firedrake/PyOP2/pyop2/base.py", line 3906, in __init__ ast_info = ast_visit(self._kernel._ast, search=ast.FunDecl) TypeError: visit() got an unexpected keyword argument 'search' If I've done everything correctly, this is with everything on master. What do you think might be wrong? cheers --cjc

2 2

Last day for abstracts!
by David Ham 30 Apr '15

30 Apr '15

Dear all, The flood of abstracts over the last 24 hours says a lot about the power of deadlines. This is the last reminder that everyone who wishes to present a poster or a talk at FEniCS '15 needs to submit an abstract today. Abstracts should be submitted on EasyChair: https://easychair.org/conferences/?conf=fenics15 I look forward to the rest of the last minute rush. David

1 0

Re: [firedrake] Nullspace & VectorSpaceBasis
by Andrew McRae 29 Apr '15

29 Apr '15

Four comments: 1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form! 2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end. 3) Could the equations given not be included in the weak form for the rest of the problem? 4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?.... Andrew On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Dear all, > > I am trying to solve a simple problem (linear water wave equations), > where I have the Laplace equation in a square domain with Neumann > boundary conditions on the three sides and on the top boundary two > equations are satisfied. > > I have the following code: > -------------------------------------- > V = FunctionSpace(mesh,"CG",1) > > eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) > phi0 = Function(V) > eta1 = Function(V) > phi1 = Function(V) > eta = TrialFunction(V) > phi = TrialFunction(V) > v = TestFunction(V) > > aeta = v*eta*ds(4) > Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx > > eta_problem = LinearVariationalProblem(aeta,Leta,eta1) > eta_solver = LinearVariationalSolver(eta_problem) > > aphi = v*phi*ds(4) > Lphi = (v*phi0 - g*dt*v*eta1)*ds(4) > > phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) > phi_solver = LinearVariationalSolver(phi_problem) > -------------------------------------- > > The error I get tells me that Object is in wrong state, and Matrix is > missing diagonal entry 39. > > I thought this had to do with the problem being singular, so I tried to > build the nullspace but I don't know how to define the VectorSpaceBasis. > Any help would be appreciated. > > Thank you very much. > > Regards, Anna. > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

3 6

Re: [firedrake] leaking memory
by David Ham 29 Apr '15

29 Apr '15

I agree about the hard limit. There is a limit to how much useful reordering we are ever going to do on the trace. Shall we set a default value for lazy_max_trace_length? Is 100 a good number? On Wed, 29 Apr 2015 at 10:00 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > On 29/04/15 08:59, Hiroe Yamazaki wrote: > > Dear all, > > > > I am experiencing a memory leak in our vertical slice model and > > wonder if there is any possibility of leaking memory associated > > with .assign(). Here I have a stand-alone code with the same > > function spaces and variable settings we use in the slice model: > > > > > https://bitbucket.org/hyamazak/test_files/src/f455e7e68722c9264845a8426ba91… > > > > What I do here is just doing > > > > new.vars += Deltavars old.vars.assign(new.vars) > > > > over and over again, and I see a constant increase of the memory > > usage both on my Linux and Macbook. I believe new.vars += Deltavars > > is not the source of leaking and wonder if there is any problem > > with doing old.vars.assign(new.vars). > > > > Do let me know if you have any ideas, > > OK, I observe the same memory leak. Here's what's going on: > > PyOP2 is running in lazy evaluation mode, which means that no > computation actually occurs until you ask to look at the resulting > data (perhaps by writing it to a file, or dropping into a solve). > Instead, we store an "execution trace" which remembers the sequence of > computations you want to do. By default this execution trace is > unbounded in length, so, since you never look at the answer in the > above code, we end up creating a trace of length 100000. This is the > cause of the memory leak. > > Arguably we should have a limit on the trace length. Can you try > setting the following parameter to see if your leak goes away? > > parameters["pyop2_options"]["lazy_max_trace_length"] = 100 > > > Cheers, > > Lawrence > > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v1 > > iQEcBAEBAgAGBQJVQJ2DAAoJECOc1kQ8PEYvr8gIALFwcyFV3XIcK+S3uzKDYkZ+ > aGR5EprrM1dvn2U4nS9eGTm52a9fxSYaZq3H++JfwdfZ0J/XWCFxfHkwAoCxmWRC > 23m+sohS5om8AjhmW4hU8S2LK6RQknDiKxt14IUE5nZXcf5aIvFF6/l3/pclptvg > OBorwxDk6LqHS17/4UOhIMjkcRLEppZ3TuH/NLqo7vftgnclf77ONjWfFP8y4DIh > 6AMPmcbx1WXLUZ+UIB7llMdaNgsH98Ei7v+hlvQnrdMxFAm4AbplAoVW6PGeTJ2z > h29W6LqVddrWaE7QMp2N+r732yY/3D7Pmjk0Iw5mO2uyYkCFrBVUUxUI+e6tqs4= > =STEw > -----END PGP SIGNATURE----- > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

2 1

leaking memory
by Hiroe Yamazaki 29 Apr '15

29 Apr '15

Dear all, I am experiencing a memory leak in our vertical slice model and wonder if there is any possibility of leaking memory associated with .assign(). Here I have a stand-alone code with the same function spaces and variable settings we use in the slice model: https://bitbucket.org/hyamazak/test_files/src/f455e7e68722c9264845a8426ba91… What I do here is just doing new.vars += Deltavars old.vars.assign(new.vars) over and over again, and I see a constant increase of the memory usage both on my Linux and Macbook. I believe new.vars += Deltavars is not the source of leaking and wonder if there is any problem with doing old.vars.assign(new.vars). Do let me know if you have any ideas, all the best, Hiroe

2 1