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)