BUG: IndexError in _pic_swarm_in_mesh on MPI ranks with empty cell partition

[hardik-corintis]

Describe the bug

_pic_swarm_in_mesh raises IndexError: index -1 is out of bounds for axis 0 with size 0
in parallel when any MPI rank owns zero cells of the parent mesh. This surfaces when
interpolating a Function defined on a Submesh back onto its parent mesh.

parent_cell_nums_local contains -1 sentinels for query points whose owning cell lives
on a different rank. On ranks with an empty cell partition, cell_closure has shape (0, k).
NumPy rejects -1 as an axis-0 index on a zero-row array.

Steps to Reproduce

1. Save the script below as mfe.py
2. mpiexec -n 2 python mfe.py

from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI

comm = COMM_WORLD

PETSc.Sys.Print("=" * 60)
PETSc.Sys.Print(f"Running on {comm.size} MPI rank(s)")
PETSc.Sys.Print("=" * 60)

PETSc.Sys.Print("\n[1/5] Building UnitCubeMesh(8, 8, 8) ...")
mesh = UnitCubeMesh(8, 8, 8)
x, y, z = SpatialCoordinate(mesh)

PETSc.Sys.Print("[2/5] Labelling x=1 face as tag 999 via RelabeledMesh ...")
V_marker = FunctionSpace(mesh, "HDiv Trace", 0)
facet_indicator = Function(V_marker).interpolate(
    conditional(x > 0.999, 1.0, 0.0)
)
mesh = RelabeledMesh(mesh, [facet_indicator], [999])
x, y, z = SpatialCoordinate(mesh)

PETSc.Sys.Print("[3/5] Creating Submesh of x=1 face ...")
surface_mesh = Submesh(mesh, mesh.topological_dimension - 1, 999)
xs, ys, zs = SpatialCoordinate(surface_mesh)

n_cells_local = mesh.cell_set.size
n_surface_cells_local = surface_mesh.cell_set.size
PETSc.Sys.Print(
    f"    rank {comm.rank}: parent cells = {n_cells_local}, "
    f"surface cells = {n_surface_cells_local}"
)

PETSc.Sys.Print("[4/5] Interpolating f_sub = y+z from surface onto parent mesh ...")
V_sub = FunctionSpace(surface_mesh, "CG", 1)
V_parent = FunctionSpace(mesh, "CG", 1)
f_sub = Function(V_sub).interpolate(ys + zs)
f_parent = Function(V_parent)
f_parent.interpolate(f_sub, allow_missing_dofs=True)  # crashes here without fix

PETSc.Sys.Print("[5/5] Sanity check ...")
max_val = comm.allreduce(f_parent.dat.data_ro.max(), op=MPI.MAX)
PETSc.Sys.Print(f"    global max |f_parent| = {max_val:.4f}  (expected ~2.0 at x=1, y=1, z=1)")

PETSc.Sys.Print("\nDONE — no IndexError raised.")

Expected behavior

The script should complete without error and print:

============================================================
Running on 2 MPI rank(s)
============================================================

[1/5] Building UnitCubeMesh(8, 8, 8) ...
[2/5] Labelling x=1 face as tag 999 via RelabeledMesh ...
[3/5] Creating Submesh of x=1 face ...
    rank 0: parent cells = 1514, surface cells = 0
[4/5] Interpolating f_sub = y+z from surface onto parent mesh ...
[5/5] Sanity check ...
    global max |f_parent| = 2.0000  (expected ~2.0 at x=1, y=1, z=1)

DONE — no IndexError raised.

Error message

Traceback (most recent call last):
  File "mfe.py", line 69, in <module>
    f_parent.interpolate(f_sub, allow_missing_dofs=True)
  ...
  File ".../firedrake/interpolation.py", line 524, in _symbolic_expressions
    vom = VertexOnlyMesh(
  File ".../firedrake/mesh.py", line 3682, in VertexOnlyMesh
    swarm, input_ordering_swarm, n_missing_points = _pic_swarm_in_mesh(
  File ".../firedrake/mesh.py", line 3928, in _pic_swarm_in_mesh
    plex_parent_cell_nums = parent_mesh.topology.cell_closure[
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index -1 is out of bounds for axis 0 with size 0

Environment:

• OS: macOS 15.6 (arm64)
• Python version: 3.11.14
• Output of pip list (relevant packages):

  firedrake        2026.5.0.dev0
  fenics-ufl       2025.3.0.dev0
  firedrake-fiat   2026.5.0.dev0
  firedrake-rtree  2026.1.0
  petsc4py         3.26.0.dev0
  petsctools       2026.1.dev0
  loopy            2025.2
  mpi4py           4.1.1
  numpy            2.4.4
  scipy            1.17.1
  h5py             3.16.0
  

[leo-collins]

I am unable to reproduce the MFE on release or main.

_pic_swarm_in_mesh raises IndexError: index -1 is out of bounds for axis 0 with size 0
in parallel when any MPI rank owns zero cells of the parent mesh.

I'm not sure this is true. For example, in test_swarm.py:

firedrake/tests/firedrake/vertexonly/test_swarm.py

Lines 408 to 413 in 547aca2

	@pytest.mark.parallel(nprocs=3) # nprocs > total number of mesh cells
	def test_pic_swarm_in_mesh_2d_3procs():
	test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=False, exclude_halos=True)
	test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=True, exclude_halos=True)
	test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=False, exclude_halos=False)
	test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=True, exclude_halos=False)

Here we call _pic_swarm_in_mesh on a UnitSquareMesh(1, 1) with three MPI ranks - this is guaranteed to have a rank which owns zero cells.

[leo-collins]

I am able to reproduce this by running the MFE on 8 ranks - see my comment here #5078 (comment)

[hardik-corintis]

Thanks a lot for having a look. I was trying to dig more into the issue, but I think the mesh partitioner is not providing the same decomposition for the test case with just 2 ranks

[leo-collins]

Fixed by #5078. Thank you @hardik-corintis !