Fault networks in parallel: keep the multigrid tail through the split, gauge the weak plane on owned probes

docs/advanced/fault-networks.md

@@ -63,7 +63,13 @@ jump between the two nodes of a cut pair, or the jump in tangential
 velocity across the layer, sampled one half-width plus a cell either
 side of the spine. Both are the fault's own throughput; a probe placed
 further out reads the surrounding flow as well and over-reads short
-strands.
+strands. The gauge is rank-local: each rank reports only the probe
+pairs it owns and omits a piece it holds no pair of, so a
+max-reduction across ranks recovers the network's answer. This
+matters because `evaluate` answers for any point it is handed,
+extrapolating from the nearest local cell when the point is not in
+the local mesh — a band-less rank would otherwise report a far-field
+extrapolation as the band's slip.
 
 `build(width=None)` keeps the older no-band path — graded refinement
 cut directly. It is split-only, and its mesh is not the one a weak
@@ -278,7 +284,16 @@ v1 scope, refused loudly outside it: planar patches (the
 near-miss — close but not crossing — is refused rather than guessed
 at). Multi-fault networks split and solve in parallel: `split_faults`
 redistributes ONCE, keyed on the union of the network's facets, and
-every split then runs with serial topology.
+every split then runs with serial topology. Either realisation runs
+its velocity block on the geometric multigrid tail the band's base
+mesh owns (`build` adopts it on the final mesh; `net.solve` says so
+when a solve falls back to algebraic multigrid). Placement is
+gather-first, so the cells gmsh fills into the carved cavity — the
+band and its graded surround — live on one rank; only the base's far
+field is balanced. On the crossing-patches fixture that is 8012 of
+8405 cells on one rank at np=4, and the solve is not faster than
+serial: parallel is a correctness mode for this path, not a speed-up,
+until the placed region is rebalanced.
 
 ## Limitations
 

src/underworld3/meshing/fault_network.py

@@ -558,7 +558,6 @@ def _build_3d_band(self, h_far=None, qdegree=2, realisation="split",
             coordinate_system_type=base.CoordinateSystem.coordinate_type,
             boundaries=Enum("boundaries", members), verbose=False)
         band = mesh.cells_labelled("Band", 71)
-        adopt_hierarchy(mesh, base, fac_zone=band)
         if realisation == "split":
             mesh = split_faults(mesh, [n for n, _P in self.prepared])
             # reduce first, then branch: the defect must raise on every
@@ -575,8 +574,15 @@ def _build_3d_band(self, h_far=None, qdegree=2, realisation="split",
                     "both sides of a pairing — the embedded "
                     "mid-surfaces are degenerate (a defect, not a "
                     "configuration error).")
-            mesh._custom_mg_fac_zone = None   # a split fault needs no patch
             band = mesh.cells_labelled("Band", 71)
+        # Adopt the base's multigrid tail on the FINAL mesh: a child of
+        # split_faults does not inherit (only Mesh.add_fault's wrapper
+        # does), and adopting first then splitting silently dropped the
+        # tail — every solver fell back to GAMG. The keying ruling
+        # stands: the split needs no FAC patch; the weak plane's patch
+        # is its band.
+        adopt_hierarchy(mesh, base,
+                        fac_zone=None if realisation == "split" else band)
         # honoured footprints: band cells within the USER patch's own
         # in-plane outline and half a width of its plane — planar patches
         # make the rule exact; the expanded margin stays unpainted. Near a
@@ -1245,9 +1251,21 @@ def damage_yield(self, velocity, dial=0.05, radius=None,
     # ------------------------------------------------------------------
     def solve(self, solver, **kwargs):
         """Solve with the fault contact imposed (thin convenience
-        wrapper over ``fault_contact.solve_with_fault``)."""
+        wrapper over ``fault_contact.solve_with_fault``). The velocity
+        preconditioner is reported when it falls back to algebraic
+        multigrid: a lost geometric tail must never decline quietly."""
+        import underworld3 as uw
         from underworld3.utilities.fault_contact import solve_with_fault
-        return solve_with_fault(solver, **kwargs)
+
+        info = solve_with_fault(solver, **kwargs)
+        if info.get("velocity_pc") == "GAMG":
+            uw.pprint(
+                "[FaultNetwork.solve] the velocity block ran on ALGEBRAIC "
+                "multigrid: this mesh owns no geometric tail, or it was "
+                "not engaged. Production runs want the FMG tail — "
+                "build(width=...) adopts it from its own base "
+                "automatically.")
+        return info
 
     # ------------------------------------------------------------------
     def slips(self, solver):
@@ -1285,11 +1303,34 @@ def slips(self, solver):
     def _slips_ti(self, solver):
         """The weak plane's slip: the tangential velocity jump across the
         band, one half-width plus a cell either side of each spine (2-D)
-        or patch (3-D, where the jump is projected onto the plane)."""
+        or patch (3-D, where the jump is projected onto the plane).
+
+        COLLECTIVE (every rank calls it), rank-local ANSWER: only probe
+        pairs whose BOTH points this rank owns count. ``evaluate``
+        answers for any point it is handed — extrapolating from the
+        nearest local cell when the point is not in the local mesh
+        (#641) — so on a distributed mesh a band-less rank would
+        otherwise report a far-field extrapolation as the band's slip
+        (measured 3x the true value at np=2). Every rank still evaluates
+        every probe — the call migrates an evaluation swarm, and a rank
+        that skipped it on a rank-local test would leave its peers
+        spinning (the conditional-collective hang) — and masks the answer
+        afterwards. A piece with no owned pair on this rank is absent
+        from the result, so a max-reduction across ranks recovers the
+        network's answer. Under co-located placement the band and its
+        skirt live on one rank, so no pair straddles a seam; a pair that
+        did would go unsampled on both sides."""
         import underworld3 as uw
 
         if self.info is None:
             raise RuntimeError("no band on this mesh: build(width=...)")
+        mesh = solver.u.mesh
+
+        def owned(points):
+            # rank-local, no communication: safe to differ across ranks
+            return mesh._robust_owning_cells(
+                np.ascontiguousarray(points)) >= 0
+
         out = {}
         for k, (name, P) in enumerate(self.prepared):
             P = np.asarray(P, dtype=float)
@@ -1302,24 +1343,34 @@ def _slips_ti(self, solver):
                 n_hat = _patch_normal(P)
                 S = np.vstack([P, 0.5 * (P + np.roll(P, -1, axis=0)),
                                P.mean(axis=0)])
+                plus, minus = S + skirt * n_hat, S - skirt * n_hat
+                # collective pair, unconditional, before any rank-local
+                # branch
                 vp = np.asarray(uw.function.evaluate(
-                    solver.u.sym, S + skirt * n_hat)).reshape(len(S), -1)
+                    solver.u.sym, plus)).reshape(len(S), -1)[:, :3]
                 vm = np.asarray(uw.function.evaluate(
-                    solver.u.sym, S - skirt * n_hat)).reshape(len(S), -1)
-                dv = (vp - vm)[:, :3]
+                    solver.u.sym, minus)).reshape(len(S), -1)[:, :3]
+                keep = owned(plus) & owned(minus)
+                if not keep.any():
+                    continue
+                dv = (vp - vm)[keep]
                 dv -= np.outer(dv @ n_hat, n_hat)
                 out[name] = float(np.linalg.norm(dv, axis=1).max())
                 continue
             P = P[:, :2]
             t = np.gradient(P, axis=0)
             t /= np.linalg.norm(t, axis=1)[:, None]
             n = np.column_stack([-t[:, 1], t[:, 0]])
+            plus, minus = P + skirt * n, P - skirt * n
             vp = np.asarray(uw.function.evaluate(
-                solver.u.sym, P + skirt * n)).reshape(len(P), -1)[:, :2]
+                solver.u.sym, plus)).reshape(len(P), -1)[:, :2]
             vm = np.asarray(uw.function.evaluate(
-                solver.u.sym, P - skirt * n)).reshape(len(P), -1)[:, :2]
-            out[name] = float(
-                np.abs(np.einsum("ij,ij->i", vp - vm, t)).max())
+                solver.u.sym, minus)).reshape(len(P), -1)[:, :2]
+            keep = owned(plus) & owned(minus)
+            if not keep.any():
+                continue
+            out[name] = float(np.abs(
+                np.einsum("ij,ij->i", (vp - vm)[keep], t[keep])).max())
         return out
 
     # ------------------------------------------------------------------

src/underworld3/utilities/fault_split.py

@@ -1357,6 +1357,20 @@ def split_faults(mesh, names, verbose=False):
             mesh._registered_children.add(out)
     for n in names:
         out = split_fault(out, n, verbose=verbose)
+    # The network's child INHERITS a mesh-owned geometric-MG tail — the
+    # same rule Mesh.add_fault applies: a cut re-represents the same
+    # grid, so the parent's coarse levels serve unchanged with the cut
+    # mesh as the finest level (#620/#629). Without this, every solver
+    # on a network split silently fell back to GAMG. The FAC zone is
+    # NOT inherited (a split fault needs no patch — the keying ruling);
+    # callers that key a zone re-adopt on the child explicitly.
+    own_tail = getattr(mesh, "_custom_mg_coarse_meshes", None)
+    if (own_tail is not None
+            and getattr(out, "_custom_mg_coarse_meshes", None) is None):
+        out._custom_mg_coarse_meshes = list(own_tail)
+        out._custom_mg_builder = getattr(mesh, "_custom_mg_builder",
+                                         "barycentric")
+        out._custom_mg_fac_zone = None
     return out
 
 

tests/parallel/ptest_0863_fault_network_3d_width_parallel.py

@@ -26,18 +26,27 @@
 # the serial run of this exact case (test_0863's end-to-end fixture):
 # peak tangential pair slip per prepared piece
 SERIAL = {"Main": 0.17567, "Cross_1": 0.00335, "Cross_2": 0.00295}
+# the weak plane's serial gauge (in-plane jump across the layer) and the
+# whole-domain integral of v.v on the same fixture, eta_1 = 0.01
+SERIAL_TI = {"Main": 0.29145, "Cross_1": 0.15717, "Cross_2": 0.13863}
+SERIAL_TI_VV = 8.35092135e-02
 
 
-def test_network_3d_width_split_solve_np2():
+def _build(realisation):
     fsA = uw.meshing.FaultSurface("Main", P_A)
     fsA.triangulate()
     fsB = uw.meshing.FaultSurface("Cross", P_B)
     fsB.triangulate()
     net = uw.meshing.FaultNetwork([fsA, fsB],
                                   hierarchy=["Main", "Cross"])
     net.prepare(h=H, ligament=1.0, verbose=False)
-    net.build(width=WIDTH, realisation="split", h_far=0.24,
+    net.build(width=WIDTH, realisation=realisation, h_far=0.24,
               margin_rings=0.5)
+    return net
+
+
+def test_network_3d_width_split_solve_np2():
+    net = _build("split")
     mesh = net.mesh
 
     # the mesh is DISTRIBUTED (the far field balanced; only the CAD
@@ -61,9 +70,46 @@ def test_network_3d_width_split_solve_np2():
     stokes.tolerance = 1e-5
     info = net.solve(stokes)
     assert info.get("converged")
+    # the geometric tail is adopted on the split child, on every rank
+    assert info.get("velocity_pc") == "custom-FMG", info
 
     slips = net.slips(stokes)                    # rank-local pairs
     for name, expected in SERIAL.items():
         peak = comm.allreduce(float(slips.get(name, 0.0)), op=max)
         assert peak == pytest.approx(expected, rel=2e-2), (
             f"{name}: parallel peak {peak:.4f} vs serial {expected}")
+
+
+def test_network_3d_width_weak_plane_solve_np2():
+    """The weak plane on the same band, distributed: the solve is the
+    serial one (a partition-independent integral says so), and the gauge
+    is reported only by the rank that owns the band's probes — a
+    band-less rank used to answer with an extrapolation 3x the slip."""
+    net = _build("ti")
+    mesh = net.mesh
+    comm = mesh.dm.comm.tompi4py()
+
+    x, y, z = mesh.X
+    v = uw.discretisation.MeshVariable("v3T", mesh, 3, degree=2)
+    p = uw.discretisation.MeshVariable("p3T", mesh, 1, degree=0,
+                                       continuous=False)
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.bodyforce = [0.0, 0.0, 0.0]
+    for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"):
+        stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall)
+    net.apply(stokes, eta_1=0.01)
+    stokes.petsc_use_pressure_nullspace = True
+    stokes.tolerance = 1e-5
+    stokes.solve()
+
+    vv = uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate()
+    assert vv == pytest.approx(SERIAL_TI_VV, rel=1e-5)
+
+    slips = net.slips(stokes)                    # rank-local, owned probes
+    band = int(np.count_nonzero(mesh.cells_labelled("Band", 71)))
+    if band == 0:
+        assert slips == {}, f"a band-less rank reported a gauge: {slips}"
+    for name, expected in SERIAL_TI.items():
+        peak = comm.allreduce(float(slips.get(name, 0.0)), op=max)
+        assert peak == pytest.approx(expected, rel=2e-2), (
+            f"{name}: parallel peak {peak:.4f} vs serial {expected}")

tests/test_0863_fault_network_3d_width.py

@@ -230,6 +230,9 @@ def test_the_realisations_solve_on_the_shared_band():
             stokes.tolerance = 1e-5
             info = net.solve(stokes)
             assert info.get("converged")
+            # the band's base owns a geometric tail; the split child must
+            # inherit it (it silently fell back to GAMG before)
+            assert info.get("velocity_pc") == "custom-FMG", info
         else:
             net.apply(stokes, eta_1=0.01)
             stokes.petsc_use_pressure_nullspace = True