From af1bec36c9df87a35768915733b2a927a6c4c5ea Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 13:04:11 -0700 Subject: [PATCH 01/19] Placement gathers the shell the seam rule needs, and nothing when the region is already interior (#670) The gather-first surgery moves base cells onto one rank so that no point it deletes or creates is shared. That needs three layers: the cells the carve drops, their vertex star (the ring the fill attaches to), and one more layer so the ring's points are unshared. _gather_region grows the star and the layer from the mark, but the mark itself reached two cell widths beyond the carve, so the margin was paid twice: on the crossing-patches fixture the gather moved 5087 of 5592 base cells for a cavity of 334, and every rank count ended with the surgery rank holding 95% of the mesh. The mark now covers what the carve drops: the victims within the clearance plus the crossed cells' vertices, one cell diameter out. Every placement path carries the same rule (lines, sheet, the 2-D and 3-D thin volumes, the ribbon). On a box large enough to hold a shell (eight times the fixture, 27,544 base cells, two crossing patches) the gather moves 5046 cells instead of 9831, identical at np=2 and np=4; the remaining extra load on the surgery rank is the shell plus the band and fill it creates, which is the design's accepted trade. A region whose star and layer already sit on one rank is placed there with no redistribution at all. The decision is collective: the per-rank star counts are gathered before it is taken. _gather_region returns the moved cell count, and the sheet and thin volume report it as info["n_gathered"]. ptest_0855 bounds it by the cells within three median cell diameters of the zone (7269 on that box): the old mask fails the test, the new one passes it at np=2, 3 and 4. The other placement paths' parallel tests and the serial thin-volume suite pass unchanged. Not changed here, recorded in #670: the network is still one region with one target rank, so two surfaces interior to different ranks are gathered together; and the moved cells stay where the surgery put them. The np=3 failures of the network solve predate this change and are #671. Underworld development team with AI support from Claude Code --- .../conforming-surfaces-and-fault-zones.md | 25 +++++++++ src/underworld3/utilities/place_surface.py | 53 ++++++++++++------- .../ptest_0855_place_thin_volume_parallel.py | 42 +++++++++++++++ 3 files changed, 102 insertions(+), 18 deletions(-) diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index 41b0813c5..c0b50096c 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -291,6 +291,31 @@ bit-identical over np = 1..5. Every refusal is collective: all ranks raise the same error, or none does. The 2-D forms are serial; a parallel call is refused rather than returning a mesh whose star-forest is silently wrong. +What the gather moves, and why, is base mesh — never the surface. The +surgery deletes base cells around the surface and creates new ones in the +cavity, and the rebuild carries the old star forest over by renumbering, so +no point the surgery deletes or creates may be shared. That needs exactly +three layers of base cells on one rank: the cells the carve drops, their +vertex star (the ring the fill attaches to), and one more layer so the +ring's own points are unshared. The mark covers the dropped cells (the +victims within the clearance plus the crossed cells' vertices, within one +cell diameter of the surface) and `_gather_region` grows the star and the +layer from it. The result is a shell about three cells thick around the +surface; a region whose shell is already interior to one rank is placed +there with nothing moved. The moved count is reported as +`info["n_gathered"]`, and `ptest_0855` bounds it by the cells within three +median cell diameters of the zone. + +The gather is one-way: the moved cells stay on the surgery rank, together +with the cells the fill creates, so that rank carries the shell plus the +band as extra load. That is the accepted trade (extra load on one rank in +exchange for no communication during the solve), and it is proportional to +the surface, not the domain. Two limits remain, recorded in #670: the whole +network is one region with one target rank, so two surfaces each interior to +a different rank are still gathered together; and a small domain cannot hold +a shell at all (three base cells reaching the walls is the whole box, which +is what the crossing-patches test fixture does). + ## The thin volume: finite-width zones, junctions in the volume `place_surface.place_thin_volume(dm, patches, width)` embeds a layer of real diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index c83b33a62..95e0b872f 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -985,7 +985,7 @@ def _place_one_parallel(dm, pts, label, label_value, clearance, spacing): comm.Allreduce(MPI.IN_PLACE, area_before, op=MPI.SUM) mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_line < (clearance + 2.0) * h_vertex) + mark[np.flatnonzero(d_line < (clearance + 1.0) * h_vertex) + vS - pStart] = 1 dm_work, moved = _gather_region(dm, mark) if moved: @@ -1606,11 +1606,14 @@ def _gather_region(dm, vertex_mark_chart, verbose=False): (feature/fault-split-node c8693579 / 1d487319; the label-driven original is measured at np=2..8 with serial-identical topology). Only the marked star moves — everything else keeps its load-balanced home — via a shell - partitioner. Returns ``(new_dm, moved)``; the input is untouched. + partitioner. Returns ``(new_dm, n_moved)``: the global count of cells + gathered, ``0`` when nothing moved (serial, or a region already + interior to one rank), so callers can both branch on it and report + it. The input is untouched. """ comm = dm.getComm().tompi4py() if comm.size == 1: - return dm, False + return dm, 0 work = dm.clone() cS, cE = work.getHeightStratum(0) @@ -1644,6 +1647,15 @@ def star_of_marked(m): counts = np.asarray(comm.allgather(len(star))) if counts.sum() == 0: raise ValueError("place_sheet: the sheet meets no cell on any rank") + if np.count_nonzero(counts) == 1: + # The region, star and layer included, is already interior to one + # rank: the seam rule holds where the mesh stands, and the surgery + # runs there with no cell moved (#670). COLLECTIVE decision — the + # counts are gathered, so every rank takes this branch together. + if verbose: + uw.pprint(f"[place_sheet] a {int(counts.sum())}-cell region is " + f"interior to rank {int(np.argmax(counts))}; no gather") + return dm, 0 target = int(np.argmax(counts)) assign = np.full(cE - cS, comm.rank, dtype=np.int32) @@ -1659,7 +1671,7 @@ def star_of_marked(m): if verbose: uw.pprint(f"[place_sheet] gathered a {int(counts.sum())}-cell region " f"onto rank {target}") - return work, True + return work, int(counts.sum()) def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, @@ -2972,14 +2984,17 @@ def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1, cells = _tet_vertices(dm) h_vertex, _h_cell = _vertex_h_3d(dm, cells, len(X)) d_sheet = _sheet_distance_within(X, sheet_pts, sheet_tris, - (clearance + 2.0) * h_vertex) - # The gather mask is a SUPERSET of everything the carve may touch: the - # victims (clearance) plus the crossed cells' vertices, which sit within - # a cell diameter of the sheet. The +2 margin covers grading between - # neighbouring cells; the carve asserts nothing shared afterwards, so an - # under-reach is loud, never silent. + (clearance + 1.0) * h_vertex) + # The gather mask covers everything the carve may DROP: the victims + # (clearance) plus the crossed cells' vertices, which sit within a cell + # diameter of the sheet. The seam rule needs the dropped cells, the ring + # (their vertex star) and one more layer so the ring's points are + # unshared; _gather_region grows the star and the layer from this mask, + # so a wider mask here only moves cells the surgery never touches (#670: + # a +2 margin gathered 91% of a fixture whose cavity was 6%). The carve + # asserts nothing shared afterwards, so an under-reach is loud. mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_sheet < (clearance + 2.0) * h_vertex) + mark[np.flatnonzero(d_sheet < (clearance + 1.0) * h_vertex) + vS - pStart] = 1 volume_before = np.array( @@ -2994,7 +3009,7 @@ def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1, cells = _tet_vertices(dm_work) h_vertex, _h_cell = _vertex_h_3d(dm_work, cells, len(X)) d_sheet = _sheet_distance_within(X, sheet_pts, sheet_tris, - (clearance + 2.0) * h_vertex) + (clearance + 1.0) * h_vertex) on_wall = _true_wall_vertex_mask(dm_work, len(X)) shared = _shared_point_flags(dm_work).astype(bool) @@ -3324,6 +3339,7 @@ def new_id(v): uw.pprint(f"[place_sheet {label!r}] placed {info['n_placed']} " f"vertices, removed {info['n_removed']}; " f"{info['n_surface_facets']} sheet faces") + info["n_gathered"] = int(moved) # cells the gather moved (#670) return new, info @@ -5836,7 +5852,7 @@ def _remove_embedded_2d(dm, label, label_value, clearance, verbose): reach_v = clearance * h_vertex mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_skin < reach_v + 2.0 * h_vertex) + mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) + vS - pStart] = 1 cS0, cE0 = dm.getHeightStratum(0) if dm.hasLabel(label): @@ -6082,7 +6098,7 @@ def remove_embedded(dm, label, label_value=1, clearance=0.6, verbose=False): d_skin = _sheet_distance(X, soup_pts, soup_tris) reach_v = clearance * h_vertex mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_skin < reach_v + 2.0 * h_vertex) + mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) + vS - pStart] = 1 # The object's own vertices must gather with it, however fat the zone. cS0, _cE0 = dm.getHeightStratum(0) @@ -6365,7 +6381,7 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, comm.Allreduce(MPI.IN_PLACE, area_before, op=MPI.SUM) mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_skin < reach_v + 2.0 * h_vertex) + mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) + vS - pStart] = 1 dm_work, moved = _gather_region(dm, mark, verbose=verbose) if moved: @@ -6955,9 +6971,9 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, h_vertex, _h_cell = _vertex_h_3d(dm, cells, len(X)) reach_v = np.maximum(clearance * h_vertex, 0.6 * width) d_skin = _sheet_distance_within(X, skin_xyz, skin_tris, - reach_v + 2.0 * h_vertex) + reach_v + 1.0 * h_vertex) mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_skin < reach_v + 2.0 * h_vertex) + mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) + vS - pStart] = 1 volume_before = np.array([_owned_cell_volume(dm)], dtype=float) @@ -6972,7 +6988,7 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, h_vertex, _h_cell = _vertex_h_3d(dm_work, cells, len(X)) reach_v = np.maximum(clearance * h_vertex, 0.6 * width) d_skin = _sheet_distance_within(X, skin_xyz, skin_tris, - reach_v + 2.0 * h_vertex) + reach_v + 1.0 * h_vertex) on_wall = _true_wall_vertex_mask(dm_work, len(X)) shared = _shared_point_flags(dm_work).astype(bool) @@ -7320,6 +7336,7 @@ def new_id(v): f"zone cells, {info['n_skin_faces']} skin faces; placed " f"{info['n_placed']} vertices, removed " f"{info['n_removed']}") + info["n_gathered"] = int(moved) # cells the gather moved (#670) return new, info diff --git a/tests/parallel/ptest_0855_place_thin_volume_parallel.py b/tests/parallel/ptest_0855_place_thin_volume_parallel.py index 9443ac1ff..64b7f29de 100644 --- a/tests/parallel/ptest_0855_place_thin_volume_parallel.py +++ b/tests/parallel/ptest_0855_place_thin_volume_parallel.py @@ -256,3 +256,45 @@ def test_2d_refusals_are_collective(): assert all(m is not None for m in messages), ( f"some rank did NOT raise: {[m is None for m in messages]}") assert len(set(messages)) == 1, "ranks raised different errors" + + +def test_the_gather_moves_only_the_shell_around_the_zone(): + """The gather exists for the seam rule: the cells the carve drops, + their vertex star, and one more layer so the ring's points are + unshared. That is a shell about three cells thick around the zone, + and nothing else may move (#670: a distance blanket two cells wide + ahead of that growth moved 91% of a fixture whose cavity was 6%, + onto one rank, for good). Measured on this box, np=2 and np=4 alike: + the shell is 5046 cells against 7269 within three median cell + diameters, and the old mask moved 9831 (this test fails on it). A + box eight times the unit fixture, so a shell fits inside it.""" + from underworld3.utilities.edge_split import cell_diameters + + comm = uw.mpi.comm + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-0.5, -0.5, -0.5), maxCoords=(1.5, 1.5, 1.5), + cellSize=0.24, refinement=1, regular=False, qdegree=2) + width = 0.045 + cells = np.asarray(mesh._cell_node_indices(1, True)) + X = np.asarray(mesh.X.coords)[:, :3] + centroid = X[cells].mean(axis=1) + diameters = np.concatenate(comm.allgather( + np.asarray(cell_diameters(mesh.dm), dtype=float))) + h_med = float(np.median(diameters)) + + def slab_distance(P): + lo, hi = P.min(axis=0) - 0.5 * width, P.max(axis=0) + 0.5 * width + return np.linalg.norm(np.maximum(np.maximum(lo - centroid, 0.0), + centroid - hi), axis=1) + + d = np.min([slab_distance(P) for P in CROSS], axis=0) + within_three = int(comm.allreduce(int((d < 3.0 * h_med).sum()), + op=MPI.SUM)) + + new, info = place_thin_volume(mesh.dm, CROSS, width=width, + label="Zone", label_value=5) + assert info["n_zone_cells"] > 0 + assert all(g == info for g in comm.allgather(info)) + assert info["n_gathered"] <= within_three, ( + f"the gather moved {info['n_gathered']} cells; the shell rule " + f"allows at most the {within_three} within three cells of the zone") From 367c3376992db5a62f8d7e2f64144f4a7101c3ea Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 13:12:49 -0700 Subject: [PATCH 02/19] The placement gather takes several regions, each to its own rank (#670) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _gather_regions generalises the one-region gather: a chart of region ids marks the vertices, each region's star and layer is claimed, regions that touch (a shared cell or vertex, on any rank) are merged by a collective union-find, and each merged region goes to the rank that already holds most of it — or stays where it is when its star and layer are already interior to one rank. One shell partition moves them all. The count returned is the regions' size (the seam rule's footprint, a function of the mesh alone) alongside the cells that actually changed rank. _gather_region is now the one-region form of it, unchanged for its callers. ptest_0857 covers the three behaviours at np=2, 3 and 4: two regions a domain apart keep their own owners, two that touch are merged, and one vertex marked deep inside every rank moves nothing. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/place_surface.py | 157 +++++++++++++----- .../ptest_0857_gather_regions_parallel.py | 105 ++++++++++++ 2 files changed, 222 insertions(+), 40 deletions(-) create mode 100644 tests/parallel/ptest_0857_gather_regions_parallel.py diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 95e0b872f..95c24fc55 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -1602,76 +1602,153 @@ def _true_wall_vertex_mask(dm, n_vertices): def _gather_region(dm, vertex_mark_chart, verbose=False): """Redistribute so every marked vertex's cell star (+1 layer) is one rank's. - A mask-driven port of the contact stream's ``_redistribute_fault_interior`` - (feature/fault-split-node c8693579 / 1d487319; the label-driven original - is measured at np=2..8 with serial-identical topology). Only the marked - star moves — everything else keeps its load-balanced home — via a shell - partitioner. Returns ``(new_dm, n_moved)``: the global count of cells + The one-region form of :func:`_gather_regions` (a 0/1 mask is one + region). Returns ``(new_dm, n_moved)``: the global count of cells gathered, ``0`` when nothing moved (serial, or a region already interior to one rank), so callers can both branch on it and report it. The input is untouched. """ + ids = (np.asarray(vertex_mark_chart) != 0).astype(np.int32) + work, n_region, _n_moved, _owner, _canon = _gather_regions( + dm, ids, verbose=verbose) + return work, n_region + + +def _gather_regions(dm, vertex_region_chart, verbose=False): + """Redistribute so each marked REGION's cell star (+1 layer) is one rank's. + + A mask-driven port of the contact stream's ``_redistribute_fault_interior`` + (feature/fault-split-node c8693579 / 1d487319; the label-driven original + is measured at np=2..8 with serial-identical topology), generalised to + several regions (#670). ``vertex_region_chart`` is chart-length: ``0`` + for an unmarked point, ``k >= 1`` the region a vertex belongs to. Each + region's star and layer go to ONE rank, chosen per region — the rank + that already holds most of it — and a region that is already interior + to one rank is left where it is. Two regions whose stars or layers + touch (share a cell or a vertex, on any rank) are merged: their + cavities would share ring points, so they must be one rank's. Only the + marked stars move — everything else keeps its load-balanced home — via + one shell partition. + + Returns ``(new_dm, n_region, n_moved, owner, canon)``: the global size + of the regions (star and layer — the seam rule's footprint, a function + of the mesh alone), the count of cells that changed rank (``0`` when + none did; the input dm is then returned untouched), ``owner`` mapping + each merged region's id to its rank, and ``canon`` mapping every input + region id to its merged id. The decisions are COLLECTIVE: overlaps and + counts are gathered before any branch. + """ comm = dm.getComm().tompi4py() + ids_in = np.asarray(vertex_region_chart, dtype=np.int32) + n_ids = int(comm.allreduce(int(ids_in.max()) if ids_in.size else 0, + op=MPI.MAX)) + if n_ids == 0: + raise ValueError("place_sheet: the sheet meets no cell on any rank") if comm.size == 1: - return dm, 0 + canon = {k: k for k in range(1, n_ids + 1)} + return dm, 0, 0, {k: 0 for k in canon}, canon work = dm.clone() cS, cE = work.getHeightStratum(0) vS, vE = work.getDepthStratum(0) pStart, pEnd = work.getChart() + pairs = set() # (a, b): regions that touch, to be merged - mark = vertex_mark_chart.astype(np.int32).copy() - mark = _propagate_vertex(work, mark, MPI.MAX, np.maximum) + def reconcile(chart): + # every rank agrees on every point it can see; where a rank's own + # id loses to a neighbour's under MAX the two regions touch + before = chart.copy() + after = _propagate_vertex(work, chart, MPI.MAX, np.maximum) + touch = (before > 0) & (after != before) + for a, b in zip(before[touch], after[touch]): + pairs.add((int(a), int(b))) + return after - def star_of_marked(m): - out = set() - for v in range(vS, vE): - if m[v - pStart]: - for q in work.getTransitiveClosure(v, useCone=False)[0]: - if cS <= int(q) < cE: - out.add(int(q)) - return out + cell_region = np.zeros(cE - cS, dtype=np.int32) - star = star_of_marked(mark) + def claim_cells(chart): + for v in range(vS, vE): + k = int(chart[v - pStart]) + if k == 0: + continue + for q in work.getTransitiveClosure(v, useCone=False)[0]: + if cS <= int(q) < cE: + c = int(q) - cS + if cell_region[c] and cell_region[c] != k: + pairs.add((int(cell_region[c]), k)) + cell_region[c] = max(cell_region[c], k) + + claim_cells(reconcile(ids_in.copy())) # One growth layer: the surgery needs every point in the closure of a # region cell unshared, and a point is unshared exactly when all its # incident cells are co-resident. - mark2 = np.zeros(pEnd - pStart, dtype=np.int32) - for c in star: - for q in work.getTransitiveClosure(c)[0]: + layer = np.zeros(pEnd - pStart, dtype=np.int32) + for c in np.flatnonzero(cell_region): + k = int(cell_region[c]) + for q in work.getTransitiveClosure(int(c) + cS)[0]: if vS <= int(q) < vE: - mark2[int(q) - pStart] = 1 - mark2 = _propagate_vertex(work, mark2, MPI.MAX, np.maximum) - star |= star_of_marked(mark2) + i = int(q) - pStart + if layer[i] and layer[i] != k: + pairs.add((int(layer[i]), k)) + layer[i] = max(layer[i], k) + claim_cells(reconcile(layer)) + + # merge touching regions: union-find on the gathered pair set, the same + # on every rank + parent = list(range(n_ids + 1)) + + def find(a): + while parent[a] != a: + parent[a] = parent[parent[a]] + a = parent[a] + return a - counts = np.asarray(comm.allgather(len(star))) + for a, b in sorted(set().union(*comm.allgather(pairs))): + ra, rb = find(a), find(b) + if ra != rb: + parent[max(ra, rb)] = min(ra, rb) + canon = {k: find(k) for k in range(1, n_ids + 1)} + for c in np.flatnonzero(cell_region): + cell_region[c] = canon[int(cell_region[c])] + + regions = sorted(set(canon.values())) + local = np.array([int((cell_region == k).sum()) for k in regions], + dtype=np.int64) + counts = np.asarray(comm.allgather(local)) # (size, n_regions) if counts.sum() == 0: raise ValueError("place_sheet: the sheet meets no cell on any rank") - if np.count_nonzero(counts) == 1: - # The region, star and layer included, is already interior to one - # rank: the seam rule holds where the mesh stands, and the surgery - # runs there with no cell moved (#670). COLLECTIVE decision — the - # counts are gathered, so every rank takes this branch together. + owner = {} + assign = np.full(cE - cS, comm.rank, dtype=np.int32) + n_moved = 0 + for j, k in enumerate(regions): + col = counts[:, j] + owner[k] = int(np.argmax(col)) + if np.count_nonzero(col) <= 1: + # star and layer already interior to one rank: the seam rule + # holds where the mesh stands, nothing moves (#670) + continue + n_moved += int(col.sum() - col[owner[k]]) + assign[cell_region == k] = owner[k] + n_region = int(counts.sum()) + if n_moved == 0: if verbose: - uw.pprint(f"[place_sheet] a {int(counts.sum())}-cell region is " - f"interior to rank {int(np.argmax(counts))}; no gather") - return dm, 0 - target = int(np.argmax(counts)) + uw.pprint(f"[place_sheet] {len(regions)} region(s), " + f"{n_region} cells, already interior to their ranks; " + f"no gather") + return dm, n_region, 0, owner, canon - assign = np.full(cE - cS, comm.rank, dtype=np.int32) - for c in star: - assign[c - cS] = target order = np.argsort(assign, kind="stable").astype(np.int32) sizes = np.bincount(assign, minlength=comm.size).astype(np.int32) - part = work.getPartitioner() part.setType(PETSc.Partitioner.Type.SHELL) part.setShellPartition(comm.size, sizes=sizes, points=order) work.distribute() if verbose: - uw.pprint(f"[place_sheet] gathered a {int(counts.sum())}-cell region " - f"onto rank {target}") - return work, int(counts.sum()) + uw.pprint(f"[place_sheet] gathered {n_moved} cells: " + + ", ".join(f"region {k} ({int(counts[:, j].sum())} " + f"cells) onto rank {owner[k]}" + for j, k in enumerate(regions))) + return work, n_region, n_moved, owner, canon def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, diff --git a/tests/parallel/ptest_0857_gather_regions_parallel.py b/tests/parallel/ptest_0857_gather_regions_parallel.py new file mode 100644 index 000000000..e8674bfd7 --- /dev/null +++ b/tests/parallel/ptest_0857_gather_regions_parallel.py @@ -0,0 +1,105 @@ +"""The region gather behind placement (:func:`place_surface._gather_regions`). + +Each marked region's star and layer go to one rank of their own; regions +that touch are merged; a region already interior to a rank is left where +it is (#670). Deterministic assertions: the region ids and owners are +gathered and identical on every rank. Run: + + mpirun -np 3 python -m pytest tests/parallel/ptest_0857_gather_regions_parallel.py --with-mpi +""" +import numpy as np +import pytest +from mpi4py import MPI + +import underworld3 as uw +from underworld3.utilities.line_cut import _coords +from underworld3.utilities.place_surface import (_gather_regions, + _shared_point_flags) + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2, + pytest.mark.tier_b, pytest.mark.timeout(300)] + + +def _box(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.12, regular=False, qdegree=1) + + +def _vertex_ids(mesh, rule): + """Chart-length region ids from a rule on vertex coordinates.""" + dm = mesh.dm + vS, vE = dm.getDepthStratum(0) + pStart, pEnd = dm.getChart() + X = _coords(dm)[: vE - vS] # DM vertex order, as the placer reads it + ids = np.zeros(pEnd - pStart, dtype=np.int32) + ids[vS - pStart: vE - pStart] = rule(X) + return ids + + +def _same_everywhere(comm, value): + return all(v == value for v in comm.allgather(value)) + + +def test_two_far_regions_keep_their_own_owners(): + comm = uw.mpi.comm + mesh = _box() + ids = _vertex_ids(mesh, lambda X: np.where( + X[:, 0] < 0.15, 1, np.where(X[:, 0] > 0.85, 2, 0))) + new, n_region, n_moved, owner, canon = _gather_regions(mesh.dm, ids) + assert _same_everywhere(comm, (n_region, n_moved, owner, canon)) + assert canon == {1: 1, 2: 2}, "two regions a domain apart were merged" + assert set(owner) == {1, 2} + assert n_region > 0 and n_moved <= n_region + # the moved cells are exactly those a region claimed away from its rank + cS, cE = new.getHeightStratum(0) + n_cells = int(comm.allreduce(cE - cS, op=MPI.SUM)) + cS0, cE0 = mesh.dm.getHeightStratum(0) + assert n_cells == int(comm.allreduce(cE0 - cS0, op=MPI.SUM)) + + +def test_touching_regions_are_merged(): + comm = uw.mpi.comm + mesh = _box() + ids = _vertex_ids(mesh, lambda X: np.where( + X[:, 0] < 0.45, 1, np.where(X[:, 0] < 0.55, 2, 0))) + _new, n_region, _n_moved, owner, canon = _gather_regions(mesh.dm, ids) + assert _same_everywhere(comm, (owner, canon)) + assert canon == {1: 1, 2: 1}, canon + assert set(owner) == {1} + + +def test_an_interior_region_is_not_moved(): + """Mark one vertex per rank, deep inside that rank's own cells: every + region's star and layer are interior, so nothing moves at all. An + elongated box, so each rank's piece is many cells deep and the vertex + farthest from any shared point has a two-cell neighbourhood of its + own.""" + comm = uw.mpi.comm + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(4.0, 1.0, 1.0), + cellSize=0.15, regular=False, qdegree=1) + dm = mesh.dm + vS, vE = dm.getDepthStratum(0) + pStart, pEnd = dm.getChart() + # shared = roots AND leaves: a seam vertex this rank owns is a root, + # absent from its own leaf list, and its star is mostly elsewhere + shared = np.asarray(_shared_point_flags(dm)).astype(bool) + shared_v = np.flatnonzero(shared[vS - pStart: vE - pStart]) + X = _coords(dm)[: vE - vS] # DM vertex order, as the placer reads it + ids = np.zeros(pEnd - pStart, dtype=np.int32) + if len(X): + if shared_v.size: + d = np.min(np.linalg.norm( + X[:, None, :] - X[shared_v][None, :, :], axis=2), axis=1) + else: + d = np.linalg.norm(X - X.mean(axis=0), axis=1) + v = int(np.argmax(d)) + ids[v + vS - pStart] = comm.rank + 1 + _new, n_region, n_moved, owner, canon = _gather_regions(dm, ids) + assert _same_everywhere(comm, (n_region, n_moved, owner, canon)) + assert n_region > 0 + assert n_moved == 0, (n_moved, owner, canon) + for k, r in owner.items(): + assert r == k - 1, (owner, "a region's owner is not the rank " + "that marked it") From 7b917e3d718a9f3d6d2463e16faa91cfb2cd957d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 13:44:16 -0700 Subject: [PATCH 03/19] Placement and split gather per region: each zone to its own rank, or not at all (#670) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The thin volume marks one region per connected component of the assembly (zones fused through shared faces are one component; zones a domain apart are separate) and gathers them with _gather_regions: regions whose shells touch are merged, each goes to the rank already holding most of it, and one whose shell is interior to a rank is not moved. The owning ranks then carve and fill their own components concurrently — each on its compacted share of the assembly and skin — and the collective rebuild sews them at once. The outcrop and ladder paths keep one region; their bowl, cap and extrusion are single-rank. The split follows the same regions. _redistribute_fault_interior takes groups of faults and delegates to _gather_regions; split_faults passes them through, and the network derives them from the placement's report of which region each embedded mid-surface's zone became (info["embedded_regions"]), so what the placement kept apart is not gathered together afterwards. The network's info now carries n_regions, n_gathered and n_moved. Measured. Two zones a domain apart on the 8x box, np=2: 304 cells moved to two different owners where one region for the pair moved 8451 to one rank; the sewn mesh has the same zone, skin and removed counts. Two faults in one network, split realisation, np=2: both regions already interior to their ranks, nothing moved, 14,339 and 14,332 cells per rank; np=4: each region to a different rank, 2507 moved. Two faults 0.8 apart on 0.126 cells merge into one region (their shells touch); 1.4 apart they do not. Tests: ptest_0855 gains the two-zone case (np=2, 4); ptest_0863 the two-fault network (np=2, 4). The placement, split and network parallel suites and the serial thin-volume and network suites pass. Not fixed here: the contact solve on the balanced two-fault layout fails in the co-located multigrid tail's setup (#671) even at np=2, where nothing moved — the band is then genuinely distributed, and the tail was built for a band on one rank. The placement and the split complete; the tail is the next item. Underworld development team with AI support from Claude Code --- .../conforming-surfaces-and-fault-zones.md | 26 ++- src/underworld3/meshing/fault_network.py | 19 +- src/underworld3/utilities/fault_split.py | 125 +++++-------- src/underworld3/utilities/place_surface.py | 164 ++++++++++++++---- .../ptest_0855_place_thin_volume_parallel.py | 27 +++ ...st_0863_fault_network_3d_width_parallel.py | 39 +++++ 6 files changed, 280 insertions(+), 120 deletions(-) diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index c0b50096c..27e4acc90 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -310,11 +310,27 @@ The gather is one-way: the moved cells stay on the surgery rank, together with the cells the fill creates, so that rank carries the shell plus the band as extra load. That is the accepted trade (extra load on one rank in exchange for no communication during the solve), and it is proportional to -the surface, not the domain. Two limits remain, recorded in #670: the whole -network is one region with one target rank, so two surfaces each interior to -a different rank are still gathered together; and a small domain cannot hold -a shell at all (three base cells reaching the walls is the whole box, which -is what the crossing-patches test fixture does). +the surface, not the domain. + +A network is gathered per region, not as a whole. The assembly's connected +components (zones fused through shared faces are one component; zones a +domain apart are separate ones) are marked as separate regions, regions +whose shells touch are merged, and each region goes to the rank that +already holds most of it, or stays where it is when its shell is already +interior to one. The surgeries then run concurrently, each owning rank +carving and filling its own components; the collective rebuild sews them +all at once. Two patches a domain apart at np=2 move 304 cells this way +where one region for the pair moved 8451, and `info["n_regions"]` and +`info["n_moved"]` report it. The split's own redistribution follows the +same regions (`split_faults(..., groups=...)`, which the network passes +from `info["embedded_regions"]`), so what the placement kept apart is not +gathered together afterwards. The outcrop and ladder paths keep one +region: their bowl, cap and extrusion machinery is single-rank. + +One limit remains: a small domain cannot hold a shell at all (three base +cells reaching the walls is the whole box, which is what the +crossing-patches test fixture does), so that fixture measures correctness +only, never balance. ## The thin volume: finite-width zones, junctions in the volume diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py index df9178066..55196e0f5 100644 --- a/src/underworld3/meshing/fault_network.py +++ b/src/underworld3/meshing/fault_network.py @@ -559,7 +559,17 @@ def _build_3d_band(self, h_far=None, qdegree=2, realisation="split", boundaries=Enum("boundaries", members), verbose=False) band = mesh.cells_labelled("Band", 71) if realisation == "split": - mesh = split_faults(mesh, [n for n, _P in self.prepared]) + # the pieces the placement kept apart stay apart in the split's + # own redistribution: one group per placed region (#670) + groups = None + regions = info.get("embedded_regions") + if regions is not None: + by_region = {} + for (name, _P), r in zip(self.prepared, regions): + by_region.setdefault(r, []).append(name) + groups = list(by_region.values()) + mesh = split_faults(mesh, [n for n, _P in self.prepared], + groups=groups) # reduce first, then branch: the defect must raise on every # rank together or not at all. A healthy pairing is a # bijection between disjoint sides: a node paired with @@ -616,7 +626,12 @@ def _build_3d_band(self, h_far=None, qdegree=2, realisation="split", "spacing": [h] * len(self.prepared), "width": float(self.width), "mesher": "network", "margin_rings": [(margin_rings, margin_rings)] - * len(self.prepared)} + * len(self.prepared), + # the placement's parallel record (#670): regions + # gathered, their size, and the cells that moved + "n_regions": int(info.get("n_regions", 1)), + "n_gathered": int(info.get("n_gathered", 0)), + "n_moved": int(info.get("n_moved", 0))} self.mesh = mesh self._make_surfaces() return self.mesh diff --git a/src/underworld3/utilities/fault_split.py b/src/underworld3/utilities/fault_split.py index 89d49646c..896b6f14a 100644 --- a/src/underworld3/utilities/fault_split.py +++ b/src/underworld3/utilities/fault_split.py @@ -1212,7 +1212,7 @@ def _fault_labels_touch_seam(dm, labels): return bool(uw.mpi.comm.allreduce(touch, op=MPI.LOR)) -def _redistribute_fault_interior(dm, labels, verbose=False): +def _redistribute_fault_interior(dm, labels, verbose=False, groups=None): """Redistribute a cut mesh so each fault's cell star is rank-interior. The default partition's balance cuts are ATTRACTED to the locally @@ -1231,98 +1231,54 @@ def _redistribute_fault_interior(dm, labels, verbose=False): keyed on the union of its faults, which is what lets every split that follows run without migrating any prior pairing. Returns a NEW dm; the input is untouched. + + ``groups`` partitions ``labels`` into the faults that must share a + rank (a junction-connected cluster); each group is one region of + :func:`place_surface._gather_regions`, moved to its own rank — or not + moved at all when its star is already interior to one — so two + faults a domain apart are never gathered together (#670). Without + ``groups`` the whole network is one region, as before. """ + from underworld3.utilities.place_surface import _gather_regions + comm = dm.getComm().tompi4py() if comm.size == 1: return dm + if groups is None: + groups = [list(labels)] - work = dm.clone() - cS, cE = work.getHeightStratum(0) - fS, fE = work.getHeightStratum(1) - vS, vE = work.getDepthStratum(0) - _pS, pEnd = work.getChart() - sf = work.getPointSF() - - def propagate(mark): - # global OR of a chart-length mark across the point SF: remote - # copies fold into the owner, the owner's verdict returns to - # every copy — after this, every rank agrees on every point it - # can see. (A cell can touch a patch vertex without holding any - # labelled face in its own closure, so purely local label - # reading under-marks near the current seam.) - tmp = mark.copy() - sf.reduceBegin(MPI.INT32_T, tmp, mark, MPI.MAX) - sf.reduceEnd(MPI.INT32_T, tmp, mark, MPI.MAX) - out = mark.copy() - sf.bcastBegin(MPI.INT32_T, mark, out, MPI.REPLACE) - sf.bcastEnd(MPI.INT32_T, mark, out, MPI.REPLACE) - return np.maximum(mark, out) - - def star_of_marked(mark): - cells = set() - for v in range(vS, vE): - if mark[v]: - for q in work.getTransitiveClosure(v, useCone=False)[0]: - if cS <= int(q) < cE: - cells.add(int(q)) - return cells - - # fault vertices (of every fault in the batch), marked globally - mark = np.zeros(pEnd, dtype=np.int32) - for name, value in labels: - if work.hasLabel(name) and \ - work.getLabel(name).getStratumSize(int(value)) > 0: - for f in work.getLabel(name).getStratumIS( - int(value)).getIndices(): + fS, fE = dm.getHeightStratum(1) + vS, vE = dm.getDepthStratum(0) + pStart, pEnd = dm.getChart() + ids = np.zeros(pEnd - pStart, dtype=np.int32) + for g, group in enumerate(groups, start=1): + for name, value in group: + if not (dm.hasLabel(name) + and dm.getLabel(name).getStratumSize(int(value)) > 0): + continue + for f in dm.getLabel(name).getStratumIS(int(value)).getIndices(): if fS <= int(f) < fE: - for q in work.getTransitiveClosure(int(f))[0]: + for q in dm.getTransitiveClosure(int(f))[0]: if vS <= int(q) < vE: - mark[int(q)] = 1 - mark = propagate(mark) - star = star_of_marked(mark) - - # one growth layer: the split's seam rule needs every point in the - # CLOSURE of a patch-vertex-star cell unshared, and a point is - # unshared exactly when all its incident cells are co-resident — - # so gather every cell touching any vertex of the star cells' - # closures. One layer is exactly sufficient (each such point's - # incident cells all touch a marked-closure vertex). - mark2 = np.zeros(pEnd, dtype=np.int32) - for c in star: - for q in work.getTransitiveClosure(c)[0]: - if vS <= int(q) < vE: - mark2[int(q)] = 1 - mark2 = propagate(mark2) - star |= star_of_marked(mark2) - - counts = np.asarray(comm.allgather(len(star))) - if counts.sum() == 0: + i = int(q) - pStart + ids[i] = max(ids[i], g) + if comm.allreduce(int(ids.max()) if ids.size else 0, op=MPI.MAX) == 0: names = sorted(name for name, _value in labels) raise RuntimeError( f"fault_split: no facets labelled {names} found on any rank.") - target = int(np.argmax(counts)) - - n_local = cE - cS - assign = np.full(n_local, comm.rank, dtype=np.int32) - for c in star: - assign[c - cS] = target - order = np.argsort(assign, kind="stable").astype(np.int32) - sizes = np.bincount(assign, minlength=comm.size).astype(np.int32) - - part = work.getPartitioner() - part.setType(PETSc.Partitioner.Type.SHELL) - part.setShellPartition(comm.size, sizes=sizes, points=order) - work.distribute() + work, n_region, n_moved, owner, _canon = _gather_regions(dm, ids) + if work is dm: + work = dm.clone() # the contract: a new dm, untouched input if verbose: names = sorted(name for name, _value in labels) - uw.pprint(f"[fault_split] {names}: star of " - f"{int(counts.sum())} cells gathered onto rank " - f"{target}; local cells now " + uw.pprint(f"[fault_split] {names}: {len(owner)} region(s) of " + f"{n_region} cells, {n_moved} moved (owners " + f"{sorted(owner.values())}); local cells now " f"{work.getHeightStratum(0)[1]}") return work -def split_faults(mesh, names, verbose=False): +def split_faults(mesh, names, verbose=False, groups=None): """Split a NETWORK of already-labelled faults, any dimension, any np. The parallel obstruction to sequential :func:`split_fault` calls is @@ -1334,6 +1290,11 @@ def split_faults(mesh, names, verbose=False): pre-pass the 2-D ``add_fault`` performs for freshly cut chains, made available for faults that already carry labels (the 3-D embedded patches, a reloaded mesh). + + ``groups``, a list of lists of names, says which faults must share a + rank (a junction-connected cluster); each group is redistributed to + its own rank, and one whose star is already interior to a rank is + not moved (#670). Faults not named in any group form one more. """ import underworld3 as uw from underworld3.discretisation import Mesh @@ -1341,9 +1302,17 @@ def split_faults(mesh, names, verbose=False): out = mesh if uw.mpi.size > 1: labels = [(n, int(mesh.boundaries[n].value)) for n in names] + label_groups = None + if groups is not None: + value = dict(labels) + label_groups = [[(n, value[n]) for n in g] for g in groups] + rest = [n for n in names if not any(n in g for g in groups)] + if rest: + label_groups.append([(n, value[n]) for n in rest]) if _fault_labels_touch_seam(mesh.dm, labels): dm = _redistribute_fault_interior(mesh.dm, labels, - verbose=verbose) + verbose=verbose, + groups=label_groups) out = Mesh(dm, simplex=mesh.dm.isSimplex(), coordinate_system_type=( mesh.CoordinateSystem.coordinate_type), diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 95c24fc55..45d97a79d 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -1599,6 +1599,39 @@ def _true_wall_vertex_mask(dm, n_vertices): return mark[vS - pStart: vS - pStart + n_vertices] == 1 +def _assembly_components(cells): + """Connected components of a standalone assembly mesh, by shared facets. + + Returns a 1-based component id per cell, the same on every rank (the + assembly is broadcast, so this is a pure function of it). Two zones of + a network that are fused touch through shared faces and are one + component; zones a domain apart are separate ones, and #670 places + each on a rank of its own. + """ + from itertools import combinations + n = len(cells) + parent = list(range(n)) + + def find(a): + while parent[a] != a: + parent[a] = parent[parent[a]] + a = parent[a] + return a + + nv = cells.shape[1] + owner_of_face = {} + for c, cell in enumerate(cells): + for f in combinations(sorted(int(v) for v in cell), nv - 1): + other = owner_of_face.setdefault(f, c) + if other != c: + ra, rb = find(other), find(c) + if ra != rb: + parent[max(ra, rb)] = min(ra, rb) + roots = np.array([find(c) for c in range(n)], dtype=np.int64) + _uniq, comp = np.unique(roots, return_inverse=True) + return (comp + 1).astype(np.int32) + + def _gather_region(dm, vertex_mark_chart, verbose=False): """Redistribute so every marked vertex's cell star (+1 layer) is one rank's. @@ -7049,15 +7082,34 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, reach_v = np.maximum(clearance * h_vertex, 0.6 * width) d_skin = _sheet_distance_within(X, skin_xyz, skin_tris, reach_v + 1.0 * h_vertex) + # One region per connected component of the assembly (#670): each + # zone's shell goes to its own rank, and a zone whose shell is already + # interior to a rank moves nothing. The outcrop and ladder paths carry + # single-rank machinery (the bowl, the cap, the extrusion) and keep one + # region. + comp_of_tet = _assembly_components(asm_tets) + n_comp = int(comp_of_tet.max()) + single = outcropping or mesher == "ladder" or n_comp == 1 mark = np.zeros(pEnd - pStart, dtype=np.int32) - mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) - + vS - pStart] = 1 + if single: + mark[np.flatnonzero(d_skin < reach_v + 1.0 * h_vertex) + + vS - pStart] = 1 + else: + for k in range(1, n_comp + 1): + xyz_k, tris_k, _ids_k = _assembly_skin( + asm_pts, asm_tets[comp_of_tet == k]) + d_k = _sheet_distance_within(X, xyz_k, tris_k, + reach_v + 1.0 * h_vertex) + near = np.flatnonzero(d_k < reach_v + 1.0 * h_vertex) + mark[near + vS - pStart] = np.maximum( + mark[near + vS - pStart], k) volume_before = np.array([_owned_cell_volume(dm)], dtype=float) comm.Allreduce(MPI.IN_PLACE, volume_before, op=MPI.SUM) - dm_work, moved = _gather_region(dm, mark, verbose=verbose) - if moved: + dm_work, moved, n_moved, owner, canon = _gather_regions( + dm, mark, verbose=verbose) + if n_moved: vS, vE = dm_work.getDepthStratum(0) pStart, pEnd = dm_work.getChart() X = _coords(dm_work)[: vE - vS] @@ -7086,16 +7138,40 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, reach_c = (np.maximum(clearance * h_cell_local, 0.6 * width) if len(cells) else np.zeros(0)) - n_region = int((d_skin < reach_v).sum()) - owners = np.asarray(comm.allgather(n_region)) - if owners.sum() == 0: - raise ValueError("the thin volume meets no cell of this mesh") - target = int(np.argmax(owners)) + mine_regions = {r for r, rk in owner.items() if rk == comm.rank} + my_comps = [k for k in range(1, n_comp + 1) + if canon[k if not single else 1] in mine_regions] + mine = bool(my_comps) + target = owner[canon[1]] # the root of the outcrop broadcasts + if mine and not single: + # this rank's share of the assembly, compacted: its components' + # cells, their nodes, and the skin of that subset (components are + # face-disjoint, so the subset's skin is its components' skins) + sel = np.isin(comp_of_tet, my_comps) + used = np.unique(asm_tets[sel]) + remap = np.full(len(asm_pts), -1, dtype=np.int64) + remap[used] = np.arange(len(used)) + asm_pts_m = asm_pts[used] + asm_tets_m = remap[asm_tets[sel]] + skin_xyz_m, skin_tris_m, skin_node_ids_m = _assembly_skin( + asm_pts_m, asm_tets_m) + skin_tris_fill_m = skin_tris_m + elif mine: + asm_pts_m, asm_tets_m = asm_pts, asm_tets + skin_xyz_m, skin_tris_m, skin_node_ids_m = ( + skin_xyz, skin_tris, skin_node_ids) + skin_tris_fill_m = skin_tris_fill + else: + asm_pts_m = np.empty((0, 3), dtype=float) + asm_tets_m = np.empty((0, 4), dtype=np.int64) + skin_xyz_m = np.empty((0, 3), dtype=float) + skin_tris_m = skin_tris_fill_m = np.empty((0, 3), dtype=np.int64) + skin_node_ids_m = np.empty(0, dtype=np.int64) failure = None victims = drop_ids = None fill = shell_vert_ids = None - if comm.rank == target: + if mine: try: deletable = near = None if outcropping: @@ -7103,9 +7179,9 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, dom_tris) deletable, near, _regions = _outcrop_frame_3d( X, on_wall, dom_verts, dom_tris, dom_region, - skin_xyz, skin_tris[band_idx]) + skin_xyz_m, skin_tris_m[band_idx]) victims, drop_ids, shell, cap_faces = _carve_around_volume_3d( - dm_work, X, cells, skin_xyz, skin_tris, reach_v, reach_c, + dm_work, X, cells, skin_xyz_m, skin_tris_m, reach_v, reach_c, held_cells, on_wall, shared, open_deletable=deletable, open_near=near) touched = set() @@ -7157,14 +7233,14 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, "boundary complex; the wall mask and the complex " "disagree") _d_band, at_band = _nearest_facet( - skin_xyz[skin_tris[band_idx]].mean(axis=1), + skin_xyz_m[skin_tris_m[band_idx]].mean(axis=1), dom_verts, dom_tris) alive = np.ones(len(X), dtype=bool) alive[np.asarray(victims, dtype=np.int64)] = False cap_nodes, hole_nodes, cap_extra, cap_tris = \ _outcrop_collar_3d( X, alive, cap_tris_mesh, dom_region[at_cap], - dom_planes, skin_xyz, skin_tris[band_idx], + dom_planes, skin_xyz_m, skin_tris_m[band_idx], dom_region[at_band], band_outline) cap_payload = { "rim_shell_local": [local[v] for v in cap_nodes], @@ -7173,8 +7249,8 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, "extra_xyz": cap_extra, } - fill = _gmsh_fill_annulus_3d(shell_xyz, shell_tris, skin_xyz, - skin_tris_fill, size_out=h, + fill = _gmsh_fill_annulus_3d(shell_xyz, shell_tris, skin_xyz_m, + skin_tris_fill_m, size_out=h, size_in=size, cap=cap_payload) (_pts, _tets, moved_nodes, skin_out, _n_shell, cap_out) = fill @@ -7182,10 +7258,10 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, raise RuntimeError( f"the gap fill moved {moved_nodes} constrained node(s); " "the cavity cannot be sewn back.") - if skin_out != len(skin_tris_fill): + if skin_out != len(skin_tris_fill_m): raise RuntimeError( f"the gap fill remeshed the skin ({skin_out} triangles " - f"for {len(skin_tris_fill)} given).") + f"for {len(skin_tris_fill_m)} given).") if cap_payload is not None and ( cap_out is None or len(cap_out) != len(cap_payload["tris"])): @@ -7207,10 +7283,10 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, # Placed rows: the assembly's nodes first, the gap fill's new points # after. A gap-fill node is a shell node (an OLD vertex), a skin node # (assembly row) or new; assembly tets reference assembly rows only. - if comm.rank == target: + if mine: fill_pts, fill_tets, _m, _s, n_shell, cap_out = fill - n_skin = len(skin_xyz) - skin_row = np.asarray(skin_node_ids, dtype=np.int64) + n_skin = len(skin_xyz_m) + skin_row = np.asarray(skin_node_ids_m, dtype=np.int64) gap_new = fill_pts[n_shell + n_skin:] def gap_code(v): @@ -7218,13 +7294,13 @@ def gap_code(v): return int(shell_vert_ids[v]) if v < n_shell + n_skin: return -(int(skin_row[v - n_shell]) + 1) - return -(len(asm_pts) + (int(v) - n_shell - n_skin) + 1) + return -(len(asm_pts_m) + (int(v) - n_shell - n_skin) + 1) made = np.array( [[gap_code(int(v)) for v in tet] for tet in fill_tets] - + [[-(int(v) + 1) for v in tet] for tet in asm_tets], + + [[-(int(v) + 1) for v in tet] for tet in asm_tets_m], dtype=np.int64) - placed = np.vstack([asm_pts, gap_new]) + placed = np.vstack([asm_pts_m, gap_new]) victims_arr = np.asarray(victims, dtype=np.int64) drop_arr = np.asarray(drop_ids, dtype=np.int64) else: @@ -7247,10 +7323,10 @@ def gap_code(v): new.createLabel(name) n_cells_local = 0 n_skin_local = 0 - if comm.rank == target: + if mine: out_label = new.getLabel(label) out_skin = new.getLabel(skin_label) - for tet in asm_tets: + for tet in asm_tets_m: joined = new.getFullJoin([int(placed_new[int(v)]) for v in tet]) if len(joined) != 1: failure = ("an assembly cell is not a cell of the sewn mesh; " @@ -7259,7 +7335,7 @@ def gap_code(v): out_label.setValue(int(joined[0]), int(label_value)) n_cells_local += 1 else: - for tri in skin_tris: + for tri in skin_tris_m: joined = new.getFullJoin( [int(placed_new[int(skin_row[int(v)])]) for v in tri]) if len(joined) != 1: @@ -7281,17 +7357,17 @@ def gap_code(v): # restores each wall's own labels. pairs = comm.bcast( sorted({p for _tri, pairs_f in removed_wall for p in pairs_f}) - if comm.rank == target else None, root=target) + if mine else None, root=target) n_wall_expect = comm.bcast( ((len(cap_out) if cap_out is not None else 0) + len(band_idx)) - if comm.rank == target else 0, root=target) + if mine else 0, root=target) for name, val in (pairs or []): if not new.hasLabel(name): new.createLabel(name) n_wall_local = 0 - if comm.rank == target and outcropping: + if mine and outcropping: n_shell_ids = np.asarray(shell_vert_ids, dtype=np.int64) - n_skin = len(skin_xyz) + n_skin = len(skin_xyz_m) def new_id(v): if v < n_shell: @@ -7300,18 +7376,18 @@ def new_id(v): return int(point_map[old_pt - pStart]) if v < n_shell + n_skin: return int(placed_new[int(skin_row[v - n_shell])]) - return int(placed_new[len(asm_pts) + (v - n_shell - n_skin)]) + return int(placed_new[len(asm_pts_m) + (v - n_shell - n_skin)]) wall_tris = ([[new_id(int(v)) for v in t] for t in cap_out] if cap_out is not None else []) n_cap_tris = len(wall_tris) wall_tris += [[int(placed_new[int(skin_row[int(v)])]) for v in t] - for t in skin_tris[band_idx]] + for t in skin_tris_m[band_idx]] centres = np.array( [fill_pts[np.asarray(t)].mean(axis=0) for t in cap_out] if cap_out is not None else np.zeros((0, 3))) if len(band_idx): - band_cen = skin_xyz[skin_tris[band_idx]].mean(axis=1) + band_cen = skin_xyz_m[skin_tris_m[band_idx]].mean(axis=1) centres = (np.vstack([centres, band_cen]) if len(centres) else band_cen) old_pts = np.vstack([tri for tri, _p in removed_wall]) @@ -7413,7 +7489,25 @@ def new_id(v): f"zone cells, {info['n_skin_faces']} skin faces; placed " f"{info['n_placed']} vertices, removed " f"{info['n_removed']}") - info["n_gathered"] = int(moved) # cells the gather moved (#670) + info["n_gathered"] = int(moved) # the regions' size, star and layer (#670) + info["n_moved"] = int(n_moved) # cells that changed rank for it + info["n_regions"] = len(set(canon.values())) + if embedded_nodes is not None: + # which region each embedded mid-surface's zone became: the split + # redistributes per region too (#670), so it must not gather what + # the placement kept apart + node_comp = np.zeros(len(asm_pts), dtype=np.int32) + for k in range(1, n_comp + 1): + node_comp[np.unique(asm_tets[comp_of_tet == k])] = k + region_of_comp = {k: (canon[k] if not single else canon[1]) + for k in range(1, n_comp + 1)} + info["embedded_regions"] = [] + for pts in embedded_nodes: + hit = np.flatnonzero( + np.all(np.isclose(asm_pts[:, None, :], pts[None, :1, :], + atol=1e-12), axis=2)) + k = int(node_comp[hit[0]]) if hit.size else 1 + info["embedded_regions"].append(region_of_comp.get(k, 1)) return new, info diff --git a/tests/parallel/ptest_0855_place_thin_volume_parallel.py b/tests/parallel/ptest_0855_place_thin_volume_parallel.py index 64b7f29de..64da736be 100644 --- a/tests/parallel/ptest_0855_place_thin_volume_parallel.py +++ b/tests/parallel/ptest_0855_place_thin_volume_parallel.py @@ -298,3 +298,30 @@ def slab_distance(P): assert info["n_gathered"] <= within_three, ( f"the gather moved {info['n_gathered']} cells; the shell rule " f"allows at most the {within_three} within three cells of the zone") + + +def test_two_zones_apart_are_two_regions(): + """Two patches a domain apart are two connected components of the + assembly, so two regions of the gather, each to its own rank + (#670): the surgeries run concurrently, and the sewn mesh is the + same one the single-region gather produced (zone, skin and removed + counts). Measured here at np=2: 304 cells moved where one region + for the pair moved 8451, and the two owners are different ranks.""" + comm = uw.mpi.comm + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-0.5, -0.5, -0.5), maxCoords=(1.5, 1.5, 1.5), + cellSize=0.24, refinement=1, regular=False, qdegree=2) + apart = [np.array([[-0.2, 0.0, 0.0], [0.2, 0.0, 0.0], + [0.2, 0.0, 0.4], [-0.2, 0.0, 0.4]]), + np.array([[0.8, 1.0, 0.6], [1.2, 1.0, 0.6], + [1.2, 1.0, 1.0], [0.8, 1.0, 1.0]])] + new, info = place_thin_volume(mesh.dm, apart, width=0.045, + label="Zone", label_value=5) + assert all(g == info for g in comm.allgather(info)) + assert info["n_regions"] == 2, info + assert info["n_zone_cells"] > 0 + assert _owned_label_count(new, "Zone", 5) == info["n_zone_cells"] + assert _owned_label_count(new, "Zone_skin", 5) == info["n_skin_faces"] + # the two regions do not share cells: what moved is at most one of + # them, never both (both to one rank would be the old behaviour) + assert info["n_moved"] < info["n_gathered"] diff --git a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py index 01da8c259..306158b41 100644 --- a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py +++ b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py @@ -113,3 +113,42 @@ def test_network_3d_width_weak_plane_solve_np2(): peak = comm.allreduce(float(slips.get(name, 0.0)), op=max) assert peak == pytest.approx(expected, rel=2e-2), ( f"{name}: parallel peak {peak:.4f} vs serial {expected}") + + +def test_two_faults_apart_are_placed_by_their_own_ranks(): + """Two faults a domain apart in one network are two regions of the + placement gather and two groups of the split's redistribution + (#670): each is placed and split by the rank that holds it, and + the network never gathers them together. On a base wide enough to + hold both shells (the width path's own box is the unit cube, so the + band builder is called with a wider one), measured at np=2: both + regions interior to their ranks, nothing moved, the mesh balanced + to within a few cells; at np=4 each region to a different rank.""" + A = np.array([[0.10, -0.20, 0.30], [0.40, -0.20, 0.30], + [0.40, -0.20, 0.70], [0.10, -0.20, 0.70]]) + B = np.array([[0.60, 1.20, 0.30], [0.90, 1.20, 0.30], + [0.90, 1.20, 0.70], [0.60, 1.20, 0.70]]) + fsA = uw.meshing.FaultSurface("West", A) + fsA.triangulate() + fsB = uw.meshing.FaultSurface("East", B) + fsB.triangulate() + net = uw.meshing.FaultNetwork([fsA, fsB], hierarchy=["West", "East"]) + net.prepare(h=H, ligament=1.0, verbose=False) + net.realisation, net.width = "split", WIDTH + net._build_3d_band(h_far=0.24, realisation="split", margin_rings=0.5, + carve_clearance=0.3, minCoords=(-0.5, -0.5, -0.5), + maxCoords=(1.5, 1.5, 1.5)) + mesh = net.mesh + comm = mesh.dm.comm.tompi4py() + info = {k: net.info[k] for k in ("n_regions", "n_gathered", "n_moved")} + assert all(g == info for g in comm.allgather(info)), info + assert info["n_regions"] == 2, info + # never both to one rank: what moved is less than the regions' size + assert info["n_moved"] < info["n_gathered"], info + # every rank holds cells, and both faults were split (their pair + # boundaries exist on the mesh) + local = int(mesh.dm.getHeightStratum(0)[1]) + assert comm.allreduce(local, op=min) > 0 + names = {b.name for b in mesh.boundaries} + for fault in ("West", "East"): + assert {f"{fault}Plus", f"{fault}Minus"} <= names, names From 0d5b8bb87f5244d924be9ab82e8c1f3080185ed3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 14:24:03 -0700 Subject: [PATCH 04/19] The geometric tail on a distributed band: zero the constrained rows by their global index, keep the better transfer (#671) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two defects in the co-located multigrid tail surfaced once a band could sit on more than one rank (the per-region placement of #670). The rotated prolongation zeroes the constrained rows — the rotated normals and the fault-pair contact rows — so they take no coarse correction. It enumerated this rank's entries of the velocity IS and passed those LOCAL block rows to zeroRows, which takes global ones: on every rank but the first the wrong rows were zeroed. Where the rows wrongly zeroed were a coarse DOF's only fine images, the Galerkin coarse operator got a zero row, and PETSc's own remedy (MatGalerkin writes an identity diagonal into zero rows) failed on a row whose diagonal was never allocated — the "New nonzero caused a malloc" in PCSetUp_MG. The rows are now the rank's ownership offset plus the local index. The "auto" transfer mode rebuilt any transfer with a zero column as a cross-partition one. On the balanced two-fault layout the co-partitioned build left 3 orphan columns (coarse DOFs under the band with no fine image) and the cross-partition build 16,791 of 96,009 at np=2 and 43,047 at np=4; the repair then injected a nearest fine DOF into every one, and the coarse operator was nonsense (the velocity KSP failed with reason -11 at its first iteration). The rebuild is now kept only where it leaves fewer orphans than the co-partitioned build, with a warning when it would have made things worse; the few genuine orphans are the repair's job, as before. Measured. Two faults a domain apart in one network, split, on a [-0.5, 1.5]^3 base: the contact solve now converges in one Newton step on custom-FMG at np=2 and np=4 with slips within 0.5% of the serial answer (the gap fill's node count varies by one or two with the partition, so the meshes are not identical). The crossing-patches fixture passes at np=2 and np=4; at np=3 it no longer fails in PCSetUp_MG but still reports the smallest piece's slip 10% low, which remains open in #671. The parallel custom-MG and rotated free-slip suites pass at np=2. ptest_0863's two-fault case now runs the solve and pins both. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/custom_mg.py | 29 ++++++++++++++--- src/underworld3/utilities/rotated_bc.py | 15 ++++++++- ...st_0863_fault_network_3d_width_parallel.py | 32 +++++++++++++++++++ 3 files changed, 71 insertions(+), 5 deletions(-) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index e463e31c8..44900f397 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -1443,10 +1443,31 @@ def build(self, solver): P = _build_parallel_transfer(*args) # "auto": a zero-column transfer means the coarse level is NOT # co-partitioned with the fine level (a fine leaf sits in an - # off-rank coarse cell). Rebuild it spanning partitions. - if (self.cross_partition == "auto" - and _count_zero_columns_parallel(P, comm) > 0): - P = _build_crosspart_transfer(*args) + # off-rank coarse cell). Rebuild it spanning partitions — + # and keep the rebuild only where it does better. On a + # placed band distributed over two ranks (#671) the + # co-partitioned build left 3 orphan columns and the + # cross-partition build 16,791 of 96,009; the repair then + # injected all of them and the coarse operator was + # nonsense (KSP reason -11 at the first iteration). The + # few genuine orphans are the repair's job below. + n_zero = _count_zero_columns_parallel(P, comm) + if self.cross_partition == "auto" and n_zero > 0: + P_x = _build_crosspart_transfer(*args) + n_x = _count_zero_columns_parallel(P_x, comm) + if n_x < n_zero: + P.destroy() + P = P_x + else: + P_x.destroy() + if n_x > n_zero: + import warnings + warnings.warn( + f"custom_mg: the cross-partition transfer " + f"{l - 1}->{l} left {n_x} coarse DOF(s) " + f"without a fine image where the " + f"co-partitioned one left {n_zero}; kept " + f"the co-partitioned transfer (#671).") # the same orphan repair the serial path has: a coarse DOF no # fine node reaches gets its nearest fine DOF as an injection if _count_zero_columns_parallel(P, comm) > 0: diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 47a06683a..2c963a305 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1566,8 +1566,21 @@ def _build_rotated_custom_Pl(solver, Q, normal_rows): vis = np.asarray(vel_is.getIndices()) g2blk = {int(g): k for k, g in enumerate(vis)} Qv = Q.createSubMatrix(vel_is, vel_is) - nrows_blk = sorted({g2blk[g] for g in normal_rows if g in g2blk}) Pfine = Qv.matMult(Ps[-1]) + # The constrained rows (rotated normals, fault-pair contact rows) take + # no coarse correction. ``g2blk`` enumerates THIS rank's entries of the + # velocity IS, so its values are local block rows; zeroRows takes + # global ones, which on the sub-matrix are the rank's ownership offset + # plus the local index. Passing the local index zeroed the wrong rows + # on every rank but the first (#671): a rank's own constrained rows + # kept their coarse correction while rows of the first rank lost + # theirs — and where those were the only fine images of a coarse + # DOF, the Galerkin coarse operator got a zero row PETSc could not + # repair (MatGalerkin writes an identity diagonal that was never + # allocated). + rstart = Pfine.getOwnershipRange()[0] + nrows_blk = sorted({rstart + g2blk[g] for g in normal_rows + if g in g2blk}) Pfine.zeroRows(nrows_blk, diag=0.0) # Remember an opportunistic mesh-owned pickup, so a later solve on this solver # (rotated or not) reuses the hierarchy instead of re-resolving it. diff --git a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py index 306158b41..f58d00d58 100644 --- a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py +++ b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py @@ -152,3 +152,35 @@ def test_two_faults_apart_are_placed_by_their_own_ranks(): names = {b.name for b in mesh.boundaries} for fault in ("West", "East"): assert {f"{fault}Plus", f"{fault}Minus"} <= names, names + + # The contact solve on the DISTRIBUTED band: the geometric tail must + # pair the base and the placed mesh across ranks. Two defects hid + # here (#671): the rotated prolongation zeroed constrained rows by + # their LOCAL index, so every rank but the first zeroed the wrong + # rows; and the "auto" transfer replaced a co-partitioned build with + # 3 orphan columns by a cross-partition one with 16,791, which the + # repair then filled with nonsense (KSP reason -11 at iteration 0). + x, y, z = mesh.X + v = uw.discretisation.MeshVariable("v2F", mesh, 3, degree=2) + p = uw.discretisation.MeshVariable("p2F", mesh, 1, degree=0, + continuous=False) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = [0.0, 0.0, 0.0] + for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"): + stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall) + net.apply(stokes) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + solve_info = net.solve(stokes) + assert solve_info.get("converged"), solve_info + assert solve_info.get("velocity_pc") == "custom-FMG", solve_info + # the serial answer on this box (np=1: West 0.13731, East 0.13653); + # the gap fill's node count varies by one or two with the partition, + # so the meshes differ slightly and the slips agree to 1%, not 5 + # digits (np=2 and np=4 give 0.13723 / 0.13712 and 0.13712 / 0.13712) + slips = net.slips(stokes) + for name, serial in (("West", 0.13731), ("East", 0.13653)): + peak = comm.allreduce(float(slips.get(name, 0.0)), op=max) + assert peak == pytest.approx(serial, rel=1e-2), (name, peak, serial) From 1047986d96f2ba3ef644de5400f2779ba670ff87 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 22:37:36 -0700 Subject: [PATCH 05/19] Fault networks in parallel: the design note, the layout throughput test, and the gather's growth as a parameter (#670, #671) docs/developer/design/fault-parallel-placement-2026-09.md is the governing note for parallel placement: the adapt-on-top frame, the three parallel strategies in the stack, what the gather moves and why, the measurements that showed it moving far too much, what PR #672 changed, the decision on long faults (faults as distributable objects, cut in the CAD along strike when one is too long for a rank, the strip along a cut plane gathered before the split), the rulings not to rebalance the child and to pre-refine the base, and what remains open. It is linked from the developer index (toctree and authority map), the conforming-surfaces subsystem note and the fault-networks page. The throughput test beside it holds everything fixed but where the faults sit relative to the seams: a 6 x 1 x 1 box at np=3, three faults one per slab, then one moved onto a seam, then all forced into one region, plus serial and GAMG baselines. Measured: per-region placement balances the mesh to 6% and gives a warm solve 1.8x serial where the old single gather gave 1.0x; one straddling fault costs 15%; the answer is layout-independent to four digits; and on this contact fixture the pressure Schur solve caps at 200 iterations in every configuration, so wall time is the pressure block's and GAMG beats the tail in wall time despite 26 velocity iterations against 4. _gather_regions takes the growth beyond the star as ``layers``. The placement was tried at ``layers=0`` (the star alone) and its own gate refused at np=2, 3 and 4: the carve drops cells beyond the marked vertices' star. Both the placement and the split keep one layer, and the note records the measurement. Underworld development team with AI support from Claude Code --- docs/advanced/fault-networks.md | 8 +- .../fault-parallel-placement-2026-09.md | 274 ++++++++++++++++++ .../design/fault_parallel_layouts.py | 104 +++++++ docs/developer/index.md | 2 + .../conforming-surfaces-and-fault-zones.md | 5 +- src/underworld3/utilities/fault_split.py | 2 +- src/underworld3/utilities/place_surface.py | 41 ++- 7 files changed, 416 insertions(+), 20 deletions(-) create mode 100644 docs/developer/design/fault-parallel-placement-2026-09.md create mode 100644 docs/developer/design/fault_parallel_layouts.py diff --git a/docs/advanced/fault-networks.md b/docs/advanced/fault-networks.md index 2424b00f7..9dbdb2cfc 100644 --- a/docs/advanced/fault-networks.md +++ b/docs/advanced/fault-networks.md @@ -292,8 +292,12 @@ gather-first, so the cells gmsh fills into the carved cavity — the band and its graded surround — live on one rank; only the base's far field is balanced. On the crossing-patches fixture that is 8012 of 8405 cells on one rank at np=4, and the solve is not faster than -serial: parallel is a correctness mode for this path, not a speed-up, -until the placed region is rebalanced. +serial. Since PR #672 the gather is per region: each zone, or +junction-connected cluster of zones, is placed and split by the rank that +holds it, and a zone whose shell is already interior to a rank is not +moved at all; a single long fault is still one rank's. The design and its +measurements are in the developer note +`fault-parallel-placement-2026-09.md`. ## Limitations diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md new file mode 100644 index 000000000..85c6d8837 --- /dev/null +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -0,0 +1,274 @@ +# Fault networks in parallel: what moves, why, and the decision (September 2026) + +This note records how the fault-network machinery behaves in parallel, +the measurements that established it, the design decision taken with +Louis on 2 September 2026, and the test that keeps the numbers honest. +It is the governing document for the parallel placement of faults; the +subsystem description is `../subsystems/conforming-surfaces-and-fault-zones.md` +and the multigrid side is `fault-patch-multigrid-2026-08.md`. Issues: +#670 (the gather), #671 (the tail), PRs #669 and #672. + +## The frame: adapt-on-top + +The base mesh is the persistent, distributed object. It holds the bulk +fields, it is refined in place by newest-vertex bisection co-partitioned +with itself, and its coarse levels form the multigrid tail. A fault +network is a *child* of that base: the band is placed into the base by +surgery, the faults are split on the placed mesh, and the child is +solved. The child is ephemeral. It is rebuilt when the base is adapted +or rebalanced, and it carries no state of its own: the fault's history +lives on the `FaultSurface` object, the bulk state on the base. + +Two consequences follow. The child does not have to inherit the base's +partition, only the two transfers must be partition-agnostic: bulk +fields base to child and back, and the multigrid tail from the base's +coarse levels to the child's finest. And the accepted trade for a fault +is **extra load on the rank that owns it, in exchange for no +communication** while it is solved. + +## Three parallel strategies coexist in the stack + +1. **The edge cut** (`line_cut`, 2-D): fully distributed. Each rank + splits its own edges where the surface crosses them; every crossing + is a pure function of the coordinates and the surface, every + tolerance is a global length, and nothing moves. +2. **Placement** (`place_sheet`, `place_thin_volume`, the 3-D band): + *gather-first*. Placement deletes the base vertices in the way and + asks gmsh to refill the cavity with the band's own mesh embedded in + it. gmsh is serial and a cavity cannot be carved or filled across a + partition seam, so the region around the band is moved onto one rank, + the surgery runs there, and every rank rebuilds its chart + collectively. The result is partition-independent by construction. +3. **The split** (contact): redistribute the fault's thin cell star onto + one rank, then split with serial topology, so every cut pair is born + on one rank. + +What the gather moves is **base mesh, never the fault**: the band is +generated from a few patch coordinates on the surgery rank. The seam +rule that forces the move is that no point the surgery deletes or +creates may be shared, because the rebuild carries the old star forest +over by renumbering. That needs exactly three layers of base cells on +one rank: the cells the carve drops, their vertex star (the ring the +fill attaches to), and one more layer so the ring's own points are +unshared. A shell about three cells thick around the band. The third +layer was tested for necessity: with the gather stopped at the star, the +placer's own gate ("the gathered region touches a shared point; the +gather mask under-reached") fires at np=2, 3 and 4 on the thin-volume +suite, so the carve drops cells beyond the marked vertices' star and the +layer stays. Two shells therefore touch at about six base cells' +separation, and two zones closer than that are one region. + +## What was wrong, measured + +All on the crossing-patches fixture of `test_0863` (unit cube, far-field +cells of about 0.126 after one refinement, band cells 0.08, width 0.04) +unless stated; the wide box is the same patches in `[-0.5, 1.5]^3`. + +**The fill is small; the gather was not.** The base has 5592 cells, the +placed mesh 8405, of which 1663 are band: the fill adds about 1150 cells, +all within a cell of the patches. But the gather moved 5087 of the 5592 +base cells onto one rank before any surgery, and at np=4 the surgery +rank ended with 8012 of 8405 cells. The mark reached two cell widths +beyond the carve *and then* `_gather_region` grew the star and the layer +from it: the margin was paid twice. + +| region on the base | cells | share | +|---|---|---| +| carve reach | 334 | 6% | +| reach + 1 cell (the dropped cells) | 1313 | 23% | +| reach + 2 cells (the old mark) | 3486 | 62% | +| old mark + star + layer | 5592 | 100% | + +The unit cube is four far-field cells across, so three layers from any +band reach the walls: **that fixture measures correctness and can never +show balance**. On the wide box the same patches gather 9831 cells under +the old mark and 5046 under the new one (CROSS patches; 8961 and 4196 for +the P_A/P_B pair), identical at np=2 and np=4. + +**One target for the whole network, and no return.** The placer took the +network as one assembly and gathered the union of every patch's region to +the single rank holding the most of it; two faults each interior to a +different rank were both moved to one of them. The moved cells stayed +there. The extra load therefore landed on the rank that was already +heaviest in that neighbourhood — anti-balanced — and was paid for with a +redistribution: neither half of the accepted trade. + +## What was done (PR #672) + +- **The mark covers only what the carve drops** (victims plus the crossed + cells' vertices, one cell diameter out), at every placement path. The + gather's own star-and-layer growth supplies the seam freedom. +- **Regions.** `_gather_regions` marks a chart of region ids, claims each + region's star and layer, merges regions whose shells touch (a collective + union-find), and sends each to the rank already holding most of it — or + leaves it where it is when it is interior to one rank. One shell + partition moves them all. +- **Per-region surgery.** One region per connected component of the + assembly; the owning ranks carve and fill their own components + concurrently, each on its compacted share of the assembly and skin, and + the collective rebuild sews them at once. Outcrop and ladder paths keep + one region (their bowl, cap and extrusion are single-rank). +- **The split follows the regions.** `split_faults(..., groups=...)` + redistributes per group through the same gather; the network passes it + the regions the placement reported (`info["embedded_regions"]`). +- **The tail on a distributed band** (#671). The rotated prolongation + zeroed its constrained rows by *local* block index where PETSc's + `zeroRows` takes global ones, so every rank but the first zeroed the + wrong rows; that is the "new nonzero caused a malloc" in `PCSetUp_MG` + (PETSc repairs a zero Galerkin row with an identity diagonal that was + never allocated). And the `cross_partition="auto"` rule replaced a + co-partitioned transfer with 3 orphan columns by a cross-partition one + with 16,791 of 96,009, which the repair then filled with nonsense. The + rows now carry the ownership offset, and the rebuild is kept only where + it leaves fewer orphans. + +Measured after: two zones a domain apart at np=2 move 304 cells to two +owners where 8451 moved to one; two faults in one network at np=2 are +both interior to their ranks, nothing moves, 14,339 / 14,332 cells; at +np=4 each region goes to a different rank. The contact solve on that +distributed band converges on the geometric tail in one Newton step at +np=2 and np=4, slips within 0.5% of serial (the gap fill's node count +varies by one or two with the partition, so the meshes are not +identical). + +## The decision + +Two designs were weighed: + +1. **Cut the faults at the partition boundary and solve locally.** Best + balance and no build communication, but every cost sits on the seam: + a fill conforming to a polyhedral seam that the band skin crosses (two + discrete surfaces to intersect consistently on two ranks), the split's + pair nodes on shared points (the one seam problem the multigrid note + lists as open), and a decomposition that changes with every rebalance. +2. **Faults as separate objects, distributed at will.** The surgery stays + serial per object, the pairs are rank-local by construction, an object + can be sent where the load is light, and the only communication is its + shell, once. Its limit is granularity: an object is a connected fault, + so one long fault is one rank's load, proportional to its area. + +**Ruling: the second is the primitive** — it is what is built — and the +first's cut moves from the partition seam into the CAD when a long fault +needs it. Cut a long fault along strike at chosen planes before meshing: +OCC cuts the band by a plane exactly, so the two skins meet on a +triangulated plane generated once, broadcast, and embedded verbatim by +both owners (the fill already gates for verbatim constraints); the base +cells straddling the plane are victims on both sides, so each rank's +cavity ring closes onto the same plane triangulation and nothing +discrete is intersected with anything discrete. The pieces are then +ordinary objects. What still crosses ranks is the split: the fault trace +runs through the cut plane, so the pair nodes there sit on shared +points, and the strip along the plane — a few cells wide — is gathered +to one of the two owners before splitting. Communication is then +proportional to the cut planes, not the fault. + +None of that is built, and it should not be until a model needs a fault +longer than one rank comfortably holds. The first thing to test when it +is: the split across one cut plane at np=3. + +Two further rulings from the same discussion: + +- **Do not rebalance the child after placement.** It fixes the solve + imbalance by destroying co-partition everywhere, and every transfer + then becomes a cross-partition search. Balance the base instead: refine + it in the broad region where faults can exist (the lithosphere, say) in + the *coarse* gmsh mesh, so the grading nests through the whole tail, + and let the partitioner balance that. +- **A rebalance of the base means a rebuild of the child.** That is + consistent with the faults being ephemeral; what must be verified is + that the base's own state (mesh variables, swarms) migrates cleanly + through a redistribution of a live mesh, which has not been exercised. + +## Still open + +- **np=3 answer on the crossing fixture** (#671): the smallest piece's + slip is 10% low. It no longer crashes; it is a wrong number. +- **The cross-partition transfer builder** loses 17 to 45 percent of the + coarse columns on a pair that is actually co-partitioned. It is only + prevented from being chosen when it is worse. +- **Slicing a long fault** (above). +- **The seam-straddling gauge**: the weak-plane slip gauge omits a probe + pair the rank does not own on both sides; unreachable while the + gathered region holds the probes, and worth a collective count if a + band is ever partitioned. + +## The throughput test + +The numbers above say what moves. The question that matters for a +time-stepping model is what the movement *costs* in the solve, and the +only fair measurement holds everything fixed except where the faults sit +relative to the seams. `fault_parallel_layouts.py` (beside this note) +does that on a 6 x 1 x 1 box whose np=3 partition is three slabs along x +with seams near x = 2 and 4, each slab about fifteen base cells long. +Four faults, 0.3 long and 0.2 tall, at least 1.5 apart so that no two +shells touch: one alone, a junction-connected pair as one cluster +(ligament 1.5; 1.0 degenerates the junction cut on this mesh), and one +more that the layout moves: + +- **local**: one fault per slab, the pair in the middle slab — nothing + should move; +- **straddle**: the fourth fault shifted onto the seam near x = 4 — it + is gathered to the majority rank, that rank carries its shell, and the + tail's transfers for the shell go cross-partition; +- **gathered**: every fault forced into one region — what the code did + before #672; +- **serial** and **GAMG instead of the tail** on the same faults, as the + baselines. + +Two fixture lessons, both found the hard way: the 6:1 box meshes with +cells up to 0.4 across, and a fault 0.4 tall in a unit-high box leaves +0.28 of clearance, so the carve reached the floor and refused ("the +cavity reached the domain wall"); and faults about eight cells apart +merge into one region, since two three-cell shells touch at six. + +The slips must agree across layouts to the fill's noise; the cost is in +the build time, the cells moved, the per-rank imbalance, the iteration +counts and the cold and warm solve times. "Warm" is a second full solve +from zero with the tail and the rotation reused — a repeat from the +converged solution takes no Newton step and measures nothing. Warm solve +time is what a time-stepping model pays per step. + +### Results (2 September 2026, 16-core workstation, 20,778-cell base) + +The junior of the crossing pair is consumed whole by the ligament cut at +1.5 on this mesh, so the fixture as run is three faults, one per slab, +and the "cluster" is a single fault. Warm is the second full solve from +zero; every solve took one Newton step. + +| layout | regions | cells moved | cells per rank (max/mean) | cold s | warm s | velocity its | pressure its | slips A / B / D | +|---|---|---|---|---|---|---|---|---| +| serial, tail | 3 | 0 | 20778 (1.00) | 97.7 | 81.0 | 3 | 200 | 0.11182 / 0.13117 / 0.11110 | +| np=3 local, tail | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 52.3 | 45.8 | 4 | 200 | 0.11181 / 0.13176 / 0.11108 | +| np=3 straddle, tail | 3 | 834 | 5624 / 6733 / 8438 (1.22) | 60.6 | 52.9 | 4 | 200 | 0.11181 / 0.13114 / 0.10926 (D moved) | +| np=3 gathered (pre-#672), tail | 1 | 5440 | 3514 / 13208 / 4054 (1.91) | 89.4 | 78.3 | 4 | 200 | 0.11182 / 0.13116 / 0.11111 | +| np=3 local, GAMG | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 24.9 | 17.6 | 26 | 200 | 0.11182 / 0.13176 / 0.11108 | + +What the table says: + +- **The answer is layout-independent.** A and D agree to four digits + across serial, local, straddle and gathered; B to 0.5%, which is the + gap fill's node-count noise between partitions. The straddle's D is a + different fault position and legitimately a different number. +- **Per-region placement is what makes np=3 worth running.** The old + single gather left one rank with 63% of the mesh and a warm solve of + 78 s against 81 s serial: no parallel gain at all. Per region, the + mesh is balanced to 6% and the warm solve is 46 s, 1.8 times serial + and 1.7 times the old gather. +- **A straddling fault costs about 15%** (53 s against 46 s): 834 cells + moved, the owner at 1.22 of the mean, and that shell's transfers + cross-partition. That is the price of one non-local fault on three + ranks, and it is the number the design decision rests on. +- **The pressure block is the wall-clock gate on this fixture, not the + velocity block.** Every configuration hits the pressure Schur solve's + 200-iteration cap, so wall time is dominated by pressure iterations, + each of which applies the velocity preconditioner. That is why GAMG, + at 26 velocity iterations against the tail's 4, is 2.6 times faster + here: its application is cheaper and the cap is the same. The tail's + advantage in velocity iterations is real and invisible in wall time + until the pressure solve converges. This is the #625 pressure cap, + measured again on a contact fixture; it is not a placement matter. + +The script beside this note (`fault_parallel_layouts.py`) regenerates +the table: + + mpirun -np 3 python -u fault_parallel_layouts.py -uw_layout local|straddle|gathered [-uw_tail 0] diff --git a/docs/developer/design/fault_parallel_layouts.py b/docs/developer/design/fault_parallel_layouts.py new file mode 100644 index 000000000..69b0e0ce6 --- /dev/null +++ b/docs/developer/design/fault_parallel_layouts.py @@ -0,0 +1,104 @@ +"""Throughput of the parallel fault network by LAYOUT, everything else +fixed: the same four faults (one alone, a junction-connected pair, one +more) on a 6 x 1 x 1 box whose np=3 partition is three slabs along x +(seams near x = 2 and 4), split realisation, contact solve. + + mpirun -np 3 python -u fault_parallel_layouts.py -uw_layout local|straddle|gathered [-uw_tail 0] + +local : one fault per slab, the pair in the middle slab -> nothing moves +straddle : the fourth fault shifted onto the seam near x = 4 -> gathered +gathered : every fault forced into ONE region (the pre-#672 behaviour) +-uw_tail 0: the geometric tail dropped -> the velocity block on GAMG + +Reports: build time, regions / cells moved, per-rank cells, Newton and +Krylov counts, cold and warm solve wall time, and the slips (which must +agree across layouts: the answer is partition-independent, the cost is +not). Design note: docs/developer/design/fault-parallel-placement-2026-09.md +""" +import time +import numpy as np +import underworld3 as uw + +params = uw.Params(layout=uw.Param("local", "local | straddle | gathered"), + tail=uw.Param(1, "1 = geometric tail (custom-FMG), 0 = GAMG"), + solve=uw.Param(1, "0 = build only")) +layout = str(params.layout) +if layout == "gathered": + # one region for the whole network: the old behaviour, for reference + from underworld3.utilities import place_surface as ps + _orig = ps._gather_regions + def one_region(dm, ids, verbose=False, layers=1): + ids = np.asarray(ids) + n_ids = dm.getComm().tompi4py().allreduce(int(ids.max()) if ids.size else 0, op=max) + work, n_region, n_moved, owner, _canon = _orig( + dm, (ids > 0).astype(np.int32), verbose=verbose, layers=layers) + return work, n_region, n_moved, owner, {k: 1 for k in range(1, n_ids + 1)} + ps._gather_regions = one_region + +from underworld3.utilities import place_surface as _ps +_g = _ps._gather_regions +_ps._gather_regions = lambda dm, ids, verbose=False, layers=1: _g(dm, ids, verbose=True, layers=layers) +H, W = 0.08, 0.04 +def patch(x0, x1, y, z0=0.3, z1=0.7): + return np.array([[x0, y, z0], [x1, y, z0], [x1, y, z1], [x0, y, z1]]) +# a 6 x 1 x 1 box: np=3 slabs of length 2 with seams near x = 2 and 4; +# faults at least 1.5 apart so no two shells touch, whichever layout +# patches 0.2 tall (z 0.4..0.6): the 6:1 box meshes with cells up to 0.4 +# across, and the carve refuses a cavity that reaches a wall +A = patch(0.40, 0.70, 0.50, 0.40, 0.60) # slab 1 +B = patch(2.20, 2.60, 0.50, 0.40, 0.60) # slab 2, senior of the pair +C = np.array([[2.20, 0.62, 0.42], [2.52, 0.30, 0.42], + [2.52, 0.30, 0.58], [2.20, 0.62, 0.58]]) # crosses B +xD = 5.20 if layout != "straddle" else 3.85 # slab 0, or ON the seam near 4 +D = patch(xD, xD + 0.30, 0.50, 0.40, 0.60) + +comm = uw.mpi.comm +t0 = time.perf_counter() +faults = [] +for name, P in (("A", A), ("B", B), ("C", C), ("D", D)): + f = uw.meshing.FaultSurface(name, P); f.triangulate(); faults.append(f) +net = uw.meshing.FaultNetwork(faults, hierarchy=["A", "B", "C", "D"]) +net.prepare(h=H, ligament=1.5, verbose=False) # 1.0 degenerates the junction cut on this mesh +net.realisation, net.width = "split", W +net._build_3d_band(h_far=0.24, realisation="split", margin_rings=0.5, + carve_clearance=0.3, minCoords=(0.0, 0.0, 0.0), + maxCoords=(6.0, 1.0, 1.0)) +t_build = time.perf_counter() - t0 +mesh = net.mesh +if not int(params.solve): + cells = comm.gather(int(mesh.dm.getHeightStratum(0)[1]), root=0) + if comm.rank == 0: + print(f"[layout] {layout} build-only np={comm.size}: regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']} cells/rank {cells}", flush=True) + raise SystemExit(0) +if not int(params.tail): + mesh._custom_mg_coarse_meshes = None +cells = comm.gather(int(mesh.dm.getHeightStratum(0)[1]), root=0) +band = comm.gather(int(np.count_nonzero(mesh.cells_labelled("Band", 71))), root=0) + +x, y, z = mesh.X +v = uw.discretisation.MeshVariable("v", mesh, 3, degree=2) +p = uw.discretisation.MeshVariable("p", mesh, 1, degree=0, continuous=False) +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 +stokes.bodyforce = [0.0, 0.0, 0.0] +for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"): + stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall) +net.apply(stokes) +stokes.petsc_use_pressure_nullspace = True +stokes.tolerance = 1e-5 +t1 = time.perf_counter(); info = net.solve(stokes); t_cold = time.perf_counter() - t1 +# warm: a full solve again, from zero, with the tail and the rotation reused +t2 = time.perf_counter(); info2 = net.solve(stokes, zero_init_guess=True); t_warm = time.perf_counter() - t2 +slips = net.slips(stokes) +peaks = {n: comm.allreduce(float(slips.get(n, 0.0)), op=max) for n in sorted(set(k for k, _ in net.prepared))} +if comm.rank == 0: + imb = max(cells) / (sum(cells) / len(cells)) + print(f"[layout] {layout} tail={int(params.tail)} np={comm.size}: build {t_build:.1f}s; " + f"regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']}; " + f"cells/rank {cells} (max/mean {imb:.2f}) band/rank {band}", flush=True) + print(f"[layout] {layout} tail={int(params.tail)}: cold {t_cold:.1f}s warm {t_warm:.1f}s; " + f"pc={info.get('velocity_pc')} newton={info.get('nonlinear_iterations')} " + f"converged={info.get('converged')} vel_its={info.get('vel_its_last')} pres_its={info.get('pres_its_last')}; " + f"warm newton={info2.get('nonlinear_iterations')} vel_its={info2.get('vel_its_last')} reused={info2.get('rotation_reused')}", flush=True) + print(f"[layout] {layout} tail={int(params.tail)}: slips " + " ".join(f"{n}={s:.5f}" for n, s in peaks.items()), flush=True) diff --git a/docs/developer/index.md b/docs/developer/index.md index c32976e1f..f7a2893b0 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -37,6 +37,7 @@ same topic are reference or historical material subordinate to the governing doc | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) | | Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) | | Rotated free-slip & wall-normal datum | [subsystems/rotated-freeslip.md](subsystems/rotated-freeslip.md) | +| Fault networks in parallel (placement, split, tail) | [design/fault-parallel-placement-2026-09.md](design/fault-parallel-placement-2026-09.md) (subsystem reference: [conforming-surfaces-and-fault-zones.md](subsystems/conforming-surfaces-and-fault-zones.md)) | | Analytic & benchmark solutions | [subsystems/analytic-solutions.md](subsystems/analytic-solutions.md) | | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) | | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) | @@ -165,6 +166,7 @@ design/TURBULENCE_MODEL_DESIGN design/declined-coord-units-proposal design/nonlinear-solver-homotopy-warmstart design/fault-zone-hybrid-architecture +design/fault-parallel-placement-2026-09 ``` ```{toctree} diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index 27e4acc90..c2d9be5c0 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -330,7 +330,10 @@ region: their bowl, cap and extrusion machinery is single-rank. One limit remains: a small domain cannot hold a shell at all (three base cells reaching the walls is the whole box, which is what the crossing-patches test fixture does), so that fixture measures correctness -only, never balance. +only, never balance. The measurements, the design decision on cutting long +faults, and the layout throughput test are recorded in +`../design/fault-parallel-placement-2026-09.md`, the governing note for +parallel placement. ## The thin volume: finite-width zones, junctions in the volume diff --git a/src/underworld3/utilities/fault_split.py b/src/underworld3/utilities/fault_split.py index 896b6f14a..1e0a69587 100644 --- a/src/underworld3/utilities/fault_split.py +++ b/src/underworld3/utilities/fault_split.py @@ -1266,7 +1266,7 @@ def _redistribute_fault_interior(dm, labels, verbose=False, groups=None): names = sorted(name for name, _value in labels) raise RuntimeError( f"fault_split: no facets labelled {names} found on any rank.") - work, n_region, n_moved, owner, _canon = _gather_regions(dm, ids) + work, n_region, n_moved, owner, _canon = _gather_regions(dm, ids, layers=1) if work is dm: work = dm.clone() # the contract: a new dm, untouched input if verbose: diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 45d97a79d..92b6d0d2e 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -1643,12 +1643,12 @@ def _gather_region(dm, vertex_mark_chart, verbose=False): """ ids = (np.asarray(vertex_mark_chart) != 0).astype(np.int32) work, n_region, _n_moved, _owner, _canon = _gather_regions( - dm, ids, verbose=verbose) + dm, ids, verbose=verbose, layers=1) return work, n_region -def _gather_regions(dm, vertex_region_chart, verbose=False): - """Redistribute so each marked REGION's cell star (+1 layer) is one rank's. +def _gather_regions(dm, vertex_region_chart, verbose=False, layers=1): + """Redistribute so each marked REGION's cell star (+ ``layers``) is one rank's. A mask-driven port of the contact stream's ``_redistribute_fault_interior`` (feature/fault-split-node c8693579 / 1d487319; the label-driven original @@ -1663,6 +1663,14 @@ def _gather_regions(dm, vertex_region_chart, verbose=False): marked stars move — everything else keeps its load-balanced home — via one shell partition. + ``layers`` is the growth beyond the star, one by default: the star is + the marked vertices' cells, and the layer makes every point in the + closure of a star cell unshared. Both the placement and the split + need it. Measured (#670): with ``layers=0`` the placement's own gate + ("the gathered region touches a shared point") fires at np=2, 3 and + 4 on the thin-volume suite — the carve drops cells beyond the marked + vertices' star, so the star alone under-reaches. + Returns ``(new_dm, n_region, n_moved, owner, canon)``: the global size of the regions (star and layer — the seam rule's footprint, a function of the mesh alone), the count of cells that changed rank (``0`` when @@ -1712,19 +1720,20 @@ def claim_cells(chart): cell_region[c] = max(cell_region[c], k) claim_cells(reconcile(ids_in.copy())) - # One growth layer: the surgery needs every point in the closure of a - # region cell unshared, and a point is unshared exactly when all its - # incident cells are co-resident. - layer = np.zeros(pEnd - pStart, dtype=np.int32) - for c in np.flatnonzero(cell_region): - k = int(cell_region[c]) - for q in work.getTransitiveClosure(int(c) + cS)[0]: - if vS <= int(q) < vE: - i = int(q) - pStart - if layer[i] and layer[i] != k: - pairs.add((int(layer[i]), k)) - layer[i] = max(layer[i], k) - claim_cells(reconcile(layer)) + # Growth layers: a point is unshared exactly when all its incident + # cells are co-resident, so each layer makes the previous cells' + # closures unshared. + for _grow in range(int(layers)): + layer = np.zeros(pEnd - pStart, dtype=np.int32) + for c in np.flatnonzero(cell_region): + k = int(cell_region[c]) + for q in work.getTransitiveClosure(int(c) + cS)[0]: + if vS <= int(q) < vE: + i = int(q) - pStart + if layer[i] and layer[i] != k: + pairs.add((int(layer[i]), k)) + layer[i] = max(layer[i], k) + claim_cells(reconcile(layer)) # merge touching regions: union-find on the gathered pair set, the same # on every rank From 1edb6f5037d6a97063b30d27902453ff9a8d5563 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 08:43:10 -0700 Subject: [PATCH 06/19] The layout table with the pressure solve converging, and the pressure tolerance as a measurement knob (#625, #670) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The first pass of the layout throughput test ran every configuration's pressure Schur solve to its 200-iteration cap. The monitor shows why: the pressure residual falls by 5e-3 in twenty iterations, then creeps to a floor of about 7e-6 relative — the inexact velocity solve inside the Schur application (3.3e-7 of its own residual) amplified by the Schur complement — while the rotated path asks the pressure solve for 1e-6. The margin rule is inverted between the two inner solves. Each wasted pressure iteration is a velocity-preconditioner apply, which is what made GAMG look 2.6 times faster than the tail. solver._rotated_pres_rtol is the measurement knob (TODO(MEASURE)); the default stays. At 1e-4 the pressure solve converges in 26 iterations and the table of record reads: serial 14.0 s warm, np=3 local 6.8 s (2.1x), one straddling fault 8.7 s (+28%), GAMG 3.3 s. The answer is unchanged to four digits. GAMG remains twice as fast as the tail on this linear, uniform-viscosity fixture; the note says so and why that is not the placement's question. Underworld development team with AI support from Claude Code --- .../fault-parallel-placement-2026-09.md | 81 +++++++++++++------ .../design/fault_parallel_layouts.py | 7 +- src/underworld3/utilities/rotated_bc.py | 11 ++- 3 files changed, 70 insertions(+), 29 deletions(-) diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md index 85c6d8837..369fbaf60 100644 --- a/docs/developer/design/fault-parallel-placement-2026-09.md +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -228,45 +228,74 @@ from zero with the tail and the rotation reused — a repeat from the converged solution takes no Newton step and measures nothing. Warm solve time is what a time-stepping model pays per step. -### Results (2 September 2026, 16-core workstation, 20,778-cell base) +### Results (2–3 September 2026, 16-core workstation, 20,778-cell base) The junior of the crossing pair is consumed whole by the ligament cut at 1.5 on this mesh, so the fixture as run is three faults, one per slab, and the "cluster" is a single fault. Warm is the second full solve from -zero; every solve took one Newton step. +zero; every solve took one Newton step; the tolerance is 1e-5. + +**The first pass was uninterpretable, and why matters.** With the rotated +path's default pressure sub-solve tolerance (a tenth of the solver +tolerance) every configuration ran the pressure Schur solve to its +200-iteration cap. The monitor showed the mechanism: the pressure +residual falls by 5e-3 in twenty iterations and then creeps — 2.6e-8 at +20, 5.5e-9 at 50, 8e-10 at 199 — against a target of 1.2e-10. The floor +is set by the inexact velocity solve inside the Schur application (the +velocity sub-solve stops at 3.3e-7 of its own residual, and the Schur +complement amplifies that to about 7e-6 relative in the pressure +residual), so the margin rule "inner converges well below outer" is +inverted between the two inner solves: the pressure cannot converge +below the velocity's floor. Each wasted pressure iteration is a +velocity-preconditioner apply, which is why in that pass GAMG (26 +velocity iterations, cheap applies) was 2.6 times faster in wall time +than the tail (4 iterations, dear applies). Recorded on #625. The +pressure tolerance is now an attribute knob (`solver._rotated_pres_rtol`, +a `TODO(MEASURE)`); the table of record uses 1e-4, ten times the solver +tolerance, at which the pressure solve converges in 26 iterations and +the outer FGMRES does the rest. The default is unchanged: what it should +be is a solver-configuration decision, not a placement one. | layout | regions | cells moved | cells per rank (max/mean) | cold s | warm s | velocity its | pressure its | slips A / B / D | |---|---|---|---|---|---|---|---|---| -| serial, tail | 3 | 0 | 20778 (1.00) | 97.7 | 81.0 | 3 | 200 | 0.11182 / 0.13117 / 0.11110 | -| np=3 local, tail | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 52.3 | 45.8 | 4 | 200 | 0.11181 / 0.13176 / 0.11108 | -| np=3 straddle, tail | 3 | 834 | 5624 / 6733 / 8438 (1.22) | 60.6 | 52.9 | 4 | 200 | 0.11181 / 0.13114 / 0.10926 (D moved) | -| np=3 gathered (pre-#672), tail | 1 | 5440 | 3514 / 13208 / 4054 (1.91) | 89.4 | 78.3 | 4 | 200 | 0.11182 / 0.13116 / 0.11111 | -| np=3 local, GAMG | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 24.9 | 17.6 | 26 | 200 | 0.11182 / 0.13176 / 0.11108 | +| serial, tail | 3 | 0 | 20778 (1.00) | 31.4 | 14.0 | 3 | 26 | 0.11182 / 0.13117 / 0.11110 | +| np=3 local, tail | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 13.5 | 6.8 | 4 | 26 | 0.11181 / 0.13176 / 0.11108 | +| np=3 straddle, tail | 3 | 834 | 5624 / 6733 / 8438 (1.22) | 16.6 | 8.7 | 4 | 26 | 0.11182 / 0.13115 / 0.10927 (D moved) | +| np=3 local, GAMG | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 9.3 | 3.3 | 26 | 26 | 0.11182 / 0.13176 / 0.11108 | +| *capped pass, for the record:* | | | | | | | | | +| serial, tail, pressure at cap | 3 | 0 | 20778 (1.00) | 97.7 | 81.0 | 3 | 200 | same | +| np=3 local, tail, pressure at cap | 3 | 122 | (1.06) | 52.3 | 45.8 | 4 | 200 | same | +| np=3 straddle, tail, pressure at cap | 3 | 834 | (1.22) | 60.6 | 52.9 | 4 | 200 | same | +| np=3 gathered (pre-#672), pressure at cap | 1 | 5440 | 3514 / 13208 / 4054 (1.91) | 89.4 | 78.3 | 4 | 200 | 0.11182 / 0.13116 / 0.11111 | +| np=3 local, GAMG, pressure at cap | 3 | 122 | (1.06) | 24.9 | 17.6 | 26 | 200 | same | What the table says: - **The answer is layout-independent.** A and D agree to four digits - across serial, local, straddle and gathered; B to 0.5%, which is the - gap fill's node-count noise between partitions. The straddle's D is a - different fault position and legitimately a different number. + across serial, local, straddle and gathered, and across both pressure + tolerances; B to 0.5%, which is the gap fill's node-count noise + between partitions. The straddle's D is a different fault position and + legitimately a different number. - **Per-region placement is what makes np=3 worth running.** The old - single gather left one rank with 63% of the mesh and a warm solve of - 78 s against 81 s serial: no parallel gain at all. Per region, the - mesh is balanced to 6% and the warm solve is 46 s, 1.8 times serial - and 1.7 times the old gather. -- **A straddling fault costs about 15%** (53 s against 46 s): 834 cells - moved, the owner at 1.22 of the mean, and that shell's transfers + single gather left one rank with 63% of the mesh and (in the capped + pass) a warm solve of 78 s against 81 s serial: no parallel gain at + all. Per region, the mesh is balanced to 6% and the warm solve is + 2.1 times serial (6.8 s against 14.0 s). +- **A straddling fault costs about 28%** (8.7 s against 6.8 s): 834 + cells moved, the owner at 1.22 of the mean, and that shell's transfers cross-partition. That is the price of one non-local fault on three - ranks, and it is the number the design decision rests on. -- **The pressure block is the wall-clock gate on this fixture, not the - velocity block.** Every configuration hits the pressure Schur solve's - 200-iteration cap, so wall time is dominated by pressure iterations, - each of which applies the velocity preconditioner. That is why GAMG, - at 26 velocity iterations against the tail's 4, is 2.6 times faster - here: its application is cheaper and the cap is the same. The tail's - advantage in velocity iterations is real and invisible in wall time - until the pressure solve converges. This is the #625 pressure cap, - measured again on a contact fixture; it is not a placement matter. + ranks, and it is the number the design decision rests on. (In the + capped pass the same difference read 15%, diluted by the wasted + pressure iterations.) +- **GAMG is still twice as fast as the tail in wall time on this + fixture** (3.3 s against 6.8 s warm) with the pressure solve fixed, + at 26 velocity iterations against 4. The tail's application is dearer + by more than the iteration ratio buys back here. This fixture is + linear viscous, uniform viscosity, one Newton step — the regime GAMG + is built for; the tail's measured advantage (#579, #576) is on banded + contrast and nonlinear solves, which this test does not exercise. It + is a red flag worth its own measurement on a contrast fixture, and it + is not a placement matter. The script beside this note (`fault_parallel_layouts.py`) regenerates the table: diff --git a/docs/developer/design/fault_parallel_layouts.py b/docs/developer/design/fault_parallel_layouts.py index 69b0e0ce6..91d37b9fd 100644 --- a/docs/developer/design/fault_parallel_layouts.py +++ b/docs/developer/design/fault_parallel_layouts.py @@ -21,7 +21,8 @@ params = uw.Params(layout=uw.Param("local", "local | straddle | gathered"), tail=uw.Param(1, "1 = geometric tail (custom-FMG), 0 = GAMG"), - solve=uw.Param(1, "0 = build only")) + solve=uw.Param(1, "0 = build only"), + pres_rtol=uw.Param(0.0, "pressure sub-solve rtol (0 = the path's default, 0.1 x tol)")) layout = str(params.layout) if layout == "gathered": # one region for the whole network: the old behaviour, for reference @@ -87,6 +88,8 @@ def patch(x0, x1, y, z0=0.3, z1=0.7): net.apply(stokes) stokes.petsc_use_pressure_nullspace = True stokes.tolerance = 1e-5 +if float(params.pres_rtol) > 0: + stokes._rotated_pres_rtol = float(params.pres_rtol) t1 = time.perf_counter(); info = net.solve(stokes); t_cold = time.perf_counter() - t1 # warm: a full solve again, from zero, with the tail and the rotation reused t2 = time.perf_counter(); info2 = net.solve(stokes, zero_init_guess=True); t_warm = time.perf_counter() - t2 @@ -97,7 +100,7 @@ def patch(x0, x1, y, z0=0.3, z1=0.7): print(f"[layout] {layout} tail={int(params.tail)} np={comm.size}: build {t_build:.1f}s; " f"regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']}; " f"cells/rank {cells} (max/mean {imb:.2f}) band/rank {band}", flush=True) - print(f"[layout] {layout} tail={int(params.tail)}: cold {t_cold:.1f}s warm {t_warm:.1f}s; " + print(f"[layout] {layout} tail={int(params.tail)} pres_rtol={params.pres_rtol}: cold {t_cold:.1f}s warm {t_warm:.1f}s; " f"pc={info.get('velocity_pc')} newton={info.get('nonlinear_iterations')} " f"converged={info.get('converged')} vel_its={info.get('vel_its_last')} pres_its={info.get('pres_its_last')}; " f"warm newton={info2.get('nonlinear_iterations')} vel_its={info2.get('vel_its_last')} reused={info2.get('rotation_reused')}", flush=True) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 2c963a305..82336b73e 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1719,7 +1719,16 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # iteration counts (measured, both velocity-block routes, isotropic and # TI); a too-loose inner solve fails by silent stagnation, not loudly. "fieldsplit_pres_ksp_type": "fgmres", - "fieldsplit_pres_ksp_rtol": str(tol * 0.1), + # TODO(MEASURE): the 0.1 margin sits BELOW the floor the inexact + # velocity solve leaves in the Schur residual (measured on the + # fault-network layouts: the pressure residual falls 5e-3 in 20 + # iterations, then creeps to the 200 cap at ~7e-6 relative, each + # iteration a velocity-preconditioner apply). The attribute is + # the knob for the measurement in + # docs/developer/design/fault-parallel-placement-2026-09.md; + # the resolved rule replaces it. + "fieldsplit_pres_ksp_rtol": str( + getattr(solver, "_rotated_pres_rtol", tol * 0.1)), "fieldsplit_pres_ksp_max_it": "200", } if Mp is not None: From 6e3362fd0c0f389dc03db7df5fe0aae505b4992e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 09:02:04 -0700 Subject: [PATCH 07/19] The layout table at the tolerance the resolution deserves; the penalty and the traction correction (#625, #670) At a solver tolerance of 1e-3 the default margins converge the pressure sub-solve in 26 iterations with no knob: the fixture's 1e-5 was one or two orders stricter than a mesh of 0.08 to 0.2 cells deserves, which is what put the pressure solve at its cap. The table carries that block too: serial 10.1 s warm, np=3 local 5.4 s, straddle 5.3 s, GAMG 2.2 s; the slips move by 0.2 to 0.7% between tolerances. The cost of one non-local fault is stated as a bound (at most a quarter of the solve, possibly nothing measurable) rather than a point value, since the three passes read 28%, 15% and none on single runs. The penalty is recorded as a recovery matter, not a discretisation one: the traction is the multiplier plus the augmentation, and the Coulomb fault law's reaction-fed normal stress needs the same correction before a penalty is used with friction (Louis's note on non-planar boundaries). Penalty 1 at 1e-5 does not lift the cap and shifts the slips by 0.6 to 0.9% at this resolution. Underworld development team with AI support from Claude Code --- .../fault-parallel-placement-2026-09.md | 46 +++++++++++++++---- .../design/fault_parallel_layouts.py | 10 ++-- 2 files changed, 43 insertions(+), 13 deletions(-) diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md index 369fbaf60..4eee1614d 100644 --- a/docs/developer/design/fault-parallel-placement-2026-09.md +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -187,6 +187,19 @@ Two further rulings from the same discussion: coarse columns on a pair that is actually co-partitioned. It is only prevented from being chosen when it is worse. - **Slicing a long fault** (above). +- **The penalty.** The runs here use `penalty = 0`, the solver default, + which the solver holds at zero because a grad-div penalty of 10 + corrupted the vertex-sampled dynamic topography recovered from the + free-slip reaction (test_1018). That is a recovery defect, not a + reason to forgo the penalty: the traction is the multiplier plus the + augmentation, r(u·n − ũ_n), and at a viscosity contrast of 1e6 the + multiplier alone carries almost none of it (Louis's note, + https://www.underworldcode.org/boundary-conditions-on-non-planar-boundaries). + The Coulomb fault law takes its normal stress from the same reaction, + so the correction belongs there too before a penalty is used with + frictional faults. Measured on this fixture: penalty 1 at tolerance + 1e-5 does not lift the pressure cap and shifts the slips by 0.6 to + 0.9% at this resolution. - **The seam-straddling gauge**: the weak-plane slip gauge omits a probe pair the rank does not own on both sides; unreachable while the gathered region holds the probes, and worth a collective count if a @@ -253,8 +266,13 @@ than the tail (4 iterations, dear applies). Recorded on #625. The pressure tolerance is now an attribute knob (`solver._rotated_pres_rtol`, a `TODO(MEASURE)`); the table of record uses 1e-4, ten times the solver tolerance, at which the pressure solve converges in 26 iterations and -the outer FGMRES does the rest. The default is unchanged: what it should -be is a solver-configuration decision, not a placement one. +the outer FGMRES does the rest. The default is unchanged, and the +better reading is Louis's: the fixture's tolerance of 1e-5 is one or two +orders stricter than a mesh of 0.08 to 0.2 cells deserves, and at 1e-3 +the default margins converge the pressure solve in the same 26 +iterations with no knob at all (the third block of the table). The +slips move by 0.2 to 0.7% between the two tolerances, which is the +loose tolerance's own noise, still below the fill's partition noise. | layout | regions | cells moved | cells per rank (max/mean) | cold s | warm s | velocity its | pressure its | slips A / B / D | |---|---|---|---|---|---|---|---|---| @@ -262,6 +280,11 @@ be is a solver-configuration decision, not a placement one. | np=3 local, tail | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 13.5 | 6.8 | 4 | 26 | 0.11181 / 0.13176 / 0.11108 | | np=3 straddle, tail | 3 | 834 | 5624 / 6733 / 8438 (1.22) | 16.6 | 8.7 | 4 | 26 | 0.11182 / 0.13115 / 0.10927 (D moved) | | np=3 local, GAMG | 3 | 122 | 6727 / 6733 / 7317 (1.06) | 9.3 | 3.3 | 26 | 26 | 0.11182 / 0.13176 / 0.11108 | +| *tolerance 1e-3, default margins (the tolerance this resolution deserves):* | | | | | | | | | +| serial, tail | 3 | 0 | 20778 (1.00) | 26.1 | 10.1 | 2 | 26 | 0.11187 / 0.13104 / 0.11116 | +| np=3 local, tail | 3 | 122 | (1.06) | 15.1 | 5.4 | 2 | 26 | 0.11187 / 0.13152 / 0.11115 | +| np=3 straddle, tail | 3 | 834 | (1.22) | 12.8 | 5.3 | 2 | 26 | 0.11197 / 0.13069 / 0.10942 (D moved) | +| np=3 local, GAMG | 3 | 122 | (1.06) | 8.3 | 2.2 | 16 | 26 | 0.11164 / 0.13162 / 0.11104 | | *capped pass, for the record:* | | | | | | | | | | serial, tail, pressure at cap | 3 | 0 | 20778 (1.00) | 97.7 | 81.0 | 3 | 200 | same | | np=3 local, tail, pressure at cap | 3 | 122 | (1.06) | 52.3 | 45.8 | 4 | 200 | same | @@ -281,15 +304,18 @@ What the table says: pass) a warm solve of 78 s against 81 s serial: no parallel gain at all. Per region, the mesh is balanced to 6% and the warm solve is 2.1 times serial (6.8 s against 14.0 s). -- **A straddling fault costs about 28%** (8.7 s against 6.8 s): 834 - cells moved, the owner at 1.22 of the mean, and that shell's transfers - cross-partition. That is the price of one non-local fault on three - ranks, and it is the number the design decision rests on. (In the - capped pass the same difference read 15%, diluted by the wasted - pressure iterations.) +- **A straddling fault costs between nothing and 28%.** At tolerance + 1e-5 it read 8.7 s against 6.8 s (28%); at 1e-3, 5.3 s against 5.4 s + (none); in the capped pass, 15%. These are single runs on a shared + workstation, so the honest statement is that one non-local fault on + three ranks costs at most a quarter of the solve and may cost nothing + measurable: 834 cells moved, the owner at 1.22 of the mean, and that + shell's transfers cross-partition. The design decision rests on that + bound, not on a point value. - **GAMG is still twice as fast as the tail in wall time on this - fixture** (3.3 s against 6.8 s warm) with the pressure solve fixed, - at 26 velocity iterations against 4. The tail's application is dearer + fixture** (3.3 s against 6.8 s warm at 1e-5; 2.2 s against 5.4 s at + 1e-3) with the pressure solve converging, at 16 to 26 velocity + iterations against 2 to 4. The tail's application is dearer by more than the iteration ratio buys back here. This fixture is linear viscous, uniform viscosity, one Newton step — the regime GAMG is built for; the tail's measured advantage (#579, #576) is on banded diff --git a/docs/developer/design/fault_parallel_layouts.py b/docs/developer/design/fault_parallel_layouts.py index 91d37b9fd..c4dd2e470 100644 --- a/docs/developer/design/fault_parallel_layouts.py +++ b/docs/developer/design/fault_parallel_layouts.py @@ -22,7 +22,9 @@ params = uw.Params(layout=uw.Param("local", "local | straddle | gathered"), tail=uw.Param(1, "1 = geometric tail (custom-FMG), 0 = GAMG"), solve=uw.Param(1, "0 = build only"), - pres_rtol=uw.Param(0.0, "pressure sub-solve rtol (0 = the path's default, 0.1 x tol)")) + pres_rtol=uw.Param(0.0, "pressure sub-solve rtol (0 = the path's default, 0.1 x tol)"), + tol=uw.Param(1e-5, "solver tolerance"), + penalty=uw.Param(0.0, "grad-div penalty (0 = the default, none)")) layout = str(params.layout) if layout == "gathered": # one region for the whole network: the old behaviour, for reference @@ -87,7 +89,9 @@ def patch(x0, x1, y, z0=0.3, z1=0.7): stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall) net.apply(stokes) stokes.petsc_use_pressure_nullspace = True -stokes.tolerance = 1e-5 +stokes.tolerance = float(params.tol) +if float(params.penalty) > 0: + stokes.penalty = float(params.penalty) if float(params.pres_rtol) > 0: stokes._rotated_pres_rtol = float(params.pres_rtol) t1 = time.perf_counter(); info = net.solve(stokes); t_cold = time.perf_counter() - t1 @@ -100,7 +104,7 @@ def patch(x0, x1, y, z0=0.3, z1=0.7): print(f"[layout] {layout} tail={int(params.tail)} np={comm.size}: build {t_build:.1f}s; " f"regions {net.info['n_regions']} gathered {net.info['n_gathered']} moved {net.info['n_moved']}; " f"cells/rank {cells} (max/mean {imb:.2f}) band/rank {band}", flush=True) - print(f"[layout] {layout} tail={int(params.tail)} pres_rtol={params.pres_rtol}: cold {t_cold:.1f}s warm {t_warm:.1f}s; " + print(f"[layout] {layout} tail={int(params.tail)} tol={params.tol} penalty={params.penalty} pres_rtol={params.pres_rtol}: cold {t_cold:.1f}s warm {t_warm:.1f}s; " f"pc={info.get('velocity_pc')} newton={info.get('nonlinear_iterations')} " f"converged={info.get('converged')} vel_its={info.get('vel_its_last')} pres_its={info.get('pres_its_last')}; " f"warm newton={info2.get('nonlinear_iterations')} vel_its={info2.get('vel_its_last')} reused={info2.get('rotation_reused')}", flush=True) From 17d67b40427a18bef608e6ddc417c6a5a9587bb4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 09:15:03 -0700 Subject: [PATCH 08/19] The hierarchy of the layout fixture, seen: one coarsening and a placed level; the preconditioner question deferred to 2-D (#670) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The renders (script beside the note) show the tail as the gmsh box at 0.24, its single bisection, and the placed mesh differing in 1,314 cells in three pockets — the finest transfer is the identity on 94% of its rows, and the band sits two cells across inside fill four to eight times larger with no grading, since the band builder hard-codes one refinement. Louis's ruling: multigrid is unsurprisingly not useful on this isoviscous two-grid fixture; the preconditioner study moves to 2-D. Recorded so the GAMG row is not read as a verdict. Underworld development team with AI support from Claude Code --- .../fault-parallel-placement-2026-09.md | 18 +++- .../design/fault_hierarchy_render.py | 101 ++++++++++++++++++ 2 files changed, 117 insertions(+), 2 deletions(-) create mode 100644 docs/developer/design/fault_hierarchy_render.py diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md index 4eee1614d..998126616 100644 --- a/docs/developer/design/fault-parallel-placement-2026-09.md +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -320,8 +320,22 @@ What the table says: linear viscous, uniform viscosity, one Newton step — the regime GAMG is built for; the tail's measured advantage (#579, #576) is on banded contrast and nonlinear solves, which this test does not exercise. It - is a red flag worth its own measurement on a contrast fixture, and it - is not a placement matter. + is not a placement matter. Seen rather than inferred (the renders in + `~/+Simulations/fault_network_3d_parallel/figures/`, 3 September): the + tail here is ONE genuine coarsening — the gmsh box at 0.24 (2433 cells, + the redundant LU coarse solve at about 13,000 velocity DOFs), its + single bisection (19,464), and the placed mesh (20,778) which differs + from the level below in 1,314 cells in three pockets, so the finest + transfer is the identity on 94% of its rows. The band is two cells + across its width inside fill cells four to eight times larger with no + grading between, because the band builder hard-codes one refinement of + the far-field box. A two-grid method with custom transfers and a large + coarse LU has nothing to gain over GAMG on an isoviscous problem, and + that is what the row shows. **Ruling (Louis, 3 September): the + preconditioner question is deferred to 2-D**, where a real hierarchy + (two or three doublings to the same finest cell) and a viscosity + contrast in the band are cheap to build and to see; this 3-D fixture + measures placement and balance, not multigrid. The script beside this note (`fault_parallel_layouts.py`) regenerates the table: diff --git a/docs/developer/design/fault_hierarchy_render.py b/docs/developer/design/fault_hierarchy_render.py new file mode 100644 index 000000000..432706b89 --- /dev/null +++ b/docs/developer/design/fault_hierarchy_render.py @@ -0,0 +1,101 @@ +"""The multigrid hierarchy of the layout fixture, seen: each level of the +tail on the z = 0.5 slice (through all three faults), a zoom on one band +with its fill, and at np=3 the partition on the same slice. + + python -u render_hierarchy.py # levels + zoom (serial) + mpirun -np 3 python -u render_hierarchy.py # partition (rank 0 renders) +""" +import numpy as np +import underworld3 as uw +import underworld3.visualisation as vis +import pyvista as pv +from mpi4py import MPI +pv.OFF_SCREEN = True + +FIG = "/Users/lmoresi/+Simulations/fault_network_3d_parallel/figures" +H, W = 0.08, 0.04 +def patch(x0, x1, y, z0=0.4, z1=0.6): + return np.array([[x0, y, z0], [x1, y, z0], [x1, y, z1], [x0, y, z1]]) +A = patch(0.40, 0.70, 0.50); B = patch(2.20, 2.60, 0.50) +C = np.array([[2.20, 0.62, 0.42], [2.52, 0.30, 0.42], [2.52, 0.30, 0.58], [2.20, 0.62, 0.58]]) +D = patch(5.20, 5.50, 0.50) +faults = [] +for name, P in (("A", A), ("B", B), ("C", C), ("D", D)): + f = uw.meshing.FaultSurface(name, P); f.triangulate(); faults.append(f) +net = uw.meshing.FaultNetwork(faults, hierarchy=["A", "B", "C", "D"]) +net.prepare(h=H, ligament=1.5, verbose=False) +net.realisation, net.width = "split", W +net._build_3d_band(h_far=0.24, realisation="split", margin_rings=0.5, carve_clearance=0.3, + minCoords=(0.0, 0.0, 0.0), maxCoords=(6.0, 1.0, 1.0)) +mesh = net.mesh +comm = uw.mpi.comm +levels = list(getattr(mesh, "_custom_mg_coarse_meshes", []) or []) + [mesh] + +def slice_edges(m, z=0.5): + grid = vis.mesh_to_pv_mesh(m) + return grid, grid.slice(normal="z", origin=(0, 0, z)) + +if comm.size == 1: + for k, lv in enumerate(levels): + n_cells = int(lv.dm.getHeightStratum(0)[1]) + n_v = int(np.diff(lv.dm.getDepthStratum(0))[0]) + print(f"[hier] level {k}: {n_cells} cells, {n_v} vertices", flush=True) + # each level on the slice, same camera + for k, lv in enumerate(levels): + grid, sl = slice_edges(lv) + pl = pv.Plotter(off_screen=True, window_size=(2400, 500)) + pl.set_background("white") + if k == len(levels) - 1: + band = np.asarray(mesh.cells_labelled("Band", 71)).astype(float) + grid.cell_data["band"] = band + sl = grid.slice(normal="z", origin=(0, 0, 0.5)) + pl.add_mesh(sl, scalars="band", cmap="RdBu_r", clim=(0, 1), show_edges=True, + edge_color="black", line_width=0.6, lighting=False, show_scalar_bar=False) + else: + pl.add_mesh(sl, color="white", show_edges=True, edge_color="black", line_width=0.6, lighting=False) + pl.view_xy(); pl.camera.parallel_projection = True + pl.camera.focal_point = (3.0, 0.5, 0.5); pl.camera.parallel_scale = 0.55 + out = f"{FIG}/hierarchy_level{k}_z05.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True) + # zoom on fault B's band, fine level: band cells red, fill/base white + grid, _ = slice_edges(mesh) + grid.cell_data["band"] = np.asarray(mesh.cells_labelled("Band", 71)).astype(float) + for tag, origin, normal, focal, scale, size in ( + ("z05", (0, 0, 0.5), "z", (2.4, 0.5, 0.5), 0.45, (1200, 1000)), + ("x24", (2.4, 0, 0), "x", (2.4, 0.5, 0.5), 0.45, (1000, 1000))): + sl = grid.slice(normal=normal, origin=origin) + pl = pv.Plotter(off_screen=True, window_size=size); pl.set_background("white") + pl.add_mesh(sl, scalars="band", cmap="RdBu_r", clim=(0, 1), show_edges=True, + edge_color="black", line_width=0.8, lighting=False, show_scalar_bar=False) + # the coarse level's edges on top, in blue, to show the nesting + _g0, sl0 = slice_edges(levels[0], z=0.5) if normal == "z" else (None, vis.mesh_to_pv_mesh(levels[0]).slice(normal="x", origin=origin)) + pl.add_mesh(sl0, color="white", opacity=0.0, show_edges=True, edge_color="blue", line_width=1.5, lighting=False) + if normal == "z": pl.view_xy() + else: pl.view_yz() + pl.camera.parallel_projection = True; pl.camera.focal_point = focal; pl.camera.parallel_scale = scale + out = f"{FIG}/hierarchy_zoomB_{tag}.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True) +else: + # partition: gather every rank's cells (P1 vertices) to rank 0 and colour by rank + from underworld3.utilities.line_cut import _coords + dm = mesh.dm; vS, vE = dm.getDepthStratum(0); X = _coords(dm)[: vE - vS] + cn = np.asarray(mesh._cell_node_indices(1, True)) + T = X[cn] # (n_cells, 4, 3) + band = np.asarray(mesh.cells_labelled("Band", 71)).astype(float) + allT = comm.gather(T, root=0); allB = comm.gather(band, root=0) + ranks = comm.gather(np.full(len(T), comm.rank, dtype=float), root=0) + if comm.rank == 0: + pts = np.vstack([t.reshape(-1, 3) for t in allT]) + n = pts.shape[0] // 4 + cells = np.hstack([np.full((n, 1), 4), np.arange(4 * n).reshape(n, 4)]).ravel() + grid = pv.UnstructuredGrid(cells, np.full(n, pv.CellType.TETRA), pts) + grid.cell_data["rank"] = np.concatenate(ranks); grid.cell_data["band"] = np.concatenate(allB) + sl = grid.slice(normal="z", origin=(0, 0, 0.5)) + pl = pv.Plotter(off_screen=True, window_size=(2400, 500)); pl.set_background("white") + pl.add_mesh(sl, scalars="rank", cmap="RdBu_r", clim=(0, comm.size - 1), show_edges=True, + edge_color="black", line_width=0.4, lighting=False, show_scalar_bar=False) + bsl = grid.threshold(0.5, scalars="band").slice(normal="z", origin=(0, 0, 0.5)) + if bsl.n_points: + pl.add_mesh(bsl, color="black", show_edges=True, edge_color="black", line_width=1.0, lighting=False) + pl.view_xy(); pl.camera.parallel_projection = True + pl.camera.focal_point = (3.0, 0.5, 0.5); pl.camera.parallel_scale = 0.55 + out = f"{FIG}/partition_np{comm.size}_z05.png"; pl.screenshot(out); pl.close(); print("[hier] wrote", out, flush=True) + print(f"[hier] np={comm.size} cells/rank {[len(t) for t in allT]} band/rank {[int(b.sum()) for b in allB]}", flush=True) From 8b1e217f303f428fc42b5e84d9bb69afe29a3e5a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 09:48:47 -0700 Subject: [PATCH 09/19] The design line after the realisations' requirements were separated: embed once, TI immediately, split locally with seam ligaments (#670) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Louis, 3 September: the embedding is the same for both realisations and only the split needs pairs on one rank; the weak and TI band need mesh continuity across a seam and nothing else. So: embed the band with its mid-surface reserved as conforming faces, cut the split only where the facet star is interior, and leave a user-set ligament at each seam crossing for the band's weak rheology to bridge — the hybrid recipe at seams, which the stepover and junction studies measured. Partition independence becomes a bound to measure. Build order: the seam-conforming fill in 2-D, the bridging split mode, the 2-D study. Underworld development team with AI support from Claude Code --- .../fault-parallel-placement-2026-09.md | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md index 998126616..01717a41f 100644 --- a/docs/developer/design/fault-parallel-placement-2026-09.md +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -166,6 +166,30 @@ None of that is built, and it should not be until a model needs a fault longer than one rank comfortably holds. The first thing to test when it is: the split across one cut plane at np=3. +**Revised on 3 September (Louis), after the hierarchy was seen and the +two realisations' requirements were separated.** The embedding is the +same for both realisations; only the split needs its pairs on one rank. +The weak and TI band need nothing beyond mesh continuity across a seam, +so with a seam-conforming fill a band can be cut anywhere, including on +a process boundary. The design line is therefore: embed the fault-zone +mesh with its mid-surface as reserved conforming faces and its cells +labelled, immediately usable by the TI formulation; then make the split +local by construction — cut only the mid-surface facets whose star is +interior to the rank and leave a ligament, a physical width the user +sets, at each seam crossing, which the band's weak rheology bridges. +That is the hybrid recipe applied at seams, and the stepover and +junction studies already measured it: a TI band bridging a cut gap +reproduces the continuous fault to 0.2%, and cuts plus a weak band cost +what the split costs and give its answer. What it gives up is +bit-identical partition independence, since the ligaments move with the +seams; that is a bound to measure at np=3, including a seam near a tip +or a junction. Build order: the seam-conforming fill in 2-D first, then +the bridging split mode (skip and report seam-touching facets rather +than refuse or redistribute), then the 2-D study of one-object cuts +against per-rank cuts with ligaments in both realisations. This +supersedes the CAD cut for the split; the per-region gather stays as +the mechanism for whole objects meanwhile. + Two further rulings from the same discussion: - **Do not rebalance the child after placement.** It fixes the solve From d9122286a00a6f1a8ed0beeb305e190b9b1f6f03 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 09:54:14 -0700 Subject: [PATCH 10/19] Interface sketch for seam-conforming placement and the bridging split (#670) The mechanism in four steps (interface planes chosen collectively, a strip alignment of cavity-adjacent cells to the plane's sides, an interface surface meshed once and broadcast, two fills and the existing rebuild extended with star-forest entries for the shared interface vertices), the signatures (place_thin_volume seams=, split_faults at_seams= and ligament=, FaultNetwork.build pass-through), the collective discipline, the tests in order (2-D serial, 2-D np=2 and 3, then 3-D), and the traps known in advance. For implementation in a fresh session. Underworld development team with AI support from Claude Code --- .../fault-parallel-placement-2026-09.md | 119 ++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/docs/developer/design/fault-parallel-placement-2026-09.md b/docs/developer/design/fault-parallel-placement-2026-09.md index 01717a41f..911da702b 100644 --- a/docs/developer/design/fault-parallel-placement-2026-09.md +++ b/docs/developer/design/fault-parallel-placement-2026-09.md @@ -365,3 +365,122 @@ The script beside this note (`fault_parallel_layouts.py`) regenerates the table: mpirun -np 3 python -u fault_parallel_layouts.py -uw_layout local|straddle|gathered [-uw_tail 0] + +## Interface sketch: seam-conforming placement and the bridging split + +Written 3 September 2026 for implementation in a fresh session; debug in +2-D serial and parallel, validate in 3-D serial and parallel. Nothing +below is built. + +### What changes and what does not + +The carve, the gmsh fill, the collective rebuild and every gate stay. +What changes is the contract on *where* a cavity may lie: today a +cavity must be interior to one rank (the gather makes it so); with seam +conformance a cavity may straddle a seam, and the two ranks fill their +own sides against an interface they both hold exactly. + +### The mechanism, in four steps + +1. **Interface planes.** For each seam a band's cavity would cross, a + plane (a line in 2-D) is chosen deterministically: through the mean + position of the seam faces inside the cavity zone, normal to the + band's strike there. Collective: every rank computes the same list. + A cavity zone touching no seam gets no plane and proceeds as now. +2. **The strip alignment.** Cells adjacent to the cavity zone are + reassigned by the side of the plane their centroid lies on, so that + within the zone the partition boundary *is* the plane, up to the + cells that straddle it, which are dropped. One shell partition, a + few cells wide along each plane. This is the only redistribution, + and it moves a strip, not a region. +3. **The interface surface.** On the lower rank of each pair, the + plane region bounded by the ring's outline (the ring faces meeting + the plane form a closed polygon) is meshed once in 2-D by gmsh with + the band's cross-section outline embedded, and broadcast to the + pair. In 2-D it is a chain of points along the line: the ring ends + and the band's two skin crossings. Both ranks embed it verbatim; + the fill's existing gating (zero moved constrained nodes, every + constraint facet present) covers it. +4. **Two fills and one rebuild.** Each rank carves its own side (the + band cells assigned by centroid side, the assembly cut by the plane + in CAD so the two halves' skins meet on the interface exactly) and + fills between its ring, its half of the skin and the interface. The + rebuild is the existing collective one with one addition: the + interface's new vertices are created on both ranks and must be the + *same* points, so `_attach_uninterp_vertex_sf` gains entries for + them — owner the lower rank, leaves on the higher, matched by the + interface's own point numbering, which both ranks hold. That is the + one deep change: today's rebuild asserts that no placed point is + ever shared. + +### Signatures + +``` +place_thin_volume(dm, patches, width, ..., seams="gather") + seams : {"gather", "conform"} + "gather" — the per-region gather (today). + "conform" — cavities may straddle seams; the four steps above. + info["seam_crossings"] : list of dicts, one per interface — + {"ranks": (r, s), "point": ..., "normal": ..., "n_strip": int} + +split_faults(mesh, names, groups=None, at_seams="gather", ligament=None) + at_seams : {"gather", "bridge"} + "gather" — redistribute per group, cut everything (today). + "bridge" — cut only facets whose closure is unshared AND farther + than ligament/2 from every interface plane; leave the + rest as ligaments for the band's rheology. + ligament : physical width (same units as the band width); required + with "bridge". + info / the child's record: "ligaments" — one entry per seam crossing + with the uncut facet count and the ligament's extent. + +FaultNetwork.build(..., seams="gather", seam_ligament=None) + passes both through; net.info gains "seam_crossings" and + "ligaments"; the band's rheology (apply / ti_fields) must cover the + ligaments, which it does already since the band is weak everywhere + it is not cut. +``` + +### Collective discipline + +Every decision above is taken from gathered data before any branch: +the plane list, the strip assignment, the interface point set, the +ligament list. A rank with no crossing still participates in each +collective. The refusals stay collective. Arm `UW_HANG_WATCHDOG` on +every run. + +### Tests, in order + +- **2-D serial**: `seams="conform"` on a single-rank mesh must reduce + to the existing path bit-for-bit (no seams, no planes). +- **2-D np=2, one band across the seam**: `conform` against `gather` + — same zone and skin counts, same volume and Euler gates, the TI + solve identical to the fill's noise; then `at_seams="bridge"` — one + ligament reported, its width the one asked for, slip along the fault + within the stepover bound (0.2% of the continuous control away from + the ligament). +- **2-D np=3**, a seam near a fault tip and one at a junction: the + bound on the ligament's effect where it is not a small perturbation. +- **3-D serial** (reduction to the existing path), then **3-D np=2 and + np=3** on the crossing fixture and the layout fixture: the same + assertions, plus the interface surface's own gates (the plane + triangulation present verbatim in both fills). +- The layout throughput test gains a `conform` column: cells moved + should be the strips only. + +### Traps known in advance + +- The rebuild currently *asserts* no placed point is shared; the + interface vertices violate that by design, and the SF extension must + come before the interpolate (the #520 ordering). +- gmsh's fill is not bit-reproducible across inputs that differ only + in ordering (the 611/612 placed-node effect); the interface surface + must be generated once and broadcast, never regenerated per rank. +- A plane chosen from the seam's mean position can cut a band near + its tip; the ligament rule must handle a crossing at a tip (the tip + itself is never split, so a seam through the tip region simply + enlarges the uncut margin). +- The 2-D thin volume has the one-region gather only; the per-region + gather (`_gather_regions`) is wired into the 3-D thin volume alone. + Wire it into 2-D first so `gather` and `conform` are compared like + for like. From 82b413e4f18b440d2a2802888a6ae268f207a7f8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 10:44:07 -0700 Subject: [PATCH 11/19] The seam ligament: place a 2-D band across a partition seam without gathering, and split it rank-local (#670) Partition-crossing structures stay transversely isotropic (the ruling of 3 September): the mesh only has to conform on each side of the seam, and the band's weak plane is the glue across it. place_thin_volume gains seams="ligament" (2-D): every rank carves its own cavities, each stops one cell short of the seam, the assembly is clipped to what the cavity holds and made manifold, and the base cells the clipped band covered are left as the ligament, labelled zone and