From df540f917f2719370550b7c6329ea863431d1978 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 14:08:35 +1100 Subject: [PATCH 01/13] test_1139: the zone-keyed FAC blocks are the same at any rank count (#925) On row 5's TI band beside a split fault (n_rungs 80) the parallel build cut the zone and the structural block at rank boundaries and held seam rows in blocks on two ranks: the velocity solve took 4 iterations serially, 118 at np 2 and more than 200 at np 4. The test asserts every patch row in one block over all ranks, the serial block count and zone block (1650 rows), at most 5 velocity iterations (4 measured with the serial blocks at np 1, 2 and 4), and an unsplit zone cut into ceil(rows / 500) pieces at fac_block_rows 500. Passes serially; fails at np 2 and 4 before the fix. Underworld development team with AI support from Claude Code --- scripts/test.sh | 2 + ...fac_zone_blocks_any_rank_count_parallel.py | 153 ++++++++++++++++++ 2 files changed, 155 insertions(+) create mode 100644 tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py diff --git a/scripts/test.sh b/scripts/test.sh index dca2b5906..89c79a062 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -331,6 +331,8 @@ if [ $PARALLEL_ONLY -eq 0 ]; then 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 + # So are the zone-keyed patch's blocks, the same at any rank count (#925). + batch 4 tests/parallel/test_1139_fac_zone_blocks_any_rank_count_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 — diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py new file mode 100644 index 000000000..c88a55889 --- /dev/null +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -0,0 +1,153 @@ +"""The finest level's FAC blocks are the same at any rank count (#925). + +With ``fac_zone`` masks keyed, the finest level's FAC patch is cut into blocks: +the zone's rows segmented along strike at ``fac_block_rows``, and, on a level +that carries a split fault, the structural rows outside the zone as one more +block (#921). Basic ASM sums the correction of every block that holds a row. + +In parallel each rank built its blocks from its own cells: its part of the zone, +with the ghost rows of its masked cells, segmented on its own, and its own +structural rows as one block. A row on a partition seam then sat in blocks on two +ranks, and rank boundaries cut both the zone and the structural block. On the +add-on matrix's row 5 (a TI band keyed beside a split fault, n_rungs 80) the +velocity solve took 4 iterations serially, 118 at np 2 and more than 200 at +np 4; at n_rungs 320 the outer solve never finished its first iteration at +np 4. Removing the rows held twice was not enough: the solve converged in +4-6 iterations at np 2 and 4 only once the blocks were the serial ones, +gathered across ranks. + +The fixture is that problem at n_rungs 80: UWTN 2026-022's problem A box, the +main fault split, and a TI band (eta_1 = 0.01 along planes parallel to it) placed +0.08 below it as a ribbon, not cut, and keyed as the zone. + +- split: every patch row is held by exactly one block over all ranks; there are + two blocks, the zone's 1650 rows under the cap and the structural rows + outside it, as in serial; and the velocity solve takes at most 5 iterations + (4 measured at np 1, 2 and 4). +- placed but not split, the zone alone at ``fac_block_rows=500``: every row is + held once, and the zone is cut into the fewest pieces of at most 500 rows, + ``ceil(rows / 500)``, whatever the rank count. + +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_1139_fac_zone_blocks_any_rank_count_parallel.py + mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.utilities.custom_mg import 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 +LENGTH = float(np.linalg.norm(B - A)) +S = LENGTH / N_RUNGS +T_HAT = (B - A) / LENGTH +N_HAT = np.array([T_HAT[1], -T_HAT[0]]) + + +def _solve_fixture(split, fac_block_rows): + """Solve the TI-band-beside-a-fault box with the band keyed as the zone; + return the solver and every rank's finest-level blocks (global rows).""" + main = A + (np.arange(N_RUNGS + 1) * S)[:, None] * T_HAT + band = A + 0.08 * N_HAT + np.linspace(0.0, 1.0, N_RUNGS + 1)[:, None] * (B - A) + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0, 0), maxCoords=(1, 1), cellSize=1 / 8, regular=False, + qdegree=2, refinement=2) + levels = base._coarse_level_meshes()[:-1] + [base] + mesh = place_fault_ribbon_2d( + base, [("Fault", main)], 2 * S, spines=[("Fault", main), ("TI", band)], + 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, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + + zone = uw.discretisation.MeshVariable("zone", mesh, 1, degree=0, continuous=False) + c = np.asarray(zone.coords) + P, Q = band[0], band[-1] + s = np.clip(((c - P) @ (Q - P)) / ((Q - P) @ (Q - P)), 0.0, 1.0) + zone.array[:, 0, 0] = (np.linalg.norm(c - (P + s[:, None] * (Q - P)), axis=1) + < S).astype(float) + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.constitutive_model.Parameters.shear_viscosity_1 = 1.0 - 0.99 * zone.sym[0] + stokes.constitutive_model.Parameters.director = sympy.Matrix([[N_HAT[0], N_HAT[1]]]) + stokes.bodyforce = [0.0, 0.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 = 1.0e-6 + + band_cells = np.asarray(mesh.cells_labelled("Band", 72), dtype=bool) + set_custom_fmg(stokes, levels, field_id=0, fac_zone=band_cells, + fac_block_rows=fac_block_rows) + stokes.solve() + entry = stokes._custom_mg["hierarchy"].level_patch_rows[-1] + pairs = [] if entry is None else entry if isinstance(entry, list) else [entry] + mine = [np.asarray(sub, dtype=np.int64) for _own, sub in pairs] + blocks = [b for rank in uw.mpi.comm.allgather(mine) for b in rank if b.size] + return stokes, blocks + + +def _times_held(blocks): + """How many blocks hold each row, over all ranks.""" + _rows, count = np.unique(np.concatenate(blocks), return_counts=True) + return count + + +@pytest.fixture(scope="module") +def split_level(): + return _solve_fixture(split=True, fac_block_rows=2000) + + +@pytest.mark.tier_b +def test_each_patch_row_is_in_one_block(split_level): + _stokes, blocks = split_level + count = _times_held(blocks) + assert count.max() == 1, ( + f"{int(np.sum(count > 1))} of {count.size} patch rows are held by two or more " + f"blocks at np {uw.mpi.size}; basic ASM adds every block's correction there (#925)") + + +@pytest.mark.tier_b +def test_the_blocks_are_the_serial_ones(split_level): + _stokes, blocks = split_level + sizes = sorted(int(b.size) for b in blocks) + # serial: the zone's 1650 rows fit under the cap as one block, and the + # structural rows outside the zone are one more + assert len(sizes) == 2 and sizes[0] == 1650, ( + f"np {uw.mpi.size}: block sizes {sizes}, serial has the zone's 1650 rows " + f"and one structural block (#925)") + + +@pytest.mark.tier_b +def test_the_velocity_solve_takes_the_serial_iterations(split_level): + stokes, _blocks = split_level + info = stokes._rotated_freeslip_info + assert info["converged"] + its = max(info["vel_its_last"]) + # 4 measured at np 1, 2 and 4; rank-cut blocks took 118 at np 2 and + # more than 200 at np 4 + assert its <= 5, f"np {uw.mpi.size}: the velocity solve took {its} iterations (#925)" + + +@pytest.mark.tier_b +def test_an_unsplit_zone_is_cut_into_the_same_pieces_at_any_rank_count(): + _stokes, blocks = _solve_fixture(split=False, fac_block_rows=500) + count = _times_held(blocks) + assert count.max() == 1, ( + f"{int(np.sum(count > 1))} of {count.size} zone rows are held by two or more " + f"blocks at np {uw.mpi.size} (#925)") + expected = int(np.ceil(count.size / 500)) + assert len(blocks) == expected, ( + f"np {uw.mpi.size}: the zone's {count.size} rows are in {len(blocks)} blocks; " + f"segmented along strike at 500 rows they are {expected} (#925)") From 32fc0814f0540cf5d994066fd7ca501dfe52988d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 14:11:28 +1100 Subject: [PATCH 02/13] FAC zone blocks are built from every rank's rows (#925) In parallel each rank cut its own part of a fac_zone mask into blocks, with the ghost rows of its masked cells, and on a split level added its own structural rows (or, holding no zone cell, its whole structural subdomain) as one more block. Seam rows sat in blocks on two ranks, and rank boundaries cut the band and the structural block. On row 5's TI band beside a split fault the velocity solve took 4 iterations serially, 118 at np 2 and more than 200 at np 4 (n_rungs 80); at n_rungs 320 the outer solve never finished its first iteration at np 4. Holding each row once was not enough; the serial blocks rebuilt at np 2 and 4 gave 4-6 iterations per solve (campaigns probe 238338f). _gathered_zone_blocks gathers each mask's rows (nodes keyed by their global row, ordered by coordinate), segments each mask along strike as the serial branch does, adds the structural rows outside every mask as one block, and gives each block whole to the rank owning most of its rows. Serially, fed the parallel branch's inputs with every row twice, it returns the serial branch's blocks exactly (test_1139's fixture, caps 2000 and 500, split and unsplit). The structural patch without a zone keeps one block per rank. Underworld development team with AI support from Claude Code --- .../subsystems/fault-split-multigrid.md | 10 ++ src/underworld3/utilities/custom_mg.py | 115 +++++++++++++----- 2 files changed, 92 insertions(+), 33 deletions(-) diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index 4f48f6fef..d365f265d 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -237,6 +237,16 @@ holds the band cells but not the fill around them. A level without a split keys 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`). +In parallel the zone's blocks are built from every rank's rows, so they are the +serial blocks at any rank count. Each mask's rows are gathered and segmented along +strike as one band, the structural rows outside every mask are one block, and each +block goes whole to the rank that owns most of its rows. Built rank by rank, a +rank's blocks held seam rows that its neighbour's blocks also held, and rank +boundaries cut the band and the structural block. On row 5 at n_rungs 80 the velocity +solve then took 4 iterations serially, 118 at np 2 and more than 200 at np 4 (#925). +Gathered, it takes 4 at np 1, 2 and 4 (`test_1139`). The structural patch without a +zone is still one block per rank, made of that rank's owned rows and halo. + ## Not yet done - A regression test for the coarse-level wall rotation itself (a rotated free-slip diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index a99497c74..5ad60d3a7 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -773,6 +773,59 @@ def _segment_block(rows, node_of_row, coords, max_rows): return pieces +def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): + """The finest level's FAC blocks in parallel, the same at any rank count. + + ``zones`` holds, per ``fac_zone`` mask, this rank's rows of the mask's cells + (ghost rows resolved to their owners') as ``(rows, node_keys, coords)``, one + entry per row; ``structural`` is this rank's structural subdomain rows, or + None. COLLECTIVE: every rank gathers every rank's rows and cuts them as the + serial branch does. Each mask is one band, segmented along its strike by + :func:`_segment_block`; the structural rows outside every mask are one more + block. Each block then goes whole to the rank that owns most of its rows. + + Built rank by rank instead, a rank's blocks held the seam rows its + neighbour's blocks also held, and rank boundaries cut the band and the + structural block, which the coarse level cannot repair: on a TI band keyed + beside a split fault the velocity solve took 4 iterations serially, 118 at + np 2 and more than 200 at np 4 (#925). A row is in one block of each mask; + masks that overlap share their common rows by design. Returns this rank's + blocks as ``(rows, rows)`` pairs.""" + everyone = comm.allgather((rstart, zones, structural)) + starts = np.array([e[0] for e in everyone]) + blocks, covered = [], [] + for k in range(len(zones)): + rows = np.concatenate([e[1][k][0] for e in everyone]).astype(np.int64) + keys = np.concatenate([e[1][k][1] for e in everyone]) + xy = np.concatenate([e[1][k][2] for e in everyone]) + rows, first = np.unique(rows, return_index=True) + keys, xy = keys[first], xy[first] + # number the nodes in coordinate order, not by their partition-dependent + # global rows, so that the segmentation sees the same points in the + # same order at any rank count + node_keys, at = np.unique(keys, return_index=True) + order = np.lexsort((node_keys,) + tuple(xy[at].T[::-1])) + position = np.empty(order.size, dtype=np.int64) + position[order] = np.arange(order.size) + node = position[np.searchsorted(node_keys, keys)] + for piece in _segment_block(rows, node, xy[at][order], max_rows): + if piece.size: + blocks.append(np.unique(piece)) + covered.append(rows) + held = [e[2] for e in everyone if e[2] is not None] + if held: + sub = np.unique(np.concatenate(held).astype(np.int64)) + extra = np.setdiff1d(sub, np.concatenate(covered)) + if extra.size: + blocks.append(extra) + mine = [] + for block in blocks: + owner = np.searchsorted(starts, block, side="right") - 1 + if np.bincount(owner, minlength=comm.size).argmax() == comm.rank: + mine.append((block, block)) + return mine + + def _append_remote_coarse_nodes(remote, coords_c, l2g_c, lay_c, nc, comm): """Append the coarse nodes of the global rows ``remote``, held on other ranks, to this rank's node table. COLLECTIVE: one scatter from the rows' @@ -1739,37 +1792,31 @@ def provenance_transfer(): f"the co-partitioned transfer (#671).") # an orphan coarse DOF (no fine image) is dropped at install # (restrict_unreached), never injected — see drop_orphans - # FAC in parallel (#670): the same blocks as the serial branch - # below, in GLOBAL row numbering — the layout's ghost-resolved - # l2g puts a masked cell's off-rank nodes into this rank's - # subdomain, which is the seam halo; the structural split - # reads this rank's owned rows of the distributed transfer; - # the cover gate is one collective verdict. + # FAC in parallel (#670): the serial branch's blocks, in GLOBAL + # row numbering. A zone's rows are gathered from every rank (the + # layout's ghost-resolved l2g names a masked cell's off-rank + # nodes by their owner's rows) and cut into blocks as one band, + # so the blocks do not depend on the partition (#925); the + # structural split reads this rank's owned rows of the + # distributed transfer; the cover gate is one collective + # verdict. if (not os.environ.get("UW_CUSTOM_MG_DISABLE_FAC") and l == nlev - 1): lay = maps[l] blocks = [] - covered = set() + zones = [] if self.fac_zone is not None: cn = self.level_meshes[l]._cell_node_indices( degree, continuous) for mk in self.fac_zone: nodes = np.unique(cn[np.asarray(mk, dtype=bool)]) - full = (nodes[:, None] * nc - + np.arange(nc)[None, :]).ravel() - rows = lay.l2g[full] - keep = rows >= 0 - rows, node_r = rows[keep], full[keep] // nc - # TODO(MEASURE): #670 knob — rows per block; - # the default comes from the size sweep. - cap = int(os.environ.get("UW_FAC_BLOCK_ROWS", - self.fac_block_rows)) - for piece in _segment_block(rows, node_r, - coords[l], cap): - piece = np.unique(piece).astype(np.int64) - if piece.size: - blocks.append((piece, piece)) - covered.update(piece.tolist()) + rows = lay.l2g[nodes[:, None] * nc + + np.arange(nc)[None, :]] + # a node is named by its largest global row, which + # every rank holding it resolves alike + key = rows.max(axis=1) + n, c = np.nonzero(rows >= 0) + zones.append((rows[n, c], key[n], coords[l][nodes[n]])) # 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 @@ -1786,8 +1833,8 @@ def provenance_transfer(): # 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) + zone_anywhere = mpi_comm.allreduce( + any(z[0].size for z in zones), op=MPI.LOR) if split_level or not zone_anywhere: gate = {} split = _fac_patch_split_parallel( @@ -1796,16 +1843,18 @@ def provenance_transfer(): cover_max=self._fac_cover_max(l), gate=gate) self.fac_gate[l] = gate or None - if split is not None: + if zone_anywhere: + # TODO(MEASURE): #670 knob — rows per block; the + # default comes from the size sweep. + cap = int(os.environ.get("UW_FAC_BLOCK_ROWS", + self.fac_block_rows)) + blocks = _gathered_zone_blocks( + zones, None if split is None else split[1], cap, + lay.rstart, mpi_comm) + elif split is not None: _own, sub = split - extra = np.asarray([g for g in sub if g not in covered], - dtype=np.int64) - if blocks: - if extra.size: - blocks.append((extra, extra)) - else: - blocks.append((np.asarray(_own, dtype=np.int64), - np.asarray(sub, dtype=np.int64))) + blocks.append((np.asarray(_own, dtype=np.int64), + np.asarray(sub, dtype=np.int64))) # every rank installs the patch, or none does: a rank with # no rows of its own still sets the smoother type n_any = comm.tompi4py().allreduce(len(blocks), op=MPI.MAX) From e67aa114bef42a61a0836fa0e51364b9daeb84e1 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 14:30:53 +1100 Subject: [PATCH 03/13] #925 review: cut the zone blocks once on rank 0, nodes by coordinate Review findings, fixed: - Every rank recomputed the segmentation (an SVD) from the gathered rows and kept the blocks it owned; a rank whose LAPACK result differed in the last bit would cut different pieces, leaving rows in no block or in two. Rank 0 now gathers, cuts and assigns, and scatters each rank its blocks, which also keeps the gathered rows off the other ranks. - Nodes keyed by global row and ordered by coordinate broke ties between the two nodes of a split pair by row number, so a piece cut fell between them according to the partition (two masks, cap 300: 18 of 29 fake partitions differed by 4-8 rows). Nodes are now taken by coordinate: a split pair stays in one piece and the zone's pieces are the same at any rank count (30 of 30 fake partitions, np 1-7, caps 2000/500/300, split, unsplit and two masks). - test_1135 compared each rank's structural rows with that rank's blocks; a zone-keyed block now goes whole to one rank, so the rows are compared over all ranks. - test_1139 asserts the 500-row cap; the docs no longer claim the serial blocks exactly (the structural halo is _fac_patch_split_parallel's, #896); dated notes in the 2026-08 design record, whose 'seam cuts are free' and 'per-rank lists with local indices' #925 contradicts. Kept: with a zone keyed, the structural block is one LU on its majority owner, as in serial. Leaving it cut per rank (probe 238338f, local mode) took 2-4x the velocity iterations at n_rungs 80 and failed every solve at n_rungs 320, np 4. Its size is #895's question. Underworld development team with AI support from Claude Code --- .../design/fault-patch-multigrid-2026-08.md | 7 ++ .../subsystems/fault-split-multigrid.md | 11 +- src/underworld3/utilities/custom_mg.py | 119 ++++++++++-------- ...1135_fac_zone_split_structural_parallel.py | 42 +++---- ...fac_zone_blocks_any_rank_count_parallel.py | 14 ++- 5 files changed, 107 insertions(+), 86 deletions(-) diff --git a/docs/developer/design/fault-patch-multigrid-2026-08.md b/docs/developer/design/fault-patch-multigrid-2026-08.md index 967ae1c76..533a52079 100644 --- a/docs/developer/design/fault-patch-multigrid-2026-08.md +++ b/docs/developer/design/fault-patch-multigrid-2026-08.md @@ -271,10 +271,17 @@ to permit it.) update the star forest consistently (the known np>=3 line-cut issue; the independent-pass design is the template). Along-strike block cuts at seams are measured free (the k-ladder). + > [2026-11-03] Not for a zone-keyed patch: cut at rank boundaries (a TI + > band beside a split fault), its blocks took the velocity solve from 4 + > iterations serially to 118 at np 2 and past 200 at np 4. The zone's blocks + > are now built from every rank's rows (#925, `test_1139`). 2. **PCASM is rank-local.** Blocks become per-rank lists with local indices; faultless ranks carry `nsd = 0`, and every IS-building path must survive empty sets (the empty-stratum `getIndices()` SEGV family, #589). + > [2026-11-03] A zone-keyed level's blocks are per-rank lists of GLOBAL + > rows: each block goes whole to the rank owning most of its rows and may + > hold other ranks' rows (#925). 3. **Pair detection must go collective or be replaced by zone keying.** Rank-local coordinate matching misses pairs that straddle a boundary and can false-match ghosts; the zone mask is rank-local by nature and diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index d365f265d..b041f0660 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -237,10 +237,13 @@ holds the band cells but not the fill around them. A level without a split keys 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`). -In parallel the zone's blocks are built from every rank's rows, so they are the -serial blocks at any rank count. Each mask's rows are gathered and segmented along -strike as one band, the structural rows outside every mask are one block, and each -block goes whole to the rank that owns most of its rows. Built rank by rank, a +In parallel the zone's blocks are built from every rank's rows. Each mask's rows are +gathered and segmented along strike as one band, with the two nodes of a split pair +kept in one piece, so the zone's pieces are the same at any rank count. The +structural rows outside every mask are one block; its halo can differ from the +serial one by a few rows where a partition boundary runs along the band (#896). Each +block goes whole to the rank that owns most of its rows, so with a zone keyed the +structural block is one LU on one rank, as in serial. Built rank by rank, a rank's blocks held seam rows that its neighbour's blocks also held, and rank boundaries cut the band and the structural block. On row 5 at n_rungs 80 the velocity solve then took 4 iterations serially, 118 at np 2 and more than 200 at np 4 (#925). diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 5ad60d3a7..3072ad7e5 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -774,56 +774,68 @@ def _segment_block(rows, node_of_row, coords, max_rows): def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): - """The finest level's FAC blocks in parallel, the same at any rank count. - - ``zones`` holds, per ``fac_zone`` mask, this rank's rows of the mask's cells - (ghost rows resolved to their owners') as ``(rows, node_keys, coords)``, one - entry per row; ``structural`` is this rank's structural subdomain rows, or - None. COLLECTIVE: every rank gathers every rank's rows and cuts them as the - serial branch does. Each mask is one band, segmented along its strike by - :func:`_segment_block`; the structural rows outside every mask are one more - block. Each block then goes whole to the rank that owns most of its rows. - - Built rank by rank instead, a rank's blocks held the seam rows its - neighbour's blocks also held, and rank boundaries cut the band and the - structural block, which the coarse level cannot repair: on a TI band keyed - beside a split fault the velocity solve took 4 iterations serially, 118 at - np 2 and more than 200 at np 4 (#925). A row is in one block of each mask; - masks that overlap share their common rows by design. Returns this rank's - blocks as ``(rows, rows)`` pairs.""" - everyone = comm.allgather((rstart, zones, structural)) - starts = np.array([e[0] for e in everyone]) - blocks, covered = [], [] - for k in range(len(zones)): - rows = np.concatenate([e[1][k][0] for e in everyone]).astype(np.int64) - keys = np.concatenate([e[1][k][1] for e in everyone]) - xy = np.concatenate([e[1][k][2] for e in everyone]) - rows, first = np.unique(rows, return_index=True) - keys, xy = keys[first], xy[first] - # number the nodes in coordinate order, not by their partition-dependent - # global rows, so that the segmentation sees the same points in the - # same order at any rank count - node_keys, at = np.unique(keys, return_index=True) - order = np.lexsort((node_keys,) + tuple(xy[at].T[::-1])) - position = np.empty(order.size, dtype=np.int64) - position[order] = np.arange(order.size) - node = position[np.searchsorted(node_keys, keys)] - for piece in _segment_block(rows, node, xy[at][order], max_rows): - if piece.size: - blocks.append(np.unique(piece)) - covered.append(rows) - held = [e[2] for e in everyone if e[2] is not None] - if held: - sub = np.unique(np.concatenate(held).astype(np.int64)) - extra = np.setdiff1d(sub, np.concatenate(covered)) - if extra.size: - blocks.append(extra) - mine = [] - for block in blocks: - owner = np.searchsorted(starts, block, side="right") - 1 - if np.bincount(owner, minlength=comm.size).argmax() == comm.rank: - mine.append((block, block)) - return mine + """The finest level's FAC blocks in parallel, built from every rank's rows. + + COLLECTIVE. Rank 0 gathers every rank's rows and cuts them as the serial + branch does: each mask is one band, segmented along its strike by + :func:`_segment_block`, and the structural rows outside every mask are one + more block. Each block goes whole to the rank that owns most of its rows. + The zone's pieces are therefore the same at any rank count; the structural + block's halo is :func:`_fac_patch_split_parallel`'s, which can differ by a + few rows where a partition boundary runs along the band. + + Cut rank by rank, a rank's blocks held the seam rows its neighbour's blocks + also held, and rank boundaries cut the band and the structural block (#925). + + Parameters + ---------- + zones : list of (ndarray, ndarray) + Per ``fac_zone`` mask, this rank's rows of the mask's cells (global, + ghost rows resolved to their owners') and each row's node coordinates. + structural : ndarray or None + This rank's structural subdomain rows (global), or None where the split + level's gate declined them or the level has no split. + max_rows : int + The cap on a zone piece's rows (``fac_block_rows``). + rstart : int + The first global row this rank owns. + comm : mpi4py communicator + + Returns + ------- + list of (ndarray, ndarray) + This rank's blocks as ``(rows, rows)`` pairs. A row is in one block of + each mask; masks that overlap share their common rows by design. + """ + everyone = comm.gather((rstart, zones, structural), root=0) + per_rank = None + if comm.rank == 0: + starts = np.array([e[0] for e in everyone]) + blocks, covered = [], [] + for k in range(len(zones)): + rows = np.concatenate([e[1][k][0] for e in everyone]).astype(np.int64) + xy = np.concatenate([e[1][k][1] for e in everyone]) + rows, first = np.unique(rows, return_index=True) + # nodes by coordinate: the segmentation sees the same points in the + # same order at any rank count, and the two nodes of a split pair + # stay in one piece + site_xy, site = np.unique(xy[first], axis=0, return_inverse=True) + for piece in _segment_block(rows, site.ravel(), site_xy, max_rows): + if piece.size: + blocks.append(np.unique(piece)) + covered.append(rows) + held = [e[2] for e in everyone if e[2] is not None] + if held: + sub = np.unique(np.concatenate(held).astype(np.int64)) + extra = np.setdiff1d(sub, np.concatenate(covered)) + if extra.size: + blocks.append(extra) + per_rank = [[] for _ in everyone] + for block in blocks: + owner = np.searchsorted(starts, block, side="right") - 1 + per_rank[int(np.bincount(owner, minlength=len(everyone)).argmax())].append( + (block, block)) + return comm.scatter(per_rank, root=0) def _append_remote_coarse_nodes(remote, coords_c, l2g_c, lay_c, nc, comm): @@ -1796,7 +1808,7 @@ def provenance_transfer(): # row numbering. A zone's rows are gathered from every rank (the # layout's ghost-resolved l2g names a masked cell's off-rank # nodes by their owner's rows) and cut into blocks as one band, - # so the blocks do not depend on the partition (#925); the + # so its pieces do not depend on the partition (#925); the # structural split reads this rank's owned rows of the # distributed transfer; the cover gate is one collective # verdict. @@ -1812,11 +1824,8 @@ def provenance_transfer(): nodes = np.unique(cn[np.asarray(mk, dtype=bool)]) rows = lay.l2g[nodes[:, None] * nc + np.arange(nc)[None, :]] - # a node is named by its largest global row, which - # every rank holding it resolves alike - key = rows.max(axis=1) n, c = np.nonzero(rows >= 0) - zones.append((rows[n, c], key[n], coords[l][nodes[n]])) + zones.append((rows[n, c], coords[l][nodes[n]])) # 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 diff --git a/tests/parallel/test_1135_fac_zone_split_structural_parallel.py b/tests/parallel/test_1135_fac_zone_split_structural_parallel.py index 689c8aaed..5360c9706 100644 --- a/tests/parallel/test_1135_fac_zone_split_structural_parallel.py +++ b/tests/parallel/test_1135_fac_zone_split_structural_parallel.py @@ -25,8 +25,10 @@ - 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 +In parallel the rows are compared over all ranks: both hierarchies are built on +the one solver, so the global numbering is the same, and a zone-keyed level's +blocks go whole to one rank each, whoever owns their rows (#925). This file +carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial branch, in the parallel batches the parallel one. @@ -50,18 +52,17 @@ 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])) + """Every row of one level's patch (owned and halo), over all ranks.""" + blocks = [] if entry is None else entry if isinstance(entry, list) else [entry] + mine = [np.asarray(sub, dtype=np.int64) for _own, sub in blocks] + held = [rows for rank in uw.mpi.comm.allgather(mine) for rows in rank] + return np.unique(np.concatenate(held)) if held else np.zeros(0, dtype=np.int64) 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.""" + zone-keyed hierarchy on the same solver; return the rows of both + finest-level patches over all ranks.""" length = float(np.linalg.norm(B - A)) spacing = length / N_RUNGS trace = A + (np.arange(N_RUNGS + 1) * spacing)[:, None] * (B - A) / length @@ -99,24 +100,23 @@ def _structural_and_zone_patches(split): @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" + assert 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)") + assert missing == 0, ( + f"{missing} of {structural.size} structural rows are not in the zone-keyed " + f"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" + assert np.setdiff1d(zoned, structural).size > 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, ( + assert 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" + assert 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)") + assert shared == 0, ( + f"{shared} structural rows are in the zone-keyed patch of a level with no " + f"split (#670's default is the zone alone)") diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index c88a55889..f95e8f43d 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -65,7 +65,7 @@ def _solve_fixture(split, fac_block_rows): base, [("Fault", main)], 2 * S, spines=[("Fault", main), ("TI", band)], margin_rings=1, clearance=0.7, split=split, mesher="ladder")[0] _x, y = mesh.X - v =uw.discretisation.MeshVariable("v", mesh, 2, degree=2) + v = uw.discretisation.MeshVariable("v", mesh, 2, degree=2) p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1, continuous=True) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) @@ -122,8 +122,8 @@ def test_each_patch_row_is_in_one_block(split_level): def test_the_blocks_are_the_serial_ones(split_level): _stokes, blocks = split_level sizes = sorted(int(b.size) for b in blocks) - # serial: the zone's 1650 rows fit under the cap as one block, and the - # structural rows outside the zone are one more + # serial: the zone's 1650 rows (macOS and Linux alike) fit under the cap as + # one block, and the structural rows outside the zone are one more assert len(sizes) == 2 and sizes[0] == 1650, ( f"np {uw.mpi.size}: block sizes {sizes}, serial has the zone's 1650 rows " f"and one structural block (#925)") @@ -148,6 +148,8 @@ def test_an_unsplit_zone_is_cut_into_the_same_pieces_at_any_rank_count(): f"{int(np.sum(count > 1))} of {count.size} zone rows are held by two or more " f"blocks at np {uw.mpi.size} (#925)") expected = int(np.ceil(count.size / 500)) - assert len(blocks) == expected, ( - f"np {uw.mpi.size}: the zone's {count.size} rows are in {len(blocks)} blocks; " - f"segmented along strike at 500 rows they are {expected} (#925)") + sizes = sorted(int(b.size) for b in blocks) + assert len(blocks) == expected and sizes[-1] <= 500, ( + f"np {uw.mpi.size}: the zone's {count.size} rows are in blocks of {sizes}; " + f"segmented along strike at 500 rows they are {expected} blocks of at most " + f"500 (#925)") From 6e7f6f4d3233c6b0b91409b0f67d9d038f2f8172 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 16:07:19 +1100 Subject: [PATCH 04/13] test_1139, test_1133: the structural patch alone is one block at any rank count (#925) With no zone keyed, the finest level's structural patch was one block per rank. On test_1139's fixture with the zone unkeyed (row 5, n_rungs 80) the velocity solve took 3 iterations serially, 4 at np 2 and 6 at np 4 (70, 96 and 143 over the solve, kaiju report 12c). test_1139 asserts one block and at most 3 iterations; test_1133 asserts one structural block and, now that the zone's blocks are gathered, the zone strip's one block of exactly 1510 rows at every rank count. Serially all pass (8 passed); at np 2 and 4 the structural assertions fail until the structural patch is gathered. Underworld development team with AI support from Claude Code --- .../test_1133_fac_patch_summary_parallel.py | 32 +++++++------- ...fac_zone_blocks_any_rank_count_parallel.py | 43 ++++++++++++++++--- 2 files changed, 53 insertions(+), 22 deletions(-) diff --git a/tests/parallel/test_1133_fac_patch_summary_parallel.py b/tests/parallel/test_1133_fac_patch_summary_parallel.py index 9757404ef..1dbf22239 100644 --- a/tests/parallel/test_1133_fac_patch_summary_parallel.py +++ b/tests/parallel/test_1133_fac_patch_summary_parallel.py @@ -18,13 +18,13 @@ partition. At every rank count they assert that every rank holds the identical summary and -that the coarse levels carry no patch. Against the serial baselines: +that the coarse levels carry no patch. Against the serial baselines, at every rank +count (the finest level's blocks are cut from the rows gathered over the ranks, +each held by one rank, #925): -- the zone strip's patch rows are exactly 1510 at np 1. At np > 1 the summed - rows can be more, because a zone block's rows include the seam halo, which two - ranks can both hold, but the union of the ranks' global rows is exactly 1510; -- the fault's structural rows are test_1106's recorded serial counts at every - rank count, and they equal the rank-local record summed over the ranks. +- the zone strip's patch is one block of exactly 1510 rows; +- the fault's structural patch is one block, its rows test_1106's recorded serial + counts, and they equal the rank-local record summed over the ranks. This file carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial baseline, in the parallel batches the partition. @@ -92,17 +92,11 @@ def test_the_zone_patch_summary_is_the_same_on_every_rank(): assert all(entry is None for entry in summary[:-1]), summary finest = summary[-1] assert finest is not None, "the zone strip installed no patch" - assert any(sum(rank) > 0 for rank in finest["block_rows_per_rank"]) - if uw.mpi.size == 1: - assert finest["owned_rows"] == ZONE_SERIAL_ROWS, finest - assert finest["blocks"] == 1, finest - else: - assert finest["owned_rows"] >= ZONE_SERIAL_ROWS, finest - entry = h.level_patch_rows[-1] - blocks = entry if isinstance(entry, list) else [entry] - local = np.concatenate([np.asarray(own, dtype=np.int64) for own, _sub in blocks]) - union = np.unique(np.concatenate(uw.mpi.comm.allgather(local))) - assert union.size == ZONE_SERIAL_ROWS, (union.size, finest) + held = [rank for rank in finest["block_rows_per_rank"] if rank] + assert held == [[ZONE_SERIAL_ROWS]], ( + f"np {uw.mpi.size}: zone blocks per rank {finest['block_rows_per_rank']}; " + f"serial has one block of {ZONE_SERIAL_ROWS} rows (#925)") + assert finest["owned_rows"] == ZONE_SERIAL_ROWS, finest @pytest.mark.tier_b @@ -134,6 +128,10 @@ def test_the_structural_patch_summary_is_the_serial_rows_at_any_rank_count(): assert finest["owned_rows"] in SERIAL_PATCH_ROWS.values(), ( f"np {uw.mpi.size}: {finest['owned_rows']} structural patch rows; serial " f"counts {SERIAL_PATCH_ROWS} (a new host records its serial count in test_1106)") + held = [rank for rank in finest["block_rows_per_rank"] if rank] + assert len(held) == 1 and len(held[0]) == 1, ( + f"np {uw.mpi.size}: structural blocks per rank {finest['block_rows_per_rank']}; " + f"serial has one block (#925)") # the summary is the rank-local record summed over the ranks entry = h.level_patch_rows[-1] blocks = entry if isinstance(entry, list) else [entry] diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index f95e8f43d..7cd55bfa5 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -16,6 +16,11 @@ 4-6 iterations at np 2 and 4 only once the blocks were the serial ones, gathered across ranks. +The structural patch alone, with no zone keyed, was one block per rank: on the +same fixture its velocity solve took 3 iterations serially, 4 at np 2 and 6 at +np 4 (70, 96 and 143 summed over the solve). It converged, but a partition +boundary cut it all the same, and it too is built from the gathered rows. + The fixture is that problem at n_rungs 80: UWTN 2026-022's problem A box, the main fault split, and a TI band (eta_1 = 0.01 along planes parallel to it) placed 0.08 below it as a ribbon, not cut, and keyed as the zone. @@ -26,7 +31,9 @@ (4 measured at np 1, 2 and 4). - placed but not split, the zone alone at ``fac_block_rows=500``: every row is held once, and the zone is cut into the fewest pieces of at most 500 rows, - ``ceil(rows / 500)``, whatever the rank count. + ``ceil(rows / 500)``, whatever the rank count; +- split with no zone keyed: the structural patch is one block, as in serial, and + the velocity solve takes at most 3 iterations (3 measured at np 1). This file carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial branch, in the parallel batches the parallel one. @@ -52,9 +59,10 @@ N_HAT = np.array([T_HAT[1], -T_HAT[0]]) -def _solve_fixture(split, fac_block_rows): - """Solve the TI-band-beside-a-fault box with the band keyed as the zone; - return the solver and every rank's finest-level blocks (global rows).""" +def _solve_fixture(split, fac_block_rows, zone_keyed=True): + """Solve the TI-band-beside-a-fault box, with the band keyed as the zone + unless ``zone_keyed`` is False; return the solver and every rank's + finest-level blocks (global rows).""" main = A + (np.arange(N_RUNGS + 1) * S)[:, None] * T_HAT band = A + 0.08 * N_HAT + np.linspace(0.0, 1.0, N_RUNGS + 1)[:, None] * (B - A) base = uw.meshing.UnstructuredSimplexBox( @@ -87,7 +95,8 @@ def _solve_fixture(split, fac_block_rows): stokes.petsc_use_pressure_nullspace = True stokes.tolerance = 1.0e-6 - band_cells = np.asarray(mesh.cells_labelled("Band", 72), dtype=bool) + band_cells = (np.asarray(mesh.cells_labelled("Band", 72), dtype=bool) + if zone_keyed else None) set_custom_fmg(stokes, levels, field_id=0, fac_zone=band_cells, fac_block_rows=fac_block_rows) stokes.solve() @@ -153,3 +162,27 @@ def test_an_unsplit_zone_is_cut_into_the_same_pieces_at_any_rank_count(): f"np {uw.mpi.size}: the zone's {count.size} rows are in blocks of {sizes}; " f"segmented along strike at 500 rows they are {expected} blocks of at most " f"500 (#925)") + + +@pytest.fixture(scope="module") +def structural_alone(): + return _solve_fixture(split=True, fac_block_rows=2000, zone_keyed=False) + + +@pytest.mark.tier_b +def test_the_structural_patch_alone_is_one_block(structural_alone): + _stokes, blocks = structural_alone + sizes = sorted(int(b.size) for b in blocks) + assert len(sizes) == 1, ( + f"np {uw.mpi.size}: the structural patch is in blocks of {sizes}; serial " + f"has one (#925)") + + +@pytest.mark.tier_b +def test_the_structural_patch_alone_takes_the_serial_iterations(structural_alone): + stokes, _blocks = structural_alone + info = stokes._rotated_freeslip_info + assert info["converged"] + its = max(info["vel_its_last"]) + # 3 measured at np 1; one block per rank took 4 at np 2 and 6 at np 4 + assert its <= 3, f"np {uw.mpi.size}: the velocity solve took {its} iterations (#925)" From fc58912c99a272c85f2bd5a69dd88d50fe581576 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 16:10:58 +1100 Subject: [PATCH 05/13] The structural FAC patch alone is gathered into one block too (#925) With no zone keyed, the finest level's structural patch was one block per rank: each rank's owned patch rows and halo. It converged, but a partition boundary cut it: on row 5 with the zone unkeyed (n_rungs 80) the velocity solve took 3 iterations serially, 4 at np 2 and 6 at np 4 (kaiju report 12c). _gathered_zone_blocks becomes _gathered_patch_blocks and takes the split's (owned, subdomain) pair: with zone rows the structural rows outside them are one block as before; without, the gathered patch and halo are one block, held by the rank that owns most of its rows, as in serial. The structural LU then sits on one rank; patch_summary's block_rows_per_rank shows it and balancing it is #895's. A serial emulation over 30 fake partitions (np 1-7, two cut directions, rows renumbered) gives the serial block exactly, owned rows included. Underworld development team with AI support from Claude Code --- .../subsystems/fault-split-multigrid.md | 37 +++++---- src/underworld3/utilities/custom_mg.py | 79 +++++++++++-------- 2 files changed, 65 insertions(+), 51 deletions(-) diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index d09e4573a..15ab20996 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -72,11 +72,11 @@ largest block, and each rank's block sizes. It is made once, with one allgather `CustomMGHierarchy.install` and the rotated path's two `_configure_pcmg` calls. Reading it makes no collective call, so any rank can read it alone. A driver that looped a reduction over its own rank's blocks hung at np 4, because the ranks made -different numbers of calls. The structural patch's owned rows are disjoint across -ranks, so their total does not depend on the partition. A zone block's rows include -the seam halo, which two ranks can both hold, so its summed rows can exceed the serial -count while the union of the ranks' rows equals it (`test_1133`). The record is each -block as handed to the installer, before PCASM adds its overlap layer at set-up. +different numbers of calls. Each block of the finest level is held whole by one rank +(see below, #925), so the summed rows are the serial count at any rank count +(`test_1133`), and `block_rows_per_rank` shows where each block sits: the structural +LU on one rank. The record is each block as handed to the installer, before PCASM +adds its overlap layer at set-up. ## Rotated walls on every level @@ -281,18 +281,21 @@ holds the band cells but not the fill around them. A level without a split keys 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`). -In parallel the zone's blocks are built from every rank's rows. Each mask's rows are -gathered and segmented along strike as one band, with the two nodes of a split pair -kept in one piece, so the zone's pieces are the same at any rank count. The -structural rows outside every mask are one block; its halo can differ from the -serial one by a few rows where a partition boundary runs along the band (#896). Each -block goes whole to the rank that owns most of its rows, so with a zone keyed the -structural block is one LU on one rank, as in serial. Built rank by rank, a -rank's blocks held seam rows that its neighbour's blocks also held, and rank -boundaries cut the band and the structural block. On row 5 at n_rungs 80 the velocity -solve then took 4 iterations serially, 118 at np 2 and more than 200 at np 4 (#925). -Gathered, it takes 4 at np 1, 2 and 4 (`test_1139`). The structural patch without a -zone is still one block per rank, made of that rank's owned rows and halo. +In parallel the finest level's blocks are built from every rank's rows. Each mask's +rows are gathered and segmented along strike as one band, with the two nodes of a +split pair kept in one piece, so the zone's pieces are the same at any rank count. +The structural rows outside every mask are one block, and with no zone keyed the +structural patch is one block; its halo can differ from the serial one by a few rows +where a partition boundary runs along the band (#896). Each block goes whole to the +rank that owns most of its rows, so the structural block is one LU on one rank, as in +serial. That rank carries the factorisation the others do not +(`patch_summary`'s `block_rows_per_rank`); balancing it is #895's question. Built +rank by rank, a rank's blocks held seam rows that its neighbour's blocks also held, +and rank boundaries cut the band and the structural block. On row 5 at n_rungs 80 +the velocity solve then took 4 iterations serially, 118 at np 2 and more than 200 at +np 4 (#925). Gathered, it takes 4 at np 1, 2 and 4 (`test_1139`). The structural +patch alone, one block per rank, converged but took 3 iterations serially, 4 at np 2 +and 6 at np 4; gathered it is one block at any rank count (`test_1139`, `test_1133`). ## Not yet done diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 6643b9d0c..d5a693555 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -622,12 +622,12 @@ def _fac_patch_summary(patch_rows, comm): list One entry per level: None where no rank has a patch, else a dict with ``blocks`` (ASM subdomains over all ranks, empty ones included), - ``owned_rows`` (the patch rows; the structural patch's are disjoint across - ranks, so its total does not depend on the partition, while a zone - block's rows include the seam halo, which two ranks can both hold), + ``owned_rows`` (the patch rows; each block is held by one rank, so + the total does not depend on the partition, #925), ``subdomain_rows`` (the patch plus its halo, summed over the ranks), ``block_rows_max``, and ``block_rows_per_rank`` (each rank's block - sizes, in rank order). + sizes, in rank order: a block goes whole to the rank that owns most + of its rows, so the structural LU sits on one rank). """ mine = [] for entry in (patch_rows or []): @@ -831,16 +831,18 @@ def _segment_block(rows, node_of_row, coords, max_rows): return pieces -def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): +def _gathered_patch_blocks(zones, structural, max_rows, rstart, comm): """The finest level's FAC blocks in parallel, built from every rank's rows. COLLECTIVE. Rank 0 gathers every rank's rows and cuts them as the serial branch does: each mask is one band, segmented along its strike by :func:`_segment_block`, and the structural rows outside every mask are one - more block. Each block goes whole to the rank that owns most of its rows. - The zone's pieces are therefore the same at any rank count; the structural - block's halo is :func:`_fac_patch_split_parallel`'s, which can differ by a - few rows where a partition boundary runs along the band. + more block; with no zone rows, the structural patch is one block. Each + block goes whole to the rank that owns most of its rows, so the structural + block is one LU on one rank, as in serial. The blocks are therefore the + same at any rank count, except that the structural halo is + :func:`_fac_patch_split_parallel`'s, which can differ by a few rows where a + partition boundary runs along the band. Cut rank by rank, a rank's blocks held the seam rows its neighbour's blocks also held, and rank boundaries cut the band and the structural block (#925). @@ -850,9 +852,11 @@ def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): zones : list of (ndarray, ndarray) Per ``fac_zone`` mask, this rank's rows of the mask's cells (global, ghost rows resolved to their owners') and each row's node coordinates. - structural : ndarray or None - This rank's structural subdomain rows (global), or None where the split - level's gate declined them or the level has no split. + Empty where no zone is keyed. + structural : (ndarray, ndarray) or None + This rank's structural patch and subdomain rows (global, owned), or + None where the split level's gate declined them or the level has no + split. max_rows : int The cap on a zone piece's rows (``fac_block_rows``). rstart : int @@ -862,7 +866,9 @@ def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): Returns ------- list of (ndarray, ndarray) - This rank's blocks as ``(rows, rows)`` pairs. A row is in one block of + This rank's blocks as ``(owned, subdomain)`` pairs: ``(rows, rows)`` + for a zone piece and for the structural rows beside a zone, the patch + and its halo for the structural patch alone. A row is in one block of each mask; masks that overlap share their common rows by design. """ everyone = comm.gather((rstart, zones, structural), root=0) @@ -880,19 +886,26 @@ def _gathered_zone_blocks(zones, structural, max_rows, rstart, comm): site_xy, site = np.unique(xy[first], axis=0, return_inverse=True) for piece in _segment_block(rows, site.ravel(), site_xy, max_rows): if piece.size: - blocks.append(np.unique(piece)) + piece = np.unique(piece) + blocks.append((piece, piece)) covered.append(rows) held = [e[2] for e in everyone if e[2] is not None] if held: - sub = np.unique(np.concatenate(held).astype(np.int64)) - extra = np.setdiff1d(sub, np.concatenate(covered)) - if extra.size: - blocks.append(extra) + sub = np.unique(np.concatenate([s for _o, s in held]).astype(np.int64)) + zone_rows = (np.concatenate(covered) if covered + else np.zeros(0, dtype=np.int64)) + if zone_rows.size: + extra = np.setdiff1d(sub, zone_rows) + if extra.size: + blocks.append((extra, extra)) + elif sub.size: + owned = np.unique(np.concatenate([o for o, _s in held]).astype(np.int64)) + blocks.append((owned, sub)) per_rank = [[] for _ in everyone] - for block in blocks: - owner = np.searchsorted(starts, block, side="right") - 1 + for owned, sub in blocks: + owner = np.searchsorted(starts, sub, side="right") - 1 per_rank[int(np.bincount(owner, minlength=len(everyone)).argmax())].append( - (block, block)) + (owned, sub)) return comm.scatter(per_rank, root=0) @@ -1842,13 +1855,13 @@ def provenance_transfer(): # an orphan coarse DOF (no fine image) is dropped at install # (restrict_unreached), never injected — see drop_orphans # FAC in parallel (#670): the serial branch's blocks, in GLOBAL - # row numbering. A zone's rows are gathered from every rank (the + # row numbering. The structural split reads this rank's owned + # rows of the distributed transfer, under one collective cover + # verdict. A zone's rows are gathered from every rank (the # layout's ghost-resolved l2g names a masked cell's off-rank - # nodes by their owner's rows) and cut into blocks as one band, - # so its pieces do not depend on the partition (#925); the - # structural split reads this rank's owned rows of the - # distributed transfer; the cover gate is one collective - # verdict. + # nodes by their owner's rows), and the zone and the structural + # rows are cut into blocks from what every rank holds, so the + # blocks do not depend on the partition (#925). if (not os.environ.get("UW_CUSTOM_MG_DISABLE_FAC") and l == nlev - 1): lay = maps[l] @@ -1889,18 +1902,16 @@ def provenance_transfer(): cover_max=self._fac_cover_max(l), gate=gate) self.fac_gate[l] = gate or None - if zone_anywhere: + # split is None on every rank or on none: the gate's + # verdict is collective + if zone_anywhere or split is not None: # TODO(MEASURE): #670 knob — rows per block; the # default comes from the size sweep. cap = int(os.environ.get("UW_FAC_BLOCK_ROWS", self.fac_block_rows)) - blocks = _gathered_zone_blocks( - zones, None if split is None else split[1], cap, + blocks = _gathered_patch_blocks( + zones if zone_anywhere else [], split, cap, lay.rstart, mpi_comm) - elif split is not None: - _own, sub = split - blocks.append((np.asarray(_own, dtype=np.int64), - np.asarray(sub, dtype=np.int64))) # every rank installs the patch, or none does: a rank with # no rows of its own still sets the smoother type n_any = comm.tompi4py().allreduce(len(blocks), op=MPI.MAX) From bed177bf5594b78702562343c700adb50f8a0a21 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 16:18:05 +1100 Subject: [PATCH 06/13] test_1133: count non-empty blocks; a rank holding none records an empty one At np 2 and 4 the zone strip is one block of 1510 rows on one rank, and the other ranks record a block of 0 rows, so comparing the per-rank lists with the serial one failed on a correct patch (kaiju report 14a). The fixture's fault band lies on one rank at np 2 and 4, so its one-block assertion guards the summary; test_1139 is where the partition cut the structural patch. Underworld development team with AI support from Claude Code --- .../test_1133_fac_patch_summary_parallel.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/tests/parallel/test_1133_fac_patch_summary_parallel.py b/tests/parallel/test_1133_fac_patch_summary_parallel.py index 1dbf22239..803632680 100644 --- a/tests/parallel/test_1133_fac_patch_summary_parallel.py +++ b/tests/parallel/test_1133_fac_patch_summary_parallel.py @@ -24,7 +24,10 @@ - the zone strip's patch is one block of exactly 1510 rows; - the fault's structural patch is one block, its rows test_1106's recorded serial - counts, and they equal the rank-local record summed over the ranks. + counts, and they equal the rank-local record summed over the ranks. (This + fixture's band lies on one rank at np 2 and 4, so the one-block assertion + guards the summary; test_1139 is the fixture whose structural patch the + partition cut.) This file carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial baseline, in the parallel batches the partition. @@ -64,6 +67,12 @@ def _stokes(mesh): return stokes +def _nonempty_blocks(level): + """The sizes of a level's non-empty blocks over all ranks; a rank holding + no block records an empty one.""" + return [rows for rank in level["block_rows_per_rank"] for rows in rank if rows] + + def _same_on_every_rank(summary): text = json.dumps(summary, sort_keys=True) return len(set(uw.mpi.comm.allgather(text))) == 1 @@ -92,8 +101,7 @@ def test_the_zone_patch_summary_is_the_same_on_every_rank(): assert all(entry is None for entry in summary[:-1]), summary finest = summary[-1] assert finest is not None, "the zone strip installed no patch" - held = [rank for rank in finest["block_rows_per_rank"] if rank] - assert held == [[ZONE_SERIAL_ROWS]], ( + assert _nonempty_blocks(finest) == [ZONE_SERIAL_ROWS], ( f"np {uw.mpi.size}: zone blocks per rank {finest['block_rows_per_rank']}; " f"serial has one block of {ZONE_SERIAL_ROWS} rows (#925)") assert finest["owned_rows"] == ZONE_SERIAL_ROWS, finest @@ -128,8 +136,7 @@ def test_the_structural_patch_summary_is_the_serial_rows_at_any_rank_count(): assert finest["owned_rows"] in SERIAL_PATCH_ROWS.values(), ( f"np {uw.mpi.size}: {finest['owned_rows']} structural patch rows; serial " f"counts {SERIAL_PATCH_ROWS} (a new host records its serial count in test_1106)") - held = [rank for rank in finest["block_rows_per_rank"] if rank] - assert len(held) == 1 and len(held[0]) == 1, ( + assert len(_nonempty_blocks(finest)) == 1, ( f"np {uw.mpi.size}: structural blocks per rank {finest['block_rows_per_rank']}; " f"serial has one block (#925)") # the summary is the rank-local record summed over the ranks From 8a11f45507027519071e5725232e165e4d630465 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 16:30:43 +1100 Subject: [PATCH 07/13] test_1139: a keyed zone with no cells installs no patch, as in serial The serial branch keys the patch on the masks whenever a zone is keyed; the parallel branch fell back to the structural split when no rank held a zone cell. On test_1139's placed but unsplit fixture with an empty mask, serial installs no patch (passes at np 1). Underworld development team with AI support from Claude Code --- ...fac_zone_blocks_any_rank_count_parallel.py | 31 ++++++++++++++----- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index 7cd55bfa5..f957e5c9b 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -33,7 +33,10 @@ held once, and the zone is cut into the fewest pieces of at most 500 rows, ``ceil(rows / 500)``, whatever the rank count; - split with no zone keyed: the structural patch is one block, as in serial, and - the velocity solve takes at most 3 iterations (3 measured at np 1). + the velocity solve takes at most 3 iterations (3 measured at np 1); +- placed but not split, a keyed mask holding no cell: no patch, as in serial. A + keyed zone beats detection; the parallel branch fell back to the structural + split when no rank held a zone cell. This file carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial branch, in the parallel batches the parallel one. @@ -59,10 +62,10 @@ N_HAT = np.array([T_HAT[1], -T_HAT[0]]) -def _solve_fixture(split, fac_block_rows, zone_keyed=True): - """Solve the TI-band-beside-a-fault box, with the band keyed as the zone - unless ``zone_keyed`` is False; return the solver and every rank's - finest-level blocks (global rows).""" +def _solve_fixture(split, fac_block_rows, keyed="band"): + """Solve the TI-band-beside-a-fault box with ``keyed`` as the zone: "band" (the TI + band's cells), "empty" (a mask holding no cell) or None; return the solver + and every rank's finest-level blocks (global rows).""" main = A + (np.arange(N_RUNGS + 1) * S)[:, None] * T_HAT band = A + 0.08 * N_HAT + np.linspace(0.0, 1.0, N_RUNGS + 1)[:, None] * (B - A) base = uw.meshing.UnstructuredSimplexBox( @@ -95,8 +98,11 @@ def _solve_fixture(split, fac_block_rows, zone_keyed=True): stokes.petsc_use_pressure_nullspace = True stokes.tolerance = 1.0e-6 - band_cells = (np.asarray(mesh.cells_labelled("Band", 72), dtype=bool) - if zone_keyed else None) + band_cells = None + if keyed is not None: + band_cells = np.asarray(mesh.cells_labelled("Band", 72), dtype=bool) + if keyed == "empty": + band_cells = np.zeros_like(band_cells) set_custom_fmg(stokes, levels, field_id=0, fac_zone=band_cells, fac_block_rows=fac_block_rows) stokes.solve() @@ -166,7 +172,7 @@ def test_an_unsplit_zone_is_cut_into_the_same_pieces_at_any_rank_count(): @pytest.fixture(scope="module") def structural_alone(): - return _solve_fixture(split=True, fac_block_rows=2000, zone_keyed=False) + return _solve_fixture(split=True, fac_block_rows=2000, keyed=None) @pytest.mark.tier_b @@ -186,3 +192,12 @@ def test_the_structural_patch_alone_takes_the_serial_iterations(structural_alone its = max(info["vel_its_last"]) # 3 measured at np 1; one block per rank took 4 at np 2 and 6 at np 4 assert its <= 3, f"np {uw.mpi.size}: the velocity solve took {its} iterations (#925)" + + +@pytest.mark.tier_b +def test_an_empty_keyed_zone_installs_no_patch(): + _stokes, blocks = _solve_fixture(split=False, fac_block_rows=2000, keyed="empty") + sizes = sorted(int(b.size) for b in blocks) + assert not sizes, ( + f"np {uw.mpi.size}: a keyed zone with no cells installed blocks of {sizes}; " + f"serial installs none (#925)") From ef7804e9a33a73e0c98af0410e73d94db38c01a6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 16:30:43 +1100 Subject: [PATCH 08/13] #925 review 2: a keyed zone beats detection in parallel too; docstring and doc wording - The parallel branch decided whether to run the structural split from whether any rank held a zone cell; it now keys on fac_zone being given, as the serial branch does. A keyed mask with no cell gives no patch on an unsplit level and the split's rows as one more block on a split one. One allreduce fewer. - _gathered_patch_blocks's docstring said it cuts the zone as the serial branch does; a piece's cut can fall a node away from the serial one at a split pair (one site here, two nodes there). - The subsystem doc's summed rows are the owned rows; the subdomain rows follow the structural halo (#896). Underworld development team with AI support from Claude Code --- .../subsystems/fault-split-multigrid.md | 7 ++-- src/underworld3/utilities/custom_mg.py | 34 +++++++++---------- 2 files changed, 21 insertions(+), 20 deletions(-) diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index 15ab20996..6f6dad2be 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -73,9 +73,10 @@ largest block, and each rank's block sizes. It is made once, with one allgather Reading it makes no collective call, so any rank can read it alone. A driver that looped a reduction over its own rank's blocks hung at np 4, because the ranks made different numbers of calls. Each block of the finest level is held whole by one rank -(see below, #925), so the summed rows are the serial count at any rank count -(`test_1133`), and `block_rows_per_rank` shows where each block sits: the structural -LU on one rank. The record is each block as handed to the installer, before PCASM +(see below, #925), so the summed patch (owned) rows are the serial count at any rank +count (`test_1133`); the subdomain rows follow the structural halo, which can be a +few rows short of serial (#896). `block_rows_per_rank` shows where each block sits: +the structural LU on one rank. The record is each block as handed to the installer, before PCASM adds its overlap layer at set-up. ## Rotated walls on every level diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index d5a693555..1ca9d6c32 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -834,15 +834,17 @@ def _segment_block(rows, node_of_row, coords, max_rows): def _gathered_patch_blocks(zones, structural, max_rows, rstart, comm): """The finest level's FAC blocks in parallel, built from every rank's rows. - COLLECTIVE. Rank 0 gathers every rank's rows and cuts them as the serial - branch does: each mask is one band, segmented along its strike by + COLLECTIVE. Rank 0 gathers every rank's rows and cuts them into the serial + branch's blocks: each mask is one band, segmented along its strike by :func:`_segment_block`, and the structural rows outside every mask are one - more block; with no zone rows, the structural patch is one block. Each + more block; with no zone keyed, the structural patch is one block. Each block goes whole to the rank that owns most of its rows, so the structural block is one LU on one rank, as in serial. The blocks are therefore the same at any rank count, except that the structural halo is :func:`_fac_patch_split_parallel`'s, which can differ by a few rows where a - partition boundary runs along the band. + partition boundary runs along the band. A segmented piece's cut can fall + a node away from the serial one: the serial branch orders a split pair's + two nodes separately, this one keeps them as one site. Cut rank by rank, a rank's blocks held the seam rows its neighbour's blocks also held, and rank boundaries cut the band and the structural block (#925). @@ -852,7 +854,7 @@ def _gathered_patch_blocks(zones, structural, max_rows, rstart, comm): zones : list of (ndarray, ndarray) Per ``fac_zone`` mask, this rank's rows of the mask's cells (global, ghost rows resolved to their owners') and each row's node coordinates. - Empty where no zone is keyed. + An empty list where no zone is keyed; the same length on every rank. structural : (ndarray, ndarray) or None This rank's structural patch and subdomain rows (global, owned), or None where the split level's gate declined them or the level has no @@ -892,10 +894,8 @@ def _gathered_patch_blocks(zones, structural, max_rows, rstart, comm): held = [e[2] for e in everyone if e[2] is not None] if held: sub = np.unique(np.concatenate([s for _o, s in held]).astype(np.int64)) - zone_rows = (np.concatenate(covered) if covered - else np.zeros(0, dtype=np.int64)) - if zone_rows.size: - extra = np.setdiff1d(sub, zone_rows) + if zones: + extra = np.setdiff1d(sub, np.concatenate(covered)) if extra.size: blocks.append((extra, extra)) elif sub.size: @@ -1887,14 +1887,15 @@ def provenance_transfer(): # 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. + # A keyed zone beats detection, as in the serial branch, + # even where its masks hold no cell on any rank. The + # split marker and the keying are the same on every rank, + # so the cover gate's collectives run on all or none. 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 + keyed = self.fac_zone is not None mpi_comm = comm.tompi4py() - zone_anywhere = mpi_comm.allreduce( - any(z[0].size for z in zones), op=MPI.LOR) - if split_level or not zone_anywhere: + if split_level or not keyed: gate = {} split = _fac_patch_split_parallel( P, coords[l - 1], coords[l], maps[l - 1], lay, nc, @@ -1904,14 +1905,13 @@ def provenance_transfer(): self.fac_gate[l] = gate or None # split is None on every rank or on none: the gate's # verdict is collective - if zone_anywhere or split is not None: + if keyed or split is not None: # TODO(MEASURE): #670 knob — rows per block; the # default comes from the size sweep. cap = int(os.environ.get("UW_FAC_BLOCK_ROWS", self.fac_block_rows)) blocks = _gathered_patch_blocks( - zones if zone_anywhere else [], split, cap, - lay.rstart, mpi_comm) + zones, split, cap, lay.rstart, mpi_comm) # every rank installs the patch, or none does: a rank with # no rows of its own still sets the smoother type n_any = comm.tompi4py().allreduce(len(blocks), op=MPI.MAX) From 9b69572c00fb67fc636d128b63820095f8f49444 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 19:02:49 +1100 Subject: [PATCH 09/13] #925 review 2: the structural block holds every rank's rows; scaling recorded The structural patch alone stays one gathered block. On row 5 with no zone at n_rungs 320 (kaiju), cut per rank the velocity solve took 67, 128, 131, 316 and 445 iterations at np 1, 2, 4, 8 and 16 (33 s to 112 s); as one block, 67 to 86 iterations and 22 to 37 s, with the owning rank at about twice the others' memory at every rank count. - test_1139: the structural block's size equals the gate's count of every rank's patch and halo rows, taken before the gather (a guard: no code lost a rank's rows); the iteration bound's comment records 3 at np 1, 2 and 4. - fault-split-multigrid.md: the structural block is one unbounded LU on one rank (#895), with the scaling numbers. - The 2026-08 design record's notes cover the structural patch as well as the zone. Underworld development team with AI support from Claude Code --- .../design/fault-patch-multigrid-2026-08.md | 12 +++++++----- .../subsystems/fault-split-multigrid.md | 9 +++++++-- ...fac_zone_blocks_any_rank_count_parallel.py | 19 ++++++++++++++----- 3 files changed, 28 insertions(+), 12 deletions(-) diff --git a/docs/developer/design/fault-patch-multigrid-2026-08.md b/docs/developer/design/fault-patch-multigrid-2026-08.md index 533a52079..541ba490e 100644 --- a/docs/developer/design/fault-patch-multigrid-2026-08.md +++ b/docs/developer/design/fault-patch-multigrid-2026-08.md @@ -273,15 +273,17 @@ to permit it.) cuts at seams are measured free (the k-ladder). > [2026-11-03] Not for a zone-keyed patch: cut at rank boundaries (a TI > band beside a split fault), its blocks took the velocity solve from 4 - > iterations serially to 118 at np 2 and past 200 at np 4. The zone's blocks - > are now built from every rank's rows (#925, `test_1139`). + > iterations serially to 118 at np 2 and past 200 at np 4. The structural + > patch alone, cut the same way, took 67 iterations serially and 445 at + > np 16 (row 5, n_rungs 320). Both are now built from every rank's rows + > (#925, `test_1139`). 2. **PCASM is rank-local.** Blocks become per-rank lists with local indices; faultless ranks carry `nsd = 0`, and every IS-building path must survive empty sets (the empty-stratum `getIndices()` SEGV family, #589). - > [2026-11-03] A zone-keyed level's blocks are per-rank lists of GLOBAL - > rows: each block goes whole to the rank owning most of its rows and may - > hold other ranks' rows (#925). + > [2026-11-03] The finest level's blocks, a keyed zone's or the structural + > patch's, are per-rank lists of GLOBAL rows: each block goes whole to the + > rank owning most of its rows and may hold other ranks' rows (#925). 3. **Pair detection must go collective or be replaced by zone keying.** Rank-local coordinate matching misses pairs that straddle a boundary and can false-match ghosts; the zone mask is rank-local by nature and diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index 6f6dad2be..b4b2ac366 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -141,8 +141,8 @@ path. | S 0.004, blocks of 2000 rows | 433 s | 1404 s | 2089 / 2058 s | 2248 / 2981 s | | S 0.015625, one block | 516 s | 766 s | 931 / 869 s | 923 / 1207 s | -The size of the LU blocks is bounded by `fac_block_rows` for keyed zones; the -structural block is one block per rank (#895). +The size of the LU blocks is bounded by `fac_block_rows` for keyed zones. The +structural block is not bounded: it is one LU, held by one rank in parallel (#895). ## Parallel @@ -297,6 +297,11 @@ the velocity solve then took 4 iterations serially, 118 at np 2 and more than 20 np 4 (#925). Gathered, it takes 4 at np 1, 2 and 4 (`test_1139`). The structural patch alone, one block per rank, converged but took 3 iterations serially, 4 at np 2 and 6 at np 4; gathered it is one block at any rank count (`test_1139`, `test_1133`). +On row 5 with no zone at n_rungs 320 (30 048 rows in the block; kaiju, 2026-11), cut +per rank the velocity solve took 67, 128, 131, 316 and 445 iterations at np 1, 2, 4, +8 and 16, and the first solve 33 s serially and 112 s at np 16. As one block it took +67 to 86 iterations and 22 to 37 s at every rank count. The owning rank's peak memory +was about twice the others' (520-540 MB against 250-290 MB) at np 2 and at np 16 alike. ## Not yet done diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index f957e5c9b..e94a870ff 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -19,7 +19,9 @@ The structural patch alone, with no zone keyed, was one block per rank: on the same fixture its velocity solve took 3 iterations serially, 4 at np 2 and 6 at np 4 (70, 96 and 143 summed over the solve). It converged, but a partition -boundary cut it all the same, and it too is built from the gathered rows. +boundary cut it all the same, and it too is built from the gathered rows. At +n_rungs 320, cut per rank, it took 67, 128, 131, 316 and 445 iterations at np 1, +2, 4, 8 and 16; as one block, 67 to 86. The fixture is that problem at n_rungs 80: UWTN 2026-022's problem A box, the main fault split, and a TI band (eta_1 = 0.01 along planes parallel to it) placed @@ -32,8 +34,9 @@ - placed but not split, the zone alone at ``fac_block_rows=500``: every row is held once, and the zone is cut into the fewest pieces of at most 500 rows, ``ceil(rows / 500)``, whatever the rank count; -- split with no zone keyed: the structural patch is one block, as in serial, and - the velocity solve takes at most 3 iterations (3 measured at np 1); +- split with no zone keyed: the structural patch is one block holding every + rank's patch and halo rows, as in serial, and the velocity solve takes at most + 3 iterations (3 measured at np 1, 2 and 4); - placed but not split, a keyed mask holding no cell: no patch, as in serial. A keyed zone beats detection; the parallel branch fell back to the structural split when no rank held a zone cell. @@ -177,11 +180,16 @@ def structural_alone(): @pytest.mark.tier_b def test_the_structural_patch_alone_is_one_block(structural_alone): - _stokes, blocks = structural_alone + stokes, blocks = structural_alone sizes = sorted(int(b.size) for b in blocks) assert len(sizes) == 1, ( f"np {uw.mpi.size}: the structural patch is in blocks of {sizes}; serial " f"has one (#925)") + # the gate counts every rank's patch and halo rows before they are gathered + gate = stokes._custom_mg["hierarchy"].fac_gate[-1] + assert sizes[0] == gate["subdomain_rows"], ( + f"np {uw.mpi.size}: the block holds {sizes[0]} rows, the ranks' patch and " + f"halo rows are {gate['subdomain_rows']} (#925)") @pytest.mark.tier_b @@ -190,7 +198,8 @@ def test_the_structural_patch_alone_takes_the_serial_iterations(structural_alone info = stokes._rotated_freeslip_info assert info["converged"] its = max(info["vel_its_last"]) - # 3 measured at np 1; one block per rank took 4 at np 2 and 6 at np 4 + # 3 measured at np 1, 2 and 4 (kaiju) and np 1 (macOS); one block per rank + # took 4 at np 2 and 6 at np 4 assert its <= 3, f"np {uw.mpi.size}: the velocity solve took {its} iterations (#925)" From 2e11a6f53ba0e8c390354d83b8fe9ac215b45955 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 19:12:38 +1100 Subject: [PATCH 10/13] #925 review 3: tests for keying that differs between ranks or holds no mask test_1139 constructs a hierarchy with fac_zone keyed on rank 0 only (others None), and with an empty list of masks; both must raise ValueError at construction. On this library neither raises: the first would hang in the build (rank 0 in the zone's gather, the others in the cover gate's allreduce), and the second is keyed with no zone. The empty-list test fails serially here (DID NOT RAISE); the per-rank test runs at np >= 2. Underworld development team with AI support from Claude Code --- ...fac_zone_blocks_any_rank_count_parallel.py | 31 ++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index e94a870ff..39e66ff33 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -40,6 +40,11 @@ - placed but not split, a keyed mask holding no cell: no patch, as in serial. A keyed zone beats detection; the parallel branch fell back to the structural split when no rank held a zone cell. +- the keying is refused at construction unless it is the same on every rank (a + mask on every rank, or None on every rank) and holds a mask: rank 0 keyed and + the others not would wait in the zone's gather on rank 0 and in the cover + gate's allreduce on the others, and an empty list of masks would be keyed but + hold no zone. This file carries no ``mpi(min_size=2)`` mark: at np 1 it checks the serial branch, in the parallel batches the parallel one. @@ -52,7 +57,7 @@ import sympy import underworld3 as uw -from underworld3.utilities.custom_mg import set_custom_fmg +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] @@ -210,3 +215,27 @@ def test_an_empty_keyed_zone_installs_no_patch(): assert not sizes, ( f"np {uw.mpi.size}: a keyed zone with no cells installed blocks of {sizes}; " f"serial installs none (#925)") + + +def _hierarchy_keyed(zone_of): + """Construct, not build, a hierarchy on a plain box with ``zone_of(n_cells)`` as + this rank's ``fac_zone``.""" + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0, 0), maxCoords=(1, 1), cellSize=1 / 4, regular=False, + qdegree=2, refinement=1) + levels = base._coarse_level_meshes()[:-1] + [base] + c_start, c_end = base.dm.getHeightStratum(0) + return CustomMGHierarchy(levels, fac_zone=zone_of(c_end - c_start)) + + +@pytest.mark.tier_b +def test_an_empty_list_of_masks_is_refused(): + with pytest.raises(ValueError, match="no mask"): + _hierarchy_keyed(lambda n: []) + + +@pytest.mark.tier_b +@pytest.mark.skipif(uw.mpi.size < 2, reason="the keying can differ only between ranks") +def test_a_zone_keyed_on_some_ranks_only_is_refused(): + with pytest.raises(ValueError, match="every rank"): + _hierarchy_keyed(lambda n: np.ones(n, dtype=bool) if uw.mpi.rank == 0 else None) From c094073dfefd69b3e492cf75efeaa70989114bbb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 19:14:43 +1100 Subject: [PATCH 11/13] #925 review 3: refuse keying that differs between ranks or holds no mask - CustomMGHierarchy checks fac_zone collectively at construction: one allgather of each rank's mask count (-1 for None). If they differ, every rank raises ValueError, where the build would otherwise wait in the zone's gather on some ranks and the cover gate's allreduce on others. This takes the place of the zone_anywhere allreduce ef7804e9 removed. - An empty list of masks raises: it was keyed but held no zone, so the structural block beside the zone was (owned, sub) in parallel and (sub, sub) in serial. - Docstrings ask for a mask on every rank (all-False where a rank holds no zone cell) or None on every rank. - fault-split-multigrid.md: owned rows equal the serial count only with no zone keyed; with a zone, both counts carry the structural halo. The n_rungs 320 scaling paragraph compares iterations only: np 1 is the same code in both series, and the per-rank series ran before the stack's development merges, so the seconds were not comparable. Serial: test_1139/1133/1135 12 passed (the per-rank test runs at np >= 2), test_1022 fac_zone passed. All library producers pass an array or None. Underworld development team with AI support from Claude Code --- .../subsystems/fault-split-multigrid.md | 20 ++++++------ src/underworld3/utilities/custom_mg.py | 31 +++++++++++++++---- ...fac_zone_blocks_any_rank_count_parallel.py | 4 +-- 3 files changed, 38 insertions(+), 17 deletions(-) diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index b4b2ac366..d9299ef9d 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -72,10 +72,12 @@ largest block, and each rank's block sizes. It is made once, with one allgather `CustomMGHierarchy.install` and the rotated path's two `_configure_pcmg` calls. Reading it makes no collective call, so any rank can read it alone. A driver that looped a reduction over its own rank's blocks hung at np 4, because the ranks made -different numbers of calls. Each block of the finest level is held whole by one rank -(see below, #925), so the summed patch (owned) rows are the serial count at any rank -count (`test_1133`); the subdomain rows follow the structural halo, which can be a -few rows short of serial (#896). `block_rows_per_rank` shows where each block sits: +different numbers of calls. With no zone keyed, the summed patch (owned) rows are the +serial count at any rank count (`test_1133`), and the subdomain rows follow the +structural halo, which can be a few rows short of serial (#896). With a zone keyed, +its blocks and the structural block beside them hold the same rows as patch and +subdomain, so both counts carry that halo. Each block of the finest level is held +whole by one rank (see below, #925), and `block_rows_per_rank` shows where each sits: the structural LU on one rank. The record is each block as handed to the installer, before PCASM adds its overlap layer at set-up. @@ -297,11 +299,11 @@ the velocity solve then took 4 iterations serially, 118 at np 2 and more than 20 np 4 (#925). Gathered, it takes 4 at np 1, 2 and 4 (`test_1139`). The structural patch alone, one block per rank, converged but took 3 iterations serially, 4 at np 2 and 6 at np 4; gathered it is one block at any rank count (`test_1139`, `test_1133`). -On row 5 with no zone at n_rungs 320 (30 048 rows in the block; kaiju, 2026-11), cut -per rank the velocity solve took 67, 128, 131, 316 and 445 iterations at np 1, 2, 4, -8 and 16, and the first solve 33 s serially and 112 s at np 16. As one block it took -67 to 86 iterations and 22 to 37 s at every rank count. The owning rank's peak memory -was about twice the others' (520-540 MB against 250-290 MB) at np 2 and at np 16 alike. +On row 5 with no zone at n_rungs 320 (30 048 rows in the block; kaiju, 2026-11), both +take 67 iterations serially, where the code is the same. Cut per rank, the velocity +solve took 128, 131, 316 and 445 iterations at np 2, 4, 8 and 16; as one block, 76 to +86. The owning rank's peak memory was about twice the others' (520-540 MB against +250-290 MB) at np 2 and at np 16 alike. ## Not yet done diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 1ca9d6c32..144ff5ae6 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -1558,6 +1558,8 @@ class CustomMGHierarchy: 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. + Pass a mask on every rank, all-False where a rank holds no zone cell, + or None on every rank. """ def __init__(self, level_meshes, builder="barycentric", field_id=None, @@ -1605,12 +1607,28 @@ def _validated_fac_zone(self, fac_zone): Validated against the FINEST level's cell count at construction — a mask of the wrong length would otherwise surface as a silently empty (or wrong) patch, which is exactly the declining-quietly - failure mode the ASM probe discipline exists to catch.""" - if fac_zone is None: - return None - masks = (list(fac_zone) if isinstance(fac_zone, (list, tuple)) - else [fac_zone]) + failure mode the ASM probe discipline exists to catch. + + Collective: the keying must be the same on every rank, a mask on + every rank (all-False where the rank holds no zone cell) or None on + every rank. The parallel finest level gathers the masks' rows, so a + rank keyed differently would wait in another collective (#925).""" + masks = None + if fac_zone is not None: + masks = (list(fac_zone) if isinstance(fac_zone, (list, tuple)) + else [fac_zone]) finest = self.level_meshes[-1] + counts = finest.dm.comm.tompi4py().allgather( + -1 if masks is None else len(masks)) + if len(set(counts)) > 1: + raise ValueError( + f"fac_zone: the ranks pass {counts} masks (-1 for None); pass " + f"a mask on every rank, all-False where a rank holds no zone " + f"cell, or None on every rank.") + if masks is None: + return None + if not masks: + raise ValueError("fac_zone holds no mask; pass None for no zone.") cS, cE = finest.dm.getHeightStratum(0) out = [] for k, mk in enumerate(masks): @@ -2201,7 +2219,8 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", 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). + factorization). In parallel pass a mask on every rank, all-False where a + rank holds no zone cell, or None on every rank. ``verbose`` prints the installed transfer sizes and the FAC cover gate's verdict per level (``hierarchy.fac_gate``, #918).""" diff --git a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py index 39e66ff33..ac925ced0 100644 --- a/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -20,8 +20,8 @@ same fixture its velocity solve took 3 iterations serially, 4 at np 2 and 6 at np 4 (70, 96 and 143 summed over the solve). It converged, but a partition boundary cut it all the same, and it too is built from the gathered rows. At -n_rungs 320, cut per rank, it took 67, 128, 131, 316 and 445 iterations at np 1, -2, 4, 8 and 16; as one block, 67 to 86. +n_rungs 320, cut per rank, it took 128, 131, 316 and 445 iterations at np 2, 4, 8 +and 16 against 67 serially; as one block, 76 to 86. The fixture is that problem at n_rungs 80: UWTN 2026-022's problem A box, the main fault split, and a TI band (eta_1 = 0.01 along planes parallel to it) placed From a2e3839612600d57928fbf2f55cbcc1afc2382a0 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 21:22:16 +1100 Subject: [PATCH 12/13] #925: the owning rank's memory, as measured The doc said the owning rank used about twice the others' memory at np 2 and np 16 alike. Kaiju report 15b: the owner peaked at 512-541 MB at np 2-16 while the others fell from 379 MB (np 2) to 248-284 MB (np 16), so the ratio is 1.35 at np 2 and 1.6-2.2 at np 4-16. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/fault-split-multigrid.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/developer/subsystems/fault-split-multigrid.md b/docs/developer/subsystems/fault-split-multigrid.md index 553bd33af..0bb75506e 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -302,8 +302,9 @@ and 6 at np 4; gathered it is one block at any rank count (`test_1139`, `test_11 On row 5 with no zone at n_rungs 320 (30 048 rows in the block; kaiju, 2026-11), both take 67 iterations serially, where the code is the same. Cut per rank, the velocity solve took 128, 131, 316 and 445 iterations at np 2, 4, 8 and 16; as one block, 76 to -86. The owning rank's peak memory was about twice the others' (520-540 MB against -250-290 MB) at np 2 and at np 16 alike. +86. The owning rank's peak memory was 512-541 MB at np 2-16; the others' +fell from 379 MB at np 2 to 248-284 MB at np 16 (1.35 times at np 2, 1.6-2.2 times +at np 4-16). ## Not yet done From e8957186c7525ab22b148c3ee85d1b8f0a900bab Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 3 Nov 2026 23:51:27 +1100 Subject: [PATCH 13/13] CI shards: move serial work off shard 4, which hit the 60-minute cap Shard 4 owns the whole --p N phase, and every new file under tests/parallel/ had also put its serial run in batch 4. It reached 55:32 on development (run 41564997729, #930) and was cancelled at 60 minutes on this PR, during test_1089 at np 2, which had not hung: development takes 48 s past the same point. Shards 2 and 3 took about 29 minutes. - The ten serial runs of tests/parallel files (test_1069 ... test_1139), about 6 minutes, go to shard 3. - test_05*/06*/07*, about 4 minutes, go to shard 2. --print-shards still gives [1,2,3,4], and check_test_coverage reports every file reachable. Underworld development team with AI support from Claude Code --- scripts/test.sh | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/scripts/test.sh b/scripts/test.sh index 6f99eaf02..aaf28de12 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -129,6 +129,12 @@ fi # Re-measure over several runs before re-assigning: runner variance here is of # the same order as the imbalance left. # +# 2026-11-03: shard 4 had grown to 55:32 on development (run 41564997729, #930) +# and ran past the 60-minute cap on PR #934, against 29 min for shards 2 and 3. +# Shard 4 owns the whole --p N phase, and every new file under tests/parallel/ +# had also put its serial run in batch 4. Those serial runs (about 6 min) moved +# to shard 3, and test_05*/06*/07* (about 4 min) to shard 2. +# # The floor is test_106*/test_107*, and it is a TAIL, not a volume: 48% of its # 149 tests finish in the first 2 minutes, then the remainder takes 31. Since # the distribution is --dist loadfile, a file is atomic on a worker, so the @@ -298,7 +304,7 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # whole band was measured 2026-09-12 and passes; a band is not a unit of trust, # and holding one back for a suspicion nobody recorded meant its level_1 files # never ran at all. - batch 4 tests/test_05*py tests/test_06*py tests/test_07*py + batch 2 tests/test_05*py tests/test_06*py tests/test_07*py # Units system tests (0800-0899: unit-aware functions, arrays, and conversions) batch 3 tests/test_08*py @@ -319,26 +325,26 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # spherical shell at np=1). Run it here so the serial job covers that path — # every other serial test of the default boundary normal is on a box, where # flat walls make the question vacuous. It also runs in the --p N batch below. - batch 4 tests/parallel/test_1069_boundary_normal_parallel.py + batch 3 tests/parallel/test_1069_boundary_normal_parallel.py # The live multigrid configuration (#892) is the same check in serial and in # parallel, so it is one file under tests/parallel/ with no min_size mark too. - batch 4 tests/parallel/test_1090_mg_live_configuration_parallel.py + batch 3 tests/parallel/test_1090_mg_live_configuration_parallel.py # So is the live finest-smoother operator of the rotated Newton loop (#894). - batch 4 tests/parallel/test_1088_rotated_interface_operator_parallel.py + batch 3 tests/parallel/test_1088_rotated_interface_operator_parallel.py # So is the band patch, the same rows at any rank count (#896). - batch 4 tests/parallel/test_1106_fac_cover_partition_parallel.py + batch 3 tests/parallel/test_1106_fac_cover_partition_parallel.py # So is the static mark of a moved base finest (#919). - batch 4 tests/parallel/test_1132_addon_partition_static_mark_parallel.py + batch 3 tests/parallel/test_1132_addon_partition_static_mark_parallel.py # So is the FAC patch summary, the same on every rank (test_1133). - batch 4 tests/parallel/test_1133_fac_patch_summary_parallel.py + batch 3 tests/parallel/test_1133_fac_patch_summary_parallel.py # So is the serial reference child's environment (#922). - batch 4 tests/parallel/test_1138_serial_reference_child_env.py + batch 3 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 + batch 3 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 + batch 3 tests/parallel/test_1135_fac_zone_split_structural_parallel.py # So are the zone-keyed patch's blocks, the same at any rank count (#925). - batch 4 tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py + batch 3 tests/parallel/test_1139_fac_zone_blocks_any_rank_count_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 —