[firedrake] Iteration space in local mass matrix assembly

[Lawrence Mitchell]

On 06/11/14 15:24, Eike Mueller wrote:
> Dear firedrakers,
> 
> when I store the local DG mass matrix in a PyOP2 dat of 3x3 matrices by
> hand, why do I need to specify the op2.i[0] iteration space as in here

This is part of the PyOP2 kernel API, you can read about it in gory
detail here http://op2.github.io/PyOP2/kernels.html#kernel-api.
Essentially the kernel that FFC generates writes to a "local iteration
space" if you don't indicate this in the map accessor then we pass a
pointer instead.

> V = FunctionSpace(mesh,'DG',1)
> V_DG0 = FunctionSpace(mesh,'DG',0)
> 
> u = TestFunction(V)
> v = TrialFunction(V)
> 
> mass = u*v*dx
> 
> mass_kernel = compile_form(mass, 'mass')[0][6]
> mass_matrix = Function(V_DG0, val=op2.Dat(V_DG0.node_set**(3*3)))
> 
> op2.par_loop(mass_kernel,
>              mass_matrix.cell_set,
>              mass_matrix.dat(op2.INC,
> mass_matrix.cell_node_map()[op2.i[0]]),
>            
>  mesh.coordinates.dat(op2.READ,mesh.coordinates.cell_node_map(),flatten=True),
>            
>  mesh.coordinates.dat(op2.READ,mesh.coordinates.cell_node_map(),flatten=True))
> 
> When I print out mass_kernel.code I get 
> 
> static inline void mass_cell_integral_0_otherwise (double A[3][3] ,
> double **vertex_coordinates , double **w0 ) { 
> […]
> }
> 
> i.e. on each cell the kernel gets passed a 3x3 matrix (= array with 9
> elements). If I leave out the "[op2.i[0]]” I only get zeros in mass_matrix.
> Also, why do the coordinates have to be passed with flatten=True keyword?

Because FFC does not distinguish between mixed and vector finite
elements, it always generates code that accesses:

x0 x1 x2 y0 y1 y2 z0 z1 z2, etc...

I.e. it expects all the x components of the vector to appear before all
the y components and so on.

However, we store the data as:

x0 y0 z0 x1 y1 z1 x2 y2 z2

If you pass flatten=False (the default)  the generated code passes a
pointer to the first "vector" entry in the Dat.  i.e. you can access

x0 y0 z0 x1 y1 z1 x2 y2 z2
^        ^        ^
|        |        |
w0[0]    w0[1]    w0[2]

and then y0 is accessed as w0[0][1]

FFC wants y0 to be at w0[1][0] so flatten=True tells PyOP2 to generate
code for pointers to every entry.  Diagram also on the referenced
webpage above.

Cheers,

Lawrence