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

[Fabio Luporini]

I recently (two-three weeks ago) fixed a bug concerning padding, packing
and unpacking. Are you sure coffee is up-to-date? If so, I'll have a look
at this in the next few days (will be back in London tomorrow) and be back
to you once fixed.

Thanks!

-- Fabio

2015-01-05 12:45 GMT+01:00 Mitchell, Lawrence <
lawrence.mitchell at imperial.ac.uk>:

>
> 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 --------------
HTML attachment scrubbed and removed