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