From 20ce01b2c22abbbb3985c76ac8198d99e2d4fcd5 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 08:39:06 -0700 Subject: [PATCH] Fault networks in parallel: the split keeps its multigrid tail, the weak-plane gauge answers only for owned probes Two defects on the finite-width network path, both silent. The split realisation ran every velocity solve on algebraic multigrid. build() adopted the base's geometric tail before splitting, and the split child did not inherit it (only Mesh.add_fault's wrapper does), so the adoption was lost and the rotated path fell back to GAMG without a word. build() now adopts the tail on the FINAL mesh, and split_faults hands a mesh-owned tail to its child the way add_fault does (the FAC zone is not inherited: a split needs no patch). net.solve() reports a fallback to algebraic multigrid so the substitution cannot decline quietly again. Serial, np=2 and np=4 all run the crossing-patches fixture on the geometric tail to the same slips (Main 0.17567 / Cross_1 0.00335 / Cross_2 0.00295). The weak-plane gauge reported 3x the true slip at np=2. Its probes sit in the band, which lives on one rank under gather-first placement; 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 the band-less rank returned a far-field extrapolation and the max across ranks took it. The gauge now masks each probe pair by ownership. Every rank still evaluates every probe: the call is collective, and skipping it on a rank-local test is the conditional-collective hang (measured: peers spin at 100%). A whole-domain integral of v.v agrees with serial to eight digits at np=2 and np=4, so the solve was never wrong; the gauge was. The previous parallel port of the FAC patch smoother is not in this change: with it on, the np=2 weak-plane velocity solve stalled at a flat residual. Whole-level smoothing on the adopted tail is the parallel configuration for now. ptest_0863 gains the weak-plane case (np=2 and np=4 checked) and pins the velocity preconditioner in both realisations; test_0863 pins it in serial. Underworld development team with AI support from Claude Code --- docs/advanced/fault-networks.md | 19 ++++- src/underworld3/meshing/fault_network.py | 75 ++++++++++++++++--- src/underworld3/utilities/fault_split.py | 14 ++++ ...st_0863_fault_network_3d_width_parallel.py | 50 ++++++++++++- tests/test_0863_fault_network_3d_width.py | 3 + 5 files changed, 145 insertions(+), 16 deletions(-) diff --git a/docs/advanced/fault-networks.md b/docs/advanced/fault-networks.md index 15609fa07..2424b00f7 100644 --- a/docs/advanced/fault-networks.md +++ b/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 diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py index 957440316..df9178066 100644 --- a/src/underworld3/meshing/fault_network.py +++ b/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,11 +1343,17 @@ 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 @@ -1314,12 +1361,16 @@ def _slips_ti(self, solver): 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 # ------------------------------------------------------------------ diff --git a/src/underworld3/utilities/fault_split.py b/src/underworld3/utilities/fault_split.py index ebb278d13..89d49646c 100644 --- a/src/underworld3/utilities/fault_split.py +++ b/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 diff --git a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py index 06749e3ce..01da8c259 100644 --- a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py +++ b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py @@ -26,9 +26,13 @@ # 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) @@ -36,8 +40,13 @@ def test_network_3d_width_split_solve_np2(): 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}") diff --git a/tests/test_0863_fault_network_3d_width.py b/tests/test_0863_fault_network_3d_width.py index 049084456..67da322c9 100644 --- a/tests/test_0863_fault_network_3d_width.py +++ b/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