From 8e93fc17d8759a7af27034dbca38ba22d84bd408 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 5 Sep 2026 23:08:01 +1000 Subject: [PATCH] fix: infer mesh cell family collectively on empty MPI ranks DMPlexIsSimplex returns false on an empty partition. Using that local value for coordinate FE construction mixed simplex and tensor bases across the same communicator, consuming different PETSc message tags. A subsequent mesh HDF5 labelsLoad then blocked in PetscSFSetUp_Basic/MPI_Waitall. Gather cell-family decisions from populated ranks before FE setup and use the same classification for element metadata. Preserve the constructor hint for an entirely empty mesh; do not change SUPG algorithms, solver tolerances or MPI providers. Add single-cell triangle/tet/quad/hex regressions that require empty ranks and subsequently load another mesh, checking P2 volume and boundary integrals. The triangle regression failed before the fix. The original SUPG migration/partition sequence plus all four regressions passed on eight Mac ranks: 15 tests in 26.05 s, 32.30 s including launcher, 2.74 GiB peak process-tree RSS. The pre-fix sequence hung and the mesh-only reproducer also hung. Gadi rerun remains pending. --- docs/developer/subsystems/meshing.md | 20 ++++++- .../discretisation/discretisation_mesh.py | 16 ++++-- .../test_0781_empty_rank_mesh_sequence.py | 53 +++++++++++++++++++ 3 files changed, 85 insertions(+), 4 deletions(-) create mode 100644 tests/parallel/test_0781_empty_rank_mesh_sequence.py diff --git a/docs/developer/subsystems/meshing.md b/docs/developer/subsystems/meshing.md index dfcc44416..be83085c1 100644 --- a/docs/developer/subsystems/meshing.md +++ b/docs/developer/subsystems/meshing.md @@ -33,6 +33,24 @@ The meshing subsystem handles computational mesh generation and manipulation for - QuadBox / HexBox # Structured meshes ``` +## Empty MPI Partitions + +Use `mesh.isSimplex` for the mesh-wide cell family. PETSc's +`mesh.dm.isSimplex()` is a rank-local query and returns `False` on a rank +with no cells, even when the distributed mesh consists of triangles or +tetrahedra. UW3 infers the family collectively from populated ranks before +constructing coordinate finite elements or element metadata. + +This is also an MPI correctness requirement: constructing simplex and +tensor-product coordinate elements on the same communicator consumes +different PETSc message tags. The first mesh may appear to construct +successfully, but a later HDF5 boundary-label load can deadlock in +`PetscSFSetUp_Basic`. This failure does not require a transport solver. + +`tests/parallel/test_0781_empty_rank_mesh_sequence.py` covers triangles, +tetrahedra, quadrilaterals and hexahedra with deliberately empty partitions, +then verifies a second mesh's volume and boundary integrals using P2 data. + ## Documentation Needs ### Critical Gaps @@ -71,4 +89,4 @@ This section needs: --- -*This document serves as a placeholder for comprehensive meshing system documentation.* \ No newline at end of file +*This document serves as a placeholder for comprehensive meshing system documentation.* diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d22d8b190..2dab069a3 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -589,9 +589,19 @@ def __init__( self._setup_symbolic_coordinates(coordinate_system_type) try: - self.isSimplex = self.dm.isSimplex() + local_simplex = self.dm.isSimplex() except: - self.isSimplex = simplex + local_simplex = simplex + + # DMPlexIsSimplex is rank-local and returns False on empty ranks. + # Coordinate FE construction must use one cell family everywhere; + # mixed simplex/tensor construction desynchronises PETSc MPI tags. + cell_start, cell_end = self.dm.getHeightStratum(0) + cell_families = self.dm.comm.tompi4py().allgather( + local_simplex if cell_end > cell_start else None + ) + populated_families = [family for family in cell_families if family is not None] + self.isSimplex = all(populated_families) if populated_families else simplex # Using WeakValueDictionary to prevent circular references self._vars = weakref.WeakValueDictionary() @@ -650,7 +660,7 @@ class ElementInfo: entities: tuple face_entities: tuple - if self.dm.isSimplex(): + if self.isSimplex: if self.dim == 2: self._element = ElementInfo("triangle", (1, 3, 3), (0, 1, 2)) else: diff --git a/tests/parallel/test_0781_empty_rank_mesh_sequence.py b/tests/parallel/test_0781_empty_rank_mesh_sequence.py new file mode 100644 index 000000000..68e7f2a54 --- /dev/null +++ b/tests/parallel/test_0781_empty_rank_mesh_sequence.py @@ -0,0 +1,53 @@ +"""Empty partitions must not change the cell family or poison later mesh loads. + +The original failure appeared in SUPG test 1077 after the empty-partition +test: rank-local DMPlexIsSimplex returned False on empty ranks, so the +coordinate FE consumed different COMM_WORLD tags. The next HDF5 label load +then hung in PetscSFSetUp_Basic. No transport solve is needed to reproduce it. +""" + +import numpy as np +import pytest +from petsc4py import PETSc + +import underworld3 as uw + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, + pytest.mark.tier_b, pytest.mark.timeout(120)] + + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("simplex", [True, False]) +def test_cell_family_and_next_mesh_load_on_empty_ranks(dim, simplex): + if simplex: + coords = np.vstack([np.zeros(dim), np.eye(dim)]) + cells = np.arange(dim + 1, dtype=PETSc.IntType).reshape(1, -1) + dm = PETSc.DMPlex().createFromCellList(dim, cells, coords) + else: + dm = PETSc.DMPlex().createBoxMesh([1] * dim, simplex=False) + first = uw.discretisation.Mesh(dm, simplex=simplex, qdegree=3) + start, end = first.dm.getHeightStratum(0) + counts = uw.mpi.comm.allgather(end - start) + assert min(counts) == 0 and sum(counts) == 1, counts + + families = uw.mpi.comm.allgather(first.isSimplex) + assert families == [simplex] * uw.mpi.size, families + expected = { (2, True): "triangle", (3, True): "tetrahedron", + (2, False): "quadrilateral", (3, False): "hexahedron" } + elements = uw.mpi.comm.allgather(first._element.type) + assert elements == [expected[dim, simplex]] * uw.mpi.size, elements + + # This public constructor exercises the HDF5 label SF exchange that hung. + second = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, + cellSize=0.25, qdegree=3, regular=False, + ) + field = uw.discretisation.MeshVariable("T", second, 1, degree=2) + field.array[:, 0, 0] = 1.0 + volume = uw.maths.Integral(second, field.sym[0]).evaluate() + assert np.isclose(volume, 1.0, rtol=1e-12, atol=1e-12), volume + for boundary in second.boundaries: + if boundary.name in ("Null_Boundary", "All_Boundaries"): + continue + area = uw.maths.BdIntegral(second, field.sym[0], boundary=boundary.name).evaluate() + assert np.isclose(area, 1.0, rtol=1e-12, atol=1e-12), (boundary.name, area)