diff --git a/docs/advanced/multigrid-preconditioning.md b/docs/advanced/multigrid-preconditioning.md index e85d9c9f9..b34a5b9c1 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,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