BUG: f.interpolate(g, subset=...) is unsafe in parallel

[ksagiyam]

f.interpolate(g, subset=...) might update some DoFs of f only on rankA, which does not own them, and not on rankB, which does own them. The values of those DoFs on rankA are then overridden with the old ones on rankB upon halo exchange.

Example

mpiexec -n 2 python ....py -dm_view ascii::ascii_info_detail

from firedrake import *
from pyop2 import op2


distribution_parameters={
    "overlap_type": (DistributedMeshOverlapType.NONE, 0),
    "partitioner_type": "simple",
}
mesh = RectangleMesh(
    2, 1, 2., 1., quadrilateral=True, distribution_parameters=distribution_parameters,
)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
g = Function(V)
f.dat.data_with_halos[:] = 1.
g.dat.data_with_halos[:] = 2.
if mesh.comm.rank == 0:
    subset = op2.Subset(mesh.topology.cell_set, [0])
else:
    subset = op2.Subset(mesh.topology.cell_set, [])
f.interpolate(g, subset=subset)
plex = mesh.topology_dm
plex.viewFromOptions("-dm_view")

Cones:
[0] Max cone size: 4
[0]: 0 <---- 5 (0)
[0]: 0 <---- 6 (0)
[0]: 0 <---- 7 (0)
[0]: 0 <---- 8 (0)
[0]: 5 <---- 1 (0)
[0]: 5 <---- 2 (0)
[0]: 6 <---- 2 (0)
[0]: 6 <---- 4 (0)
[0]: 7 <---- 4 (0)
[0]: 7 <---- 3 (0)
[0]: 8 <---- 3 (0)
[0]: 8 <---- 1 (0)
[1] Max cone size: 4
[1]: 0 <---- 5 (-1)
[1]: 0 <---- 6 (0)
[1]: 0 <---- 7 (0)
[1]: 0 <---- 8 (0)
[1]: 5 <---- 2 (0)
[1]: 5 <---- 1 (0)
[1]: 6 <---- 2 (0)
[1]: 6 <---- 4 (0)
[1]: 7 <---- 4 (0)
[1]: 7 <---- 3 (0)
[1]: 8 <---- 3 (0)
[1]: 8 <---- 1 (0)
coordinates with 1 fields
  field 0 with 2 components
Process 0:
  (   1) dof  2 offset   0 0. 0.
  (   2) dof  2 offset   2 0. 1.
  (   3) dof  2 offset   4 1. 0.
  (   4) dof  2 offset   6 1. 1.
Process 1:
  (   1) dof  2 offset   0 1. 0.
  (   2) dof  2 offset   2 1. 1.
  (   3) dof  2 offset   4 2. 0.
  (   4) dof  2 offset   6 2. 1.

PetscSF Object: 2 MPI processes
  type: basic
  [0] Number of roots=9, leaves=3, remote ranks=1
  [0] 3 <- (1,1)
  [0] 4 <- (1,2)
  [0] 7 <- (1,5)
  [1] Number of roots=9, leaves=0, remote ranks=0

plex points:

rank 0:
2----6---(4)
|         |
5    0   (7)
|         |
1----8---(3)

rank 1:
          2----6----4
          |         |
          5    0    7
          |         |
          1----8----3    () = ghost

mesh.cell_closure:

0
[[1 2 3 4 5 7 8 6 0]]
1
[[2 1 4 3 5 7 6 8 0]]

mesh.coordinates.dat.data_with_halos:

0
[[0. 0.]
 [0. 1.]
 [1. 1.]
 [1. 0.]]
1
[[1. 0.]
 [1. 1.]
 [2. 1.]
 [2. 0.]]

Result of f.interpolate(g, subset=subset):

Before halo exchange (f.dat._data):

0
[2. 2. 2. 2.]
1
[1. 1. 1. 1.]

After halo exchange (f.dat.data_with_halos):

0
[2. 2. 1. 1.]  # <- 2nd and 3rd entries have been overridden with the old values.
1
[1. 1. 1. 1.]

[ksagiyam]

It looks that the workaround is to assemble the interpolation matrix A. assemble(action(A, g)) then computes owned DoFs first.

from firedrake import *
from pyop2 import op2


distribution_parameters={
    "overlap_type": (DistributedMeshOverlapType.NONE, 0),
    "partitioner_type": "simple",
}
mesh = RectangleMesh(
    2, 1, 2., 1., quadrilateral=True, distribution_parameters=distribution_parameters,
)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
g = Function(V)
f.dat.data_with_halos[:] = 1.
g.dat.data_with_halos[:] = 2.
if mesh.comm.rank == 0:
    subset = op2.Subset(mesh.topology.cell_set, [0])
else:
    subset = op2.Subset(mesh.topology.cell_set, [])
v0 = Coargument(V.dual(), 0)
v1 = TrialFunction(V)
interp = interpolate(v1, v0, subset=subset)
A = assemble(interp)
f = assemble(action(A, g))

f.dat._data:

0
[2. 2. 0. 0.]
1
[2. 2. 0. 0.]

f.dat.data_with_halos:

0
[2. 2. 2. 2.]
1
[2. 2. 0. 0.]

[connorjward]

f.interpolate(g, subset=...) might update some DoFs of f only on rankA, which does not own them, and not on rankB, which does own them. The values of those DoFs on rankA are then overridden with the old ones on rankB upon halo exchange.

Is this not expected behaviour? If you are writing to ghost values I wouldn't expect them to propagate unless you are INC-ing into the result. I would think that subset should be 'collective' in that if ghost DoFs are addressed then the equivalent owned DoFs on another process should be too.

[ksagiyam]

Ranks 0 and 1 agree on what subset is in the above; it can be obtained from the global expression conditional(0 < x < 1, 1, 0).

[connorjward]

Ahh because this is the cell subset not the node subset.

I can see now why assembling the matrix does the right thing. The long term proper fix would be to do interpolation point/node-wise I expect. David has mentioned in the past that the current approach is quite unpleasant.