diff --git a/docs/developer/design/fault-patch-multigrid-2026-08.md b/docs/developer/design/fault-patch-multigrid-2026-08.md index 967ae1c76..541ba490e 100644 --- a/docs/developer/design/fault-patch-multigrid-2026-08.md +++ b/docs/developer/design/fault-patch-multigrid-2026-08.md @@ -271,10 +271,19 @@ 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 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] 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 8653bd7d2..075322adc 100644 --- a/docs/developer/subsystems/fault-split-multigrid.md +++ b/docs/developer/subsystems/fault-split-multigrid.md @@ -72,11 +72,14 @@ 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. 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. ## Rotated walls on every level @@ -140,8 +143,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 @@ -296,6 +299,28 @@ 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 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`). +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 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 - A regression test for the coarse-level wall rotation itself (a rotated free-slip diff --git a/scripts/test.sh b/scripts/test.sh index 2304e88a9..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,24 +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 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 — diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 813125e5d..dd5081cd6 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -625,12 +625,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 []): @@ -834,6 +834,84 @@ def _segment_block(rows, node_of_row, coords, max_rows): return pieces +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 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 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. 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). + + 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. + 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 + 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 ``(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) + 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: + 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([s for _o, s in held]).astype(np.int64)) + if zones: + extra = np.setdiff1d(sub, np.concatenate(covered)) + 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 owned, sub in blocks: + owner = np.searchsorted(starts, sub, side="right") - 1 + per_rank[int(np.bincount(owner, minlength=len(everyone)).argmax())].append( + (owned, sub)) + return comm.scatter(per_rank, root=0) + + 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' @@ -1483,6 +1561,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, @@ -1530,12 +1610,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): @@ -1779,37 +1875,28 @@ 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. 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 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] 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, :]] + n, c = np.nonzero(rows >= 0) + 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 @@ -1821,14 +1908,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(len(blocks) > 0, - 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, @@ -1836,16 +1924,15 @@ def provenance_transfer(): cover_max=self._fac_cover_max(l), gate=gate) self.fac_gate[l] = gate or None - if 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))) + # split is None on every rank or on none: the gate's + # verdict is collective + 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, 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) @@ -2139,7 +2226,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_1133_fac_patch_summary_parallel.py b/tests/parallel/test_1133_fac_patch_summary_parallel.py index 9757404ef..803632680 100644 --- a/tests/parallel/test_1133_fac_patch_summary_parallel.py +++ b/tests/parallel/test_1133_fac_patch_summary_parallel.py @@ -18,13 +18,16 @@ 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 + 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,17 +101,10 @@ 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) + 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 @pytest.mark.tier_b @@ -134,6 +136,9 @@ 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)") + 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 entry = h.level_patch_rows[-1] blocks = entry if isinstance(entry, list) else [entry] 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 new file mode 100644 index 000000000..ac925ced0 --- /dev/null +++ b/tests/parallel/test_1139_fac_zone_blocks_any_rank_count_parallel.py @@ -0,0 +1,241 @@ +"""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 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. At +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 +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; +- 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. +- 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. + + 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 CustomMGHierarchy, set_custom_fmg +from underworld3.utilities.place_surface import place_fault_ribbon_2d + +pytestmark = [pytest.mark.level_2] + +A, B = np.array([0.30, 0.45]), np.array([0.70, 0.55]) +N_RUNGS = 80 +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, 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( + 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 = 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() + 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 (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)") + + +@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)) + 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)") + + +@pytest.fixture(scope="module") +def structural_alone(): + return _solve_fixture(split=True, fac_block_rows=2000, keyed=None) + + +@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)") + # 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 +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, 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)" + + +@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)") + + +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)