diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 42afa698..06374ab1 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -8226,6 +8226,34 @@ def eval_metric(centroids): def eval_metric(centroids): return numpy.asarray(_sampler(metric_sym, centroids)).reshape(-1) + def marking_metric(centroids): + """``eval_metric``, clipped, with "nobody has cells" decided together. + + Returns ``None`` when NO rank owns cells, which the marking loops + read as "mark nothing". The point of routing the emptiness question + through a reduction is that the loops below then have a collective + fact to stop on, instead of each rank deciding for itself whether + to leave — which is the defect this replaced (#512). + + Note on what this does NOT fix. #512 describes the rank-local + ``if cur_h.size:`` guards as skipping a collective, on the grounds + that ``eval_metric`` is ``global_evaluate`` for a field or + expression metric. Measured at np=2, that is not so: with one rank + asleep for 10 s, the other's ``global_evaluate`` returned in 0.06 s, + both for in-domain points and for points that strand outside the + mesh entirely. It does not wait for its peers on these shapes, so a + cell-less rank skipping it does not hang. The guards are routed + through here for uniformity with the stop below, not because they + were hanging. + """ + + # mpi4py allreduce defaults to MPI.SUM, as at the collective stop + # further down. + if uw.mpi.comm.allreduce(int(centroids.shape[0])) == 0: + return None + + return numpy.clip(eval_metric(centroids), 1e-30, None) + def cell_geometry(dm): """Per-cell centroid and size for every cell of a simplicial DM. @@ -8458,8 +8486,8 @@ def _relax_generation(engine_obj, carry, rcarry): n_gen = dim * max_levels for level in range(n_gen): centroids, cur_h, cs = cell_geometry(current_dm) - if cur_h.size: - M = numpy.clip(eval_metric(centroids), 1e-30, None) + M = marking_metric(centroids) + if M is not None and cur_h.size: h_target = 1.0 / numpy.sqrt(M) sel = numpy.where(cur_h > h_target)[0] if node_budget is not None and sel.size > node_budget: @@ -8543,8 +8571,8 @@ def _relax_generation(engine_obj, carry, rcarry): current_dm = base_finest for level in range(n_pass): centroids, _proxy_h, cs = cell_geometry(current_dm) - if centroids.shape[0]: - M = numpy.clip(eval_metric(centroids), 1e-30, None) + M = marking_metric(centroids) + if M is not None and centroids.shape[0]: h_target = 1.0 / numpy.sqrt(M) diameter = edge_split.cell_diameters(current_dm) sel = numpy.where(diameter > h_target)[0] @@ -8639,6 +8667,13 @@ def _relax_generation(engine_obj, carry, rcarry): regions=rcarry) n_gen = dim * max_levels for level in range(n_gen): + # Rank-local marking and a rank-local break, deliberately: this + # is the pure-Python cell-list engine, reached only when the + # native uwnvb transform is absent, and that combination raises + # NotImplementedError at np>1 further up. It is serial by + # construction, so the collective discipline the other engines + # need here (#512) has nothing to protect. Port it before + # porting this loop. centroids, cur_h, cids = nvb.centroids_h() M = numpy.clip(eval_metric(centroids), 1e-30, None) h_target = 1.0 / numpy.sqrt(M) @@ -8696,17 +8731,30 @@ def _relax_generation(engine_obj, carry, rcarry): current_dm = base_finest # base finest, current geometry for level in range(max_levels): centroids, cur_h, cs = cell_geometry(current_dm) - if cur_h.size == 0: - break # Metric M = 1/h_target² at the cell centroids (parent field). # A callable is evaluated directly on the centroids; a field/expr # goes through global_evaluate (parallel) or evaluate (serial). - M = numpy.clip(eval_metric(centroids), 1e-30, None) - h_target = 1.0 / numpy.sqrt(M) + M = marking_metric(centroids) + if M is None: + break # no rank has cells: leaving together - refine = numpy.where(cur_h > h_target)[0] - if refine.size == 0: + if cur_h.size: + h_target = 1.0 / numpy.sqrt(M) + refine = numpy.where(cur_h > h_target)[0] + else: + refine = numpy.empty(0, dtype=int) + + # Collective stop. `refine.size == 0` is a rank-local fact: this + # rank's cells are all fine enough, which says nothing about its + # peers'. Breaking on it alone takes the rank out of the loop for + # good while the others go round again into the DM refinement, + # which IS collective — so the job hangs, and not latently. + # Measured at np=2 with a callable metric demanding h=0.01 inside + # r<0.15 of one corner and nothing elsewhere: the rank holding no + # corner cells left at the first level and `mpirun` reached its + # timeout with neither rank returning. See test_0873. + if uw.mpi.comm.allreduce(int(refine.size)) == 0: if verbose: uw.pprint(0, f"[adapt] level {level}: no cells need refinement") break diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py index af1b9189..427bc34f 100644 --- a/src/underworld3/utilities/reconnect.py +++ b/src/underworld3/utilities/reconnect.py @@ -833,20 +833,37 @@ def _rebuild_point_sf(new, dm, point_map, nroots): # petsc4py will not narrow an index array for us — the graph must arrive as # PETSc's own integer type or `setGraph` raises an unsafe-cast TypeError. new_sf = PETSc.SF().create(comm=dm.comm) - if ilocal is None or not len(ilocal): + + # A rank that shares nothing has no leaves to renumber, but it still has to + # reach the verdict below with its peers — which is why the emptiness test + # is read here and acted on after, rather than returning straight out of it. + shares_nothing = ilocal is None or not len(ilocal) + if shares_nothing: + leaves = local = remote_index = None + deleted_here = False + else: + leaves = np.asarray(ilocal, dtype=np.int64) + local = point_map[leaves - pStart] + remote_index = leaf_new[leaves - pStart] + deleted_here = bool((local < 0).any() or (remote_index < 0).any()) + + # Assert-class: fires only on an internal contract violation. It is still + # raised collectively — one rank raising while its peers carry on into the + # next collective is a hang on top of the error, and the error is the thing + # worth seeing (#512). + offenders = [rank for rank, bad + in enumerate(dm.comm.tompi4py().allgather(deleted_here)) if bad] + if offenders: + raise RuntimeError( + f"reconnect: a shared point was deleted on rank(s) {offenders}. " + "The removal pass must freeze the seam; see remove_vertices.") + + if shares_nothing: new_sf.setGraph(nroots, np.zeros(0, dtype=PETSc.IntType), np.zeros(0, dtype=PETSc.IntType)) _install_point_sf(new, new_sf) return - leaves = np.asarray(ilocal, dtype=np.int64) - local = point_map[leaves - pStart] - remote_index = leaf_new[leaves - pStart] - if (local < 0).any() or (remote_index < 0).any(): - raise RuntimeError( - "reconnect: a shared point was deleted. The removal pass must " - "freeze the seam; see remove_vertices.") - remote = np.empty((len(leaves), 2), dtype=PETSc.IntType) remote[:, 0] = np.asarray(iremote).reshape(-1, 2)[:, 0] remote[:, 1] = remote_index diff --git a/tests/parallel/ptest_0845_relax_pinned_band_parallel.py b/tests/parallel/ptest_0845_relax_pinned_band_parallel.py index ed37c053..097bfa6c 100644 --- a/tests/parallel/ptest_0845_relax_pinned_band_parallel.py +++ b/tests/parallel/ptest_0845_relax_pinned_band_parallel.py @@ -61,7 +61,21 @@ def _pinned_indices(mesh, name): return np.asarray(iset.getIndices(), dtype=np.int64) - vS -def test_pinned_set_is_partition_independent(): +def test_every_rank_pins_every_band_vertex_it_holds(): + """A vertex in the band is pinned on EVERY rank that holds it. + + The previous form of this test allgathered the pinned coordinates, took the + union, and compared ``len(union)`` with ``allreduce(len(union), MAX)``. The + union is built from the same gathered list on every rank, so those two + numbers are the same by construction and the assertion could not fail + whatever the pinning did (#512). + + What can fail is this: a rank-local pin misses a vertex its peer pinned, so + the owner is free to move a vertex the seam expects to stay. Every rank + therefore has to pin every band vertex present in its own coordinate array, + not merely some of them. + """ + mesh, surface = _fixture() name = mesh.label_interface_band(surface, offset=0.0, halo=1) X = _coords(mesh) @@ -70,13 +84,21 @@ def test_pinned_set_is_partition_independent(): # Compare the pinned COORDINATES, not counts: a shared vertex is held by # every rank on the seam, so a count double-counts it and would mask exactly # the defect this test exists to catch. - local = {(round(float(x), 12), round(float(y), 12)) for x, y in X[idx]} - gathered = uw.mpi.comm.allgather(local) - union = set().union(*gathered) + def key(x, y): + return (round(float(x), 12), round(float(y), 12)) + + local = {key(x, y) for x, y in X[idx]} + union = set().union(*uw.mpi.comm.allgather(local)) + + assert union, "nothing pinned; the fixture is not exercising the band" + + held_here = {key(x, y) for x, y in X} + missing = (union & held_here) - local - total = uw.mpi.comm.allreduce(len(union), op=MPI.MAX) - assert len(union) == total - assert total > 0, "nothing pinned; the fixture is not exercising the band" + assert not missing, ( + f"rank {uw.mpi.rank} holds {len(missing)} band vertex/vertices that a " + f"peer pinned and it did not: {sorted(missing)[:5]} — a rank-local pin " + "lets the owner move a vertex the seam expects to stay") def test_pinned_vertices_including_shared_ones_do_not_move(): diff --git a/tests/parallel/test_0873_adapt_collective_stop_mpi.py b/tests/parallel/test_0873_adapt_collective_stop_mpi.py new file mode 100644 index 00000000..95cf494e --- /dev/null +++ b/tests/parallel/test_0873_adapt_collective_stop_mpi.py @@ -0,0 +1,104 @@ +"""The adapt level loop must stop on a collective fact, not a rank-local one (#512). + +Each level of the marking loop ends with a DM refinement, which is collective. +Deciding to leave that loop is therefore a decision every rank has to take +together. The sbr path used to break on ``refine.size == 0`` — "this rank's own +cells are all fine enough" — which says nothing about anyone else's. The rank +that finished first left for good and its peers went round again into the +refinement and waited. + +It is not latent. Measured at np=2 on `development` before the fix: a callable +metric demanding h = 0.01 within r < 0.15 of one corner and nothing elsewhere, +on a 132-cell-per-rank base, hung with neither rank returning; `mpirun` reached +its timeout. + +Run under MPI:: + + mpirun -np 2 python -m pytest --with-mpi \ + tests/parallel/test_0873_adapt_collective_stop_mpi.py + +A callable metric is used deliberately, so that nothing here depends on how +``global_evaluate`` behaves — it is evaluated rank-locally. #512 attributes the +hazard to ``eval_metric`` being collective; measured at np=2 that is not so, +with one rank asleep for 10 s the other's ``global_evaluate`` returned in +0.06 s, for in-domain and for stranded points alike. The collective in this +loop is the refinement, not the metric. +""" + +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [ + pytest.mark.mpi(min_size=2), + pytest.mark.timeout(300), + pytest.mark.level_2, + pytest.mark.tier_a, +] + +ENGINES = ["sbr", "edge_split", "nvb"] + + +def _corner_metric(points): + """Fine in one corner, already satisfied everywhere else. + + This is what splits the ranks: whichever rank holds no corner cells has + nothing to refine at the first level, while its peer has work for several. + """ + + p = np.asarray(points) + r = np.sqrt((p[:, 0] - 0.02) ** 2 + (p[:, 1] - 0.02) ** 2) + + return 1.0 / np.where(r < 0.15, 0.01, 1.0) ** 2 + + +def _base_mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + cellSize=0.2, + regular=False, + qdegree=2, + refinement=1, + ) + + +def test_premise_the_metric_splits_the_ranks(): + """Test premise: without a rank that runs out of work, nothing here is tested.""" + + mesh = _base_mesh() + cs, ce = mesh.dm.getHeightStratum(0) + + centroid = np.asarray( + [mesh.dm.computeCellGeometryFVM(c)[1] for c in range(cs, ce)] + ).reshape(ce - cs, -1)[:, : mesh.dim] + + wants_refining = int((_corner_metric(centroid) > 1.0).sum()) + per_rank = uw.mpi.comm.allgather(wants_refining) + + assert sum(per_rank) > 0, f"no rank has anything to refine: {per_rank}" + assert min(per_rank) == 0, ( + f"every rank has work at np={uw.mpi.size} ({per_rank}) — this test " + "cannot detect the rank-local stop it exists for" + ) + + +@pytest.mark.parametrize("engine", ENGINES) +def test_adapt_returns_when_one_rank_runs_out_of_work(engine): + """The whole assertion is that this returns. The failure mode is a hang.""" + + mesh = _base_mesh() + + child = mesh.adapt(_corner_metric, max_levels=4, engine=engine) + + c0, c1 = child.dm.getHeightStratum(0) + total = uw.mpi.comm.allreduce(c1 - c0) + base = uw.mpi.comm.allreduce( + mesh.dm.getHeightStratum(0)[1] - mesh.dm.getHeightStratum(0)[0] + ) + + assert total > base, ( + f"{engine}: the corner metric demands refinement, but the child has " + f"{total} cells against the base's {base}" + )