Incorrect assembly of surface integrals when diagonal=True

[pbrubeck]

Consider a form that only depends on the jump of CG(1) test and trial functions.
Calling assemble(a, diagonal=False) gives a matrix of zeros, which is the expected output.
On the other hand, assemble(a, diagonal=True) gives non-zeros. For DG(1), the diagonal is also incorrectly computed.

Here is the minimal failing example:

from firedrake import *

mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
outer_jump = lambda w, n: outer(w('+'), n('+')) + outer(w('-'), n('-'))
a = inner(outer_jump(v, n), outer_jump(u, n)) * dS
A = assemble(a)
print(A.petscmat.view())

diag = assemble(a, diagonal=True)
with diag.dat.vec as x:
    print(x.view())

Giving the output:

Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 0.)  (1, 0.)  (2, 0.) 
row 1: (0, 0.)  (1, 0.)  (2, 0.) 
row 2: (0, 0.)  (1, 0.)  (2, 0.) 
None
Vec Object: 1 MPI processes
  type: seq
0.
2.
0.
None

[wence-]

TSFC produces a kernel to compute the diagonal with this:

        # Diagonal assembly is obtained by using the test indices for
        # the trial space as well.
        a, _ = argument_multiindices
        argument_multiindices = (a, a)

I think this isn't right for interior facets (I am not sure there is something that is correct).

If I just use jump(u), jump(v) for your example, the local tensor is 4x4

In [91]: A
Out[91]: 
array([[ 1.,  0.,  0., -1.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [-1.,  0.,  0.,  1.]])

In [92]: V.interior_facet_node_map().values
Out[92]: array([[1, 0, 2, 1]], dtype=int32)

So you can see that the +1 cancels with the -1 on insertion.

The 4x4 matrix is logically split into blocks:

+/+ | +/-
----+----
-/+ | -/-

I think for the +/+ and -/- one should take the diagonals, but one also needs to take bits from +/- and -/+ and I don't see if that is something obvious.

In the interim, I think we should just raise a NotImplementedError inside TSFC if requesting the diagonal of an interior facet term (so that the error is loud rather than quiet).

[pbrubeck]

Shouldn't the current implementation give the right result for V = FunctionSpace(mesh, "DG", 1)? This is not the case:

Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 0.133975)  (1, -0.5)  (2, 0.5)  (3, -0.133975) 
row 1: (0, -0.5)  (1, 1.86603)  (2, -1.86603)  (3, 0.5) 
row 2: (0, 0.5)  (1, -1.86603)  (2, 1.86603)  (3, -0.5) 
row 3: (0, -0.133975)  (1, 0.5)  (2, -0.5)  (3, 0.133975) 
None
Vec Object: 1 MPI processes
  type: seq
0.633975
2.36603
2.36603
0.633975
None

[wence-]

I thought so, but seemingly not. I got lost debugging the kernels yesterday, but I think the naive "just replace the index that ranges over the trial space with that of the test space" doesn't work right.

[sv2518]

I have just had a go at this too see what happens if we use Slate. I adjusted the script to

from firedrake import *

mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
outer_jump = lambda w, n: outer(w('+'), n('+')) + outer(w('-'), n('-'))
a = inner(outer_jump(v, n), outer_jump(u, n)) * dS

A = assemble(Tensor(a))
A2 = assemble(a)
print(A.petscmat.view())
print(A2.petscmat.view())

diag = assemble(Tensor(a), diagonal=True)
with diag.dat.vec as x:
    print(x.view())

and I get

Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 0.)  (1, 0.) 
row 1: (0, 0.)  (1, 2.)  (2, 0.) 
row 2: (1, 0.)  (2, 0.) 
None
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 0.)  (1, 0.)  (2, 0.) 
row 1: (0, 0.)  (1, 0.)  (2, 0.) 
row 2: (0, 0.)  (1, 0.)  (2, 0.) 
None
Vec Object: 1 MPI processes
  type: seq
0.
2.
0.

I ran this on the branch of this PR here #2292.
So it looks like the assembly of a Slate Tensor on this form for a CG discretisation would be wrong too.
If I do the same for DG and I assemble the diagonal on the Slate tensor of the form I get the correct answer for the diagonal but TSFC and Slate don't assemble the same matrix.

row 0: (0, 0.133975)  (1, -0.5) 
row 1: (0, -0.5)  (1, 1.86603) 
row 2: (2, 1.86603)  (3, -0.5) 
row 3: (2, -0.5)  (3, 0.133975) 
None
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 0.133975)  (1, -0.5)  (2, 0.5)  (3, -0.133975) 
row 1: (0, -0.5)  (1, 1.86603)  (2, -1.86603)  (3, 0.5) 
row 2: (0, 0.5)  (1, -1.86603)  (2, 1.86603)  (3, -0.5) 
row 3: (0, -0.133975)  (1, 0.5)  (2, -0.5)  (3, 0.133975) 
None