Fault networks in parallel: co-located multigrid tail, gather-safe pair diagnostics

src/underworld3/function/_dminterp_wrapper.pyx

@@ -178,7 +178,18 @@ cdef class CachedDMInterpolationInfo:
                                            <size_t*> &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}")

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).

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"))

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()

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}")