From 948c1f5666f9286da96dc7d749c1163bcfc7d631 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 15 Jul 2026 23:45:25 -1000 Subject: [PATCH 1/2] refactor(custom-mg): D-79..D-87 readability rows + minimal-control multigrid design note Wave D behaviour-neutral readability pass over utilities/custom_mg.py (docs/reviews/2026-07 REMEDIATION-WORKLIST, D-custom_mg table; all nine rows re-verified live at the post-#369 base): - D-79 LevelLayout NamedTuple replaces the magic lay[3] tuple indexing - D-80 _clone_dm_with_solver_discretisation extracted (copyDS trick documented once, on the helper) - D-81 _assert_no_zero_columns_serial extracted beside the parallel guard - D-82 getNumFields except narrowed to PETSc.Error + sanctioned comment - D-83 parallel-path parameters renamed to the serial vocabulary (coarse_coords/fine_coords/coarse_layout/fine_layout) - D-84 semicolon-packed lines split (6 sites) - D-85 'sub' prefix-string collision renamed vel_prefix, moved into the verbose block - D-86 RBF r=0 clamp commented (r^2 log r -> 0; clamp keeps log finite) - D-87 TODO(deprecate) lifespan marker on the legacy finest-only path, stating the one behavioural difference (no BC-per-level reduction) Plus docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md (Status: PROPOSED): survey of the five existing MG paths, a proposed minimal-control auto ladder, and seven numbered maintainer questions. No implementation. Bit-identical check: custom-MG suites identical pre/post edits (serial 1015+1016+1017: 14 passed; np2 1017: 6 passed). Serial 'level_1 and tier_a' gate identical to the development baseline (386 passed / 12 skipped / 2 xfailed / 1 xpassed). Scanner green; docs build green. Underworld development team with AI support from Claude Code --- .../MULTIGRID_MINIMAL_CONTROL_2026-07.md | 159 +++++++++++++++++ docs/developer/index.md | 1 + src/underworld3/utilities/custom_mg.py | 164 +++++++++++------- 3 files changed, 266 insertions(+), 58 deletions(-) create mode 100644 docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md diff --git a/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md b/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md new file mode 100644 index 00000000..e2454cad --- /dev/null +++ b/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md @@ -0,0 +1,159 @@ +# Multigrid with minimal user control — design note + +**Status**: PROPOSED (2026-07). No implementation in this PR — this note +surveys what exists after #369 and asks the maintainer to settle the +default-behaviour questions at the end. + +**Maintainer direction** (2026-07): *"the laplacian solvers all benefit +from MG; the real question is how much control is handed to / required +from the user? I would prefer that this is, by default, minimal."* + +--- + +## 1. What exists today (post-#369) + +Five distinct ways a UW3 solver ends up with (or without) multigrid: + +### 1.1 `solver.preconditioner` — the one documented knob + +`"auto"` (default) / `"fmg"` / `"gamg"` on `SolverBaseClass` +(`cython/petsc_generic_snes_solvers.pyx`). Resolution happens at +`_build` time in `_apply_preconditioner_options()`: + +- `"auto"`: **native PETSc geometric FMG** when the mesh carries a + genuine refinement hierarchy (`len(mesh.dm_hierarchy) > 1`, i.e. + built with `refinement >= 1`), otherwise GAMG. Deliberately + conservative: it only ever *adds* geometric MG, never rewrites a PC + the user (or an internal utility solver) configured — an explicit + `pc_type` in the options DB latches `_pc_user_override` and wins + forever after. +- Native FMG is **locked out for single-field solvers** (#276): PETSc's + `DMCreateInjection` fails unpredictably (geometry × degree × + refinement) on refined DMPlex for scalar/vector discretisations, so + `"auto"`/`"fmg"` on a scalar solver falls back to GAMG with a + warning. The Stokes **velocity block** (`fieldsplit_velocity_`) is + the validated native-FMG path. + +### 1.2 Raw PETSc options — the expert escape hatch + +`solver.petsc_options[...]` reaches everything (`pc_type=mg`, SNESFAS, +smoother choices, level counts). The `_pc_user_override` latch above is +what keeps 1.1 from fighting this. This is also how internal utility +solvers (projections, smoothers) carry their own tuned PC bundles. + +### 1.3 `solver.set_custom_fmg(coarse_meshes, ...)` — custom-P geometric MG + +The generalized custom-prolongation hierarchy +(`utilities/custom_mg.py`, `CustomMGHierarchy`): barycentric/RBF +node-level transfers, BC-per-level reduction, Galerkin (RAP) coarse +operators. Works for **any coarse/fine mesh pair** (nested or not), +**any single field** (`field_id=0` targets the Stokes velocity block), +and — since #369 — **in parallel** (co-partitioned rank-local fast path +plus a cross-partition fallback selected automatically on the +zero-column signature). Notably it does *not* need `DMCreateInjection`, +so it does not suffer #276's single-field lock-out. + +### 1.4 Mesh-owned auto-injection on `adapt()` children (#369) + +A refinement child from `mesh.adapt(...)` carries +`mesh._custom_mg_coarse_meshes` (the static coarse tail). The first +`solve()` of *any* solver on such a mesh lazily builds and installs a +custom-P hierarchy (`auto_inject_custom_mg`) — geometric MG on the +refined operator with **no user call at all**. It is opportunistic and +fail-safe: any build/validation failure warns and falls back to the +solver's default PC (a solve must never crash because an optional +preconditioner upgrade failed). + +### 1.5 Legacy `solver.set_custom_mg(...)` — deprecated + +Serial-only, single-field, finest-only reduction (no BC-per-level). +Deprecated in favour of `set_custom_fmg`; carries a +`TODO(deprecate)` lifespan marker (this PR, D-87). + +**Adjacent fact**: SNESFAS (nonlinear multigrid) is reachable through +1.2 alone today (validated feasible via options; no UW3 surface). + +## 2. The problem + +The capability matrix is excellent but the *decision burden* sits with +the user: whether their Poisson solve gets MG depends on how the mesh +was built (`refinement=`), which solver class they used (velocity block +vs scalar), whether the mesh is an adapt child, and whether they know +`set_custom_fmg` exists. The #276 lock-out means the most common case — +**a scalar Laplacian on an ordinary refined mesh** — silently gets +GAMG even though a robust geometric path (1.3) exists one call away. + +## 3. Proposal — minimal-control default + +One principle: **solvers auto-configure the best available MG; the user +supplies nothing.** Expert control remains via exactly two escape +hatches: `solver.petsc_options` (always wins, latched) and +`set_custom_fmg` (explicit hierarchy). + +Proposed `"auto"` resolution ladder, per solver, at `_build` time: + +1. **User options present** (`_pc_user_override`) → hands off. +2. **Explicit `set_custom_fmg` hierarchy** → build + install it (loud + failure — the user asked for it). +3. **Adapt-child mesh-owned hierarchy** → auto-inject custom-P + (fail-safe, as today). +4. **Refinement hierarchy + multi-field velocity block** → native + PETSc FMG (as today). +5. **Refinement hierarchy + single-field solver** → *NEW*: + auto-build a custom-P FMG from `mesh.dm_hierarchy`'s coarse levels + (closing the #276 gap with the injection-free path) instead of + falling back to GAMG. +6. **No hierarchy** → GAMG (as today). + +Steps 1–4 and 6 are current behaviour; only step 5 is new machinery +(and it reuses 1.3 wholesale — the coarse meshes already exist as +`mesh._coarse_level_meshes()`). + +Under this ladder `solver.preconditioner` keeps its three values but +`"fmg"` on a single-field solver stops meaning "warn and use GAMG" and +starts meaning "custom-P geometric FMG". + +## 4. Questions for the maintainer + +1. **Adopt step 5?** Should `"auto"` route single-field solvers with a + refinement hierarchy to custom-P FMG (closing the #276 GAMG + fallback), or does custom-P need more production mileage before it + is a silent default? (An intermediate: `"fmg"` explicitly requests + it, `"auto"` keeps GAMG.) +2. **Failure loudness.** Auto-injection today warns-and-falls-back + (1.4). Is that the right contract for *every* automatic MG decision, + or should an explicit `preconditioner="fmg"` request fail hard when + the hierarchy cannot be built? +3. **Utility solvers.** Internal solvers (projections, mesh-smoothing + Poisson solves) are currently excluded via their own tuned options + + the override latch. Keep the latch as the exclusion mechanism, or + add an explicit `wants_mg = False` class attribute so exclusion is + visible in the solver class rather than in options state? +4. **Naming.** `solver.preconditioner` takes algorithm names + (`"fmg"`/`"gamg"`). Under the purposeful-naming ruling, should the + user surface become capability-shaped — + e.g. `solver.multigrid = "auto" | "geometric" | "algebraic" | "off"` + — with the current spellings kept as accepted aliases? Or is + `preconditioner` (an established PETSc term) exempt? +5. **Legacy `set_custom_mg` removal.** It is deprecated with a + lifespan marker. Remove in the next release cycle (with test_1015's + legacy cases), or keep one more cycle? +6. **SNESFAS scope.** Nonlinear MG stays options-only (out of the + minimal-control surface), or should a future `"auto"` consider it + for strongly nonlinear Laplacian-family solves? +7. **`mesh.dm_hierarchy` as the trigger.** The whole ladder keys off + `refinement >= 1` at mesh build. Is it acceptable that a mesh built + without `refinement=` never gets geometric MG silently (GAMG only), + or should UW3 build a default hierarchy (e.g. `refinement=1`) for + the common constructors so step 5 usually has something to work + with? + +## 5. References + +- `docs/developer/design/GENERALIZED_FMG_HIERARCHY_AND_ADAPT.md` — + the custom-P hierarchy design (BC-per-level invariant). +- `docs/advanced/multigrid-preconditioning.md` — user-facing guide. +- `utilities/custom_mg.py` — implementation (readability pass D-79..87 + in the PR that carries this note). +- #276 (single-field native-FMG lock-out), #369 (parallel custom-MG + + auto-injection), #290 (custom prolongation). diff --git a/docs/developer/index.md b/docs/developer/index.md index a720ba89..eed57fbf 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -149,6 +149,7 @@ design/ADAPTIVE_MESHING_DESIGN design/mesh-adaptation-formulation design/LAYER2_SBR_ADAPT_ON_TOP design/NVB_GRADED_ADAPT +design/MULTIGRID_MINIMAL_CONTROL_2026-07 design/ARCHITECTURE_ANALYSIS design/MATHEMATICAL_MIXIN_DESIGN design/COORDINATE_MIGRATION_GUIDE diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 9d4a252d..9023f610 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -38,6 +38,8 @@ operators from the DM hierarchy, so our explicit ``P`` is used. """ +from typing import NamedTuple + import numpy as np from petsc4py import PETSc @@ -64,14 +66,18 @@ def barycentric_prolongation(coarse_coords, fine_coords): for i in range(fine_coords.shape[0]): s = simp[i] if s < 0: # outside coarse hull → nearest coarse DOF - rows.append(i); cols.append(int(tree.query(fine_coords[i])[1])); vals.append(1.0) + rows.append(i) + cols.append(int(tree.query(fine_coords[i])[1])) + vals.append(1.0) continue verts = tri.simplices[s] T = tri.transform[s] bary = T[:dim].dot(fine_coords[i] - T[dim]) w = np.append(bary, 1.0 - bary.sum()) for k in range(dim + 1): - rows.append(i); cols.append(int(verts[k])); vals.append(float(w[k])) + rows.append(i) + cols.append(int(verts[k])) + vals.append(float(w[k])) return sp.csr_matrix((vals, (rows, cols)), shape=(fine_coords.shape[0], coarse_coords.shape[0])) @@ -85,6 +91,8 @@ def rbf_prolongation(coarse_coords, fine_coords, smooth=0.0): from scipy.spatial.distance import cdist def phi(r): + # r² log r → 0 as r → 0; the clamp only keeps log(0) finite at + # coincident points — it does not perturb the kernel value. r = np.where(r == 0.0, 1e-30, r) return r ** 2 * np.log(r) @@ -220,7 +228,8 @@ def _reduce_to_global(dm, full_coords): gvec = dm.getGlobalVec() cdim = full_coords.shape[1] if lvec.getLocalSize() != full_coords.shape[0]: - dm.restoreLocalVec(lvec); dm.restoreGlobalVec(gvec) + dm.restoreLocalVec(lvec) + dm.restoreGlobalVec(gvec) raise RuntimeError( f"custom_mg: local DOF count {lvec.getLocalSize()} != coord count " f"{full_coords.shape[0]} — degree/continuity mismatch") @@ -252,7 +261,10 @@ def _field_subdm(dm, field_id): try: if dm.getNumFields() <= 1: return dm - except Exception: + except PETSc.Error: + # Sanctioned: a DM whose DS is not created yet cannot report its + # field count; treat it as single-field rather than mask a real + # multi-field mistake with a broad swallow. return dm _iset, sub = dm.createSubDM(field_id) return sub @@ -278,24 +290,31 @@ def _reduced_map(dm, field_id=None): return r2f, n_full -def _coarse_reduced_map(solver, coarse_mesh, field_id=None): - """A COARSE level's BC-constrained reduced map, with NO throwaway solver. - - Clone the coarse mesh DM and copy the (built) finest ``solver``'s fields + DS - onto it. The DS carries UW's exact essential-BC definitions (the custom - essential-field boundaries), and is a topology-independent discretisation spec, - so ``createDS`` constrains the matching boundary DOFs on ANY coarse mesh that - carries the same boundary labels (nested or not). The resulting global section - gives the same reduced map a full same-discretisation solver would — validated - identical to the old ``level_solver_factory`` path. Leak-free: DM ops only, no - SNES / JIT. +def _clone_dm_with_solver_discretisation(solver, coarse_mesh): + """Clone a coarse mesh DM carrying the finest ``solver``'s discretisation. - NOTE: serial, like ``_reduced_map`` (uses local indices); parallel correctness - is a Phase-3 item.""" + Copy the (built) solver's fields + DS onto the clone. The DS carries UW's + exact essential-BC definitions (the custom essential-field boundaries), and + is a topology-independent discretisation spec, so ``createDS`` constrains + the matching boundary DOFs on ANY coarse mesh that carries the same + boundary labels (nested or not). The resulting global section gives the + same reduced map a full same-discretisation solver would — validated + identical to the old ``level_solver_factory`` path. Leak-free: DM ops + only, no SNES / JIT.""" cdm = coarse_mesh.dm.clone() solver.dm.copyFields(cdm) solver.dm.copyDS(cdm) cdm.createDS() + return cdm + + +def _coarse_reduced_map(solver, coarse_mesh, field_id=None): + """A COARSE level's BC-constrained reduced map, with NO throwaway solver + (see :func:`_clone_dm_with_solver_discretisation` for the copyDS trick). + + NOTE: serial, like ``_reduced_map`` (uses local indices); parallel correctness + is a Phase-3 item.""" + cdm = _clone_dm_with_solver_discretisation(solver, coarse_mesh) return _reduced_map(cdm, field_id) @@ -317,14 +336,26 @@ def _reduced_transfer(coarse_coords, fine_coords, r2f_c, r2f_f, ncomp, builder): # its block of P from its LOCAL (ghost-inclusive) coarse coords — point-location # is rank-local. The reduced global numbering rides the DM global section. # --------------------------------------------------------------------------- # +class LevelLayout(NamedTuple): + """Parallel DOF layout of one MG level. + + ``l2g[i]`` is the GLOBAL reduced index of local DOF ``i`` — ghost-resolved + to the owner's global index, and ``-1`` for a BC-constrained DOF. + ``[rstart, rend)`` is this rank's owned global range; ``n_full`` is the + local (ghost-inclusive) full DOF count.""" + + l2g: np.ndarray + rstart: int + rend: int + n_full: int + + def _level_dof_layout(dm, field_id=None): - """Parallel DOF layout for one level: ``(l2g, rstart, rend, n_full)``. + """Parallel DOF layout for one level (a :class:`LevelLayout`). - ``l2g[i]`` is the GLOBAL reduced index of local DOF ``i`` — ghost-resolved to - the owner's global index, and ``-1`` for a BC-constrained DOF. Built by - scattering each owned global index out to the local (incl. ghost) layout via - ``globalToLocal`` (constrained local DOFs have no global source, so they keep - the pre-set ``-1``). ``[rstart, rend)`` is this rank's owned global range.""" + ``l2g`` is built by scattering each owned global index out to the local + (incl. ghost) layout via ``globalToLocal`` (constrained local DOFs have no + global source, so they keep the pre-set ``-1``).""" sub = _field_subdm(dm, field_id) gv = sub.getGlobalVec() lv = sub.getLocalVec() @@ -336,20 +367,18 @@ def _level_dof_layout(dm, field_id=None): n_full = lv.getLocalSize() sub.restoreGlobalVec(gv) sub.restoreLocalVec(lv) - return l2g, rstart, rend, n_full + return LevelLayout(l2g, rstart, rend, n_full) def _coarse_dof_layout(solver, coarse_mesh, field_id=None): """Parallel coarse-level DOF layout, no throwaway solver — same copyDS trick - as :func:`_coarse_reduced_map` but returning the parallel ``l2g`` layout.""" - cdm = coarse_mesh.dm.clone() - solver.dm.copyFields(cdm) - solver.dm.copyDS(cdm) - cdm.createDS() + as :func:`_coarse_reduced_map` but returning the parallel layout.""" + cdm = _clone_dm_with_solver_discretisation(solver, coarse_mesh) return _level_dof_layout(cdm, field_id) -def _build_parallel_transfer(cc, fc, lay_c, lay_f, ncomp, builder, comm): +def _build_parallel_transfer(coarse_coords, fine_coords, coarse_layout, + fine_layout, ncomp, builder, comm): """One reduced->reduced prolongation as an MPIAIJ matrix. Node-level barycentric/RBF weights are built rank-locally (coarse LOCAL coords @@ -357,11 +386,11 @@ def _build_parallel_transfer(cc, fc, lay_c, lay_f, ncomp, builder, comm): fine OWNED DOF becomes a global row; its coarse contributions map through the coarse ``l2g`` to global columns (off-rank columns are fine for MPIAIJ). Constrained coarse DOFs (``l2g == -1``) drop out -> reduced->reduced.""" - l2g_c, cstart, cend, _ = lay_c - l2g_f, fstart, fend, _ = lay_f - Pn = builder(cc, fc).tocsr() # (n_f_nodes, n_c_nodes), local + l2g_c = coarse_layout.l2g + l2g_f, fstart, fend = fine_layout.l2g, fine_layout.rstart, fine_layout.rend + Pn = builder(coarse_coords, fine_coords).tocsr() # (n_f_nodes, n_c_nodes), local nloc_f = fend - fstart - nloc_c = cend - cstart + nloc_c = coarse_layout.rend - coarse_layout.rstart P = PETSc.Mat().create(comm=comm) P.setSizes(((nloc_f, None), (nloc_c, None))) @@ -386,7 +415,7 @@ def _build_parallel_transfer(cc, fc, lay_c, lay_f, ncomp, builder, comm): return P -def _gather_coarse_cloud(cc, lay_c, ncomp, comm): +def _gather_coarse_cloud(coarse_coords, coarse_layout, ncomp, comm): """All-gather the (small) coarse node cloud across ranks and deduplicate. Returns ``(coords_u, cols_u)``: ``coords_u`` is the FULL coarse node cloud @@ -396,16 +425,18 @@ def _gather_coarse_cloud(cc, lay_c, ncomp, comm): coordinate; constrained nodes stay in the cloud as barycentric vertices but carry ``-1`` so they are dropped from the transfer columns.""" m4 = comm.tompi4py() - ncn = cc.shape[0] - cols_local = np.asarray(lay_c[0], dtype=np.int64).reshape(ncn, ncomp) - cc_all = np.vstack(m4.allgather(np.ascontiguousarray(cc, dtype=float))) + ncn = coarse_coords.shape[0] + cols_local = np.asarray(coarse_layout.l2g, dtype=np.int64).reshape(ncn, ncomp) + cc_all = np.vstack(m4.allgather( + np.ascontiguousarray(coarse_coords, dtype=float))) cols_all = np.vstack(m4.allgather(np.ascontiguousarray(cols_local))) _key, uidx = np.unique(np.round(cc_all, 9), axis=0, return_index=True) uidx = np.sort(uidx) return cc_all[uidx], cols_all[uidx] -def _build_crosspart_transfer(cc, fc, lay_c, lay_f, ncomp, builder, comm): +def _build_crosspart_transfer(coarse_coords, fine_coords, coarse_layout, + fine_layout, ncomp, builder, comm): """One reduced->reduced prolongation for a NON-NESTED (independently partitioned) coarse level. @@ -417,13 +448,13 @@ def _build_crosspart_transfer(cc, fc, lay_c, lay_f, ncomp, builder, comm): small — that is what makes a coarse MG level coarse), so point location spans partitions. Columns are the coarse GLOBAL reduced indices (off-rank columns are fine for MPIAIJ); constrained coarse DOFs (col < 0) drop out.""" - _l2g_c, cstart, cend, _ = lay_c - l2g_f, fstart, fend, _ = lay_f - coords_u, cols_u = _gather_coarse_cloud(cc, lay_c, ncomp, comm) + l2g_f, fstart, fend = fine_layout.l2g, fine_layout.rstart, fine_layout.rend + coords_u, cols_u = _gather_coarse_cloud(coarse_coords, coarse_layout, + ncomp, comm) - Pn = builder(coords_u, fc).tocsr() # (n_f_nodes_local, Nu) + Pn = builder(coords_u, fine_coords).tocsr() # (n_f_nodes_local, Nu) nloc_f = fend - fstart - nloc_c = cend - cstart + nloc_c = coarse_layout.rend - coarse_layout.rstart P = PETSc.Mat().create(comm=comm) P.setSizes(((nloc_f, None), (nloc_c, None))) @@ -453,12 +484,14 @@ def _count_zero_columns_parallel(P, comm): image). Column sums via P^T·1 (barycentric/RBF weights are positive, so a zero sum is an empty column). Used to auto-detect a cross-partition point-location miss (non-nested coarse level) and, separately, as the hard guard below.""" - ones_f = P.createVecLeft(); ones_f.set(1.0) + ones_f = P.createVecLeft() + ones_f.set(1.0) colsum = P.createVecRight() P.multTranspose(ones_f, colsum) nzero_local = int((colsum.array == 0.0).sum()) nzero = comm.tompi4py().allreduce(nzero_local) - ones_f.destroy(); colsum.destroy() + ones_f.destroy() + colsum.destroy() return nzero @@ -472,6 +505,17 @@ def _assert_no_zero_columns_parallel(P, comm): f"image) — BC-per-level reduction failed; coarse operator would be singular.") +def _assert_no_zero_columns_serial(P_csr, level): + """Serial zero-column guard — same physics as the parallel guard above: a + coarse DOF with no fine image makes the Galerkin coarse operator singular.""" + nzero = int((np.asarray((P_csr != 0).sum(axis=0)).ravel() == 0).sum()) + if nzero: + raise RuntimeError( + f"transfer {level - 1}->{level} has {nzero} zero columns (coarse DOFs " + f"with no fine image) — BC-per-level reduction failed; coarse " + f"operator would be singular.") + + def _configure_pcmg(pc, Ps): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -583,7 +627,6 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): outer_pc.setUp() vel_ksp = outer_pc.getFieldSplitSubKSP()[0] vel_pc = vel_ksp.getPC() - sub = vel_pc.getOptionsPrefix() or "fieldsplit_velocity_" A_vv, P_vv = vel_pc.getOperators() # capture before reset (reset drops them) # 3. fresh PCMG on the velocity sub-block from our Ps @@ -598,8 +641,9 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): if verbose: from underworld3 import mpi + vel_prefix = vel_pc.getOptionsPrefix() or "fieldsplit_velocity_" mpi.pprint(f"[{solver.name}] custom FMG installed on velocity block: " - f"{len(Ps) + 1} levels, sub-prefix {sub!r}, " + f"{len(Ps) + 1} levels, sub-prefix {vel_prefix!r}, " f"P sizes {[tuple(P.getSize()) for P in Ps]}") @@ -680,9 +724,11 @@ def build(self, solver): c = np.asarray(mesh._get_coords_for_basis(degree, continuous)) finest = (k == nlev - 1) if parallel: + # ``maps`` holds a LevelLayout per level in parallel, and a + # bare reduced->full index array per level in serial. lay = (_level_dof_layout(solver.dm, self.field_id) if finest else _coarse_dof_layout(solver, mesh, self.field_id)) - nfull = lay[3] + nfull = lay.n_full maps.append(lay) else: rmap, nfull = (_reduced_map(solver.dm, self.field_id) if finest @@ -692,7 +738,8 @@ def build(self, solver): if nfull % c.shape[0] != 0: raise RuntimeError( f"level {k}: full DOFs {nfull} not divisible by nodes {c.shape[0]}") - coords.append(c); ncomp.append(nc) + coords.append(c) + ncomp.append(nc) if len(set(ncomp)) != 1: raise RuntimeError(f"inconsistent component counts across levels: {ncomp}") @@ -727,12 +774,7 @@ def build(self, solver): else: Pr = _reduced_transfer(coords[l - 1], coords[l], maps[l - 1], maps[l], nc, self.builder) - zc = int((np.asarray((Pr != 0).sum(axis=0)).ravel() == 0).sum()) - if zc: - raise RuntimeError( - f"transfer {l-1}->{l} has {zc} zero columns (coarse DOFs with " - f"no fine image) — BC-per-level reduction failed; coarse " - f"operator would be singular.") + _assert_no_zero_columns_serial(Pr, l) Ps.append(_to_petsc_aij(Pr)) self.transfers = Ps return Ps @@ -753,8 +795,8 @@ def _assert_finest_matches_operator(solver, finest_map, parallel): if op_n <= 0: return if parallel: - _l2g, rstart, rend, _n = finest_map - red_n = int(solver.dm.comm.tompi4py().allreduce(int(rend - rstart))) + red_n = int(solver.dm.comm.tompi4py().allreduce( + int(finest_map.rend - finest_map.rstart))) else: red_n = int(len(finest_map)) # r2f length = reduced global size if red_n != op_n: @@ -876,6 +918,12 @@ def inject_custom_mg(solver): return # ---- legacy finest-only path (back-compat, serial only) ----------------- + # TODO(deprecate): remove together with SolverBaseClass.set_custom_mg and + # test_1015's legacy cases. The one behavioural difference from the + # hierarchy path: NO BC-per-level reduction — transfers are built on the + # full node cloud, which is only valid when the coarse levels are + # non-nested / unconstrained (a nested coarse boundary DOF coinciding with + # a BC-removed fine DOF would give a zero column -> singular coarse op). _require_serial("legacy custom_mg (set_custom_mg)") coarse_meshes = cfg["coarse_meshes"] builder = _BUILDERS[cfg["kind"]] From 115126426455dd31a2f186b79d5022fd9ccbd75d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 16 Jul 2026 10:03:35 -1000 Subject: [PATCH 2/2] docs(design): record the maintainer's seven multigrid rulings (2026-07-17); SNESFAS not-viable notes The minimal-control design note moves PROPOSED -> RULED: GAMG stays the automatic default (custom-P FMG opt-in); two-tier failure loudness; wants_mg class attribute; solver.multigrid capability spelling; set_custom_mg removed next cycle; no default hierarchies. SNESFAS is ruled NOT viable at present (no good preconditioners; loses the robust linear-solver path) - noted in the solver-strategies catalogue and the plasticity-solvers skill as future investigation only. Underworld development team with AI support from Claude Code --- .claude/skills/plasticity-solvers/SKILL.md | 10 +++++++ .../MULTIGRID_MINIMAL_CONTROL_2026-07.md | 30 +++++++++++++++++-- .../design/solver-strategies-catalogue.md | 10 +++++++ 3 files changed, 47 insertions(+), 3 deletions(-) diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md index 9c0c243b..5b468ef6 100644 --- a/.claude/skills/plasticity-solvers/SKILL.md +++ b/.claude/skills/plasticity-solvers/SKILL.md @@ -129,3 +129,13 @@ Override with `consistent_tangent=True` / `False` / `"continuation"`. Footnote: before this work UW3 differentiated the flux with the viscosity still wrapped, so `∂η/∂(grad v)` was dropped and viscoplastic solves silently ran the Picard tangent — the origin of the "~20 iterations is intrinsic" folklore. + +## SNESFAS — do not reach for it + +Nonlinear multigrid (SNESFAS) looks tempting for hard viscoplastic solves but is +**not a viable option** at present (maintainer ruling 2026-07-17): there are no +good preconditioners for the nonlinear hierarchy, and it abandons the robust +linear-solver path (consistent tangent / continuation + fieldsplit + MG) that +this skill is built around. It stays options-only for experiments; treat it as a +future investigation. See `docs/developer/design/solver-strategies-catalogue.md` +and `MULTIGRID_MINIMAL_CONTROL_2026-07.md` (ruling 6). diff --git a/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md b/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md index e2454cad..30e9d1a1 100644 --- a/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md +++ b/docs/developer/design/MULTIGRID_MINIMAL_CONTROL_2026-07.md @@ -1,8 +1,32 @@ # Multigrid with minimal user control — design note -**Status**: PROPOSED (2026-07). No implementation in this PR — this note -surveys what exists after #369 and asks the maintainer to settle the -default-behaviour questions at the end. +**Status**: RULED (2026-07-17) — the maintainer settled all seven questions; +rulings recorded below. No implementation in this PR; the rulings govern the +implementation PR that follows. + +## Maintainer rulings (2026-07-17) + +1. **Automatic routing**: GAMG remains the automatic default; custom-P + geometric FMG is **opt-in**, not auto-selected — it needs more production + mileage before it becomes a silent default (#276's fallback stands). +2. **Failure loudness — two-tier**: automatic/default MG decisions warn and + fall back to a working solve; an EXPLICIT user request (e.g. geometric FMG) + that cannot be honoured fails hard. +3. **Utility-solver exclusion**: an explicit `wants_mg = False` class + attribute — the solver class declares its MG policy readably; the options + latch is retired when this lands. +4. **User spelling**: `solver.multigrid = "auto" | "geometric" | "algebraic" + | "off"` — capability words per the purposeful-naming ruling. + `solver.preconditioner` remains the expert/PETSc-level knob. +5. **`set_custom_mg`**: removed next release cycle (one full deprecation + cycle, the Wave C pattern). +6. **SNESFAS**: options-only, and explicitly **not a viable option** at + present — "we don't have good preconditioners for this one and we seem to + lose the robust linear solver path if we adopt this strategy." Noted here + and in the solver skills as a future investigation only. +7. **Default hierarchies**: NO — meshes do not build a refinement hierarchy + by default. Hierarchy cost stays explicit and predictable; users who want + geometric MG ask for `refinement>=1`. **Maintainer direction** (2026-07): *"the laplacian solvers all benefit from MG; the real question is how much control is handed to / required diff --git a/docs/developer/design/solver-strategies-catalogue.md b/docs/developer/design/solver-strategies-catalogue.md index f23a338f..5294b9dc 100644 --- a/docs/developer/design/solver-strategies-catalogue.md +++ b/docs/developer/design/solver-strategies-catalogue.md @@ -355,6 +355,16 @@ degrades on poor cells (the GAMG anisotropy section above) — mesh-quality monitoring is therefore not aesthetic, it directly predicts solver robustness. +## SNESFAS (nonlinear multigrid) — NOT a viable option (ruled 2026-07-17) + +SNESFAS remains reachable via PETSc options for experimentation, but it is +**not part of any recommended strategy**: we do not have good preconditioners +for the nonlinear hierarchy, and adopting it means losing the robust +linear-solver path (Newton/Picard + fieldsplit + the MG ladder above) that the +rest of this catalogue is built on. Treat it as a future investigation, not a +tool to reach for when a solve struggles — the failure-class map below never +routes to it. See `MULTIGRID_MINIMAL_CONTROL_2026-07.md`, ruling 6. + ## Failure-class → strategy map (the picking guide) ```