Document why the multigrid coarse level must be small (#644) (#666)

* 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

* 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

Author: lmoresi

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,165 @@ 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
+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 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:
 
 ```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"
+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"
+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.
+
+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.
+
+### When the coarse operator has a null space: SVD
+
+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.
 ```
 
-Because `"auto"` does not overwrite options you set yourself, you can keep
-`preconditioner = "auto"` and just override the coarse solver as above.
+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
+
+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.preconditioner = "auto"
+stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"
+```
+
+```{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