Placement and split gather per region; the geometric tail solves on a distributed band (#670, #671)

docs/advanced/fault-networks.md

@@ -292,8 +292,12 @@ 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.
+serial. Since PR #672 the gather is per region: each zone, or
+junction-connected cluster of zones, is placed and split by the rank that
+holds it, and a zone whose shell is already interior to a rank is not
+moved at all; a single long fault is still one rank's. The design and its
+measurements are in the developer note
+`fault-parallel-placement-2026-09.md`.
 
 ## Limitations
 

docs/developer/design/fault_hierarchy_render.py

@@ -0,0 +1,101 @@
+"""The multigrid hierarchy of the layout fixture, seen: each level of the
+tail on the z = 0.5 slice (through all three faults), a zoom on one band
+with its fill, and at np=3 the partition on the same slice.
+
+    python -u render_hierarchy.py                 # levels + zoom (serial)
+    mpirun -np 3 python -u render_hierarchy.py    # partition (rank 0 renders)
+"""
+import numpy as np
+import underworld3 as uw
+import underworld3.visualisation as vis
+import pyvista as pv
+from mpi4py import MPI
+pv.OFF_SCREEN = True
+
+FIG = "/Users/lmoresi/+Simulations/fault_network_3d_parallel/figures"
+H, W = 0.08, 0.04
+def patch(x0, x1, y, z0=0.4, z1=0.6):
+    return np.array([[x0, y, z0], [x1, y, z0], [x1, y, z1], [x0, y, z1]])
+A = patch(0.40, 0.70, 0.50); B = patch(2.20, 2.60, 0.50)
+C = np.array([[2.20, 0.62, 0.42], [2.52, 0.30, 0.42], [2.52, 0.30, 0.58], [2.20, 0.62, 0.58]])
+D = patch(5.20, 5.50, 0.50)
+faults = []
+for name, P in (("A", A), ("B", B), ("C", C), ("D", D)):
+    f = uw.meshing.FaultSurface(name, P); f.triangulate(); faults.append(f)
+net = uw.meshing.FaultNetwork(faults, hierarchy=["A", "B", "C", "D"])
+net.prepare(h=H, ligament=1.5, verbose=False)
+net.realisation, net.width = "split", W
+net._build_3d_band(h_far=0.24, realisation="split", margin_rings=0.5, carve_clearance=0.3,
+                   minCoords=(0.0, 0.0, 0.0), maxCoords=(6.0, 1.0, 1.0))
+mesh = net.mesh
+comm = uw.mpi.comm
+levels = list(getattr(mesh, "_custom_mg_coarse_meshes", []) or []) + [mesh]
+
+def slice_edges(m, z=0.5):
+    grid = vis.mesh_to_pv_mesh(m)
+    return grid, grid.slice(normal="z", origin=(0, 0, z))
+
+if comm.size == 1:
+    for k, lv in enumerate(levels):
+        n_cells = int(lv.dm.getHeightStratum(0)[1])
+        n_v = int(np.diff(lv.dm.getDepthStratum(0))[0])
+        print(f"[hier] level {k}: {n_cells} cells, {n_v} vertices", flush=True)
+    # each level on the slice, same camera
+    for k, lv in enumerate(levels):
+        grid, sl = slice_edges(lv)
+        pl = pv.Plotter(off_screen=True, window_size=(2400, 500))
+        pl.set_background("white")
+        if k == len(levels) - 1:
+            band = np.asarray(mesh.cells_labelled("Band", 71)).astype(float)
+            grid.cell_data["band"] = band
+            sl = grid.slice(normal="z", origin=(0, 0, 0.5))
+            pl.add_mesh(sl, scalars="band", cmap="RdBu_r", clim=(0, 1), show_edges=True,
+                        edge_color="black", line_width=0.6, lighting=False, show_scalar_bar=False)
+        else:
+            pl.add_mesh(sl, color="white", show_edges=True, edge_color="black", line_width=0.6, lighting=False)
+        pl.view_xy(); pl.camera.parallel_projection = True
+        pl.camera.focal_point = (3.0, 0.5, 0.5); pl.camera.parallel_scale = 0.55
+        out = f"{FIG}/hierarchy_level{k}_z05.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True)
+    # zoom on fault B's band, fine level: band cells red, fill/base white
+    grid, _ = slice_edges(mesh)
+    grid.cell_data["band"] = np.asarray(mesh.cells_labelled("Band", 71)).astype(float)
+    for tag, origin, normal, focal, scale, size in (
+            ("z05", (0, 0, 0.5), "z", (2.4, 0.5, 0.5), 0.45, (1200, 1000)),
+            ("x24", (2.4, 0, 0), "x", (2.4, 0.5, 0.5), 0.45, (1000, 1000))):
+        sl = grid.slice(normal=normal, origin=origin)
+        pl = pv.Plotter(off_screen=True, window_size=size); pl.set_background("white")
+        pl.add_mesh(sl, scalars="band", cmap="RdBu_r", clim=(0, 1), show_edges=True,
+                    edge_color="black", line_width=0.8, lighting=False, show_scalar_bar=False)
+        # the coarse level's edges on top, in blue, to show the nesting
+        _g0, sl0 = slice_edges(levels[0], z=0.5) if normal == "z" else (None, vis.mesh_to_pv_mesh(levels[0]).slice(normal="x", origin=origin))
+        pl.add_mesh(sl0, color="white", opacity=0.0, show_edges=True, edge_color="blue", line_width=1.5, lighting=False)
+        if normal == "z": pl.view_xy()
+        else: pl.view_yz()
+        pl.camera.parallel_projection = True; pl.camera.focal_point = focal; pl.camera.parallel_scale = scale
+        out = f"{FIG}/hierarchy_zoomB_{tag}.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True)
+else:
+    # partition: gather every rank's cells (P1 vertices) to rank 0 and colour by rank
+    from underworld3.utilities.line_cut import _coords
+    dm = mesh.dm; vS, vE = dm.getDepthStratum(0); X = _coords(dm)[: vE - vS]
+    cn = np.asarray(mesh._cell_node_indices(1, True))
+    T = X[cn]                                          # (n_cells, 4, 3)
+    band = np.asarray(mesh.cells_labelled("Band", 71)).astype(float)
+    allT = comm.gather(T, root=0); allB = comm.gather(band, root=0)
+    ranks = comm.gather(np.full(len(T), comm.rank, dtype=float), root=0)
+    if comm.rank == 0:
+        pts = np.vstack([t.reshape(-1, 3) for t in allT])
+        n = pts.shape[0] // 4
+        cells = np.hstack([np.full((n, 1), 4), np.arange(4 * n).reshape(n, 4)]).ravel()
+        grid = pv.UnstructuredGrid(cells, np.full(n, pv.CellType.TETRA), pts)
+        grid.cell_data["rank"] = np.concatenate(ranks); grid.cell_data["band"] = np.concatenate(allB)
+        sl = grid.slice(normal="z", origin=(0, 0, 0.5))
+        pl = pv.Plotter(off_screen=True, window_size=(2400, 500)); pl.set_background("white")
+        pl.add_mesh(sl, scalars="rank", cmap="RdBu_r", clim=(0, comm.size - 1), show_edges=True,
+                    edge_color="black", line_width=0.4, lighting=False, show_scalar_bar=False)
+        bsl = grid.threshold(0.5, scalars="band").slice(normal="z", origin=(0, 0, 0.5))
+        if bsl.n_points:
+            pl.add_mesh(bsl, color="black", show_edges=True, edge_color="black", line_width=1.0, lighting=False)
+        pl.view_xy(); pl.camera.parallel_projection = True
+        pl.camera.focal_point = (3.0, 0.5, 0.5); pl.camera.parallel_scale = 0.55
+        out = f"{FIG}/partition_np{comm.size}_z05.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True)
+        print(f"[hier] np={comm.size} cells/rank {[len(t) for t in allT]} band/rank {[int(b.sum()) for b in allB]}", flush=True)

docs/developer/design/fault_parallel_layouts.py

@@ -0,0 +1,111 @@
+"""Throughput of the parallel fault network by LAYOUT, everything else
+fixed: the same four faults (one alone, a junction-connected pair, one
+more) on a 6 x 1 x 1 box whose np=3 partition is three slabs along x
+(seams near x = 2 and 4), split realisation, contact solve.
+
+    mpirun -np 3 python -u fault_parallel_layouts.py -uw_layout local|straddle|gathered [-uw_tail 0]
+
+local     : one fault per slab, the pair in the middle slab -> nothing moves
+straddle  : the fourth fault shifted onto the seam near x = 4 -> gathered
+gathered  : every fault forced into ONE region (the pre-#672 behaviour)
+-uw_tail 0: the geometric tail dropped -> the velocity block on GAMG
+
+Reports: build time, regions / cells moved, per-rank cells, Newton and
+Krylov counts, cold and warm solve wall time, and the slips (which must
+agree across layouts: the answer is partition-independent, the cost is
+not). Design note: docs/developer/design/fault-parallel-placement-2026-09.md
+"""
+import time
+import numpy as np
+import underworld3 as uw
+
+params = uw.Params(layout=uw.Param("local", "local | straddle | gathered"),
+                   tail=uw.Param(1, "1 = geometric tail (custom-FMG), 0 = GAMG"),
+                   solve=uw.Param(1, "0 = build only"),
+                   pres_rtol=uw.Param(0.0, "pressure sub-solve rtol (0 = the path's default, 0.1 x tol)"),
+                   tol=uw.Param(1e-5, "solver tolerance"),
+                   penalty=uw.Param(0.0, "grad-div penalty (0 = the default, none)"))
+layout = str(params.layout)
+if layout == "gathered":
+    # one region for the whole network: the old behaviour, for reference
+    from underworld3.utilities import place_surface as ps
+    _orig = ps._gather_regions
+    def one_region(dm, ids, verbose=False, layers=1):
+        ids = np.asarray(ids)
+        n_ids = dm.getComm().tompi4py().allreduce(int(ids.max()) if ids.size else 0, op=max)
+        work, n_region, n_moved, owner, _canon = _orig(
+            dm, (ids > 0).astype(np.int32), verbose=verbose, layers=layers)
+        return work, n_region, n_moved, owner, {k: 1 for k in range(1, n_ids + 1)}
+    ps._gather_regions = one_region
+
+from underworld3.utilities import place_surface as _ps
+_g = _ps._gather_regions
+_ps._gather_regions = lambda dm, ids, verbose=False, layers=1: _g(dm, ids, verbose=True, layers=layers)
+H, W = 0.08, 0.04
+def patch(x0, x1, y, z0=0.3, z1=0.7):
+    return np.array([[x0, y, z0], [x1, y, z0], [x1, y, z1], [x0, y, z1]])
+# a 6 x 1 x 1 box: np=3 slabs of length 2 with seams near x = 2 and 4;
+# faults at least 1.5 apart so no two shells touch, whichever layout
+# patches 0.2 tall (z 0.4..0.6): the 6:1 box meshes with cells up to 0.4
+# across, and the carve refuses a cavity that reaches a wall
+A = patch(0.40, 0.70, 0.50, 0.40, 0.60)           # slab 1
+B = patch(2.20, 2.60, 0.50, 0.40, 0.60)           # slab 2, senior of the pair
+C = np.array([[2.20, 0.62, 0.42], [2.52, 0.30, 0.42],
+              [2.52, 0.30, 0.58], [2.20, 0.62, 0.58]])   # crosses B
+xD = 5.20 if layout != "straddle" else 3.85       # slab 0, or ON the seam near 4
+D = patch(xD, xD + 0.30, 0.50, 0.40, 0.60)
+
+comm = uw.mpi.comm
+t0 = time.perf_counter()
+faults = []
+for name, P in (("A", A), ("B", B), ("C", C), ("D", D)):
+    f = uw.meshing.FaultSurface(name, P); f.triangulate(); faults.append(f)
+net = uw.meshing.FaultNetwork(faults, hierarchy=["A", "B", "C", "D"])
+net.prepare(h=H, ligament=1.5, verbose=False)   # 1.0 degenerates the junction cut on this mesh
+net.realisation, net.width = "split", W
+net._build_3d_band(h_far=0.24, realisation="split", margin_rings=0.5,
+                   carve_clearance=0.3, minCoords=(0.0, 0.0, 0.0),
+                   maxCoords=(6.0, 1.0, 1.0))
+t_build = time.perf_counter() - t0
+mesh = net.mesh
+if not int(params.solve):
+    cells = comm.gather(int(mesh.dm.getHeightStratum(0)[1]), root=0)
+    if comm.rank == 0:
+        print(f"[layout] {layout} build-only np={comm.size}: regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']} cells/rank {cells}", flush=True)
+    raise SystemExit(0)
+if not int(params.tail):
+    mesh._custom_mg_coarse_meshes = None
+cells = comm.gather(int(mesh.dm.getHeightStratum(0)[1]), root=0)
+band = comm.gather(int(np.count_nonzero(mesh.cells_labelled("Band", 71))), root=0)
+
+x, y, z = mesh.X
+v = uw.discretisation.MeshVariable("v", mesh, 3, degree=2)
+p = uw.discretisation.MeshVariable("p", mesh, 1, degree=0, continuous=False)
+stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+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)
+stokes.petsc_use_pressure_nullspace = True
+stokes.tolerance = float(params.tol)
+if float(params.penalty) > 0:
+    stokes.penalty = float(params.penalty)
+if float(params.pres_rtol) > 0:
+    stokes._rotated_pres_rtol = float(params.pres_rtol)
+t1 = time.perf_counter(); info = net.solve(stokes); t_cold = time.perf_counter() - t1
+# warm: a full solve again, from zero, with the tail and the rotation reused
+t2 = time.perf_counter(); info2 = net.solve(stokes, zero_init_guess=True); t_warm = time.perf_counter() - t2
+slips = net.slips(stokes)
+peaks = {n: comm.allreduce(float(slips.get(n, 0.0)), op=max) for n in sorted(set(k for k, _ in net.prepared))}
+if comm.rank == 0:
+    imb = max(cells) / (sum(cells) / len(cells))
+    print(f"[layout] {layout} tail={int(params.tail)} np={comm.size}: build {t_build:.1f}s; "
+          f"regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']}; "
+          f"cells/rank {cells} (max/mean {imb:.2f}) band/rank {band}", flush=True)
+    print(f"[layout] {layout} tail={int(params.tail)} tol={params.tol} penalty={params.penalty} pres_rtol={params.pres_rtol}: cold {t_cold:.1f}s warm {t_warm:.1f}s; "
+          f"pc={info.get('velocity_pc')} newton={info.get('nonlinear_iterations')} "
+          f"converged={info.get('converged')} vel_its={info.get('vel_its_last')} pres_its={info.get('pres_its_last')}; "
+          f"warm newton={info2.get('nonlinear_iterations')} vel_its={info2.get('vel_its_last')} reused={info2.get('rotation_reused')}", flush=True)
+    print(f"[layout] {layout} tail={int(params.tail)}: slips " + " ".join(f"{n}={s:.5f}" for n, s in peaks.items()), flush=True)

docs/developer/design/figures/seam-ligament/dump_fault_mesh.py

@@ -0,0 +1,55 @@
+"""Dump the split long-fault mesh (serial or ligament) for the diagram."""
+import os, numpy as np
+from mpi4py import MPI
+import underworld3 as uw
+from underworld3.utilities.place_surface import _shared_point_flags, _cells_anticlockwise, _cell_centroids_of
+comm = MPI.COMM_WORLD
+params = uw.Params(uw_seams="ligament")
+SEAMS = str(params.uw_seams)
+H = 0.02
+main = np.column_stack([np.full(36, 0.5), np.linspace(0.15, 0.85, 36)])
+base = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=4 * H, regular=False, refinement=1, qdegree=2)
+net = uw.meshing.FaultNetwork([("Main", main)])
+net.prepare(h=H, verbose=False)
+net.build(base=base, width=2 * H, realisation="split", max_levels=1, seams=SEAMS, band=2*H, ramp=6*H)
+mesh = net.mesh; dm = mesh.dm
+lig = net.ligament_cells()
+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 = 1.0
+if lig is not None:
+    net.apply(stokes, eta_1=0.01)
+    eta = np.asarray(net.ti["eta_1"].array).ravel()
+else:
+    net.apply(stokes)
+    eta = np.ones(int(dm.getHeightStratum(0)[1]))
+vS, vE = dm.getDepthStratum(0); pStart, _ = dm.getChart()
+X = np.asarray(mesh.X.coords)[:vE - vS]
+cells = _cells_anticlockwise(dm, X)
+shared = _shared_point_flags(dm).astype(bool)[vS - pStart: vE - pStart]
+cat = np.zeros(len(cells), dtype=int); cat[np.asarray(net.info["band"])] = 1
+if lig is not None: cat[lig] = 2
+edges = {}
+for side in ("Plus", "Minus"):
+    lbl = "Main" + side; out = []
+    if dm.hasLabel(lbl):
+        val = int(mesh.boundaries[lbl].value)
+        if dm.getLabel(lbl).getStratumSize(val):
+            for e in dm.getLabel(lbl).getStratumIS(val).getIndices():
+                a, b = (int(q) - vS for q in dm.getCone(int(e)))
+                out.append((a, b))
+    edges[side] = np.asarray(out, dtype=int).reshape(-1, 2)
+# tips: vertices of degree 1 in the Plus-edge graph
+from collections import Counter
+deg = Counter(edges["Plus"].ravel().tolist())
+tips = np.array([vtx for vtx, d in deg.items() if d == 1], dtype=int)
+# replica pairs (minus -> plus vertex ids) for drawing the duplicated nodes
+pairs = mesh._fault_point_pairs["Main"]
+rep = np.array([[qm - vS, qp - vS] for qm, qp in pairs.items() if vS <= qm < vE and vS <= qp < vE], dtype=int).reshape(-1, 2)
+np.savez(os.path.join(os.environ.get("SP", "."), f"faultmesh_{SEAMS}_np{comm.size}_rank{comm.rank}.npz"),
+         X=X, cells=cells, cat=cat, eta=eta, shared=shared, plus=edges["Plus"], minus=edges["Minus"], tips=tips, rep=rep)
+if comm.rank == 0:
+    print(f"[{SEAMS} np{comm.size}] dumped", flush=True)

docs/developer/design/figures/seam-ligament/plot_along_strike.py

@@ -0,0 +1,24 @@
+import os, sys, numpy as np
+import pyvista as pv
+pv.OFF_SCREEN = True
+SP = os.environ.get("SP", ".")
+pl = pv.Plotter(off_screen=True, window_size=(1600, 1600))
+colors = {0: "lightgrey", 1: "gold", 2: "tomato"}
+edge_colors = ["blue", "green"]
+for r in range(2):
+    d = np.load(os.path.join(SP, f"np2_rank{r}.npz"))
+    X = np.column_stack([d["X"], np.zeros(len(d["X"]))])
+    cells = d["cells"]
+    faces = np.column_stack([np.full(len(cells), 3), cells]).ravel()
+    grid = pv.PolyData(X, faces)
+    grid.cell_data["cat"] = d["cat"] + 3 * r
+    pl.add_mesh(grid, scalars="cat", show_edges=True, cmap=["#dddddd", "gold", "tomato", "#aaaaaa", "orange", "red"], clim=(0, 5), show_scalar_bar=False, line_width=0.5)
+    sh = d["shared"]
+    if sh.any():
+        pl.add_points(X[sh], color="black", point_size=6)
+    for a, b in d["edges"]:
+        pl.add_lines(np.array([X[a], X[b]]), color=edge_colors[r], width=4)
+pl.view_xy()
+pl.camera.zoom(1.4)
+pl.screenshot(os.path.join(SP, "np2_partition.png"))
+print("saved")

docs/developer/design/figures/seam-ligament/plot_cross.py

@@ -0,0 +1,23 @@
+import os, numpy as np
+import pyvista as pv
+pv.OFF_SCREEN = True
+SP = os.environ.get("SP", ".")
+pl = pv.Plotter(off_screen=True, shape=(1, 2), window_size=(2400, 1200))
+for k, mode in enumerate(("gather", "ligament")):
+    pl.subplot(0, k)
+    for r in range(2):
+        d = np.load(os.path.join(SP, f"cross_{mode}_rank{r}.npz"))
+        X = np.column_stack([d["X"], np.zeros(len(d["X"]))])
+        cells = d["cells"]
+        faces = np.column_stack([np.full(len(cells), 3), cells]).ravel()
+        grid = pv.PolyData(X, faces)
+        grid.cell_data["cat"] = d["cat"] + 3 * r
+        pl.add_mesh(grid, scalars="cat", show_edges=True, cmap=["#dddddd", "gold", "tomato", "#aaaaaa", "orange", "red"], clim=(0, 5), show_scalar_bar=False, line_width=0.5)
+        if d["shared"].any():
+            pl.add_points(X[d["shared"]], color="black", point_size=7)
+        for a, b in d["edges"]:
+            pl.add_lines(np.array([X[a], X[b]]), color=["blue", "green"][r], width=4)
+    pl.add_text(f"seams='{mode}'", font_size=14)
+    pl.view_xy(); pl.camera.zoom(1.5)
+pl.screenshot(os.path.join(SP, "cross_modes.png"))
+print("saved")