June 2016 - firedrake - imperial.ac.uk

Naming components of mixed function
by Anna Kalogirou 14 Jun '16

14 Jun '16

Dear all, How can I name the components of a function in a mixed function space, so that I can then easily visualize the results based on their names? For example, say that I have V = FunctionSpace(mesh, "CG", 1) W = V*V w = Function(W) w1, w2 = split(w) output_file = File('results.pvd') output_file.write(w1, w2, time=t) I would like to give names to w1 and w2. Regards, Anna. -- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

2 2

Re: [firedrake] Naming components of mixed function
by Colin Cotter 14 Jun '16

14 Jun '16

What kind of visualisation do you want to do? In Paraview, you can just to visualise a single component of a vector field. all the best --cjc On 14 June 2016 at 14:00, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Dear all, > > How can I name the components of a function in a mixed function space, > so that I can then easily visualize the results based on their names? > > For example, say that I have > > V = FunctionSpace(mesh, "CG", 1) > W = V*V > > w = Function(W) > w1, w2 = split(w) > > output_file = File('results.pvd') > output_file.write(w1, w2, time=t) > > I would like to give names to w1 and w2. > > Regards, > > Anna. > > -- > > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] Wrong matrix assembly
by David Ham 14 Jun '16

14 Jun '16

On Mon, 13 Jun 2016 at 12:40 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Dear David, > > Thank you for sharing the detailed calculation. I am trying to find out > what is causing the discrepancy. > > I have a few of questions: > 1. Why did you multiply (mu_bar[2] * phi_2(x))? > 2a. Does firedrake pick up that mu_bar is a function of x (the exact > expression in file parameters.py > <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…>, > line 88)? > parameters.py:88 interpolates the analytic expression into the P1 Function mu_bar. This means that the full form of the Function is: mu_bar(x) = sum_i mu_bar[i] * phi_i(x) Because mu_bar[3] and mu_bar[4] are 0, the only term in this sum which is non-zero for the [3, 3] entry is: mu_bar[2] * phi_2(x) Note that mu_bar is not the analytic expression itself, it is the P1 interpolant of the analytic expression. > 2b. In my Matlab code I use a 2-point Gauss quadrature. What quadrature is > Firedrake using? > Firedrake calculates the required quadrature degree to achieve complete quadrature and uses that. In this case, because the degree of the integrand is 3, a two point quadrature will be selected. Firedrake employs Gauss-Lobatto quadrature. Regards, David > > It's a bit confusing that the results "almost" agree on the first nodes, > but are quite different closer to the "waterline" point (where mu_bar > becomes zero). > > Thank you. > > Regards, > > Anna. > > > > On 08/06/16 15:38, David Ham wrote: > > OK. Let's look at entry 3,3 in the matrix (ie the last nonzero entry - > what matlab would call the 4,4 entry). > > Note that the last two entries in mu_bar are zero. > > This means that the only nonzero integral in the calculation of this entry > is: > > \int_{cell 3} 2/dt * phi_3(x) * phi_3(x) * (mu_bar[2] * phi_2(x)) dx > > changing coordinates to local we have: > > int_0^1 1/det(J) * 2/dt * X * X * (1-X) * mu_bar[2] dX > > det(J) = 0.25 > dt = L/(Nx*2*pi*sqrt(g*H0)) = 0.04017220189012352 > > int_0^1 X*X*(1-X) dX = int_0^1 X^2 -X^4 dX = 1/12 > > 1/12 * 2 / dt * 0.25 * mu_bar[2] = 0.20775215537495284 > > Note that this is the same as the Firedrake result and not the same as the > Matlab code. I therefore claim that the problem is in the Matlab code, and > not a bug in Firedrake > > Regards, > > David > > > On Wed, 8 Jun 2016 at 14:44 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> > wrote: > >> Yes, it is correct. In that code I define a very small mesh just to >> inspect the values of the matrix. The same variables in Matlab are: >> >> mu_bar = >> >> 0.346930948820870 >> 0.283267933530252 >> 0.200300676691942 >> 0 >> 0 >> >> B_mu = >> >> (1,1) 1.377548049348976 >> (2,1) 0.655870455292249 >> (1,2) 0.655870455292249 >> (2,2) 2.344153409095711 >> (3,2) 0.506503005236205 >> (2,3) 0.506503005236205 >> (3,3) 1.641943147430029 >> (4,3) 0.280003094180926 >> (3,4) 0.280003094180926 >> (4,4) 0.405861612962805 >> >> You can see that even though mu_bar is the same as the one obtained from >> firedrake, the matrix entries are different. I don't understand why.. The >> rest of the matrices, i.e. M, A and C_lam, have the same entries as in >> firedrake. >> >> Thanks, >> >> Anna. >> >> >> >> On 08/06/16 14:38, Fabio Luporini wrote: >> >> Hello Anna, >> >> under the Inequality constraint folder, I'm running "python buoy-swe.py" >> and, besides the TSFC/COFFEE prints, I'm obtaining: >> >> mu_bar >> ------ >> [ 0.34693095 0.28326793 0.20030068 0. -0. ] >> >> Matrix b_mu >> ----------- >> [[ 1.37331778 0.6536432 0. 0. 0. ] >> [ 0.6536432 2.33042494 0.50155807 0. 0.] >> [ 0. 0.50155807 1.54031885 0.20775216 0.] >> [ 0. 0. 0.20775216 0.20775216 0. ] >> [ 0. 0. 0. 0. 0. ]] >> >> What values do you expect for b_mu ? Am I running the right thing? >> >> I tried switching off some optimisations but I still get the same output. >> >> Thanks! >> >> -- Fabio >> >> 2016-06-08 15:16 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: >> >>> Hi Fabio, >>> >>> Sorry for that. It should be public now. >>> >>> Anna. >>> >>> >>> >>> On 08/06/16 13:56, Fabio Luporini wrote: >>> >>> Hello Anna, >>> >>> the code is not accessible as it appears to be in a private repository. >>> Could you fix this? Thanks >>> >>> -- Fabio >>> >>> 2016-06-08 14:10 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: >>> >>>> Dear all, >>>> >>>> I have a simple question regarding a matrix assembly in Firedrake which >>>> I believe is not computed correctly. >>>> >>>> The relevant code can be found here >>>> <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…> >>>> and the issue is with the B_mu matrix in solvers.py, line 51. I compare the >>>> result with my own Matlab code and the two don't agree. Also, the rest of >>>> the matrices M, A, C_lam agree in both sets of code. >>>> >>>> Can someone help with what the problem might be? >>>> >>>> Thanks, >>>> >>>> Anna. >>>> >>>> -- >>>> >>>> >>>> Dr Anna Kalogirou >>>> Research Fellow >>>> School of Mathematics >>>> University of Leeds >>>> >>>> http://www1.maths.leeds.ac.uk/~matak/ >>>> >>>> >>> >>> >>> _______________________________________________ >>> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

1 0

Re: [firedrake] Wrong matrix assembly
by David Ham 13 Jun '16

13 Jun '16

OK. Let's look at entry 3,3 in the matrix (ie the last nonzero entry - what matlab would call the 4,4 entry). Note that the last two entries in mu_bar are zero. This means that the only nonzero integral in the calculation of this entry is: \int_{cell 3} 2/dt * phi_3(x) * phi_3(x) * (mu_bar[2] * phi_2(x)) dx changing coordinates to local we have: int_0^1 1/det(J) * 2/dt * X * X * (1-X) * mu_bar[2] dX det(J) = 0.25 dt = L/(Nx*2*pi*sqrt(g*H0)) = 0.04017220189012352 int_0^1 X*X*(1-X) dX = int_0^1 X^2 -X^4 dX = 1/12 1/12 * 2 / dt * 0.25 * mu_bar[2] = 0.20775215537495284 Note that this is the same as the Firedrake result and not the same as the Matlab code. I therefore claim that the problem is in the Matlab code, and not a bug in Firedrake Regards, David On Wed, 8 Jun 2016 at 14:44 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Yes, it is correct. In that code I define a very small mesh just to > inspect the values of the matrix. The same variables in Matlab are: > > mu_bar = > > 0.346930948820870 > 0.283267933530252 > 0.200300676691942 > 0 > 0 > > B_mu = > > (1,1) 1.377548049348976 > (2,1) 0.655870455292249 > (1,2) 0.655870455292249 > (2,2) 2.344153409095711 > (3,2) 0.506503005236205 > (2,3) 0.506503005236205 > (3,3) 1.641943147430029 > (4,3) 0.280003094180926 > (3,4) 0.280003094180926 > (4,4) 0.405861612962805 > > You can see that even though mu_bar is the same as the one obtained from > firedrake, the matrix entries are different. I don't understand why.. The > rest of the matrices, i.e. M, A and C_lam, have the same entries as in > firedrake. > > Thanks, > > Anna. > > > > On 08/06/16 14:38, Fabio Luporini wrote: > > Hello Anna, > > under the Inequality constraint folder, I'm running "python buoy-swe.py" > and, besides the TSFC/COFFEE prints, I'm obtaining: > > mu_bar > ------ > [ 0.34693095 0.28326793 0.20030068 0. -0. ] > > Matrix b_mu > ----------- > [[ 1.37331778 0.6536432 0. 0. 0. ] > [ 0.6536432 2.33042494 0.50155807 0. 0.] > [ 0. 0.50155807 1.54031885 0.20775216 0.] > [ 0. 0. 0.20775216 0.20775216 0. ] > [ 0. 0. 0. 0. 0. ]] > > What values do you expect for b_mu ? Am I running the right thing? > > I tried switching off some optimisations but I still get the same output. > > Thanks! > > -- Fabio > > 2016-06-08 15:16 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: > >> Hi Fabio, >> >> Sorry for that. It should be public now. >> >> Anna. >> >> >> >> On 08/06/16 13:56, Fabio Luporini wrote: >> >> Hello Anna, >> >> the code is not accessible as it appears to be in a private repository. >> Could you fix this? Thanks >> >> -- Fabio >> >> 2016-06-08 14:10 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: >> >>> Dear all, >>> >>> I have a simple question regarding a matrix assembly in Firedrake which >>> I believe is not computed correctly. >>> >>> The relevant code can be found here >>> <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…> >>> and the issue is with the B_mu matrix in solvers.py, line 51. I compare the >>> result with my own Matlab code and the two don't agree. Also, the rest of >>> the matrices M, A, C_lam agree in both sets of code. >>> >>> Can someone help with what the problem might be? >>> >>> Thanks, >>> >>> Anna. >>> >>> -- >>> >>> >>> Dr Anna Kalogirou >>> Research Fellow >>> School of Mathematics >>> University of Leeds >>> >>> http://www1.maths.leeds.ac.uk/~matak/ >>> >>> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

2 1

Setting function value on elements manually
by George Ovchinnikov 10 Jun '16

10 Jun '16

Dear all, I have a simple function W = FunctionSpace(mesh, "DG", 0) w = Function(W) and I want to set i-th element of it manually. For this I do: for j in range((n+1)**2): w.interpolate(Expression("1")) l = [int(x) for x in bin(j)[2:]] l = ((n+1)**2 - len(l))*[0] + l w_vec = w.vector().array() for i in range(n*n): if l[i] == 1: w_vec[i] = 1 else: w_vec[i] = 0.001 File("w" + str(j) + ".pvd") << w but from the last line I see the the same picture every time. What am I doing wrong? Sincerely, George

2 1

Re: [firedrake] Valgrind and Firedrake
by Lawrence Mitchell 09 Jun '16

09 Jun '16

Please send emails about firedrake issues to the firedrake list, or submit github issues. The intersection of the expertise is larger than that of any one individual, moreover you're much more likely to get an answer quickly. > On 9 Jun 2016, at 14:51, Nicolas Barral <n.barral(a)imperial.ac.uk> wrote: > > Hi Lawrence, > > I am still trying to find the bug I mentioned on IRC two weeks ago (point evaluation raising a not in domain error for points in the domain), and I have been totally unsuccessful so I need your help. > > You suggested using valgrind to find a memory issue, which I tried with the regular build of Firedrake (not in debug mode). Even with the suppression file suggested on Python's website, I get a many thousand line log, most errors being unusable because the debug symbols were not exported. Do any of the lines look like they are memory errors in likely candidates. That is, either code compiled by us (library has some kind of md5 hash in its name) or else inside libspatialindex? > So I thought I should try to compile python with the adequate debug options, and use this python with Firedrake. But how do I do that ? I tried to run firedrake-install with a PATH pointing to my custom build of Python (Miklos's suggestion), but it fails when checking/installing virtualenv (it keeps stopping at "Requirement already satisfied Virtual env installed. Please run firedrake-install again."). It sees the virtualenv from the standard python, which might be a problem ? So I'm stuck at that point. I do not know how to install multiple versions of python and get them to pick up the correct virtualenv and similar. I would have thought that your hand-installed python would not find the globally installed virtualenv stuff, but maybe not. Can you just make a virtualenv "by hand". We effectively just run "python -c 'import virtualenv'" and use that. So if that doesn't work that's somewhere to start. Normally if you're on linux you can just install the package that adds debug symbols for python (on ubuntu this is called python-debug I think). > My plan was then to compile the modules/dependencies in debug mode (notably the libspatialindex). Would that be really helpful ? And how do I do that since I can't just run configure/make in the right directory ? firedrake-install runs configure/make for you. You can determine how to installed the shared libraries in debug mode and add that as an option to firedrake-install, we will gladly accept patches for this. > Another idea came to me, and I searched for the cells of the mesh surrounding the points "not in domain" and computed their barycentric coordinates in the corresponding cell. It turns out that the point is always on a vertex or an edge (one or two barycentric coordinates ~ 1e-16). Could there be a bug in libspatialindex or the at function in these cases ? So what you're saying is that the point is right on the boundary between two cells? I guess there are two ways the point location can fail. Either the point is not found in any bounding box. Or, the point is found by libspatialindex in a bounding box. But that bucket does not actually contain the cell in question, so then the linear search for point location may fail due to floating point rounding. Which of these two are occurring? All of this code is compiled by us and is pretty straight line C, so you can just compile with debugging and step through it in the failing case and figure out what's going on. > I wanted to talk to you about all that in person, but I'm not sure we'll see each other before a while, since I'll be in Bath next week until wednesday... The faulty piece of code is likely to be removed when we will have a good Firedrake implementation of the interpolation and the metric computation, but maybe there's something deeper I'm not seeing which is worth being understood, so any help would be helpful :) I am at Imperial on Thursday next week. Lawrence

3 3

Re: [firedrake] Wrong matrix assembly
by Fabio Luporini 08 Jun '16

08 Jun '16

Hello Anna, under the Inequality constraint folder, I'm running "python buoy-swe.py" and, besides the TSFC/COFFEE prints, I'm obtaining: mu_bar ------ [ 0.34693095 0.28326793 0.20030068 0. -0. ] Matrix b_mu ----------- [[ 1.37331778 0.6536432 0. 0. 0. ] [ 0.6536432 2.33042494 0.50155807 0. 0.] [ 0. 0.50155807 1.54031885 0.20775216 0.] [ 0. 0. 0.20775216 0.20775216 0. ] [ 0. 0. 0. 0. 0. ]] What values do you expect for b_mu ? Am I running the right thing? I tried switching off some optimisations but I still get the same output. Thanks! -- Fabio 2016-06-08 15:16 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: > Hi Fabio, > > Sorry for that. It should be public now. > > Anna. > > > > On 08/06/16 13:56, Fabio Luporini wrote: > > Hello Anna, > > the code is not accessible as it appears to be in a private repository. > Could you fix this? Thanks > > -- Fabio > > 2016-06-08 14:10 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: > >> Dear all, >> >> I have a simple question regarding a matrix assembly in Firedrake which I >> believe is not computed correctly. >> >> The relevant code can be found here >> <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…> >> and the issue is with the B_mu matrix in solvers.py, line 51. I compare the >> result with my own Matlab code and the two don't agree. Also, the rest of >> the matrices M, A, C_lam agree in both sets of code. >> >> Can someone help with what the problem might be? >> >> Thanks, >> >> Anna. >> >> -- >> >> >> Dr Anna Kalogirou >> Research Fellow >> School of Mathematics >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~matak/ >> >> > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

2 1

Hadamard (element wise) product
by Floriane Gidel [RPG] 08 Jun '16

08 Jun '16

Dear all, Some Hadamard products are involved in my equations, and I did not find a way to implement it. The idea is to apply a element by element product between two matrices A and B, where : - A and B have the same size, let say tensor of dim (N,) - M = A*B is also a tensor of dim (N,) where each elements M[i] = A[i]*B[i] (not summed!) ; (same with rank 2 tensors) So I would like to use the same kind of product as the dot product, but without summing the elements. Is that possible ? Thanks, Floriane

2 2

Re: [firedrake] Wrong matrix assembly
by Fabio Luporini 08 Jun '16

08 Jun '16

Hello Anna, the code is not accessible as it appears to be in a private repository. Could you fix this? Thanks -- Fabio 2016-06-08 14:10 GMT+02:00 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk>: > Dear all, > > I have a simple question regarding a matrix assembly in Firedrake which I > believe is not computed correctly. > > The relevant code can be found here > <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…> > and the issue is with the B_mu matrix in solvers.py, line 51. I compare the > result with my own Matlab code and the two don't agree. Also, the rest of > the matrices M, A, C_lam agree in both sets of code. > > Can someone help with what the problem might be? > > Thanks, > > Anna. > > -- > > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > >

2 1

Wrong matrix assembly
by Anna Kalogirou 08 Jun '16

08 Jun '16

Dear all, I have a simple question regarding a matrix assembly in Firedrake which I believe is not computed correctly. The relevant code can be found here <https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6…> and the issue is with the B_mu matrix in solvers.py, line 51. I compare the result with my own Matlab code and the two don't agree. Also, the rest of the matrices M, A, C_lam agree in both sets of code. Can someone help with what the problem might be? Thanks, Anna. -- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

1 0