A split fault level keeps its structural rows in a zone-keyed FAC patch (#921)

docs/developer/design/fault-parallel-placement-2026-09.md

@@ -834,7 +834,8 @@ overlap layer.
 The last row is the default now when a zone is keyed: `fac_zone` masks
 give the blocks, the structural split is off, the sub-solve is LU with a
 nonzero shift, the overlap one layer (`UW_FAC_MODE`, `UW_FAC_SUB_PC`,
-`UW_FAC_OVERLAP`, `UW_FAC_STRUCTURAL` remain as measurement knobs). One
+`UW_FAC_OVERLAP` remain as measurement knobs; `UW_FAC_STRUCTURAL`, listed here
+at the time, and `UW_FAC_COVER` were removed in 2026-11 by #921 and #918). One
 parallel trap on the way: the structural split's cover gate is a
 collective, and a rank with no band cells decided alone to take it while
 the others skipped it, a hang at np=4; the decision is now made from

docs/developer/subsystems/fault-split-multigrid.md

@@ -282,6 +282,20 @@ mesh-owned and auto-injected hierarchies do not print, and their record is at
 runs on, which leaves out pairs that touch a split coarse level. Parallel builds gate
 only the finest. Tests: `test_1134` in `tests/` and `tests/parallel/`.
 
+With `fac_zone` masks keyed, the masks' cells give the patch's blocks. On a finest
+level that carries a split fault, the structural rows are added as one more block,
+whatever the masks cover, under the split level's 98 % gate. On a zone-keyed level,
+the gate's verdict in `fac_gate` is the structural block's alone: the zone's blocks are
+installed either way. On a placed split fault the structural rows are the band, the
+split duplicates and the fill around them, which the coarse level cannot represent
+and whole-level SOR converges slowly. On the add-on matrix's row 5 (a TI band keyed
+beside a split fault at h_b/S = 24), the velocity solve took 820 iterations serially, summed
+over 17 applications, with the TI blocks alone and 94 with the structural block added
+(#921). Keying the fault's band mask as well is not enough (344): the mask
+holds the band cells but not the fill around them. A level without a split keys its
+patch on the masks alone. That is #670's default, measured on a placed weak plane
+when the patch replaced the level smoother (`test_1135`).
+
 ## Not yet done
 
 - A regression test for the coarse-level wall rotation itself (a rotated free-slip

scripts/test.sh

@@ -335,6 +335,8 @@ if [ $PARALLEL_ONLY -eq 0 ]; then
   batch 4 tests/parallel/test_1138_serial_reference_child_env.py
   # So is the split level's FAC patch, taken over the old cover gate (#918).
   batch 4 tests/parallel/test_1134_fac_cover_gate_split_parallel.py
+  # So is the zone-keyed patch of a split level, which keeps its structural rows (#921).
+  batch 4 tests/parallel/test_1135_fac_zone_split_structural_parallel.py
   # test_106*/test_107* were held back pending "a triage/deselect decision". That
   # decision is made: the band is not the unit of exclusion. Of the 243 tests in
   # the deferred bands only 26 are level_3 or tier_c; 55 are level_1 AND tier_a —

src/underworld3/utilities/custom_mg.py

@@ -686,9 +686,8 @@ def _fac_patch_split(P_csr, coords_c, coords_f, map_c, map_f, nc,
     than ``cover_max`` of the level: a uniform pair refines everywhere, and
     whole-level smoothing is then the right configuration. The hierarchy
     passes 0.75, or 0.98 on a 2-D level that carries a split fault
-    (``_FAC_COVER_MAX``, ``_FAC_COVER_MAX_SPLIT``); ``cover_max=None`` takes
-    the split whatever its size. A dict passed as ``gate`` receives the
-    verdict: ``rows``, ``subdomain_rows`` and ``taken``.
+    (``_FAC_COVER_MAX``, ``_FAC_COVER_MAX_SPLIT``). A dict passed as
+    ``gate`` receives the verdict: ``rows``, ``subdomain_rows`` and ``taken``.
     """
     P = P_csr.tocsr()
     n = P.shape[0]
@@ -788,7 +787,7 @@ def _fac_patch_split(P_csr, coords_c, coords_f, map_c, map_f, nc,
     node_in_sub[node_f[halo_rows]] = True
 
     sub = np.flatnonzero(node_in_sub[node_f])
-    taken = cover_max is None or sub.size <= cover_max * n
+    taken = sub.size <= cover_max * n
     if gate is not None:
         gate.update(rows=int(n), subdomain_rows=int(sub.size), taken=bool(taken))
     if not taken:
@@ -959,7 +958,7 @@ def _fac_patch_split_parallel(P, coords_c, coords_f, lay_c, lay_f, nc, comm,
     sub_local = np.flatnonzero(node_in_sub[node_of_row]) if n else np.zeros(0, dtype=np.int64)
     n_sub = comm.allreduce(int(sub_local.size), op=MPI.SUM)
     n_all = comm.allreduce(int(n), op=MPI.SUM)
-    taken = cover_max is None or n_sub <= cover_max * n_all
+    taken = n_sub <= cover_max * n_all
     if gate is not None:
         gate.update(rows=int(n_all), subdomain_rows=int(n_sub), taken=bool(taken))
     if not taken:
@@ -1479,9 +1478,11 @@ class CustomMGHierarchy:
         smoother (#629): a boolean mask over the finest mesh's cells — the
         painted / placed band the modeler authored — or a LIST of masks for
         one ASM block per fault segment (masks may overlap; a junction cell
-        in two masks sits inside both blocks). The structural patch (the
-        cut/split-inserted rows the transfer cannot represent) is unioned in
-        automatically. ``None`` keys the patch on the structural rows alone.
+        in two masks sits inside both blocks). When the finest level carries
+        a split fault, its structural patch (the rows the transfer cannot
+        represent) is added as one more block, unless the split level's cover
+        gate declines it; a level without a split keys the patch on the masks
+        alone (#921). ``None`` keys the patch on the structural rows alone.
     """
 
     def __init__(self, level_meshes, builder="barycentric", field_id=None,
@@ -1809,28 +1810,30 @@ def provenance_transfer():
                                 if piece.size:
                                     blocks.append((piece, piece))
                                     covered.update(piece.tolist())
-                    # TODO(MEASURE): #670 knob — UW_FAC_STRUCTURAL=0 keeps the
-                    # fault-zone blocks alone, without the structural
-                    # (non-nested transfer) rows; remove when settled.
+                    # A placed split fault's structural rows (its band, the
+                    # split duplicates and the fill around them) are rows the
+                    # coarse level cannot represent, and whole-level SOR
+                    # converges them slowly: with a TI zone keyed beside a
+                    # split fault at h_b/S = 24, leaving them out took 820
+                    # velocity iterations over 17 applications against 94
+                    # (#921). They join the zone's blocks whatever the zone
+                    # covers, under the split level's gate (#918). A level
+                    # without a split keeps the zone's blocks alone, the
+                    # default #670 measured on a placed weak plane. The split
+                    # marker is the same on every rank.
                     split = None
+                    split_level = fault_mg.is_split_level(self.level_meshes[l])
                     # the split's cover gate is collective: decide whether
                     # to take it from what ANY rank holds, not this rank
                     mpi_comm = comm.tompi4py()
                     zone_anywhere = mpi_comm.allreduce(len(blocks) > 0,
                                                        op=MPI.LOR)
-                    if (os.environ.get("UW_FAC_STRUCTURAL", "0") != "0"
-                            or not zone_anywhere):
-                        # rows added to a zone keyed on a level without a
-                        # split are taken whatever their size; otherwise the
-                        # level's gate (#918). The split marker is the same on
-                        # every rank.
-                        split_level = fault_mg.is_split_level(self.level_meshes[l])
+                    if split_level or not zone_anywhere:
                         gate = {}
                         split = _fac_patch_split_parallel(
                             P, coords[l - 1], coords[l], maps[l - 1], lay, nc,
                             mpi_comm,
-                            cover_max=(None if zone_anywhere and not split_level
-                                       else self._fac_cover_max(l)),
+                            cover_max=self._fac_cover_max(l),
                             gate=gate)
                         self.fac_gate[l] = gate or None
                     if split is not None:
@@ -1920,23 +1923,17 @@ def provenance_transfer():
                                 if piece.size:
                                     blocks.append((piece, piece))
                                     covered[piece] = True
-                        # The patch smoother REPLACES whole-level smoothing,
-                        # so it inherits every row the coarse level cannot
-                        # represent — the STRUCTURAL (non-identity) transfer
-                        # patch — whether or not the physics zone covers it.
-                        # A zone away from the cut/split leaves those rows
-                        # smoothed nowhere and the velocity solve stalls at
-                        # its cap (measured: any off-cut zone, #629). Add
-                        # the uncovered structural rows as one more block.
+                        # A split level's structural rows join the zone's
+                        # blocks as one more block under the split level's
+                        # gate, as in the parallel branch (#921, #918); a
+                        # level without a split keeps the zone's blocks alone
+                        # (#670).
                         split = None
-                        if os.environ.get("UW_FAC_STRUCTURAL", "0") != "0":
+                        if fault_mg.is_split_level(self.level_meshes[l]):
                             gate = {}
                             split = _fac_patch_split(
                                 Pr, coords[l - 1], coords[l], maps[l - 1],
-                                maps[l], nc,
-                                cover_max=(_FAC_COVER_MAX_SPLIT
-                                           if fault_mg.is_split_level(
-                                               self.level_meshes[l]) else None),
+                                maps[l], nc, cover_max=self._fac_cover_max(l),
                                 gate=gate)
                             self.fac_gate[l] = gate or None
                         if split is not None:
@@ -2138,9 +2135,11 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric",
     zone (#629): a boolean mask over the solver mesh's cells — e.g.
     ``mesh.cells_labelled("Band")`` for a placed ribbon, or
     ``mesh.cells_supporting("Fault")`` for a split surface's support — or a
-    list of masks for per-segment ASM blocks. The structural (cut/split)
-    rows are unioned in automatically; the patch is solved strongly once
-    per smoother application (basic ASM, shifted factorization).
+    list of masks for per-segment ASM blocks. On a split finest level the
+    structural rows are added as one more block, unless its gate declines
+    them; without a split the masks alone key the patch (#921). The patch is
+    solved strongly once per smoother application (basic ASM, shifted
+    factorization).
 
     ``verbose`` prints the installed transfer sizes and the FAC cover gate's
     verdict per level (``hierarchy.fac_gate``, #918)."""

tests/parallel/test_1135_fac_zone_split_structural_parallel.py

@@ -0,0 +1,122 @@
+"""A split finest level keeps its structural rows in the FAC patch when a zone is keyed (#921).
+
+The finest level's FAC patch is keyed either on the structural rows (the rows
+its transfer cannot represent, found from the transfer) or, when ``fac_zone``
+masks are passed, on the masks' cells. With a zone keyed, the structural rows
+joined the zone's blocks only under ``UW_FAC_STRUCTURAL=1``, against the
+documentation, which said they were always added. On a placed split fault those
+rows are the band, the split duplicates and the fill around them, which the
+coarse level cannot represent, and whole-level SOR converges them slowly: on the
+add-on matrix's row 5 (a TI band keyed beside a split fault, h_b/S = 24) the
+velocity solve took 820 iterations, summed over 17 applications, against 94 with
+the union.
+
+Now a level that carries a split fault always has its structural rows in the
+patch, whatever the zone. A level without a split (a painted or placed band)
+keeps #670's default, the zone's blocks alone; that default was measured on a
+placed weak plane with no split.
+
+Both tests key a zone strip near the left wall, away from the fault, on UWTN
+2026-022's linear fault fixture at 80 rungs (the fixture of test_1106), and
+compare with the structural patch the same level gets without a zone:
+
+- split: every structural row is in the zone-keyed patch, and so are the
+  strip's rows;
+- placed but not split: no structural row is in it, and the structural patch
+  is not empty, so a union would have shown.
+
+In parallel each rank compares its own rows: both hierarchies are built on the
+one solver, so the global numbering is the same. This file carries no
+``mpi(min_size=2)`` mark: at np 1 it checks the serial branch, in the parallel
+batches the parallel one.
+
+    python -m pytest tests/parallel/test_1135_fac_zone_split_structural_parallel.py
+    mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1135_fac_zone_split_structural_parallel.py
+"""
+import numpy as np
+import pytest
+
+import underworld3 as uw
+from underworld3.utilities.custom_mg import CustomMGHierarchy, set_custom_fmg
+from underworld3.utilities.place_surface import place_fault_ribbon_2d
+
+pytestmark = [pytest.mark.level_2]
+
+A, B = np.array([0.30, 0.45]), np.array([0.70, 0.55])
+N_RUNGS = 80
+# the zone: cells whose centroid lies in this strip, clear of the fault's
+# refinement and its halo
+STRIP = (0.03, 0.12)
+
+
+def _patch_rows(entry):
+    """Every row of one level's patch on this rank (owned and halo)."""
+    if entry is None:
+        return np.zeros(0, dtype=np.int64)
+    blocks = entry if isinstance(entry, list) else [entry]
+    return np.unique(np.concatenate([np.asarray(sub, dtype=np.int64)
+                                     for _own, sub in blocks]))
+
+
+def _structural_and_zone_patches(split):
+    """Solve the fixture with the structural patch (no zone), then build the
+    zone-keyed hierarchy on the same solver; return this rank's rows of both
+    finest-level patches."""
+    length = float(np.linalg.norm(B - A))
+    spacing = length / N_RUNGS
+    trace = A + (np.arange(N_RUNGS + 1) * spacing)[:, None] * (B - A) / length
+    base = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0, 0), maxCoords=(1, 1), cellSize=1 / 8, regular=True,
+        qdegree=2, refinement=2)
+    levels = base._coarse_level_meshes()[:-1] + [base]
+    mesh = place_fault_ribbon_2d(base, [("Fault", trace)], 2 * spacing, margin_rings=1,
+                                 clearance=0.7, split=split, mesher="ladder")[0]
+    x, y = mesh.X
+    v = uw.discretisation.MeshVariable("v", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1)
+    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
+    for wall in ("Bottom", "Top", "Left", "Right"):
+        stokes.add_dirichlet_bc((y - 0.5, 0.0), wall)
+    if split:
+        stokes.add_fault_bc(0, boundary="Fault")
+    stokes.petsc_use_pressure_nullspace = True
+    stokes.tolerance = 1e-6
+    set_custom_fmg(stokes, levels, field_id=0)
+    stokes.solve()
+    structural = _patch_rows(stokes._custom_mg["hierarchy"].level_patch_rows[-1])
+
+    cS, cE = mesh.dm.getHeightStratum(0)
+    centroids = np.array([mesh.dm.computeCellGeometryFVM(c)[1][:2]
+                          for c in range(cS, cE)]).reshape(-1, 2)
+    strip = (centroids[:, 0] > STRIP[0]) & (centroids[:, 0] < STRIP[1])
+    zoned = CustomMGHierarchy(levels + [mesh], field_id=0, fac_zone=strip)
+    zoned.build(stokes)
+    return structural, _patch_rows(zoned.level_patch_rows[-1])
+
+
+@pytest.mark.tier_b
+def test_a_split_level_keeps_its_structural_rows_with_a_zone_keyed():
+    structural, zoned = _structural_and_zone_patches(split=True)
+    assert uw.mpi.comm.allreduce(structural.size) > 0, "the split level has no structural patch"
+    missing = np.setdiff1d(structural, zoned).size
+    assert uw.mpi.comm.allreduce(missing) == 0, (
+        f"rank {uw.mpi.rank}: {missing} of {structural.size} structural rows are not in "
+        f"the zone-keyed patch (#921)")
+    # and the zone's own rows are there too: a dropped zone would leave the
+    # structural patch alone, which also contains every structural row
+    zone_only = np.setdiff1d(zoned, structural).size
+    assert uw.mpi.comm.allreduce(zone_only) > 0, "the zone strip is not in the patch"
+
+
+@pytest.mark.tier_b
+def test_an_unsplit_level_keys_its_patch_on_the_zone_alone():
+    structural, zoned = _structural_and_zone_patches(split=False)
+    assert uw.mpi.comm.allreduce(structural.size) > 0, (
+        "the placed level has no structural patch, so a union could not show")
+    assert uw.mpi.comm.allreduce(zoned.size) > 0, "the zone strip keyed no patch"
+    shared = np.intersect1d(structural, zoned).size
+    assert uw.mpi.comm.allreduce(shared) == 0, (
+        f"rank {uw.mpi.rank}: {shared} structural rows are in the zone-keyed patch of "
+        f"a level with no split (#670's default is the zone alone)")