November 2016 - firedrake - imperial.ac.uk

Best way to detect memory leaks
by Justin Chang 01 Dec '16

01 Dec '16

Hi all, What's the best way or tool to detect for potential memory leaks for Firedrake? Normally I use the PETSc log functions but in my Firedrake code I am doing lots of extra curricular programming which isn't seen by PETSc. Thanks, Justin

2 3

Re: [firedrake] Mixed system
by David Ham 01 Dec '16

01 Dec '16

Yes. This is because trace elements were very recently (after my email) merged into other components of Firedrake and the globals branch needs to catch up. I've just pushed the relevant commits to the real-function-space branch so that might work now. Regards, David On Mon, 7 Nov 2016 at 17:15 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: Dear David, I have just updated firedrake using the branches you mention in the email below but I now get some error related to the mesh (?), when running various codes which were previously working. Any idea what this is? Thanks, Anna. File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/functionspace.py", line 42, in make_scalar_element mesh.init() File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/mesh.py", line 898, in init self._callback(self) File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/mesh.py", line 1245, in callback dim=geometric_dim) File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/functionspace.py", line 174, in VectorFunctionSpace return FunctionSpace(mesh, element, name=name) File "<decorator-gen-280>", line 2, in FunctionSpace File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/PyOP2/pyop2/profiling.py", line 59, in wrapper return f(*args, **kwargs) File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/functionspace.py", line 127, in FunctionSpace check_element(element) File "/e/data1/users/matak/opt/firedrake-installation/firedrake/src/firedrake/firedrake/functionspace.py", line 85, in check_element ufl.TraceElement, ufl.HDivElement, ufl.HCurlElement): AttributeError: 'module' object has no attribute 'TraceElement' On 04/11/16 16:16, David Ham wrote: A brief status update on this: the real-function-spaces branch of Firedrake taken together with the globals_petsc_changes branch of PyOP2 no longer exhibits the direct crash observed. However I am also yet to actually successfully solve using a block preconditioner. I need to consult with Lawrence to establish whether the problem is another bug or pilot error. Regards, David On Mon, 24 Oct 2016 at 11:37 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > On 24 Oct 2016, at 11:13, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> wrote: > > Hi Lawrence, > > Any update regarding the bug reported below? This appears to be more involved than I had first thought. So no luck yet, sorry! David and I need to find some time together to work out the differences/corner cases. Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake(a)imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake

2 1

eigenvalue problems with SLEPc
by Francis Poulin 01 Dec '16

01 Dec '16

Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Are there any other examples that someone might be able to share that might to able to do this? Cheers, 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

3 12

Re: [firedrake] no normal flow Boundary Conditions
by David Ham 29 Nov '16

29 Nov '16

Hi Francis, On Tue, 29 Nov 2016 at 04:29 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > Some progress since last I wrote: > > 1) I was able to get the QG linear solver working. I do have a couple of > questions though. > > > https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_qg.… > > First, the solver that Jemma and I used for the nonlinear QG code seems to > yield the trivial solution. Any idea why that is? > I don't think it's finding a trivial solution. I think that the issue is that you never solve the system. Having set up your linear solver, you are missing: psi_solver.solve() However, the linear solver in your code has another problem. The equation isn't symmetric (the beta*phi*psi.dx(0)*dx term) so CG is unlikely to work. For 2D problems, using a direct solver is usually the right start: solver_parameters={ 'ksp_type':'preonly', 'pc_type':'lu' } > Second, should we be able to use matplotlib to plot the solutions? I was > able to get a demo to work but I can't do it for this script. > There are two problems here. The first is that you are calling the plot command from Matplotlib. Unsurprisingly, this does not know about Firedrake Functions. You need to call the Firedrake plot command. p = plot(psi_soln) p.show() However, more importantly, the matplotlib functionality is mostly useful for debugging small systems. Your 50 x 50 mesh is already going to be very slow via this route. > > 2) The SW problem is not working much better I'm afraid. It does not > converge to a solution, which is I suppose better than converging to an > incorrect solution. I have a couple of concerns. > > > https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw.… > > First, i am trying to impose no-normal flow BCs. I am trying to do what I > understood Colin suggested but since it's not converging, maybe I > misunderstood. > > Second, I presume there is a better solver that I should be using? > Same problem as above. You are attempting to use CG for a non-symmetric system. In this case it's a saddle point system. Try using a direct solver instead. Regards, David > > Thanks in advance for your help. > > > Cheers, 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 <(519)%20888-4567> > > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Francis Poulin [fpoulin(a)uwaterloo.ca] > *Sent:* Saturday, November 26, 2016 6:54 PM > > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] no normal flow Boundary Conditions > Hello, > > Minor progress. > > I found out that the problem with my QG code was how I imposed the > boundary conditions. Now I am happy to say that they both run. > Unfortunately, they both give me a zero solution even though the winds are > not zero. I plotting the winds and so I know that they're not zero. Maybe > I"ll have more luck tomorrow. > > Cheers, > 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 <(519)%20888-4567> > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Francis Poulin [fpoulin(a)uwaterloo.ca] > *Sent:* Saturday, November 26, 2016 12:58 PM > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] no normal flow Boundary Conditions > > Hello again, > > Minor update on my wind-driven gyre problem. > > I just put together a SW version of the code that solve the Linear Stommel > problem. You can find the code in the same github repo, > > > https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_sw.… > > To my surprise the solver does not complain, and it is exactly the same > thing as in the QG version so maybe that part of it is ok. Unfortunately, > the solution seems to be zero. So this code is at least giving me an > answer, but a wrong one. > > Cheers, 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 <(519)%20888-4567> > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Francis Poulin [fpoulin(a)uwaterloo.ca] > *Sent:* Saturday, November 26, 2016 7:36 AM > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] no normal flow Boundary Conditions > > Hello, > > Thank you Colin and everyone else for your responses. That was very > helpful. > > We are working on a bunch of problems right now, both QG and SW. We have > many different codes working in FEniCS but I am going to start translating > them into Firedrake. Maybe make a demo if people seem interested? > > The first code is to solve the Linear Stommel problem in one-layer QG. I > was able to translate most of it ok but I seem to be having problems with > the solver for some reason. > If anyone had any advice on what I am doing wrong that would be greatly > appreciated. You can find a link on github > > > https://github.com/francispoulin/firedrakeQG/blob/master/linear_stommel_qg.… > > After I get this working I will like to solve a variety of other > problems. One is the Munk problem, then look at non-linear gyre solutions, > then find these solutions in SW. > > Cheers, 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 <(519)%20888-4567> > > ------------------------------ > *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:* Friday, November 25, 2016 4:15 AM > *To:* firedrake > *Subject:* Re: [firedrake] no normal flow Boundary Conditions > > Hi Francis, > That depends on the formulation of the equations. If you are using a > streamfunction-vorticity formulation as for QG, then you need to set a zero > boundary condition for the streamfunction, as a scalar. If you are using a > velocity formulation based on H(div) spaces (BDM, RT etc.) then there > are only u.n DOFs on the boundary, so setting: > > bcs = [DirichletBC(Z.sub(0), Constant((0, 0)), (blah,)), ... > > actually only sets the normal component. That's one of the reasons that I like these spaces for GFD. > > > cheers > > --cjc > > > On 25 November 2016 at 01:17, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > Hello, > > Sorry for the bother but I have a simple question to ask. > > I am currently developing some code to solve for wind-driven gyres and > basin modes using QG and SW. My question is how is it best to impose > no-normal flow BCs, > > $$\vec u \cdot \hat n = 0$$ > > for any boundary. I found this in the Navier-Stokes example which says > that we should impose zero on particular boundaries. I guess that would > work but if one had a more complicated geometry I imagine that would have > it's problems. > > bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), > DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))] > > > Any suggestions would be greatly appreciated. > > Cheers, > 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 > > > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

3 7

Re: [firedrake] [petsc-users] pcfieldsplitting not working for VI
by Justin Chang 29 Nov '16

29 Nov '16

Okay, switched over to the firedrake list On Tue, Nov 29, 2016 at 2:29 AM, Lawrence Mitchell < lawrence.mitchell(a)imperial.ac.uk> wrote: > I suspect this is not a petsc issue per se. Shall we migrate this to the > firedrake list? > > Lawrence > > On 29 Nov 2016, at 07:40, Justin Chang <jychang48(a)gmail.com> wrote: > > FWIW I ended up doing this: > > rs_solver.setDM(W._dm()) > > and was left with this error: > > Traceback (most recent call last): > File "sphere.py", line 180, in <module> > virs() > File "sphere.py", line 163, in virs > rs_solver.solve(None,sol_vec) > File "PETSc/SNES.pyx", line 537, in petsc4py.PETSc.SNES.solve > (src/petsc4py.PETSc.c:170012) > petsc4py.PETSc.Error: error code 75 > [0] SNESSolve() line 4060 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/snes/ > interface/snes.c > [0] SNESSolve_VINEWTONRSLS() line 424 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/snes/ > impls/vi/rs/virs.c > [0] MatGetSubMatrix() line 7929 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/mat/ > interface/matrix.c > [0] MatGetSubMatrix_Nest() line 451 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/mat/ > impls/nest/matnest.c > [0] MatNestFindSubMat() line 424 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/mat/ > impls/nest/matnest.c > [0] MatNestFindIS() line 363 in /private/var/folders/92/ > 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-gYEj0p-build/src/mat/ > impls/nest/matnest.c > [0] Arguments are incompatible > [0] Could not find index set > > I must be doing something wrong here. > > Justin > > On Tue, Nov 29, 2016 at 1:24 AM, Justin Chang <jychang48(a)gmail.com> wrote: > >> Lawrence, >> >> I added the following line: >> >> rs_solver.setDM(a._dm) >> >> And it gives this error: >> >> Traceback (most recent call last): >> File "sphere.py", line 128, in <module> >> rs_solver.setDM(a._dm) >> AttributeError: 'Form' object has no attribute '_dm' >> >> Thanks, >> Justin >> >> On Tue, Nov 29, 2016 at 1:12 AM, Lawrence Mitchell < >> lawrence.mitchell(a)imperial.ac.uk> wrote: >> >>> >>> >>> > On 29 Nov 2016, at 02:02, Justin Chang <jychang48(a)gmail.com> wrote: >>> > >>> > Why is this happening? It seems to me VINEWTONRSLS is not >>> understanding that I have a two-field formulation. >>> >>> For the built in solvers, we set up a dmshell on the snes that knows how >>> to do field splitting. Because you set your snes up by hand, it isn't >>> available. You'll need to call snes.setDM. The mixed finctionspace provides >>> the correct dm via a ._dm property I think. >>> >>> Lawrence >>> >> >> >

2 3

Re: [firedrake] no normal flow Boundary Conditions
by Colin Cotter 29 Nov '16

29 Nov '16

Hi Francis, That depends on the formulation of the equations. If you are using a streamfunction-vorticity formulation as for QG, then you need to set a zero boundary condition for the streamfunction, as a scalar. If you are using a velocity formulation based on H(div) spaces (BDM, RT etc.) then there are only u.n DOFs on the boundary, so setting: bcs = [DirichletBC(Z.sub(0), Constant((0, 0)), (blah,)), ... actually only sets the normal component. That's one of the reasons that I like these spaces for GFD. cheers --cjc On 25 November 2016 at 01:17, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > Sorry for the bother but I have a simple question to ask. > > I am currently developing some code to solve for wind-driven gyres and > basin modes using QG and SW. My question is how is it best to impose > no-normal flow BCs, > > $$\vec u \cdot \hat n = 0$$ > > for any boundary. I found this in the Navier-Stokes example which says > that we should impose zero on particular boundaries. I guess that would > work but if one had a more complicated geometry I imagine that would have > it's problems. > > bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), > DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))] > > > Any suggestions would be greatly appreciated. > > Cheers, > 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 > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 4

Re: [firedrake] no normal flow Boundary Conditions
by Andrew McRae 25 Nov '16

25 Nov '16

On 25 November 2016 at 01:41, Justin Chang <jychang48(a)gmail.com> wrote: > For rectangular domains and vector Lagrange elements you just need to > constraint one component for the normal. For example > > bc1 = DirichletBC(Z.sub(0).sub(0),value,(1,2)) #x component > bc2 = DirichletBC(Z.sub(0).sub(1),value,(3,4)) #y component > Does this actually work? I know doing setting individual components for VFS inside Mixed is flaky, per https://github.com/ firedrakeproject/firedrake/issues/920. However, the snippet in that bug report sets multiple components one-by-one, which might be what leads to the breakage. > > Not sure how you would do this for a more complicated geometry > > Justin > > On Thu, Nov 24, 2016 at 7:17 PM Francis Poulin <fpoulin(a)uwaterloo.ca> > wrote: > >> Hello, >> >> Sorry for the bother but I have a simple question to ask. >> >> I am currently developing some code to solve for wind-driven gyres and >> basin modes using QG and SW. My question is how is it best to impose >> no-normal flow BCs, >> >> $$\vec u \cdot \hat n = 0$$ >> >> for any boundary. I found this in the Navier-Stokes example which says >> that we should impose zero on particular boundaries. I guess that would >> work but if one had a more complicated geometry I imagine that would have >> it's problems. >> >> bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), >> DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))] >> >> >> Any suggestions would be greatly appreciated. >> >> Cheers, >> 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 >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >

2 1

no normal flow Boundary Conditions
by Francis Poulin 25 Nov '16

25 Nov '16

Hello, Sorry for the bother but I have a simple question to ask. I am currently developing some code to solve for wind-driven gyres and basin modes using QG and SW. My question is how is it best to impose no-normal flow BCs, $$\vec u \cdot \hat n = 0$$ for any boundary. I found this in the Navier-Stokes example which says that we should impose zero on particular boundaries. I guess that would work but if one had a more complicated geometry I imagine that would have it's problems. bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))] Any suggestions would be greatly appreciated. Cheers, 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

2 1

Re: [firedrake] no normal flow Boundary Conditions
by Andrew McRae 25 Nov '16

25 Nov '16

The answer depends somewhat on what finite element space you are using for velocity. --- If you are using a CG/DG VectorFunctionSpace, it's possible to impose BCs strongly that set one cartesian component to 0. I think you access this by doing, e.g., DirichletBC(V.sub(0), 0, (1,)). This has the limitation that you mentioned (i.e. it only works if that boundary is axis-aligned). (Note: The code snippet you posted is setting both components of velocity, not just the normal component) Unfortunately setting just one component of velocity isn't implemented for mixed function spaces, see https://github.com/firedrakeproject/firedrake/ issues/920 --- An alternative is to impose the boundary condition weakly. I think this is as simple as adding a dot(v, n)*dot(u, n)*ds term to your equation, where v is a test function from the same space as u. This is equivalent to saying "integral of gamma * u.n over the boundary = 0" for all appropriate test functions gamma, which is clearly a weak version of u.n = 0. This works for arbitrary geometry. (And if your velocity field is DG then you should be imposing the BC weakly anyway, according to the FE experts on the list!) --- Finally, if you are using a vector-valued finite element space, e.g. an RT or BDM FunctionSpace, you can set no-normal flow BCs strongly with DirichletBC(V, 0, (1, 2, 3, 4)). This sets the boundary degrees of freedom to 0, and the boundary degrees of freedom precisely represent the normal flow velocity due to the way the FE space is defined [actually, it's scaled so I think it's a flux rather than a velocity, but scaled 0 is still 0]. Again, this works for arbitrary geometry. --- Hopefully one of these is suitable (imposing the BC weakly should always work). Best, Andrew On 25 November 2016 at 01:17, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > Sorry for the bother but I have a simple question to ask. > > I am currently developing some code to solve for wind-driven gyres and > basin modes using QG and SW. My question is how is it best to impose > no-normal flow BCs, > > $$\vec u \cdot \hat n = 0$$ > > for any boundary. I found this in the Navier-Stokes example which says > that we should impose zero on particular boundaries. I guess that would > work but if one had a more complicated geometry I imagine that would have > it's problems. > > bcs = [DirichletBC(Z.sub(0), Constant((1, 0)), (4,)), > DirichletBC(Z.sub(0), Constant((0, 0)), (1, 2, 3))] > > > Any suggestions would be greatly appreciated. > > Cheers, > 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 > >

1 0

Re: [firedrake] Getting centre mass of the cell
by David Ham 18 Nov '16

18 Nov '16

Dear George, Well one way to do something like this is to observe that the node of the DG0 space is at the centre of mass. Accordingly if you project mesh.coordinates into a VectorFunctionSpace with element DG 0 then you get a Function whose n-th value is the centre of mass of cell n. Regards, David On Fri, 18 Nov 2016 at 22:21 George Ovchinnikov <lives9(a)gmail.com> wrote: > Hello, > > There is a nice function, which gives me cell by its coordinates > "mesh.locate_cell(x)" > > I wonder, does it have some kind of inverse analogue, so I could get > coordinates for a cell (given by id of cell)? > > Sincerely, > George > > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

1 0