BUG: interpolating from a Submesh onto its parent zeroes interface dofs in parallel

Aperta
  • Qualcosa non va?
#5,520 2 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
4/5
Tempo stimato
3-5 giorni
Idoneità per principianti
15/100
Tipo di issue
Bug
Chiarezza
Abbastanza chiara
Stato di attività
Attiva
Stack tecnologico
hpc, python
Ambito
distributed-systems

Direzione di ricerca

Start from the reproduction script submesh_interp_repro.py and run it under mpiexec -n 2, where 239 submesh dofs on the x = 0.5 interface read back as 0. The cause is in how the interpolation writes into ghost entries of the parent space when submesh extraction changes dof ownership, so the entry points are the subset= and allow_missing_dofs=True paths of interpolate. Done means the script reports 0 wrong dofs on 2, 3, 4 and 6 ranks, with the result independent of partition.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

bug difficult

Describe the bug
Interpolating a Function from a Submesh onto a space on the parent mesh, with interpolate(..., subset=parent.cell_subset(...)) or with allow_missing_dofs=True, is exact in serial. In parallel, some dofs on the boundary of the submesh come back as 0. Which dofs, and how many, depends on the number of ranks and the partition.

As far as we can tell, the cause is that extracting the submesh reassigns dof ownership: a dof on the submesh boundary can be owned by one rank on the submesh and by another on the parent. The rank that owns a submesh cell around the dof writes the correct value, but on the parent that entry may be a ghost. The parent's owner, which only has cells outside the submesh around that dof, writes nothing. Since the interpolation writes rather than accumulates, the owner's 0 then overwrites the ghost value in the halo exchange.

Steps to Reproduce

  1. Script (submesh_interp_repro.py): a 12x12x12 unit cube, submesh = the half x < 0.5, P2 field 1 + x + y + z.
import firedrake as fd
from firedrake.petsc import PETSc

mesh = fd.UnitCubeMesh(12, 12, 12)
x = fd.SpatialCoordinate(mesh)
indicator = fd.Function(fd.FunctionSpace(mesh, "DG", 0)).interpolate(fd.conditional(x[0] < 0.5, 1.0, 0.0))
parent = fd.RelabeledMesh(mesh, [indicator], [999])
submesh = fd.Submesh(parent, 3, 999)

V_sub = fd.FunctionSpace(submesh, "CG", 2)
V_parent = fd.FunctionSpace(parent, "CG", 2)
xs = fd.SpatialCoordinate(submesh)
f_sub = fd.Function(V_sub).interpolate(1 + xs[0] + xs[1] + xs[2])

for label, kwargs in [
    ("subset=", {"subset": parent.cell_subset(999)}),
    ("allow_missing_dofs=True", {"allow_missing_dofs": True}),
]:
    f_parent = fd.assemble(fd.interpolate(f_sub, V_parent, **kwargs))
    # Read the parent values back onto the submesh and count the nodes that changed.
    round_trip = fd.assemble(fd.interpolate(f_parent, V_sub))
    wrong = fd.Function(V_sub).interpolate(fd.conditional(abs(round_trip - f_sub) > 1e-10, 1.0, 0.0))
    with wrong.dat.vec_ro as v:
        n_wrong = int(v.sum())
    PETSc.Sys.Print(f"{mesh.comm.size} ranks, {label}: {n_wrong} of {V_sub.dim()} submesh dofs wrong")
  1. Command executed:
for n in 1 2 3 4 6; do mpiexec -n $n python submesh_interp_repro.py; done

Output:

1 ranks, subset=: 0 of 8125 submesh dofs wrong
1 ranks, allow_missing_dofs=True: 0 of 8125 submesh dofs wrong
2 ranks, subset=: 239 of 8125 submesh dofs wrong
2 ranks, allow_missing_dofs=True: 239 of 8125 submesh dofs wrong
3 ranks, subset=: 0 of 8125 submesh dofs wrong
3 ranks, allow_missing_dofs=True: 0 of 8125 submesh dofs wrong
4 ranks, subset=: 0 of 8125 submesh dofs wrong
4 ranks, allow_missing_dofs=True: 0 of 8125 submesh dofs wrong
6 ranks, subset=: 7 of 8125 submesh dofs wrong
6 ranks, allow_missing_dofs=True: 7 of 8125 submesh dofs wrong

On 2 ranks, all 239 wrong dofs read back as exactly 0, and they all lie on the x = 0.5 interface. The read-back step is not the culprit: after the accumulate-and-average transfer described below, the same read-back reports 0 wrong dofs.

Expected behavior
The result should not depend on the partition: every dof in the closure of the submesh cells should carry the submesh value on any number of ranks, as it does in serial.

Error message
None: the wrong values are silent.

Environment:

  • OS: macOS (Darwin 25.3.0)
  • Python version: 3.14.6
  • Firedrake 2025.10.4, petsc4py 3.24.3. Also reproduced on Firedrake 2026.4.0.
  • MPI: Open MPI 5.0.9
  • Relevant packages from pip list:
firedrake        2025.10.4
firedrake-fiat   2025.10.1
fenics-ufl       2025.2.1
loopy            2025.2
petsc4py         3.24.3
mpi4py           4.1.2
numpy            2.4.4
  • No relevant environment variables or modifications

Additional Info
Our workaround: interpolate as usual into a zeroed target, interpolate a field of ones the same way to count how many times each dof is written, sum both onto the owners with dat.local_to_global_begin/end(Access.INC), then divide the values by the counts wherever the count is positive. Every contribution to a given dof carries the same value, so this gives the exact result on any partition. It would be good if the interpolation onto the parent were partition-independent on its own, or, if writing into ghost entries is expected here, if this were documented.

Lingua principale
Python
Stelle
684
Fork
204
Merge medio
4g 14h
PR unite (30g)
46

Preparare l'ambiente

  • Nessun Dockerfile né file Docker Compose
  • Ha un modello di pull request
  • Nessuna guida per i contributori

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di firedrakeproject/firedrake

Tutte le issue di firedrakeproject/firedrake

Issue simili

Altre issue su Python

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.