diff --git a/src/underworld3/function/_dminterp_wrapper.pyx b/src/underworld3/function/_dminterp_wrapper.pyx index 6c10c7ae..966dc0e0 100644 --- a/src/underworld3/function/_dminterp_wrapper.pyx +++ b/src/underworld3/function/_dminterp_wrapper.pyx @@ -178,7 +178,18 @@ cdef class CachedDMInterpolationInfo: &cells_view[0], 1 if hint_authoritative else 0) else: - ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, NULL, 0) + # No hint array to pass (no points, or no cells) — but the POLICY + # still has to be forwarded. Hardcoding 0 here made a rank with zero + # local points disagree with its peers about whether the hint is + # authoritative, and petsc_tools.c takes the DMLocatePoints branch + # when it is not. DMLocatePoints is COLLECTIVE on the mesh DM, so + # that rank blocked inside DMGetBoundingBox -> MPI_Allreduce while + # the others bypassed and ran on to DMSwarmMigrate -> MPI_Comm_dup: + # a deadlock whenever a query set leaves some rank empty (#611). + # The policy is a mesh capability and already agrees across ranks; + # it just has to survive the trip. + ierr = DMInterpolationSetUp_UW(self._ipInfo, dm, 0, 1, NULL, + 1 if hint_authoritative else 0) if ierr != 0: DMInterpolationDestroy(&self._ipInfo) raise RuntimeError(f"DMInterpolationSetUp_UW failed with error {ierr}") diff --git a/src/underworld3/function/petsc_tools.c b/src/underworld3/function/petsc_tools.c index 96936d38..7bec52e7 100644 --- a/src/underworld3/function/petsc_tools.c +++ b/src/underworld3/function/petsc_tools.c @@ -69,7 +69,17 @@ PetscErrorCode DMInterpolationSetUp_UW(DMInterpolationInfo ctx, DM dm, PetscBool PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs)); for (p = 0; p < N; ++p) foundProcs[p] = size; cellSF = NULL; - if (owning_cell && hintAuthoritative) { + /* N == 0 is included deliberately. An empty hint array reaches C as a NULL + `owning_cell`, which flipped this test and sent a rank with no local points + down the DMLocatePoints branch ALONE -- and DMLocatePoints is collective on + the mesh DM communicator (the comment below is right that the + Allreduce(foundProcs) is a COMM_SELF no-op, but DMLocatePoints itself is + not). Three ranks bypassed while the empty one blocked inside + DMGetBoundingBox -> MPI_Allreduce, deadlocking the job (#611). With no + points there is nothing to locate, so bypassing is trivially correct and + the branch now depends only on `hintAuthoritative`, which is a mesh + capability and agrees across ranks. */ + if ((owning_cell || N == 0) && hintAuthoritative) { /* Bypass DMLocatePoints when the caller supplies an AUTHORITATIVE hint (ported from feature/dminterp-bypass-element-check, 17a5a8d). diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 3c184d4a..e463e31c 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -867,6 +867,65 @@ def _count_zero_columns_parallel(P, comm): return nzero +def _repair_zero_columns_parallel(P, coords_f, fine_layout, coords_u, cols_u, + ncomp, comm): + """Give each unreached coarse DOF a nearest-fine-DOF entry (parallel). + + The counterpart of :func:`_repair_zero_columns_serial` for the + cross-partition build: a coarse DOF no fine node reaches (the fine + mesh gathered onto a surgery rank, a placed band, a relaxed child) + leaves an empty column and a singular Galerkin coarse operator. Each + rank reads the empty columns it OWNS off P^T.1, the orphan list is + all-gathered (small), every rank offers its nearest OWNED fine DOF + of the same component, and the rank holding the nearest one sets the + injection entry. Returns ``(P, n_repaired)``; ``P`` is re-assembled. + """ + m4 = comm.tompi4py() + ones_f = P.createVecLeft() + ones_f.set(1.0) + colsum = P.createVecRight() + P.multTranspose(ones_f, colsum) + cstart, _cend = colsum.getOwnershipRange() + zero_local = np.flatnonzero(colsum.array == 0.0) + cstart + ones_f.destroy() + colsum.destroy() + zero = np.concatenate(m4.allgather(zero_local.astype(np.int64))) + if zero.size == 0: + return P, 0 + # global column -> (coarse node coordinate, component) via the gathered cloud + where = {int(cols_u[n, c]): (n, c) + for n in range(cols_u.shape[0]) for c in range(ncomp) + if cols_u[n, c] >= 0} + from scipy.spatial import cKDTree + l2g_f, fstart, fend = fine_layout.l2g, fine_layout.rstart, fine_layout.rend + tree = cKDTree(coords_f) if coords_f.shape[0] else None + offers = [] # (distance, global fine row) + for col in zero.tolist(): + node_c, comp = where[int(col)] + best = (np.inf, -1) + if tree is not None: + k = min(8, coords_f.shape[0]) + d, idx = tree.query(coords_u[node_c], k=k) + for dd, i in zip(np.atleast_1d(d), np.atleast_1d(idx)): + grow = int(l2g_f[int(i) * ncomp + comp]) + if fstart <= grow < fend: # an OWNED fine row + best = (float(dd), grow) + break + offers.append(best) + # the nearest offer across ranks wins each orphan + dist = np.array([o[0] for o in offers]) + rows = np.array([o[1] for o in offers], dtype=np.int64) + all_dist = np.vstack(m4.allgather(dist)) + all_rows = np.vstack(m4.allgather(rows)) + winner = np.argmin(all_dist, axis=0) + for j, col in enumerate(zero.tolist()): + if winner[j] == m4.rank and np.isfinite(all_dist[winner[j], j]): + P.setValues([int(all_rows[winner[j], j])], [int(col)], [1.0], + addv=PETSc.InsertMode.INSERT_VALUES) + P.assemble() + return P, int(zero.size) + + def _assert_no_zero_columns_parallel(P, comm): """Parallel zero-column guard: a coarse DOF with no fine image -> singular Galerkin coarse operator.""" @@ -1388,6 +1447,20 @@ def build(self, solver): if (self.cross_partition == "auto" and _count_zero_columns_parallel(P, comm) > 0): P = _build_crosspart_transfer(*args) + # the same orphan repair the serial path has: a coarse DOF no + # fine node reaches gets its nearest fine DOF as an injection + if _count_zero_columns_parallel(P, comm) > 0: + coords_u, cols_u = _gather_coarse_cloud( + coords[l - 1], maps[l - 1], nc, comm) + P, n_rep = _repair_zero_columns_parallel( + P, coords[l], maps[l], coords_u, cols_u, nc, comm) + if n_rep: + import warnings + warnings.warn( + f"custom_mg: parallel transfer {l - 1}->{l} had " + f"{n_rep} coarse DOF(s) with no fine image " + f"(non-nested levels); repaired by " + f"nearest-fine-DOF injection.") _assert_no_zero_columns_parallel(P, comm) Ps.append(P) else: @@ -1775,20 +1848,19 @@ def build_transfers(solver, field_id=None): solver._record_pc_fallback( "custom_mg.transfer_builder", requested=_b, - installed=f"{_attempts[_i + 1]} (DENSE transfer)", + installed=f"{_attempts[_i + 1]} (local kd-tree RBF)", reason="build_failed", - detail=f"{exc}; the RBF rescue is a performance cliff — " - f"its transfer is dense (nnz/row == n_coarse), see #424") + detail=f"{exc}; the local RBF stencils are wider than the " + f"barycentric ones, so the Galerkin coarse " + f"operators fatten (#429)") warnings.warn( f"custom_mg: {_b} transfer build failed ({exc}); " - f"retrying with the '{_attempts[_i + 1]}' builder, which " - f"has global support and cannot leave a coarse DOF " - f"without a fine image. NOTE the RBF transfer is DENSE " - f"(nnz/row == n_coarse), so the Galerkin coarse operators " - f"are dense too — this rescues correctness but does not " - f"scale. If it fires on a production-sized problem, treat " - f"it as a performance cliff and fix the cause, not the " - f"symptom (#424).") + f"retrying with the '{_attempts[_i + 1]}' builder — the " + f"sparse, linear-exact local kd-tree RBF (#429), whose " + f"kNN stencils reach coarse DOFs the barycentric simplex " + f"does not. Its wider stencils fatten the Galerkin coarse " + f"operators; if this fires routinely, fix the level " + f"geometry rather than live with the fallback.") continue solver._record_pc_fallback( "custom_mg.build", @@ -1957,6 +2029,63 @@ def inject_custom_mg(solver): _install_transfers(solver, Ps, verbose=cfg.get("verbose", False)) +def _colocate_level(coarse_mesh, fine_mesh): + """Redistribute one coarse level so each coarse cell lives on the rank + that holds the fine cells over it (nearest owned fine centroid, by a + global minimum). A placed fine mesh is gathered onto its surgery rank + while the tail stays load-balanced; the transfer then pairs a fine + node with a coarse cell on another rank and coarse DOFs lose every + fine image (measured: 488 of 5614 on the S-fault rig at np=2, and the + repaired transfer does not precondition). Co-resident levels are the + same construction ptest_0004 uses for a reloaded hierarchy. Returns a + new Mesh, or ``coarse_mesh`` itself when nothing moves.""" + import underworld3 as uw + from scipy.spatial import cKDTree + + dm = coarse_mesh.dm + comm = dm.getComm().tompi4py() + if comm.size == 1: + return coarse_mesh + cS, cE = dm.getHeightStratum(0) + cen_c = np.array([dm.computeCellGeometryFVM(c)[1] for c in range(cS, cE)]) + fdm = fine_mesh.dm + fS, fE = fdm.getHeightStratum(0) + # OWNED fine cells only: a ghost cell belongs to another rank + fsf = fdm.getPointSF() + try: + _n, ileaf, _r = fsf.getGraph() + ghost = set(int(q) for q in ileaf) + except (ValueError, TypeError): + ghost = set() + owned_f = [c for c in range(fS, fE) if c not in ghost] + cen_f = (np.array([fdm.computeCellGeometryFVM(c)[1] for c in owned_f]) + if owned_f else np.zeros((0, cen_c.shape[1]))) + # every rank offers its nearest owned fine cell to EVERY coarse centroid + # in the mesh (all-gathered: the coarse levels are small) + cen_all = np.vstack(comm.allgather(cen_c)) + if cen_f.shape[0]: + d_local = cKDTree(cen_f).query(cen_all)[0] + else: + d_local = np.full(len(cen_all), np.inf) + d_all = np.vstack(comm.allgather(d_local)) + owner_all = np.argmin(d_all, axis=0).astype(np.int32) + off = np.cumsum([0] + comm.allgather(len(cen_c))) + assign = owner_all[off[comm.rank]:off[comm.rank + 1]] + if not comm.allreduce(int((assign != comm.rank).sum())): + return coarse_mesh + work = dm.clone() + part = work.getPartitioner() + part.setType(PETSc.Partitioner.Type.SHELL) + order = np.argsort(assign, kind="stable").astype(np.int32) + sizes = np.bincount(assign, minlength=comm.size).astype(np.int32) + part.setShellPartition(comm.size, sizes=sizes, points=order) + work.distribute() + return uw.discretisation.Mesh( + work, simplex=coarse_mesh.dm.isSimplex(), qdegree=coarse_mesh.qdegree, + coordinate_system_type=coarse_mesh.CoordinateSystem.coordinate_type, + boundaries=coarse_mesh.boundaries, verbose=False) + + def adopt_hierarchy(mesh, base_mesh, fac_zone=None, builder=None): """Make ``mesh`` OWN the multigrid hierarchy of ``base_mesh`` — the static coarse tail every solver built on ``mesh`` then drives @@ -1977,8 +2106,15 @@ def adopt_hierarchy(mesh, base_mesh, fac_zone=None, builder=None): # Mesh._adopt_cut_child applies; a plain refined base contributes its # static level wraps (coarsest .. base-finest) own = getattr(base_mesh, "_custom_mg_coarse_meshes", None) - mesh._custom_mg_coarse_meshes = (list(own) + [base_mesh] if own is not None - else list(base_mesh._coarse_level_meshes())) + tail = (list(own) + [base_mesh] if own is not None + else list(base_mesh._coarse_level_meshes())) + # In parallel the placed mesh is gathered onto its surgery rank while + # the tail is load-balanced: co-locate every level with the finest so + # the transfers pair rank-locally (the coarse levels are small; the + # fine mesh and its FAC patch never move). + if mesh.dm.getComm().getSize() > 1: + tail = [_colocate_level(level, mesh) for level in tail] + mesh._custom_mg_coarse_meshes = tail mesh._custom_mg_builder = (builder if builder is not None else getattr(base_mesh, "_custom_mg_builder", "barycentric")) diff --git a/src/underworld3/utilities/fault_contact.py b/src/underworld3/utilities/fault_contact.py index 951b4ad2..73814c39 100644 --- a/src/underworld3/utilities/fault_contact.py +++ b/src/underworld3/utilities/fault_contact.py @@ -1058,18 +1058,40 @@ def fault_normal_traction(solver, boundary, solve_result): return s_coord[order], sig[order] -def fault_pair_jumps(solver, boundary, solve_result): +def fault_pair_jumps(solver, boundary, solve_result, gather=False): """The velocity jump at every coincident pair, from the solve. - Returns ``(coords, jumps, normals)`` on this rank — the pair position, - the full jump vector :math:`v^+ - v^-`, and the fault unit normal — - in any dimension. Reads the composite solution ``solve_result["U"]`` + Returns ``(coords, jumps, normals)`` — the pair position, the full + jump vector :math:`v^+ - v^-`, and the fault unit normal — in any + dimension. Reads the composite solution ``solve_result["U"]`` through the pairing, which is the only correct route: the pair coordinates are identical, so field queries by position see one side only. The tangential part of the jump is the slip (a scalar against the in-fault tangent in 2-D, an in-plane vector in 3-D); the normal part is the leak, held at machine zero by the strong constraint. + + The pairs are rank-local (the split keeps a fault rank-interior), so + a rank without the fault returns empty arrays. ``gather=True`` + all-gathers the three arrays so every rank holds the whole fault — + the form a diagnostic that goes on to make collective calls + (``evaluate``, a write) must use, or the ranks diverge and hang. """ + coords, jumps, normals = _fault_pair_jumps_local(solver, boundary, + solve_result) + if not gather: + return coords, jumps, normals + comm = solver.mesh.dm.comm.tompi4py() + if comm.size == 1: + return coords, jumps, normals + dim = solver.mesh.dim + parts = comm.allgather((np.asarray(coords, dtype=float).reshape(-1, dim), + np.asarray(jumps, dtype=float).reshape(-1, dim), + np.asarray(normals, dtype=float).reshape(-1, dim))) + return tuple(np.vstack([p[k] for p in parts]) for k in range(3)) + + +def _fault_pair_jumps_local(solver, boundary, solve_result): + """The rank-local half of :func:`fault_pair_jumps`.""" dm = solver.dm dim = solver.mesh.dim lsec = dm.getLocalSection() diff --git a/tests/parallel/ptest_0859_fault_network_parallel.py b/tests/parallel/ptest_0859_fault_network_parallel.py new file mode 100644 index 00000000..50881785 --- /dev/null +++ b/tests/parallel/ptest_0859_fault_network_parallel.py @@ -0,0 +1,75 @@ +"""The fault network in parallel: build, junction glue, contact solve at +np=2 — the serial answer, and geometric FMG with NO preconditioner +fallback. + +The placed mesh is gathered onto its surgery rank while the multigrid +tail stays load-balanced; unless the tail is co-located with the finest +level (custom_mg.adopt_hierarchy), the transfer pairs a fine node with a +coarse cell on another rank, coarse DOFs lose every fine image and the +build degrades to the local-RBF rescue (measured on the S-fault rig: +488 orphan coarse DOFs on one level at np=2). This test is the guard. +Run with: + mpirun -np 2 python -m pytest tests/parallel/\ +ptest_0859_fault_network_parallel.py --with-mpi +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.parallel_safe, pytest.mark.level_2, + pytest.mark.tier_b] + +H = 0.03 +WIDTH = 0.04 +# the serial answer of the same network (tests/test_0859, glued): the +# peak tangential jump per piece, read on every rank after an all-gather +SERIAL = {"Main": 0.4373, "Cont": 0.3603, "Splay": 0.1044} + + +def _pieces(): + main = np.column_stack([np.linspace(0.25, 0.50, 12), np.full(12, 0.5)]) + cont = np.column_stack([np.linspace(0.55, 0.75, 9), np.full(9, 0.5)]) + s = np.linspace(0.0, 1.0, 8) + splay = np.column_stack([0.38 + 0.12 * s, 0.5 + 0.18 * s]) + return [("Main", main), ("Cont", cont), ("Splay", splay)] + + +def test_network_glue_solve_np2(): + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=8 * H, + regular=False, refinement=1, qdegree=2) + pieces = _pieces() + net = uw.meshing.FaultNetwork(pieces, hierarchy=[n for n, _p in pieces]) + net.prepare(h=H, ligament=1.0, verbose=False) + net.build(base=base, width=WIDTH, realisation="split", max_levels=1) + + mesh = net.mesh + x, y = mesh.X + v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, + continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = \ + net.junction_patch(eta_0=1.0) + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((2.0 * (y - 0.5), 0.0), wall) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + net.apply(stokes) + info = net.solve(stokes) + assert info.get("converged"), "the contact solve did not converge" + + # geometric FMG survived the partition: nothing was swapped for a + # rescue builder or the default preconditioner + fallbacks = getattr(stokes, "pc_fallbacks", {}) or {} + assert not fallbacks, f"preconditioner fallback recorded: {fallbacks}" + + # the pairs are rank-local; the peak per piece is a global max + local = net.slips(stokes) + comm = mesh.dm.comm.tompi4py() + for name, expected in SERIAL.items(): + peak = comm.allreduce(float(local.get(name, 0.0)), op=max) + assert peak == pytest.approx(expected, rel=2e-2), ( + f"{name}: parallel peak slip {peak:.4f} vs serial {expected}") diff --git a/tests/parallel/test_1076_global_evaluate_empty_rank.py b/tests/parallel/test_1076_global_evaluate_empty_rank.py new file mode 100644 index 00000000..e8ec8a41 --- /dev/null +++ b/tests/parallel/test_1076_global_evaluate_empty_rank.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +"""global_evaluate must not deadlock when a rank receives no query points (#611). + +The cell-location policy (`mesh._hint_is_authoritative`) decides whether the +barycentric hint may bypass PETSc's `DMLocatePoints`. `DMLocatePoints` is +COLLECTIVE on the mesh DM communicator, so every rank has to reach the same +verdict — and the policy is a mesh capability, so they do. + +It was then thrown away in transit. `CachedDMInterpolationInfo.create_structure` +had no hint array to pass when a rank held zero points, and passed +`hintAuthoritative = 0` hardcoded along with the NULL. That rank alone took the +`DMLocatePoints` branch and blocked inside `DMGetBoundingBox -> MPI_Allreduce`, +while its peers bypassed and ran on to `DMSwarmMigrate -> MPI_Comm_dup`. + +Measured before the fix, at np=4: three ranks in `MPI_Comm_dup`, one in +`MPI_Allreduce`, no progress in 900 s on a 300-point query. + +The trigger is a query set that leaves some rank empty. Biasing every point into +x > 0.5 does it at np=4 — one rank owns only x < 0.5 — but not at np=2, where +the coarser partition leaves both ranks straddling the split. That is exactly +why this hid at np=2 and appeared at np=4. +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _mesh_and_field(tag): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)) + field = uw.discretisation.MeshVariable(tag, mesh, mesh.dim, degree=2) + return mesh, field + + +def test_global_evaluate_with_points_only_on_some_ranks(): + """The regression. Every point in x > 0.5, so some rank owns none of them.""" + _mesh, field = _mesh_and_field("u1076a") + rng = np.random.default_rng(42) + coords = rng.random((300, 2)) + coords[:, 0] = 0.5 + 0.5 * coords[:, 0] + + result = uw.function.global_evaluate(field.sym, coords) + assert result.shape[0] == 300, ( + f"rank {uw.mpi.rank}: expected 300 results, got {result.shape[0]}") + + +def test_global_evaluate_with_points_everywhere_still_works(): + """Negative control. + + The unbiased query already passed before the fix, so if this ever fails the + change has broken the ordinary path rather than repaired the empty-rank one. + """ + _mesh, field = _mesh_and_field("u1076b") + rng = np.random.default_rng(42) + coords = rng.random((300, 2)) + + result = uw.function.global_evaluate(field.sym, coords) + assert result.shape[0] == 300, ( + f"rank {uw.mpi.rank}: expected 300 results, got {result.shape[0]}")