On Mon, 13 Jun 2016 at 12:40 Anna Kalogirou <a.kalogirou@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, 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@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 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 list
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
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake