From 9ce33f9d3d2af5d193dd0d0cecd45a7813a78e2d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 1 Sep 2026 18:20:26 -0700 Subject: [PATCH 1/2] Document why the MG coarse level must be small (#644, #643) bknight1 hit an OOM at Pawsey: the FMG coarse solve is redundant+LU, which replicates the coarse factorization on every rank, and his coarsest level still carried 14k velocity DOFs - ~300 MB/rank replicated, killed at the 900 MB cap. Agreed on the issue that this is not a bug but a documentation gap. The premise the docs never stated: multigrid is cheap because the bottom of the hierarchy is solved exactly at negligible cost, and that is only true if the bottom is genuinely tiny. Nothing in the method makes it so and nothing warns when it is not. A 14k-DOF coarsest level means the hierarchy was built from too fine a base with too few levels, so the premise never held - the fix is another level underneath, not a different bottom solve. Doubling the base cell size and adding one refinement level lands on the same finest grid with every coarser level ~4x smaller in 2-D (8x in 3-D); bknight1 measured that as 683 -> 416 MB/rank for +11% linear-solve time, with the finest grid and the SNES iteration sequence unchanged. Also records that this is PETSc's own choice, not ours: PCSetUp_MG selects PCREDUNDANT above one rank and PCLU at one rank, its source comment reading "coarse solve is (redundant) LU by default" (src/ksp/pc/impls/mg/mg.c), which I confirmed by inspection at np = 1, 2 and 4. We set it explicitly only so a bundle cannot inherit a sibling's leftover options (#468). Deleting the keys and letting PETSc decide would change nothing - which is the point, because it puts the attention back on the hierarchy. Two corrections to what was there. The page claimed the default coarse solve was a bare `lu` in both the bundle listing and the parallel section, and told users to set redundant themselves for large partitions - it has been the default all along. And the SVD path is documented as what it actually is: keyed on VERIFIED rotation null modes, not on a problem being "rotated", because a rotated problem with Dirichlet walls has none and a blanket dense SVD coarse solve cost most of ~0.8 s per V-cycle application in #622. Adds a snippet to check the hierarchy you actually built (verified to run) and a rule of thumb: low thousands of DOFs at the bottom, and if it reaches five figures add a level. docs-build: succeeded, no warnings on this page. Underworld development team with AI support from Claude Code --- docs/advanced/multigrid-preconditioning.md | 138 +++++++++++++++++++-- 1 file changed, 127 insertions(+), 11 deletions(-) diff --git a/docs/advanced/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md index e85d9c9f9..3c713e4ab 100644 --- a/docs/advanced/multigrid-preconditioning.md +++ b/docs/advanced/multigrid-preconditioning.md @@ -71,7 +71,8 @@ pc_mg_galerkin = "both" # RAP coarse operators (see below) mg_levels_ksp_type = "chebyshev" mg_levels_pc_type = "sor" mg_levels_ksp_max_it = 4 -mg_coarse_pc_type = "lu" # direct coarse solve +mg_coarse_pc_type = "redundant" # direct coarse solve, on every rank +mg_coarse_redundant_pc_type = "lu" ``` **Why Galerkin coarse operators?** Underworld3 does not install @@ -92,23 +93,138 @@ This is the key reason to prefer FMG when you adapt: GAMG on the next solve, so nothing breaks — you simply lose the geometric path until a hierarchy is available again. -## Parallel coarse solve +## The coarse level must be small, or the bottom solve stops being free -The default coarse solver is a serial direct solve (`mg_coarse_pc_type = "lu"`), -which is the fast, simple choice in serial and for modest core counts. For large -parallel partitions, replicate or gather the (small) coarse grid instead: +Multigrid is cheap because each level costs a fraction of the one above it, and +the bottom of the hierarchy is solved **exactly** at negligible cost. The word +doing the work in that sentence is *negligible*. A direct (or SVD) coarse solve +is only free if the coarsest level is genuinely tiny; nothing in the method makes +it so, and nothing warns you when it is not. + +### What the default does, and why + +The bundle sets a redundant direct solve: ```python -# redundant: copy the coarse grid to every rank and solve it there -stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "redundant" -stokes.petsc_options["fieldsplit_velocity_mg_coarse_redundant_pc_type"] = "lu" +mg_coarse_pc_type = "redundant" +mg_coarse_redundant_pc_type = "lu" +``` + +`redundant` copies the whole coarse system to **every rank** and factors it +there. For a genuinely small system that is the right trade: communicating +across ranks to solve a few thousand unknowns collectively costs more than +solving the same thing everywhere. It is also np-safe, which a bare `lu` is not +— a serial LU cannot factor a distributed matrix and fails at `np > 1` with +`DIVERGED_LINEAR_SOLVE` after zero iterations. In serial, `redundant` + `lu` is +identical to `lu`. + +This is not an Underworld idiosyncrasy, and it is worth being clear about that +before blaming the choice for a memory problem. **PETSc makes the same choice on +its own.** With nothing set under `mg_coarse_*`, `PCSetUp_MG` selects +`PCREDUNDANT` when the communicator has more than one rank and `PCLU` when it +does not — its own source comment reads "coarse solve is (redundant) LU by +default" (`src/ksp/pc/impls/mg/mg.c`). Confirmed by inspection at `np` = 1, 2 +and 4 on PETSc 3.25. + +Underworld sets it explicitly rather than inheriting it, so a bundle applied +after another cannot pick up a sibling's leftover options (#468) — but the value +is the one PETSc would have picked anyway. Deleting these keys and letting PETSc +decide would change nothing here; the leverage is in the hierarchy, not the +bottom solve. + +### The failure mode: replicated memory + +Because the factorization is replicated, its memory is paid **per rank**, and +the total scales with the coarse system size *times* the core count. That is +harmless for a small bottom level and fatal for a large one. + +Measured on Setonix (issue #644), a production model at `np = 30` under a +900 MB/rank cgroup cap: + +| Base mesh | Refinement | Coarsest velocity DOFs | Peak RSS/rank | Linear solve | +|---|---|---|---|---| +| 10 km | 2 | 14,195 | **683 MB** — OOM killed | ~100 s | +| 20 km | 3 | ~1,900 | **416 MB** | 111 s (+11%) | + +Both configurations have essentially the same *finest* grid (233k vs 227k +velocity DOFs) and converge identically — the same SNES iteration sequence. The +only real difference is how far down the hierarchy goes. Shrinking the coarsest +level alone removed 261 MB/rank, which matches the gather-plus-fill arithmetic +(~45 MB for the gathered copy, ~232 MB for the factor at a 5.16x fill ratio). + +### The fix is the hierarchy, not the coarse solver + +The instinct on an OOM is to reach for a different bottom solve. That treats the +symptom. A 14,000-DOF coarsest level means the hierarchy was built from too fine +a base with too few levels — the bottom of the ladder never got small, so the +"solve the coarse problem exactly, for free" premise never held. + +Keep the finest resolution you need and add levels underneath it: + +```python +# OOM: fine base, few levels -> the bottom is still large +mesh = uw.meshing.UnstructuredSimplexBox(..., cellSize=0.10, refinement=2) + +# same finest grid, one more level, a much smaller bottom +mesh = uw.meshing.UnstructuredSimplexBox(..., cellSize=0.20, refinement=3) +``` + +Each refinement level divides cell size by two, so in 2-D the level above has +about 4x the cells (8x in 3-D). Doubling the base cell size and adding one level +therefore lands on the same finest grid while making every coarser level — the +bottom one included — about 4x (or 8x) smaller. + +Check what you actually built rather than assuming; `dm_hierarchy[0]` is the +coarsest level and `[-1]` is the mesh you solve on: + +```python +for i, dm in enumerate(mesh.dm_hierarchy): + cells = uw.mpi.comm.allreduce(dm.getStratumSize("depth", mesh.dim)) + uw.pprint(f"level {i}: {cells} cells") +``` -# or telescope onto a sub-communicator for very large runs -# stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "telescope" +As a rule of thumb, aim for a coarsest level in the low thousands of DOFs or +fewer. If it runs to five figures, add a level rather than changing the bottom +solve. + +### When the coarse operator has a null space: SVD + +One case does require overriding the bottom solve. A **rotated** velocity block +whose domain leaves rigid rotations free — a closed circle (one mode), a +spherical shell (three) — hands those modes to the Galerkin-coarsened coarse +operator, because the rotations survive the constraint. `redundant`/LU hits a +zero pivot on a singular coarse operator and fails with `SUBPC_ERROR` (outer +reason `-11`). SVD is null-space robust: + +```python +from underworld3.utilities.multigrid_options import geometric_mg_bundle + +geometric_mg_bundle(coarse="svd") # null-space-robust bottom solve +``` + +The solver already does this for you, and it is deliberately **not** a blanket +choice for every rotated problem. It keys on the count of *verified rotation* +modes, because a rotated problem with Dirichlet walls, or a split-fault contact +box, has none — and there an SVD is pure cost. A 3-D P2 coarse level is a +**dense** factorization: measured, it accounted for most of ~0.8 s per V-cycle +application at healthy iteration counts (#622). + +Two things follow. Do not reach for `svd` because a problem is "rotated"; reach +for it when the bottom solve reports a zero pivot. And because SVD is dense, it +is affordable only on a genuinely small coarse level — which makes the sizing +guidance above a hard requirement on this path rather than a preference. + +### Very large partitions: telescope + +For core counts where even a small coarse system is replicated too many times, +gather it onto a sub-communicator instead of onto every rank: + +```python +stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "telescope" ``` Because `"auto"` does not overwrite options you set yourself, you can keep -`preconditioner = "auto"` and just override the coarse solver as above. +`preconditioner = "auto"` and override just the coarse solver. ## Benchmark: FMG vs GAMG on a deforming adaptive mesh From ffd8920f3db1aca0766db02566af055ef942092f Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 1 Sep 2026 20:18:27 -0700 Subject: [PATCH 2/2] Adversarial review of my own #666: five defects, two of them wrong advice Reviewed the docs change by running its claims rather than re-reading them. Five findings, fixed here. 1. TELESCOPE CANNOT WORK HERE. The page offered `mg_coarse_pc_type = "telescope"` for large partitions, and my rewrite promoted it to its own subsection. PCTELESCOPE has no DMPlex support and every UW3 mesh is a DMPlex, so the snippet aborts the solve: "Support for DMPLEX is currently not available", PETSc error 56, reproduced at np=2. It is now a warning that it is unavailable - which strengthens the point, since it means there is NO bottom-solve escape hatch from an over-large coarse level. 2. THE SVD SNIPPET WAS INERT. `geometric_mg_bundle(coarse="svd")` builds an MGSettings and discards it; nothing is applied. It read as a user action and did nothing. Removed - the solver already selects SVD from the verified rotation-mode count. 3. THE OVERRIDE PROMISE IS FALSE ON THE ROTATED PATH - the one path that wants SVD. That route writes its bundle under the velocity sub-PC's own prefix, applies it, then deletes the keys; `solver.petsc_options[...]` sits under the SNES prefix and never reaches it. Verified the override DOES survive on the standard path (fieldsplit_velocity_mg_coarse_pc_type = svd after solve), so the promise is now scoped to that path and the rotated exception is a warning. 4. "LANDS ON THE SAME FINEST GRID" OVERCLAIMED. Measured, cellSize 0.10/ref2 gives 3,872 finest cells against 0.20/ref3 at 4,224 - 9% MORE, not equal, while the coarsest drops to 27%. #644's own numbers go the other way (227k vs 233k). Now stated as a few percent either side, with both measurements. 5. THE RULE OF THUMB HAD NO np OR MEMORY BUDGET IN IT, yet those are what set the threshold - the replicated factor is paid per rank against a per-rank cap. Rewritten around that arithmetic, calibrated on #644's 14k DOFs -> ~300 MB/rank at 5.16x fill. Claims that survived the review, all verified rather than asserted: the dm_hierarchy snippet works on a bare mesh (no solver needed); the per-level growth is exactly 4x in 2-D and 8.00x in 3-D (66/264/1056, 184/1472/11776); PETSc's own source picks redundant above one rank and lu at one. docs-build: succeeded, no warnings on this page. Underworld development team with AI support from Claude Code --- docs/advanced/multigrid-preconditioning.md | 97 ++++++++++++++-------- 1 file changed, 62 insertions(+), 35 deletions(-) diff --git a/docs/advanced/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md index 3c713e4ab..b34a5b9c1 100644 --- a/docs/advanced/multigrid-preconditioning.md +++ b/docs/advanced/multigrid-preconditioning.md @@ -169,10 +169,19 @@ mesh = uw.meshing.UnstructuredSimplexBox(..., cellSize=0.10, refinement=2) mesh = uw.meshing.UnstructuredSimplexBox(..., cellSize=0.20, refinement=3) ``` -Each refinement level divides cell size by two, so in 2-D the level above has -about 4x the cells (8x in 3-D). Doubling the base cell size and adding one level -therefore lands on the same finest grid while making every coarser level — the -bottom one included — about 4x (or 8x) smaller. +Each refinement level divides cell size by two, so the level above has 4x the +cells in 2-D and 8x in 3-D — measured exactly, 66 / 264 / 1056 and +184 / 1472 / 11776 on a `refinement=2` box. Doubling the base cell size and +adding one level therefore lands on the same finest grid while making every +coarser level — the bottom one included — about 4x (or 8x) smaller. + +"The same finest grid" is approximate, because the base mesh is unstructured and +gmsh does not produce exactly 4x the cells for half the target size. Measured on +a unit box, `cellSize=0.10, refinement=2` gives 3,872 finest cells against +`cellSize=0.20, refinement=3` at 4,224 — 9% more, while the coarsest level drops +to 27% of its former size. The production model in #644 landed the other way: +227k against 233k. Expect a few percent either side, and check rather than +assume. Check what you actually built rather than assuming; `dm_hierarchy[0]` is the coarsest level and `[-1]` is the mesh you solve on: @@ -183,48 +192,66 @@ for i, dm in enumerate(mesh.dm_hierarchy): uw.pprint(f"level {i}: {cells} cells") ``` -As a rule of thumb, aim for a coarsest level in the low thousands of DOFs or -fewer. If it runs to five figures, add a level rather than changing the bottom -solve. +There is no single DOF number that is "small enough", because the quantity that +has to fit is the replicated factor **per rank**, against your per-rank memory +budget — the core count is in the arithmetic, not just the coarse size. The #644 +measurement calibrates it: ~14k velocity DOFs cost ~300 MB/rank at a 5.16x LU +fill, so a 900 MB/rank cap is already lost before the rest of the solver is +counted, while ~1.9k DOFs cost little enough to disappear. -### When the coarse operator has a null space: SVD - -One case does require overriding the bottom solve. A **rotated** velocity block -whose domain leaves rigid rotations free — a closed circle (one mode), a -spherical shell (three) — hands those modes to the Galerkin-coarsened coarse -operator, because the rotations survive the constraint. `redundant`/LU hits a -zero pivot on a singular coarse operator and fails with `SUBPC_ERROR` (outer -reason `-11`). SVD is null-space robust: +As a working rule, keep the coarsest level in the low thousands of DOFs, and +treat five figures as a bug in the hierarchy rather than a tuning question. If +you are near the line, estimate it rather than guessing: the factor grows faster +than the DOF count, and you pay it on every rank at once. -```python -from underworld3.utilities.multigrid_options import geometric_mg_bundle +### When the coarse operator has a null space: SVD -geometric_mg_bundle(coarse="svd") # null-space-robust bottom solve +One case does require a different bottom solve, and the solver handles it for +you. A **rotated** velocity block whose domain leaves rigid rotations free — a +closed circle (one mode), a spherical shell (three) — hands those modes to the +Galerkin-coarsened coarse operator, because the rotations survive the +constraint. `redundant`/LU hits a zero pivot on a singular coarse operator and +fails with `SUBPC_ERROR` (outer reason `-11`), so the rotated path selects SVD +instead. + +It is deliberately **not** a blanket choice for every rotated problem. It keys +on the count of *verified rotation* modes, because a rotated problem with +Dirichlet walls, or a split-fault contact box, has none — and there an SVD is +pure cost. A 3-D P2 coarse level is a **dense** factorization: measured, it +accounted for most of ~0.8 s per V-cycle application at healthy iteration counts +(#622). + +```{warning} +There is no `petsc_options` override for the coarse solve **on the rotated +path**. That route writes its multigrid bundle under the velocity sub-PC's own +prefix, applies it, and deletes the keys again; anything you set through +`solver.petsc_options[...]` lives under the solver's SNES prefix and never +reaches it. The override below works on the standard path only. ``` -The solver already does this for you, and it is deliberately **not** a blanket -choice for every rotated problem. It keys on the count of *verified rotation* -modes, because a rotated problem with Dirichlet walls, or a split-fault contact -box, has none — and there an SVD is pure cost. A 3-D P2 coarse level is a -**dense** factorization: measured, it accounted for most of ~0.8 s per V-cycle -application at healthy iteration counts (#622). +Because SVD is dense, it is affordable only on a genuinely small coarse level — +which makes the sizing guidance above a hard requirement on this path rather +than a preference. -Two things follow. Do not reach for `svd` because a problem is "rotated"; reach -for it when the bottom solve reports a zero pivot. And because SVD is dense, it -is affordable only on a genuinely small coarse level — which makes the sizing -guidance above a hard requirement on this path rather than a preference. +### Overriding the coarse solve yourself -### Very large partitions: telescope - -For core counts where even a small coarse system is replicated too many times, -gather it onto a sub-communicator instead of onto every rank: +On the standard (non-rotated) path, `"auto"` does not overwrite options you set +yourself, so you can keep `preconditioner = "auto"` and change just the bottom +solve: ```python -stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "telescope" +stokes.preconditioner = "auto" +stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" ``` -Because `"auto"` does not overwrite options you set yourself, you can keep -`preconditioner = "auto"` and override just the coarse solver. +```{note} +`PCTELESCOPE` — the usual PETSc answer to "replicating onto every rank is too +many copies", which gathers the coarse system onto a sub-communicator instead — +**is not available here**. It has no DMPlex support, and every Underworld mesh +is a DMPlex, so requesting it aborts the solve with `Support for DMPLEX is +currently not available` (PETSc error 56). There is no bottom-solve escape hatch +from an over-large coarse level: the hierarchy is the lever. +``` ## Benchmark: FMG vs GAMG on a deforming adaptive mesh