On 22 Dec 2014, at 14:15, Eike Mueller <eike.h.mueller@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