[firedrake] Issues with LMA of vertical derivative on extruded meshes

[Lawrence Mitchell]

On 22 Dec 2014, at 14:15, Eike Mueller <eike.h.mueller at googlemail.com> wrote:

> Dear firedrakers,
> 
> I came across a problem when doing a local matrix assembly of the weak derivative matrix from a DG space (W3) to a vertical-only HDiv space (W2_vert), i.e.:
> 
> order_vertical = 0
> U2 = FiniteElement('DG',interval,0)
> V0 = FiniteElement('CG',interval,1)
> V1 = FiniteElement('DG',interval,order_vertical)
> 
> W2_vert = FunctionSpace(mesh,HDiv(OuterProductElement(U2,V0)))
> W3 = FunctionSpace(mesh,OuterProductElement(U2,V1))
> 
> phi = TestFunction(W3)
> psi = TrialFunction(W3)
> u = TestFunction(W2_vert)
> w = TrialFunction(W2_vert)
> 
> form_D = phi*div(w)*dx
> form_DT = div(u)*psi*dx
> 
> I work with an extruded mesh with one vertical layer, for details, see the attached code which exposes the issue.
> 
> When I do a local matrix assembly for form_D and form_DT, I get the following entries of the LMA's (the grid has three cells, and I attach a 2x1 or 1x2 matrix representing the LMA at each cell, hence the three entries in the arrays below).
> 
> *** D:  W2_vert -> W3 ***
> [[-1.  1.]
>  [-1.  1.]
>  [ 1. -1.]]
> *** DT: W3 -> W2_vert ***
> [[-1. -1.]
>  [-1. -1.]
>  [ 1.  1.]]
> 
> i.e. it looks fine for D, but the result is clearly wrong for DT, which should be the transpose of D. It does not even contain the same number of +1's and -1's.
> I can't really see what I'm doing wrong here. Furthermore, it does seem to work if I change set order_vertical from 0 to 1.
> 
> I have to admit that I copied the code in build_lma() method from Lawrence, without trying to understand it in detail, but I thought it was sufficiently generic to deal with any UFL forms?
> 
> Thanks a lot for any help, I'm lost as to what is going wrong here,

This is a side-effect of switching coffee optimisations on by default.  The code to build the LMA matrix is slightly hacky.  Basically it takes the kernel used to generate a local element tensor a matrix, and relies on being able to type pun to insert into a Dat of appropriate size.

When running with coffee optimisations on, the matrix assembly kernel is modified to pad the local element tensor to a multiple of the vector length (for aligned accesses).

But the code generation doesn't then do the right thing for inserting into the Dat (there needs to be a pack-unpack stage).

You can get what plausibly looks correct by turning coffee optimisations off (at least the number of 1s and -1s is correct).  I note that there's no way to do this on a per assemble call.

It happens to work for order_vertical=1 because in that case the element tensor is /already/ a multiple of the vector length in size and hence does not get padded.

The right way to do this is probably to have an assemble LMA operator I suppose.


Lawrence
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20150105/9c524b22/attachment.sig>