Generalized FMG hierarchy: custom-P transfers for arbitrary nested/non-nested grids (Layer 1) by lmoresi · Pull Request #290 · underworldcode/underworld3

lmoresi · 2026-06-28T20:37:24Z

What

Generalize UW3's geometric FMG so the multigrid hierarchy can use arbitrary coarse grids (can be nested or non-nested, uniform or locally (SBR) refined) by building the prolongations in house (barycentric / local-RBF) and installing them via PC.setMGInterpolation, with Galerkin RAP forming the coarse operators. This decouples geometric multigrid from uniform refine() nesting, enabling locally-adapted hierarchies while keeping FMG.

Independent of any adapter ("adapt-on-top" is a separate follow-on). The existing Mesh(refinement=N) uniform FMG path is unchanged.

Status — Layer 1, Phase-1 hardening complete

Home: utilities/custom_mg.py. Entry point: set_custom_fmg(solver, coarse_meshes, builder=, field_id=).

Working + tested — serial test_1014/1015/1016/1017 (20) + parallel tests/parallel/test_1017_custom_mg_parallel_mpi.py (np=2):

CustomMGHierarchy, set_custom_fmg, sbr_refine / sbr_refine_where (local Skeleton-Based Refinement; no MMG; on-rank; conforming).
BC applied at every level (transfers reduced→reduced) with a zero-column guard — the load-bearing invariant (PETSc native FMG does this via constrained sections; custom-P must replicate it).
Scalar / single-field-vector (top-level PC) and Stokes velocity block (field_id=0).
Serial and parallel (np>1), nested co-partitioned.

Hardening steps (all done)

Stokes / saddle-point velocity block. The velocity sub-PC is unreachable until the monolithic Jacobian is assembled (PCFieldSplit forms A_vv via MatCreateSubMatrix; snes.setUp builds structure only). So _install_velocity_block_transfers forces a Jacobian assembly (computeFunction+computeJacobian at the zero guess; max_it=0 fallback), reaches the velocity sub-PC, resets it and rebuilds a fresh PCMG from our P (mirrors the proven standalone recipe; sidesteps the MatProductReplaceMats live-swap bug), and re-attaches the coupled Stokes nullspace. Wired into SNES_Stokes_SaddlePt.solve as a guarded no-op.
Dropped the throwaway-solver factory. Each coarse level's BC-constrained reduced map is derived directly from its DM via copyFields + copyDS + createDS (_coarse_reduced_map) — the DS carries UW's exact essential-BC definitions and is topology-independent, so it constrains any coarse mesh sharing the solver's boundary labels. Validated byte-identical to the old factory path; leak-free (no SNES / JIT). set_custom_fmg / build no longer take a level_solver_factory.
Parallel (np>1). Nested co-partitioned transfers: each rank builds its block of P rank-locally (ghost-inclusive coarse coords → every owned fine node lands in a local coarse simplex), reduced global numbering rides the DM global section (_level_dof_layout), transfers assemble as MPIAIJ (constrained coarse DOFs drop → reduced→reduced), parallel zero-column guard via Pᵀ·1. The legacy finest-only set_custom_mg path stays serial-only and raises loudly at np>1.

Validation

	iters	vs GAMG	notes
Scalar jump-coeff Poisson, 5-level (3 uniform + 2 SBR)	3 FMG	46
SolCx velocity block (η-jump 1e6), 3-level	6	~198	matches GAMG solution to 1.5e-9
SolCx velocity block vs native FMG	6	(native 5)	ties native FMG; value = works where native FMG cannot (non-nested / adapted)
Parallel np=1/2/4 (Poisson + Stokes)	4 / 6		matches GAMG and each other across rank counts

Out of scope (later)

Non-nested custom-P in parallel (cross-rank point location) — serial-only. Layer-2 adapt-on-top adapter is a separate follow-on.

Design: docs/developer/design/GENERALIZED_FMG_HIERARCHY_AND_ADAPT.md

Underworld development team with AI support from Claude Code

Add set_custom_mg() + utilities/custom_mg.py: drive PCMG geometric multigrid with a prolongation we build ourselves (barycentric or RBF) from independent, possibly non-nested coarse meshes. Coarse operators come from Galerkin RAP (P^T A P), so only the prolongation list is supplied. This decouples geometric multigrid from a nested refine() hierarchy. Mechanism: inject the custom P before the first PCSetUp (avoids the MatProductReplaceMats shape error a live matrix swap triggers under Galerkin), and KSPSetDMActive(OPERATOR, False) so PETSc uses our explicit P instead of re-deriving DM-hierarchy interpolation. Finest level reduced to the BC-eliminated global ordering; coarse levels use full DOF coords. Validated (scalar Poisson, refinement=2 box): custom barycentric P from independent coarse meshes reaches FMG iteration counts (1-3 vs FMG 2) and the same solution. On a mesh with NO refine hierarchy (FMG unavailable -> GAMG, 13 iters) custom-P geometric MG converges in ~9 iters. test_1015 (4 tests); test_1014 (11) still pass. Single-field (scalar/vector) only; Stokes velocity-block path is a follow-up. Underworld development team with AI support from Claude Code

…coped SBR Phase 1 of hardening the custom-prolongation work into a generalized FMG hierarchy. - CustomMGHierarchy: adapter-agnostic multi-level hierarchy that builds custom-P transfers (barycentric/rbf) with essential BCs applied at EVERY level (transfers map reduced->reduced). This is the load-bearing invariant: on an exactly-nested hierarchy, omitting per-level reduction makes a coarse boundary DOF coincide with a BC-removed fine DOF -> zero column -> singular Galerkin coarse operator. build() guards against zero columns. - set_custom_fmg(): register a hierarchy + per-level solver factory; built and installed at solve() time (build-time injection, DMActive(OPERATOR,False), Galerkin RAP). Scalar / single-field-vector (top-level PC); Stokes velocity block is Phase 2. - sbr_refine / sbr_refine_where: local Skeleton-Based Refinement (no MMG, on-rank, conforming) with SCOPED dm_plex_transform_type (leaking it globally breaks UW's uniform refine() with err73). - Legacy finest-only path kept for back-compat (test_1015). Validated: scalar jump-coefficient Poisson on a 5-level (3 uniform + 2 SBR) hierarchy converges in 3 FMG iters vs GAMG 46, solution matches to 2.4e-8. test_1016 (3 tests); test_1015 (4) + test_1014 (11) still pass. Design: docs/developer/design/GENERALIZED_FMG_HIERARCHY_AND_ADAPT.md Underworld development team with AI support from Claude Code

Layer 1 (generalized FMG hierarchy) is an independent, working capability (arbitrary nested/non-nested coarse grids -> geometric FMG, BC-per-level). Current scope experimental: scalar/single-field, serial, factory API. Remaining before merge-ready general feature: Stokes velocity-block, drop the factory (label-based BC reduction), parallel. Underworld development team with AI support from Claude Code

Integrate the generalized FMG hierarchy into the saddle-point solver's velocity sub-block. set_custom_fmg(..., field_id=0) now drives geometric multigrid on the velocity block with our supplied (barycentric / RBF) prolongations + Galerkin RAP. The velocity sub-PC is unreachable until the monolithic Jacobian is assembled (PCFieldSplit forms A_vv via MatCreateSubMatrix; snes.setUp builds structure only -> err73). So _install_velocity_block_transfers: forces a Jacobian assembly (computeFunction + computeJacobian at the zero guess; throwaway max_it=0 fallback), reaches the velocity sub-PC, reset()s it and rebuilds a FRESH PCMG from our P (mirrors the proven standalone recipe; sidesteps the MatProductReplaceMats live-swap bug), re-attaches the coupled Stokes nullspace. _configure_pcmg now derives the options prefix from pc.getOptionsPrefix() so both the scalar top-level PC and the velocity sub-PC are configured correctly. inject_custom_mg is wired into SNES_Stokes_SaddlePt.solve (guarded no-op unless set_custom_fmg was called), injected after setFromOptions / nullspace, before the real solve. Validated (SolCx eta_B=1e6, 3-level nested hierarchy, in-solver via solver.solve()): velocity block converges in 6 MG iters vs GAMG 198, solution matches the GAMG reference to 1.5e-9. test_1017 (barycentric + rbf). test_1014/1015/1016 unchanged (20 pass total). Underworld development team with AI support from Claude Code

…, Step 1) Underworld development team with AI support from Claude Code

Each coarse level's BC-constrained reduced map is now derived directly from the coarse mesh DM via _coarse_reduced_map: clone the DM, copy the (built) finest solver's fields + DS onto it, createDS, read the global section. The DS carries UW's exact essential-BC definitions and is a topology-independent discretisation spec, so it constrains the matching boundary DOFs on ANY coarse mesh that shares the solver's boundary labels (nested or not). This removes the throwaway-solver factory (heavy, leak-prone, and an API wart): set_custom_fmg and CustomMGHierarchy.build no longer take a level_solver_factory. Leak-free — DM ops only, no SNES / JIT. Validated identical to the old factory path (reduced maps byte-equal: 126 + 510 DOFs on the SolCx 3-level hierarchy). test_1014/1015/1016/1017 all green (20). custom-P vs native FMG unchanged (6 vs 5 iters, sol match 3.7e-10). Underworld development team with AI support from Claude Code

Underworld development team with AI support from Claude Code

…tial) The reduced maps use rank-local DOF indices and the prolongations assemble as serial AIJ, so at np>1 custom-P would silently build wrong P (or hit a cryptic broadcast error). Add _require_serial: set_custom_fmg and inject_custom_mg now raise a clear NotImplementedError on >1 ranks, pointing to preconditioner='fmg'/'gamg'. Serial behaviour unchanged. Parallel test tests/parallel/test_1017_custom_mg_serial_guard_mpi.py (np=2) asserts both the set-time and legacy solve-time guards fire. Full parallel custom-P (nested co-partitioned, rank-local point location + ghosting, MPIAIJ assembly with global-section reduction) remains a designed fast-follow. Underworld development team with AI support from Claude Code

Underworld development team with AI support from Claude Code

…(Phase 1, Step 3) The hierarchy path now builds parallel-correct transfers. Each rank builds its block of P rank-locally: ghost-inclusive coarse coords mean every owned fine node lands in a local coarse simplex (verified 0 misses np=2/4), and the reduced global numbering rides the DM global section via a clean local->global-reduced map (_level_dof_layout scatters owned global indices out through globalToLocal; constrained DOFs stay -1, ghosts resolve to the owner's global index). Transfers assemble as MPIAIJ (owned fine rows, global coarse cols incl. off-rank); constrained coarse DOFs drop -> reduced->reduced. Parallel zero-column guard via P^T·1 + allreduce. _coarse_dof_layout reuses the leak-free copyDS trick for coarse levels. CustomMGHierarchy.build branches serial (scipy CSR) vs parallel (MPIAIJ). The guard now blocks only the legacy finest-only path (still serial). Validated np=1/2/4: scalar Poisson custom-P 4 iters and Stokes SolCx velocity block 6 iters, both matching the GAMG reference (and each other across rank counts). New tests/parallel/test_1017_custom_mg_parallel_mpi.py (scalar + Stokes correctness + legacy serial-guard). Serial 20 unchanged. Underworld development team with AI support from Claude Code

Underworld development team with AI support from Claude Code

Copilot

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

…+ Constrained tests) Verify custom-P works for every solver family that consumes the mesh, so a mesh-owned adapted hierarchy drives MG uniformly: - SNES_Vector.solve was MISSING the inject_custom_mg hook (only SNES_Scalar and SNES_Stokes_SaddlePt had it) — so Vector_Projection silently stayed on GAMG. Add the hook (mirrors the scalar one). All 8 solver-subclass solve() overrides delegate to a base solve(), so this completes coverage: scalar branch -> Poisson/Darcy/Projection/AdvDiff/Diffusion; vector branch -> Vector_Projection; velocity-block -> Stokes/VE/Constrained/NS. - test_1017: add SNES_Vector (Vector_Projection, top-level vector PC, field_id=None) and Stokes_Constrained (free-slip multipliers, grouped [p,h] split, field_id=0) serial tests — both drive custom-P 'mg', constrained matches analytic SolCx to 4.4e-4. - parallel test: add Vector (passes np=2) and Constrained. The Constrained parallel test is SKIPPED: Stokes_Constrained is not parallel-safe yet — it segfaults at np>1 independently of custom-P (canonical test_1062_constrained_solcx also segfaults at np=2 under plain GAMG). It auto-enables when the constrained solver becomes parallel-ready. Serial: custom_mg 22 + projections + constrained green. Parallel np=2: scalar + Stokes-velocity + vector + legacy-guard pass. Underworld development team with AI support from Claude Code

… limitation Underworld development team with AI support from Claude Code

Underworld development team with AI support from Claude Code

…k fieldsplit The Stokes velocity-block injection forces computeJacobian, then reaches the fieldsplit BEFORE the SNES solve. SNESSolve normally wires the freshly-assembled Jacobian into the KSP/outer-PC; doing it ourselves was missing, so for some configurations (e.g. Stokes on an SBR-refined child mesh with a mesh-variable bodyforce) the outer PC carried an UNASSEMBLED operator and PCSetUp failed with "Matrix must be set first" (err73). Fix: snes.getKSP().setOperators(J, Pmat) after forcing assembly, before outer_pc.setUp(). Found via the Layer-2 adapt-on-top prototype (Stokes on a fault-refined child). Diagnosed: ksp.A set but asm=False while J.assembled=True -> operator not wired. New regression test_custom_fmg_stokes_on_sbr_child (Stokes velocity-block custom-P on an SBR child with a mesh-variable bodyforce) — fails pre-fix, passes after. test_1014/1015/1016/1017 green (serial); parallel np2 unchanged. Underworld development team with AI support from Claude Code

Copilot

Pull request overview

Copilot reviewed 7 out of 7 changed files in this pull request and generated 4 comments.

+    csr = csr.tocsr()
+    M = PETSc.Mat().createAIJ(
+        size=csr.shape,
+        csr=(csr.indptr.astype("int32"), csr.indices.astype("int32"), csr.data),
+    )


+    M = np.block([[Acc, Pc], [Pc.T, np.zeros((dim + 1, dim + 1))]])
+    Minv = np.linalg.inv(M)
+    B = np.hstack([phi(cdist(fine_coords, coarse_coords)),
+                   np.ones((fine_coords.shape[0], 1)), fine_coords])
+    Praw = (B @ Minv)[:, :nc]


+    # geometric MG crushes GAMG on the eta-jump velocity block
+    assert vksp.getIterationNumber() <= 15
+    assert vksp.getIterationNumber() < iters_g


+    assert vksp.getPC().getType() == "mg"
+    assert vksp.getPC().getMGLevels() == len(coarse) + 1
+    assert vksp.getIterationNumber() <= 25
+    assert sol.velocity_error(s.u) < 1e-2


…olve not inv Applied on this branch (feature/rotated-freeslip-bc, stacked on #290) so the fixes flow through with #293; the #290 review is answered with this reference. - _to_petsc_aij: cast CSR indptr/indices to PETSc.IntType instead of a hard int32. On a 64-bit PetscInt build a large mesh can exceed int32 range → overflow / mis- addressed entries. - rbf_prolongation: compute the prolongation via np.linalg.solve(M, B.T).T rather than forming np.linalg.inv(M) explicitly (faster, more numerically stable; M is symmetric so B M^-1 = solve(M, B^T)^T). The two test_1017 iteration-count bounds Copilot flagged are kept: the margins are large (custom MG ~6 iters vs GAMG ~198) and pair with the "beats GAMG" assertion, so they are low flake-risk and preserve a useful regression signal. Verified: test_1017 serial 5/5, test_1017 parallel np=2 (4 passed, 1 skipped=#291), test_1018 rotated 6/6. Underworld development team with AI support from Claude Code

lmoresi · 2026-07-02T01:26:23Z

Thanks — Copilot's 4 comments are addressed on the stacked branch feature/rotated-freeslip-bc (#293), commit 3dc7b893, since that branch carries custom_mg.py and is where the follow-on work lives:

custom_mg.py:179 (int32 CSR indices) — fixed: _to_petsc_aij now casts indptr/indices to PETSc.IntType, matching the PETSc build's integer width (avoids int32 overflow / mis-addressing on a 64-bit PetscInt build with large meshes).
custom_mg.py:97 (explicit inverse) — fixed: rbf_prolongation now uses np.linalg.solve(M, B.T).T instead of forming np.linalg.inv(M) (faster, more stable; M is symmetric so B M⁻¹ = solve(M, Bᵀ)ᵀ).
test_1017:91 and test_1017:110 (absolute iteration bounds) — kept intentionally. The margins are large (custom MG ~6 iters vs GAMG ~198) and each bound is paired with the "custom MG beats GAMG" assertion, so they are low flake-risk and preserve a useful regression signal. For tests, these tight margin cases are actually better because they will trigger if (say) gamg improves and that's sometimes what we want to know.

Verified after the change: test_1017 serial 5/5, test_1017 parallel np=2 (4 passed, 1 skipped = the unrelated #291 constrained-parallel), and the dependent test_1018 rotated 6/6.

…not inv - _to_petsc_aij: cast CSR indptr/indices to PETSc.IntType instead of a hard int32. On a 64-bit PetscInt build a large mesh can exceed int32 range → overflow / mis- addressed entries. - rbf_prolongation: compute the prolongation via np.linalg.solve(M, B.T).T rather than forming np.linalg.inv(M) explicitly (faster, more numerically stable; M is symmetric so B M^-1 = solve(M, B^T)^T). (The two test_1017 iteration-count bounds Copilot flagged are kept intentionally: the margins are large — custom MG ~6 iters vs GAMG ~198 — and pair with the "beats GAMG" assertion, so they are low flake-risk and preserve a useful regression signal.) Underworld development team with AI support from Claude Code

lmoresi · 2026-07-02T01:28:51Z

Follow-up: these two fixes are now applied directly on this branch (feature/custom-mg-prolongation, commit b7d1c4fd), not only via the stacked #293 — so they merge with #290 itself. (The identical patch on #293 will drop as a no-op when #293 rebases onto the merged #290.) Validated identically: test_1017 serial 5/5, parallel np=2 (4 passed, 1 skipped=#291), test_1018 6/6.

Copilot review Rebased onto development (which now has #290 + the merged boundary_flux primitive, #294). The rotated σ_nn recovery is a rotated-frame reading of the same Consistent Boundary Flux, so collapse the duplication and answer the Copilot review: - Copilot pyx:7597 (nonlinear/preamble): the rotated path runs the standard pre-solve preamble (DM time, auxiliary vector, _update_constants) before delegating, and now GUARDS the unsupported cases — it is a single LINEAR solve, so picard!=0 / warm-start raise NotImplementedError rather than silently returning a single-linearisation answer. - Copilot rotated_bc:65 & :577 (boundary label): use the consolidated "UW_Boundaries" label via the shared boundary_flux._boundary_stratum_is (survives mesh adaptation; clear error for unknown names) — the per-boundary label lookups are gone. - Unification (my own review note): delete rotated_bc._recover_sigma_nn_2d and reuse boundary_flux._desmear (σ_nn = de-smear of −reaction, mean-removed); the field hand-off reuses boundary_flux.write_boundary_scalar_field. Net −133 lines in rotated_bc. _desmear gains partial_reaction: rotated passes False (its reaction is the ASSEMBLED operator Q(A·u−b), complete at every node → OVERWRITE across ranks) vs boundary_flux's default True (raw per-rank residual → SUM). This is the fix for the np=4 double-count. Verified: test_1018 serial 6/6; test_1064 parallel 5/5 at -n 2 and -n 4. Underworld development team with AI support from Claude Code

…e-surface hand-off (#293) * feat(stokes): rotated strong free-slip BC with σ_nn reaction (increment 1) Add `add_rotated_freeslip_bc(boundary, normal=None)` on the Stokes saddle: rotate the boundary velocity DOFs into a per-node (normal, tangential) frame and impose the rotated normal component as an exact Dirichlet constraint (v·n̂=0 to machine precision), with the constraint reaction giving the consistent boundary normal traction σ_nn (`boundary_normal_traction`) — no augmented-Lagrangian splitting. New module `utilities/rotated_bc.py` (pure Python; solve() delegates when rotated BCs are registered, no-op otherwise): - dimension-general SVD constraint frame: a node's accumulated face normals span a rank-r subspace; the r rows spanning it are constrained, the (dim-r) tangential complement is free. Handles 2D face/corner, 3D face/edge/corner, and regional spherical boxes (curved caps + planar sides + their edges/corners). - per-boundary normal source: geometric facet normal (computeCellGeometryFVM, 2D+3D) or an analytic sympy normal (exact radial X/|X| on caps) or a constant. - rigid-rotation gauge auto-removed only when it is a genuine null space of the constrained problem (circular/spherical free-slip), never on straight walls. Validated (test_1018): box rotated free-slip on 4 walls reproduces native essential free-slip to 4e-6; annulus per-node radial free-slip gives machine-zero radial leakage with the gauge removed. Direct LU solve here; geometric-FMG path (the two custom_mg velocity-block injections) is the next increment. Also fixes the stale set_custom_mg docstring (documents the Stokes velocity-block path, which is wired via inject_custom_mg). Underworld development team with AI support from Claude Code * feat(stokes): rotated free-slip solved by geometric FMG (iterative by default) Replace the direct-LU rotated solve with a self-contained fieldsplit-Schur KSP on the rotated operator (LU is almost never the right option; kept only behind solver._rotated_use_lu). A plain rotated Mat carries no DM field info, so UW3's DM-coupled fieldsplit cannot split it — the split is built from explicit velocity/pressure index sets, and the velocity sub-PC is driven by: - geometric FMG on the custom (PR#290) prolongation, rotated P̂ = Q_v·P via setMGInterpolation (needs no DM), when a hierarchy is registered (set_custom_fmg); the rotated block A_vv = Q_v A_vv Q_vᵀ is formed from the rotated operator automatically and only the FINE prolongation is rotated (Galerkin coarse operators auto-correct); - GAMG otherwise. selfp Schur + jacobi pressure PC → the SolCx velocity block converges in ~6 outer iterations. The coupled Stokes null space (constant pressure ⊕ rotated rigid rotation) is attached to the rotated operator. Two fixes so the strong BC is EXACT under an iterative solve (independent of tolerance): - zero the RHS at constrained rows (zeroRowsColumns leaves b untouched, so a nonzero b there leaks straight into the solution); - zero the (fully decoupled) constrained DOFs post-solve — an identity constraint row only drives its residual below tolerance in a Krylov solve, so the DOF is ~tol, not 0, until pinned exactly here. Null-space vectors are also zeroed at the constrained rows for compatibility. test_1018 gains a geometric-FMG velocity-block case; annulus leakage is now machine-zero at the default (loose) tolerance. No regressions in the custom-mg Stokes suite. Underworld development team with AI support from Claude Code * fix(rotated-bc): parallel-safety pass (distributed Q, collective-safe, guards) Move the rotated free-slip construction off serial-only constructs so it no longer segfaults at np>1: - build_rotation now assembles a DISTRIBUTED PETSc Mat with the operator's row layout, each rank setting only its owned rows (the per-node dim×dim block is node-local — a node's velocity components live on one DMPlex point owned by one rank — so Q is block-diagonal with no off-rank columns), instead of a global scipy matrix replicated on every rank; - assemble the Jacobian before build_rotation so its parallel layout is final; - post-solve zeroing of the constrained DOFs and the null-space-vector zeroing use ownership-relative local indices, not global indices into a local array; - the rigid-rotation gauge is removed on the GLOBAL vector with PETSc dots (a local nodal sum double-counts shared nodes); new _rigid_rotation_global helper; - _rotation_is_nullspace uses collective PETSc norms and no longer early-returns on a per-rank empty normal_rows (that desynced the collective and could deadlock); - guard dm.getStratumIS against a NULL IS on ranks that own no part of a boundary (the actual np>1 segfault: getIndices() on a null IS); - evaluate an analytic (sympy) boundary normal via a single lambdify, not per-node .subs() (orders of magnitude faster; also removes a per-rank serialisation). Serial (np=1) behaviour unchanged: test_1018 still passes. The core rotated solve runs correctly and fast at np=2 (verified by staged bisection). A hang in the full solve() wrapper at np>1 is still under investigation and tracked separately. Underworld development team with AI support from Claude Code * fix(rotated-bc): parallel np>1 wrapper "hang" was a global-into-local index crash The rotated free-slip s.solve() appeared to hang at np=2/4 (one rank ~100% CPU, the other idle). It was not a hang: the RHS constrained-row zeroing indexed the LOCAL vector slice with GLOBAL row indices, ba = bhat.getArray(); ba[normal_rows] = 0.0 # normal_rows are GLOBAL which is in-range on rank 0 (ownership starts at 0, so global == local) but overflows on every other rank (ownership starts at rstart > 0). Rank 1 therefore raised an IndexError and left the collective while rank 0 proceeded into the iterative solve and blocked in the next collective (getGlobalVec) — the classic asymmetric-crash-that-looks-like-a-hang. This is the same global-into-local bug class already fixed for the post-solve zeroing; line 222 (and the LU-path pin write) were missed. Fixes, in rotated_bc.py: - solve_rotated_freeslip: zero bhat at the constrained rows using ownership- relative local indices (g - rstart), matching the post-solve zeroing; same for the opt-in LU pressure-pin write. - _solve_rotated_iterative: give each solve a UNIQUE PETSc options prefix and drop the keys afterwards, so sequential rotated solves (multiple solvers / time steps) do not share or accumulate global-options state. - solve_rotated_freeslip: after the field scatter, refresh the enhanced-variable gvec cache and clear canonical-data / mark the mesh lvec stale, mirroring the normal solve so downstream reads (var.data/array, checkpoint, stats) aren't stale. Validated np={1,2,4}: box velocity L2 matches serial to ~1e-10, annulus velocity L2 to ~1e-8 (iterative-solver tolerance), radial leakage byte-identical. New tests/parallel/test_1064_rotated_freeslip_parallel.py (box + annulus, MPI-safe Integral/BdIntegral diagnostics) passes at -n 2 and -n 4; serial test_1018 still 3/3. Underworld development team with AI support from Claude Code * test(rotated-bc): cover custom-FMG rotated annulus free-slip in parallel Add test_rotated_freeslip_annulus_fmg_partition_independent to the parallel suite: rotated radial free-slip on a nested annulus hierarchy (coarse -> refine -> refine) whose velocity block is CUSTOM GEOMETRIC FMG via set_custom_fmg (no GAMG, no direct solve). Asserts the solve converges and that velocity L2 + radial leakage reproduce the serial reference in parallel. Verified: 3 passed at -n 2 and -n 4. This closes the coverage gap on the full stack FMG x rotated x annulus x np>1 (previously only GAMG rotated box/annulus and serial FMG box were tested). Underworld development team with AI support from Claude Code * feat(rotated-bc): parallel-safe, corner-correct sigma_nn (CBF topography) boundary_normal_traction was serial-only and rang at corners. Two changes: 1. Corner-correct recovery. Read the CARTESIAN nodal reaction r_c = A u - b and use R_i = n_hat_i . r_c(node_i) as the sigma_nn nodal load, instead of the rotated frame's "normal row". At a node shared with another rotated-free-slip boundary the rotated frame's first SVD row is a MIX of both walls' normals, so reading it gave a wrong stress component and large corner spikes; n_hat . r_c is the true normal traction for this boundary. Whole-boundary SolCx relL2 vs analytic sigma_yy drops from ~0.43 to 0.042 at res48 (interior unchanged ~0.046), converging to 0.030 at res96, corr 0.999 -- matching the standalone strong-Dirichlet reaction method. 2. Parallel safety. The reaction is scattered to a local vector (ghosts included) and read by LOCAL section offset (no global-into-local overflow -- the same crash class as the solve-path hang). The consistent P2 line mass is assembled GLOBALLY by a coordinate-keyed allgather of the boundary elements (shared facets de-duplicated), so every rank solves the identical 1D system and the mean-removal gauge is global; each rank returns sigma_nn at its own local nodes. Verified partition-independent: res48 relL2/corr byte-identical at np=1/2/4. Underworld development team with AI support from Claude Code * test(rotated-bc): cover sigma_nn recovery accuracy (serial) + partition-independence - test_1018: sigma_nn from boundary_normal_traction on the SolCx top boundary matches the analytic sigma_yy (whole boundary relL2 < 0.08, corr > 0.99) -- guards the Cartesian-reaction + n_hat projection (corner-correct) + consistent P2 line mass. - test_1064: the same recovery is partition-independent -- whole-boundary relL2/|corr| (gathered + de-duplicated + broadcast) match the serial reference at np=2 and np=4, and stay accurate. Guards the parallel reaction read + allgather'd consistent mass. Verified: test_1018 4 passed serial; test_1064 4 passed at -n 2 and -n 4. Underworld development team with AI support from Claude Code * feat(rotated-bc): lumped-mass sigma_nn de-smear (monotone, no Gibbs overshoot) The consistent P2 line mass overshoots where the boundary traction jumps (e.g. across a viscosity contrast) — a Gibbs wiggle that, for a free surface, injects a spurious surface-velocity pulse. The lever is mass LUMPING, not element order: the lumped (diagonal, row-sum) boundary mass is an M-matrix, so its de-smear cannot overshoot. boundary_normal_traction gains mass={"lumped"(default),"consistent"}: - "lumped": sigma = -R / m_lumped (m_lumped = row sums h*[1/6,2/3,1/6] of the P2 line mass). Monotone at discontinuities, a purely local division (no global mass solve), and marginally MORE accurate on SolCx than consistent (whole-boundary relL2 0.0399 vs 0.0419 at res48, converging 0.0555 -> 0.0399 -> 0.0284 at res 24/48/96). - "consistent": the previous full-mass solve, kept for smooth tractions. Prototype (P1-consistent/P1-lumped/P2-lumped) showed P2-lumped is best: it keeps the P2 node density AND is monotone. Assembled from the same coordinate-keyed element allgather, so both modes stay partition-independent (byte-identical np=1/2/4). Tests: test_1018 adds a total-variation no-overshoot guard (lumped TV ~ analytic TV, and < consistent TV); test_1064 golden updated to the lumped default. test_1018 5/5 serial; test_1064 4/4 at -n 2 and -n 4. Underworld development team with AI support from Claude Code * feat(rotated-bc): dynamic_topography surface field — free-surface hand-off Expose the rotated-free-slip sigma_nn to the free-surface machinery as a MeshVariable. solver.dynamic_topography(boundary, field, buoyancy_scale=1, mass="lumped") writes h = -(sigma_nn - mean)/(rho g) onto a scalar surface field (P1 recommended) at the boundary nodes, from the last solve's constraint reaction. The 3-number topography integrator drives node motion from a surface field, so this is the form it consumes: the field is usable symbolically (BdIntegral) and the interior nodes are untouched. Parallel gotcha fixed: the field must be written ONCE from a local numpy copy, not per-node. A per-element write to var.data fires the variable write-callback each time, and the number of boundary nodes differs per rank (a rank may own none of the boundary), so per-element writes desync the callback's collective and deadlock. Verified: the field BdIntegral over the boundary is byte-identical at np=1/2/4. Tests: test_1018 (field reproduces analytic SolCx topography at the top vertices, corr>0.99, and is BdIntegral-usable); test_1064 (BdIntegral L2 partition-independent). test_1018 6/6 serial; test_1064 5/5 at -n 2 and -n 4. Underworld development team with AI support from Claude Code * refactor(rotated-bc): unify σ_nn onto shared boundary_flux; address #293 Copilot review Rebased onto development (which now has #290 + the merged boundary_flux primitive, #294). The rotated σ_nn recovery is a rotated-frame reading of the same Consistent Boundary Flux, so collapse the duplication and answer the Copilot review: - Copilot pyx:7597 (nonlinear/preamble): the rotated path runs the standard pre-solve preamble (DM time, auxiliary vector, _update_constants) before delegating, and now GUARDS the unsupported cases — it is a single LINEAR solve, so picard!=0 / warm-start raise NotImplementedError rather than silently returning a single-linearisation answer. - Copilot rotated_bc:65 & :577 (boundary label): use the consolidated "UW_Boundaries" label via the shared boundary_flux._boundary_stratum_is (survives mesh adaptation; clear error for unknown names) — the per-boundary label lookups are gone. - Unification (my own review note): delete rotated_bc._recover_sigma_nn_2d and reuse boundary_flux._desmear (σ_nn = de-smear of −reaction, mean-removed); the field hand-off reuses boundary_flux.write_boundary_scalar_field. Net −133 lines in rotated_bc. _desmear gains partial_reaction: rotated passes False (its reaction is the ASSEMBLED operator Q(A·u−b), complete at every node → OVERWRITE across ranks) vs boundary_flux's default True (raw per-rank residual → SUM). This is the fix for the np=4 double-count. Verified: test_1018 serial 6/6; test_1064 parallel 5/5 at -n 2 and -n 4. Underworld development team with AI support from Claude Code * fix(rotated-bc): address #293 re-review — DM vector-pool leaks + stale docstring Copilot re-review (7 comments) — resource-management + docs: - DM getGlobalVec() pool leaks: vectors that PERSIST beyond solve_rotated_freeslip (returned in the info dict: U, the LU-branch Uhat; the pressure null-space vector pv) now use dm.createGlobalVec(); the borrowed temporaries (U0, F0) are returned with restoreGlobalVec(); the transient rigid-rotation vector from _rigid_rotation_global is restored at every call site (gauge removal, _rotated_nullspace, _rotation_is_nullspace), and the transient duplicates in _rotation_is_nullspace are destroyed. Prevents pool exhaustion / memory growth over repeated (time-stepping) solves. - Module docstring corrected: the solve is iterative fieldsplit-Schur (custom-FMG / GAMG) by default with LU opt-in, and σ_nn reuses the shared boundary_flux de-smear — not the "direct LU, FMG later" of the original increment-1 note. Deferred (noted on the PR): building Q's identity via a per-row Mat.setValue loop is a one-time-per-solve setup that only matters at very large scale; left as a follow-up to avoid perturbing the well-tested Q construction. Verified: test_1018 serial 6/6; test_1064 parallel 5/5 at -n 2 and -n 4. Underworld development team with AI support from Claude Code

…og, value-first call-site sweep (WE-01..03,05,06,08,09,10) (#338) * docs(WE-01): adopt the one-governing-doc-per-topic authority map Repoint CLAUDE.md's Data Access 'Authoritative Reference' from the stale UW3_Style_and_Patterns_Guide.md to subsystems/data-access.md (the guide it crowned teaches patterns the code deprecates at runtime — DOC-04), and record the Style Charter §10 authority table in docs/developer/index.md as the master authority index. The Charter is added to the Getting Started toctree (removes a baseline 'not included in any toctree' warning). Finding: DOC-04 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-02): de-drift the Style Guide's four stale normative sections Rewrites the sections DOC-01 verified as contradicting the settled standards: - Docstring format: the 'Markdown Docstrings for pdoc/pdoc3' section is replaced by the NumPy/Sphinx RST standard (worked example with :math: and Parameters/Returns/Examples/Notes; conversion tracked in docs/plans/docstring-conversion-plan.md), per Style Charter section 6. - Doc file format: Quarto .qmd prescription (zero .qmd files exist in the repo) replaced by MyST .md/Sphinx guidance matching CLAUDE.md; migration table row updated. - Data access examples: 'Preferred' coordinate examples now use the real, runnable API — mesh.X.coords (read), mesh.deform() (coordinate changes), and the swarm.coords getter/setter for particle positions. The previous 'Preferred' example swarm.data += displacement raises AttributeError (getter-only property — SWARM-13 evidence); mesh.data warns at runtime. The private-attribute migration advice (swarm._particle_coordinates, mesh._deform_mesh presented as the NEW pattern) is deleted. - Front matter: the 21-line Quarto YAML header is replaced by a minimal MyST title block, and the guide now states that the UW3 Style Charter is the normative contract and wins on conflict. All replacement examples verified against current source: Swarm.coords setter (swarm.py), Mesh.deform (discretisation_mesh.py:3133), uw.synchronised_array_update / NDArray_With_Callback.delay_callbacks_global. Findings: DOC-01, SWARM-13 (style-guide part). Underworld development team with AI support from Claude Code * docs(WE-03): regenerate the docstring review queue; add the sweep to the release checklist The queue (last generated 2026-01-13, cdf5bb2) misrepresented the codebase both ways: it flagged now-complete items (solve, SNES_Scalar) as missing and contained zero entries for the June 2026 API (DOC-02). Regenerated over src/underworld3/**/*.py + **/*.pyx at the current tip. Two bugs in scripts/docstring_sweep.py's regex-based Cython parser made the regenerated queue lie about .pyx docstrings and are fixed as part of making the regeneration meaningful: - the indent group '(\s*)' with re.MULTILINE consumed preceding blank lines, shifting the computed definition line so the docstring search started ON the def/class line and always missed; - the docstring search started at the definition line rather than after the (possibly multi-line) signature, so long signatures hid their docstrings; - raw-string docstrings (r""", the norm in the solver .pyx) were not recognised. DOC-02 cross-validation on the regenerated queue now passes: solve / SNES_Scalar in the solver pyx are no longer flagged 'none'; the queue contains the June API (add_nitsche_bc, add_rotated_freeslip_bc, boundary_flux, set_custom_fmg, consistent_jacobian: 13 mentions) and flags the DOC-05 targets (Swarm.advection x2, read_timestep, write_proxy) as undocumented. Also adds the sweep to the quarterly release checklist (guides/release-process.md) so the queue cannot go stale unnoticed again. Findings: DOC-02 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-05): backfill the changelog for May - early July 2026; add the changelog sweep to the release checklist The changelog (the quarterly CIG/stakeholder record) ended in April 2026 while ~117 first-parent commits landed May through early July (DOC-03). Backfilled at the existing conceptual granularity — 14 grouped entries, grouped by subsystem rather than by PR, matching the established format (### Title (Month Year), bold lead sentence, hyphen bullets, inline PR references): - New '2026 Q3 (July - September)' section: the July 2026 quality campaign (#309-#313, #317, #322-#326, #329, #334 as grouped entries), rotated strong free-slip / boundary traction / dynamic topography (#293, #294, #298, #306), generalized geometric multigrid via custom prolongation (#290, #297), consistent Jacobian tangent (#258), swarm correctness (#216, #313, #323, #329), numpy 2 support (#301, #305). - Extended '2026 Q2' section with the May-June entries: mesh adaptation movers (#190, #209, #213, #228, #259, #264, #266), moving-mesh field transfer / deform() (#246, #249, #251), semi-Lagrangian accuracy controls (#164, #183, #185-#189, #208, #220), snapshot/checkpoint toolkit (#146, #195, #196, #198), Stokes_Constrained (#224, #229, #240, #265), local-h Nitsche + boundary-slip surfaces (#225, #241, #275), units interoperability (#277, #278, #283, #284), memory/evaluation/solver infrastructure (#161, #177-#179, #181, #182, #222, #237, #250, ...). Every entry is backed by a merged commit on development (verified against git log --first-parent aed517f..3184a40). Also adds a quarterly-changelog sweep step beside the docstring sweep in the release checklist (guides/release-process.md) per DOC-03's proposed fix. Findings: DOC-03 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-06): status headers on the unmarked design docs (per-doc git verification) Adds one-to-three-line Status markers to the 13 design docs that lacked one, following the directory's existing conventions (**Status**: line under the title; status: key inside existing YAML frontmatter for the three frontmatter-only docs), and corrects the stale 'Design Phase' marker on MATHEMATICAL_MIXIN_DESIGN.md (the mixin ships in utilities/mathematical_mixin.py). Every stamp was verified against git history (git log --follow dates) and the current source tree before writing: - Implemented: jacobian-consistent-tangent (PR #258, c63cd70), fmg-checkpoint-hierarchy (3cd73cd), petsc-dmplex-checkpoint-reload-plan (PR #146, write_timestep(petsc_reload=True) in tree), fault-refinement-simplification (smooth_mesh_interior / metric_density_from_gradient / fault_comb_metric all in tree), MATHEMATICAL_MIXIN_DESIGN. - Current reference/contract: mesh-adaptation-formulation, ND_UNITS_BOUNDARY_CONTRACT (PR #278, e0ece9a). - Investigation records (preserved via PR #245, 34a9dd4; production geometric-MG is custom prolongation, PR #290): snesfas-feasibility, snesfas-vanka-feasibility-study. - Design notes / prototypes with honest gaps: in_memory_checkpoint_design (not implemented, per its own trailing Status section), submesh-solver-architecture (extract_region/extract_surface exist; coarsened_companion does not). - Historical: ARCHITECTURE_ANALYSIS (persistence.py layout superseded), COORDINATE_MIGRATION_GUIDE (transition shipped), WHY_UNITS_NOT_DIMENSIONALITY (decision record). The audit's ~16 estimate over-counted: re-derived at this tip, 13 docs were unmarked plus one marked-but-stale (DOC-07). Findings: DOC-07 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-08): convert units.py public docstrings Google -> NumPy style Docstring-only conversion of the 18 public module-level functions that carried Google-style Args:/Returns:/Raises:/Examples: labels (check_units_consistency, get_dimensionality, get_units, non_dimensionalise, show_nondimensional_form, simplify_units, create_quantity, convert_units, to_base_units, to_reduced_units, to_compact, get_scaling_coefficients, set_scaling_coefficients, validate_expression_units, assert_dimensionality, validate_coordinates_dimensionality, enforce_units_consistency, require_units_if_active, convert_angle_to_degrees) to the NumPy/Sphinx standard (Style Charter section 6). dimensionalise was already NumPy style; one-line docstrings and private helpers are untouched. No code, signature, or behaviour changes (verified: every diff hunk is inside a docstring; ast.parse clean). Finding: API-12 (docs/reviews/2026-07/API-CONSISTENCY-REVIEW.md). Underworld development team with AI support from Claude Code * docs(WE-09): sweep call sites of the newer BC methods to value-first (conds, boundary, ...) order Wave C (#334) made the ORIGINAL value-first order canonical for add_nitsche_bc / add_rotated_freeslip_bc / add_constraint_bc (maintainer decisions D2/D3; Style Charter section 6) with deprecation shims for the legacy boundary-first and g= spellings. This sweep updates every call site of those THREE methods to the canonical order so nothing in the repository exercises the shims — 74 sites total: - tests/: 63 call sites across 12 files (test_1017, test_1018, test_1060, test_1061, test_1062, test_1064, test_1065 x2 serial; parallel test_1017, test_1062, test_1063, test_1064). tests/test_0641_wave_c_api_shims.py is deliberately untouched — its legacy-order calls ARE the deprecation contract. - docs/: 7 sites (curved-boundary-conditions.md x4, CONSTRAINED_FREESLIP_MULTIPLIER.md call + signature line, examples/submesh_investigation/test_region_ds_nitsche.py). - .claude/skills/: 3 sites (adapt-on-top-faults x2, free-surface-convection x1). - CLAUDE.md: 1 signature reference (free-slip BC preference section). The ~1,370 legacy-trio (add_dirichlet_bc/add_natural_bc/add_essential_bc) sites already conform and are untouched per the D2 decision. The audit review documents under docs/reviews/2026-07/ record the pre-decision state as evidence and are not swept. Discovered while verifying the swept tests run warning-free: the Wave C zero-datum guard in add_rotated_freeslip_bc rejects FLOAT zero (sympy.sympify(0.0) != 0 is structurally True), so the canonical add_rotated_freeslip_bc(0.0, boundary) raises NotImplementedError while conds=0 works. Filed as issue #336 with a TODO(BUG) marker at the guard (comment-only src touch); the swept call sites use the working integer form add_rotated_freeslip_bc(0, boundary). No fix applied here (Charter section 9 scope discipline). Findings: API-01/API-02 sweep (WE-09, REMEDIATION-WORKLIST.md). Underworld development team with AI support from Claude Code

lmoresi added 11 commits June 28, 2026 10:31

docs(custom_mg): mark Stokes velocity-block integration done (Phase 1…

cb77289

…, Step 1) Underworld development team with AI support from Claude Code

docs(custom_mg): mark factory removal done (Phase 1, Step 2)

2e2c4f9

Underworld development team with AI support from Claude Code

docs(custom_mg): record np>1 serial guard (Phase 1, Step 3)

caff202

Underworld development team with AI support from Claude Code

docs(custom_mg): record parallel np>1 done (Phase 1, Step 3)

af8e734

Underworld development team with AI support from Claude Code

lmoresi marked this pull request as ready for review June 29, 2026 04:20

Copilot AI review requested due to automatic review settings June 29, 2026 04:20

Copilot AI reviewed Jun 29, 2026

lmoresi added 2 commits June 29, 2026 14:43

docs(custom_mg): solver-family coverage matrix + constrained-parallel…

32fec5a

… limitation Underworld development team with AI support from Claude Code

lmoresi mentioned this pull request Jun 29, 2026

Stokes_Constrained segfaults at np>1 in the interior-multiplier section reduction #291

Open

lmoresi added 2 commits June 29, 2026 15:01

docs(custom_mg): link constrained-parallel segfault to issue #291

d7a4ca1

Underworld development team with AI support from Claude Code

lmoresi mentioned this pull request Jul 1, 2026

Rotated strong free-slip BC + σ_nn dynamic-topography (parallel), free-surface hand-off #293

Merged

lmoresi requested a review from Copilot July 2, 2026 01:12

Copilot started reviewing on behalf of lmoresi July 2, 2026 01:12 View session

Copilot AI reviewed Jul 2, 2026

View reviewed changes

lmoresi merged commit f6a75ef into development Jul 2, 2026
1 check passed

lmoresi mentioned this pull request Jul 3, 2026

custom_mg prolongation builders should use the coarse DMPlex + native point location, not scipy Delaunay over the DOF cloud #307

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Generalized FMG hierarchy: custom-P transfers for arbitrary nested/non-nested grids (Layer 1)#290

Generalized FMG hierarchy: custom-P transfers for arbitrary nested/non-nested grids (Layer 1)#290
lmoresi merged 16 commits into
developmentfrom
feature/custom-mg-prolongation

lmoresi commented Jun 28, 2026 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI left a comment

Uh oh!

lmoresi commented Jul 2, 2026 •

edited

Loading

Uh oh!

lmoresi commented Jul 2, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

lmoresi commented Jun 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What

Status — Layer 1, Phase-1 hardening complete

Hardening steps (all done)

Validation

Out of scope (later)

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

lmoresi commented Jul 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lmoresi commented Jul 2, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

lmoresi commented Jun 28, 2026 •

edited

Loading

lmoresi commented Jul 2, 2026 •

edited

Loading