BUG: Parallel assembly of expressions involving SpatialCoordinates, on submeshes, can give the wrong answer

[ABaierReinio]

There seems to be a bug where submesh assembly is not done correctly in parallel. See the MFE below.

I compute two quantities, qty_1 and qty_2, by assembling two UFL forms over the submesh. The UFL forms are the same, but for qty_1 it uses the parent SpatialCoordinates, and for qty_2 the submesh SpatialCoordinates.

In serial this gives the same answer, the script output is: qty_1 is 2.83297e+00, qty_2 is 2.83297e+00
In parallel (8 cores) it gives: qty_1 is 2.61235e+00, qty_2 is 2.83297e+00

So it seems that assembling parent mesh SpatialCoordinates may be bugged in parallel.

MFE:

from firedrake import *

# Create the unlabelled mesh (unit square)
mesh = UnitSquareMesh(4, 4)
x_sc, y_sc = SpatialCoordinate(mesh)

# Label the mesh, creating the left-half, right-half and interface
DG0 = FunctionSpace(mesh, "DG", 0)
DGT0 = FunctionSpace(mesh, "DGT", 0)

tag_L = Function(DG0).interpolate(conditional(lt(x_sc, 0.5), 1, 0))                     # Left-half, of dimension 2
tag_R = Function(DG0).interpolate(conditional(gt(x_sc, 0.5), 1, 0))                     # Right-half, of dimension 2
tag_G = Function(DGT0).interpolate(conditional(lt(abs(x_sc - 0.5), 1.0e-8), 1, 0))      # Interface, of dimension 1

label_L = 1
label_R = 2
label_G = 3

mesh = RelabeledMesh(mesh, [tag_L, tag_R, tag_G], [label_L, label_R, label_G])
x_sc, y_sc = SpatialCoordinate(mesh)

# Create the submeshes
mesh_L = Submesh(mesh, 2, label_L)      # Left-half mesh, of dimension 2
mesh_R = Submesh(mesh, 2, label_R)      # Right-half mesh, of dimension 2
mesh_G = Submesh(mesh, 1, label_G)      # Interface mesh, of dimension 1

# Create the measure on the interface
metadata = {"quadrature_degree": 15}
ds_G = Measure("dx", domain=mesh_G, intersect_measures=[ds(mesh_L), ds(mesh_R), dx(domain=mesh)], metadata=metadata)

# Create a UFL expression using the parent SpatialCoordinates
ufl_expr_1 = exp(x_sc + y_sc)

# Same expression, replacing the parent SpatialCoordinates with the interface SpatialCoordinates
ufl_expr_2 = replace(ufl_expr_1, {SpatialCoordinate(mesh) : SpatialCoordinate(mesh_G)})

# Integrate over the interface, and print the result
qty_1 = assemble(ufl_expr_1 * ds_G)
qty_2 = assemble(ufl_expr_2 * ds_G)

PETSc.Sys.Print("qty_1 is %.5e, qty_2 is %.5e" % (qty_1, qty_2))

# Running in parallel (8 cores) gives:  qty_1 is 2.61235e+00, qty_2 is 2.83297e+00
# Running in serial gives:              qty_1 is 2.83297e+00, qty_2 is 2.83297e+00

[ABaierReinio]

Update: The correct answer is obtained in parallel by using (as suggested by @pbrubeck )

ds_G = Measure("dx", domain=mesh_G, intersect_measures=[ds(mesh_L), ds(mesh_R), dS(domain=mesh)], metadata=metadata)

instead of

ds_G = Measure("dx", domain=mesh_G, intersect_measures=[ds(mesh_L), ds(mesh_R), dx(domain=mesh)], metadata=metadata).

Maybe we should throw an error if the user tries to do the latter, since it involves the intersection of different-dimensional measures?