Re: [firedrake] Wrong matrix assembly
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@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@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@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/14334d3c20b9f10ed7c1246cde9e3cb60b1c75e4/Inequality%20constraint/?at=master> 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
_______________________________________________ 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
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/14334d3c20b9f10ed7c1246cde9e3cb60b1c75e4/Inequality%20constraint/parameters.py?at=master&fileviewer=file-view-default>, line 88)? 2b. In my Matlab code I use a 2-point Gauss quadrature. What quadrature is Firedrake using? 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@leeds.ac.uk <mailto:a.kalogirou@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@leeds.ac.uk <mailto:a.kalogirou@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@leeds.ac.uk <mailto:a.kalogirou@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/14334d3c20b9f10ed7c1246cde9e3cb60b1c75e4/Inequality%20constraint/?at=master> 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/ <http://www1.maths.leeds.ac.uk/%7Ematak/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Anna Kalogirou
-
David Ham