FMG geometric-MG coarse solve uses PCREDUNDANT: coarse LU replicated on every rank — OOM at scale when the base mesh is fine

[bknight1]

Problem

When using FMG with the default redundant coarse solver (pc_type=redundant + pc_factor_mat_solver_type=lu), the coarse-grid LU factorization is replicated on every MPI rank. For models with large base meshes (e.g., 10 km cells, 800x230 km domain), the coarsest level still has ~14k velocity DOFs, causing:

• ~300 MB/rank replicated memory from the coarse LU factorization (fill ratio 5.16× on 14k DOFs)
• Total ~9 GB equivalent replicated across 30 ranks
• OOM kills at 900 MB/rank cap on Setonix/Pawsey

Root Cause

The multigrid_options.py default geometric_mg_bundle() sets mg_coarse_pc_type=redundant (see multigrid_options.py:262). This is documented as "(small) coarse system" but fails when base mesh resolution is high.

Related: PR #213 (_use_iterative_solver) already avoids MUMPS in parallel for the same reason, using redundant+svd instead — but svd is infeasible at this scale (dense 14k×14k).

Solution Found

Using a finer base mesh with more refinement levels keeps the finest grid identical while shrinking the coarsest level:

# Original (OOM)
-uw_csize_km 10 -uw_refinement 2   # 14k coarse DOFs → ~300 MB/rank replicated

# Fixed (fits in 900 MB/rank cap)
-uw_csize_km 20 -uw_refinement 3   # ~1.9k coarse DOFs → divided / minimal replication

Config	Coarsest vDOFs	Linear solve time (serial)	Peak RSS/rank (np=30)
csize10 ref2	14,195	~100 s	683 MB ❌ OOM at 900 MB cap
csize20 ref3	~1,900	111 s (+11%)	416 MB ✅

• Finest grid essentially unchanged (233k vs 227k velocity DOFs)
• Trade-off: +11 % linear-solve time buys −39 % per-rank peak RSS
• Causal check: shrinking only the coarsest level removed 261 MB/rank, matching the redundant gather + LU-fill arithmetic (~45 MB gathered copy + 232 MB factor ≈ 280 MB)

Proposed Fix Options

1. Size-aware default in multigrid_options.py (preferred):
   • Query coarsest-level DOF count from mesh.dm_hierarchy
   • Keep redundant+lu only below a threshold (≈3–5k coarse DOFs)
   • Above threshold use a non-replicating bottom solve (gamg, or asm) — note distributed MUMPS must stay excluded per PR Parallel adaptive seam-spike fix: mover heap corruption + locator hardening + remesh field-transfer redesign (Phase 1 + 2) #213's finding that MUMPS-in-parallel corrupts in this build
2. Document the workaround: publish the coarser-base + extra-refinement recipe and warn on large-coarse redundant
3. Expose the threshold as a solver option for user control

Related

• PR Parallel adaptive seam-spike fix: mover heap corruption + locator hardening + remesh field-transfer redesign (Phase 1 + 2) #213 / docs/developer/design/parallel-repeated-solve-corruption.md: MUMPS-in-parallel heap corruption in this build (why mumps-based coarse solves are excluded here)

Verification

Setonix (Pawsey), production-geometry model, np=30, per-rank RSS sampled every 2 s:

• Control (csize10 ref2): post-solve 683 MB/rank uniform across all 30 ranks (spread 13 MB) → step killed at the 900 MB/task cgroup cap
• Treatment (csize20 ref3): post-solve 398–427 MB/rank (median 416), window max 576 MB, clean completion
• Convergence identical between arms (same SNES iteration sequence)

[lmoresi]

I don't this is a bug. Redundant solve means we do that everywhere the same. Why would we do that ? Because the problem is so tiny that it's not worth communicating across cores to solve it collectively. Make sure you are using a very coarse mesh as your base or you will run into trouble.

[bknight1]

I agree not a bug per se, but probably needs documenting in the fmg docs that it can cause an OOM error

[lmoresi]

Documented in #666.

Two things came out of writing it that are worth having on this thread.

The choice of redundant is PETSc's, not ours. PCSetUp_MG selects PCREDUNDANT above one rank and PCLU at one rank — its 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. We set the keys explicitly only so a bundle cannot inherit a sibling's leftover options (#468); deleting them and letting PETSc decide would give the same solver. So option 1's alternative bottom solves would be overriding PETSc's default rather than correcting an Underworld one — worth knowing before spending effort there.

The page was actively misleading, in the direction that made this harder to diagnose. It stated the default coarse solve was a bare lu, in both the bundle listing and the parallel section, and then told users to set redundant themselves for large partitions. Anyone reading it would have concluded redundant was opt-in and not suspected replicated memory. That is corrected.

The framing the docs now lead with is your own finding rather than a workaround: multigrid is only cheap because the bottom level is solved exactly at negligible cost, which holds only if the bottom is genuinely small. A 14k-DOF coarsest level means the hierarchy was built from too fine a base with too few levels, so the premise never held — and csize20 ref3 is not a workaround, it is the correctly-built hierarchy for that finest grid. Your measurement makes the point cleanly: same finest grid, same SNES iteration sequence, 683 → 416 MB/rank for +11%.

The size-aware warning from option 1 is not in the PR, since the call here was documentation. It still looks worth doing as a follow-up — the check is cheap, and the failure currently presents as an OOM kill with no hint of the cause. Happy to open that separately if you want it.

Underworld development team with AI support from Claude Code