- firedrake - imperial.ac.uk

Boundary IDs on extruded meshes
by Hanie Joodat 01 Aug '16

01 Aug '16

Dear all, I want to enforce DirichletBCs on an extruded mesh. For a normal unit cube tetrahedron mesh I could use the following face ID's: 1: x ==0 2: x==1 3: y==0 4: y==1 5: z==0 6: z==1 However, when I do the following: meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = ExtrudedMesh(meshbase, Nz, Dz) I get this error: Traceback (most recent call last): File "VMS_patch_test_3D_Double_firedrake.py", line 51, in <module> FrontBackBC1 = DirichletBC(wSpace.sub(0).sub(1),Constant(0.0),(3,4)) File "/data/Software/firedrake/src/firedrake/firedrake/bcs.py", line 68, in __init__ raise NotImplementedError("Indexed VFS bcs not implemented on extruded meshes") NotImplementedError: Indexed VFS bcs not implemented on extruded meshes What should I do to properly identify my boundary faces? Regards, Hanie

1 0

Re: [firedrake] Unexpected asymmetry of solution of Poisson equation
by Andrew McRae 01 Aug '16

01 Aug '16

Miklos is exactly right. What you see in the very corners is the following: if you have a quad with 3 zeroes and a nonzero (let's say a 1), if the quad is split along a 0-0 diagonal the you get a 0-0-0 triangle and a 0-1-0 triangle, but if the quad is split a 0-1 diagonal then you get two 0-0-1 triangles. It's even worse if the quad is +1 -1 -1 +1 in which case you get a 'valley' or a 'ridge', even though the actual function is symmetric about 0. To answer your last question, the triangular mesh isn't symmetric in that way. On 1 August 2016 at 11:08, George Ovchinnikov <lives9(a)gmail.com> wrote: > Dear all, > > I'm solving Poisson equation and I found some irregularity. Here the code: > > from firedrake import * > > mesh = UnitSquareMesh(8, 8, quadrilateral=True) > > V = FunctionSpace(mesh, "CG", 1) > W = FunctionSpace(mesh, "DG", 0) > w = Function(W) > w.interpolate(Expression(1)) #Initial material distribution > > bc0 = DirichletBC(V, Expression(("0")), 1) > bc1 = DirichletBC(V, Expression(("0")), 2) > bc2 = DirichletBC(V, Expression(("0")), 3) > bc3 = DirichletBC(V, Expression(("0")), 4) > bcs = [bc0, bc1, bc2, bc3] > > u = TrialFunction(V) > v = TestFunction(V) > > f = Function(V) > f.interpolate(Expression("1")) > L = f * v * dx > rhs = assemble(L) > > a = (w * dot(grad(u), grad(v))) * dx > frm = assemble(a, bcs = bcs) > sol = Function(V) > solve(frm, sol, rhs, solver_parameters = {'ksp_type': 'preonly', > 'pc_type': 'lu'}) > File("u.pvd") << sol > > The asymmetry of solution (see the figure attached) is clearly visible: > for example, upper left element is not the same as upper right. > > I'm wondering what causes this and can this be remedied? > > Also, I'm getting exactly same result (with the same asymmetry) on > triangular, not quadrilateral mesh. Why solutions of triangular and > quadrilateral meshes are the same? > > Sincerely, > George >

2 1

Unexpected asymmetry of solution of Poisson equation
by George Ovchinnikov 01 Aug '16

01 Aug '16

Dear all, I'm solving Poisson equation and I found some irregularity. Here the code: from firedrake import * mesh = UnitSquareMesh(8, 8, quadrilateral=True) V = FunctionSpace(mesh, "CG", 1) W = FunctionSpace(mesh, "DG", 0) w = Function(W) w.interpolate(Expression(1)) #Initial material distribution bc0 = DirichletBC(V, Expression(("0")), 1) bc1 = DirichletBC(V, Expression(("0")), 2) bc2 = DirichletBC(V, Expression(("0")), 3) bc3 = DirichletBC(V, Expression(("0")), 4) bcs = [bc0, bc1, bc2, bc3] u = TrialFunction(V) v = TestFunction(V) f = Function(V) f.interpolate(Expression("1")) L = f * v * dx rhs = assemble(L) a = (w * dot(grad(u), grad(v))) * dx frm = assemble(a, bcs = bcs) sol = Function(V) solve(frm, sol, rhs, solver_parameters = {'ksp_type': 'preonly', 'pc_type': 'lu'}) File("u.pvd") << sol The asymmetry of solution (see the figure attached) is clearly visible: for example, upper left element is not the same as upper right. I'm wondering what causes this and can this be remedied? Also, I'm getting exactly same result (with the same asymmetry) on triangular, not quadrilateral mesh. Why solutions of triangular and quadrilateral meshes are the same? Sincerely, George

2 1

Re: [firedrake] Product in an expression
by David Ham 01 Aug '16

01 Aug '16

Hi Floriane, It looks rather like you are trying to express the basis function expansion at the Firedrake level. This is something which is likely to be quite hard because you're operating at the wrong abstraction level. Can you post at a weak form (i.e. PDEs and function spaces) level what it is you are trying to achieve? Many thanks, David On Mon, 1 Aug 2016 at 16:29 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > Dear all, > > > I would like to compute the following sum over i : phi(x,y,z,t) = > psi_i(x,y,t)varphi_i(z), where : > > - each psi_i (x,y,t) is a Firedrake function defined on a horizontal plane, > > -varphi_i is defined by a lagrange expansion of order n, that is: > varphi_i(z) = product( (z-z_k)/(z_i - z_k) ) for k = 0:n , and k different > from i. Each z_k can be computed analytically and denotes the discretised z > coordinate. For instance, z_k can be defined as z_k = (n-k)H0/n , for z in > [0,H0]. > > > How can I define varphi_i so that I can compute the dot product > dot(psi,varphi) on Firedrake ? I tried two options: > > - use the symbpy package to compute the expression of varphi_i with the > 'product' function. But in that case the expression of varphi_i depends on > z as a symbol, and I did not manage to use it correctly on Firedrake > afterwards, > > - define varphi_i(z) as a Firedrake function in z, but in that case, how > can I compute the product ? I would normally use an expression, but is > there a 'product expression' that I could interpolate ? Something like : > varphi[i].interpolate(Expression( "product((x[2]-z[k]) / (z[i]-z[k]) )", > z=z, i=i, for k=0:n)) > > > Thank you for your help, > > Floriane >

2 1

Product in an expression
by Floriane Gidel [RPG] 01 Aug '16

01 Aug '16

Dear all, I would like to compute the following sum over i : phi(x,y,z,t) = psi_i(x,y,t)varphi_i(z), where : - each psi_i (x,y,t) is a Firedrake function defined on a horizontal plane, -varphi_i is defined by a lagrange expansion of order n, that is: varphi_i(z) = product( (z-z_k)/(z_i - z_k) ) for k = 0:n , and k different from i. Each z_k can be computed analytically and denotes the discretised z coordinate. For instance, z_k can be defined as z_k = (n-k)H0/n , for z in [0,H0]. How can I define varphi_i so that I can compute the dot product dot(psi,varphi) on Firedrake ? I tried two options: - use the symbpy package to compute the expression of varphi_i with the 'product' function. But in that case the expression of varphi_i depends on z as a symbol, and I did not manage to use it correctly on Firedrake afterwards, - define varphi_i(z) as a Firedrake function in z, but in that case, how can I compute the product ? I would normally use an expression, but is there a 'product expression' that I could interpolate ? Something like : varphi[i].interpolate(Expression( "product((x[2]-z[k]) / (z[i]-z[k]) )", z=z, i=i, for k=0:n)) Thank you for your help, Floriane

1 0

DirichletBC not working anymore
by Justin Chang 31 Jul '16

31 Jul '16

Hi all, Last year (~8 months ago) the following worked: mesh = UnitSquareMesh(10,10) V = VectorFunctionSpace(mesh, "CG", 1) P = FunctionSpace(mesh, "CG", 1) W = V*Q bc1 = DirichletBC(W.sub(0).sub(0), Expression("cos(pi*x[0])*sin(pi*x[1])"), (1,2)) Today, I get these errors: Traceback (most recent call last): File "Darcy_LS.py", line 32, in <module> bc1 = DirichletBC(W.sub(0).sub(0), Expression("cos(pi*x[0])*sin(pi*x[1])"), (1,2)) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/bcs.py", line 54, in __init__ self.function_arg = g File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/bcs.py", line 100, in function_arg g = function.Function(self._function_space).interpolate(g) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/function.py", line 243, in __init__ val=val, name=name, dtype=dtype) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/function.py", line 79, in __init__ self.dat = function_space.make_dat(val, dtype, self.name(), uid=self.uid) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/functionspaceimpl.py", line 798, in make_dat raise ValueError("Can't build Function on %s function space" % self.typ) AttributeError: 'ProxyFunctionSpace' object has no attribute 'typ' What's going on, and how do I work around this? Also, is it suppose to say 'self.type' instead of 'self.typ'? Thanks, Justin

2 5

Re: [firedrake] Runtime visualization
by David Ham 28 Jul '16

28 Jul '16

Hi Justin, That's a very well timed question. Happily Fangyi (Thomas) Zhou is with us as a summer student working on this precise problem. He's started out by implementing some functionality for 1 and 2 dimensions using matplotlib. If you're feeling brave you can try the plotting and plotting-1d branches but they will hopefully land in master in the next week or so. Meanwhile I'm at WCCM and was in a visualisation session this morning. It appears that what both PyFR and GMSH are now doing is using catalyst, which is a runtime paraview interface (http://www.paraview.org/in-situ/) I think it would be really cool if Thomas has a look at that next (although I haven't yet had a chance to discuss that with him since we're not in the same location or timezone right now). Regards, David On Wed, 27 Jul 2016 at 13:56 Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > One of the beautiful things about FEniCS is its ability to produce runtime > plots of the solution. Is there (or will there be) a way to do this in > Firedrake? It is a little bit cumbersome to dump the solution into a pvd > file and (re)open it in paraview/VisIt each time you want to tweak > something. > > Thanks, > Justin >

2 1

Re: [firedrake] different SNES convergence on different machine?
by Colin Cotter 27 Jul '16

27 Jul '16

Hi TJ, We are currently discussing this in Slack if you would like to pop in? On 27 July 2016 at 13:54, Sun, Tianjiao <tianjiao.sun14(a)imperial.ac.uk> wrote: > Hi Lawrence, > > 20x20 is the smallest I tried. For smaller mesh the problem is not > well-conditioned. I can think of some other test cases for smaller mesh. > > I tried to run with snes_fd, and I got much similar but not still not > identical results on both machines, so it does look like something to do > with the Jacobian. I will try comparing the Jacobians and report back. > > Thanks, > > -TJ > > On 27/07/16 13:36, Lawrence Mitchell wrote: > > > On 27/07/16 13:09, Tianjiao Sun wrote: > > Hi team, > > The attached script is a trimmed down version of the sort of problem I'm > trying to run. Running it on most machines (I've tried 4 machines), I > got this output: > > $ python mytest.py > 0 SNES Function norm 0.873814 > 1 SNES Function norm 0.245545 > 2 SNES Function norm 0.0178611 > 3 SNES Function norm 0.000156473 > 4 SNES Function norm 1.26023e-09 > > However, on my workstation in college (where I do most work > unfortunately), I got: > > 0 SNES Function norm 0.873814 > 1 SNES Function norm 0.432782 > 2 SNES Function norm 0.0478636 > 3 SNES Function norm 0.0405004 > 4 SNES Function norm 0.0405004 > 5 SNES Function norm 0.0405004 > 6 SNES Function norm 0.0405004 > 7 SNES Function norm 0.0405004 > 8 SNES Function norm 0.0405004 > 9 SNES Function norm 0.0405004 > 10 SNES Function norm 0.0405004 > > I have tried to update and also reinstall firedrake, and pip list shows > same version for FIAT, petsc, petsc4py, UFL (firedrake version is shown > as 'unknown' though). The snes objects are identical as well. Is there > anything else that I can check to see where the discrepancy occurs? > > Does this problem exhibit on a small mesh? > > Can you try comparing the jacobians on a machine that works and one > that fails? > > Lawrence > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

different SNES convergence on different machine?
by Tianjiao Sun 27 Jul '16

27 Jul '16

Hi team, The attached script is a trimmed down version of the sort of problem I'm trying to run. Running it on most machines (I've tried 4 machines), I got this output: $ python mytest.py 0 SNES Function norm 0.873814 1 SNES Function norm 0.245545 2 SNES Function norm 0.0178611 3 SNES Function norm 0.000156473 4 SNES Function norm 1.26023e-09 However, on my workstation in college (where I do most work unfortunately), I got: 0 SNES Function norm 0.873814 1 SNES Function norm 0.432782 2 SNES Function norm 0.0478636 3 SNES Function norm 0.0405004 4 SNES Function norm 0.0405004 5 SNES Function norm 0.0405004 6 SNES Function norm 0.0405004 7 SNES Function norm 0.0405004 8 SNES Function norm 0.0405004 9 SNES Function norm 0.0405004 10 SNES Function norm 0.0405004 I have tried to update and also reinstall firedrake, and pip list shows same version for FIAT, petsc, petsc4py, UFL (firedrake version is shown as 'unknown' though). The snes objects are identical as well. Is there anything else that I can check to see where the discrepancy occurs? Many thanks, -TJ

2 2

segfaults on Xenial
by Stephan Kramer 27 Jul '16

27 Jul '16

Dear firedrake devs After upgrading my workstation to Ubuntu Xenial, and rebuilding firedrake from scratch I get several segfaults in the running of firedrake tests. The ones I've seen so far were segfaulting in petsc: kernel -> MatSetValuesLocal() -> ISLocalToGlobalMappingApply() If I switch on the debug PyOP2 option and run it in gdb, it seems the values passed to MatSetValuesLocal are all sensible. Unfortunately if I rebuild firedrake with a debugging-petsc the tests all pass (i.e. the segfault goes away). Any further hints how to debug this? Cheers Stephan

2 2