No link between mesh and DM coordinates breaks VertexOnlyMesh

[ReubenHill]

Failing Script

Consider the below (as reported by @JDBetteridge )

from firedrake import *

comm = COMM_WORLD
pcg = PCG64(seed=333)
rng = RandomGenerator(pcg)
dim = 2
deg = 1
Npoints = 20
samples = 10

# Pick a bunch of random points (same set on each rank)
if comm.rank == 0:
    points = rng.uniform(size=(Npoints, dim))
else:
    points = np.zeros((Npoints, dim))
comm.Bcast(points, root=0)

mesh = SquareMesh(32, 32, 2)
# Changing the coordinates breaks things!
mesh.coordinates.dat.data[:, :] -= 0.5

vom = VertexOnlyMesh(mesh, points)

if run on more than one rank, this fails with

Traceback (most recent call last):
  File "test3.py", line 22, in <module>
    vom = VertexOnlyMesh(mesh, points)
  File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/mesh.py", line 2002, in VertexOnlyMesh
    dmcommon.label_pic_parent_cell_info(swarm, mesh)
  File "firedrake/cython/dmcommon.pyx", line 2710, in firedrake.cython.dmcommon.label_pic_parent_cell_info
    parent_cell_num, reference_coord = parentmesh.locate_cell_and_reference_coordinate(swarm_coords[i])
TypeError: an integer is required
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1

Thoughts

My belief is that the line

mesh.coordinates.dat.data[:, :] -= 0.5

is naughty.

The mesh is sitting on top of a DMPlex which has an implicit coordinates field (made when _from_cell_list in firedrake/mesh.py is called) which I don’t think is updated.

This break things like creating a VertexOnlyMesh which relies on that underlying-plex coordinate field to decide on parallel distribution of the DMSwarm and whether-or-not the DMSwarm coordinates (the VertexOnlyMesh vertices) are within the coordinates field of the mesh:

firedrake/firedrake/mesh.py

Line 2119 in 9e7c358

swarm.setPointCoordinates(coords, redundant=False, mode=PETSc.InsertMode.INSERT_VALUES)

Specific cause of the observed error

When calling

firedrake/firedrake/cython/dmcommon.pyx

Line 2710 in 9e7c358

parent_cell_num, reference_coord = parentmesh.locate_cell_and_reference_coordinate(swarm_coords[i])

at

firedrake/firedrake/mesh.py

Line 2002 in 9e7c358

dmcommon.label_pic_parent_cell_info(swarm, mesh)

some of the swarm coordinates are no longer inside the mesh as reported by the mesh's coordinates field (mesh.dat.coordinates) despite the DMPlex happily reporting that they are inside the mesh and therefore not filtering them out.
At this point in the code it’s assumed that you’ve already filtered out all out-of-(local)-mesh points at the previously mentioned

firedrake/firedrake/mesh.py

Line 2119 in 9e7c358

swarm.setPointCoordinates(coords, redundant=False, mode=PETSc.InsertMode.INSERT_VALUES)

line.
Therefore

firedrake/firedrake/cython/dmcommon.pyx

Line 2710 in 9e7c358

parent_cell_num, reference_coord = parentmesh.locate_cell_and_reference_coordinate(swarm_coords[i])

returns aNone which is not the PetscInt type that cython is expecting.

Fixes

Quick (specific to this issue)

We could switch to filtering points using locate_cell_and_reference_coordinate before rank locating them relative to the DMPlex

Harder (specific to this issue)

We could go “full firedrake” and, rather than using PETSc’s in built swarm functionality, abandon use of the DMPlex coordinates field altogether.

Ideal (more general)

Find a way to maintain the link between the DMPlex coordinates field and firedrake mesh.coordinates field.

[ReubenHill]

No longer failing