Scalable linear solvers: plan, decisions and traceability #21

Tracking issue for the scalable linear-solver block (part of the master plan, #22; driver companion #23). Full design, measurements and per-stage records in the project documents design/solvers.md and design/solvers-*-status.md (kept outside git). This issue records the decisions, the stages and their PRs so that every PR of the block can reference it.

Last update: 2026-11-13 — S0 to S5 done (S5: HyperFEM #283, HyperFEMKernels #52 merged). Pending administrative items: HyperFEM release with S4-B/S5, README recipe of the hybrid Newton for electro-mechanics. MPI: R5 proposal (tracked in #23).

Goal and acceptance

Replace LU as the linear solver inside Newton so that the full solve! scales:

• Time: solve! ≥ ×5 vs Gridap+LU from ~200k DOFs (measured on EM, D4).
• Algorithmic scalability: Krylov iterations nearly mesh-independent (±30 % from 50k to 1M DOFs).
• Newton unchanged: same Newton iterations and solution as LU.
• GPU/MPI-compatible design; every feature with tests.

Decisions

• D1 Libraries (2026-10-23): PETSc via GridapPETSc in production (GAMG with rigid-body modes, BoomerAMG, fieldsplit, MUMPS, PCMG); AlgebraicMultigrid.jl for prototyping and CI; Krylov and matrix-free multigrid in Julia. The installed binary (GridapPETSc 0.5.8 → PETSc_jll 3.22.1) has everything the plan needs (checked in S0). No first-class pure-Julia solver backend.
• D2 Code location (2026-10-23): the contract is Gridap's LinearSolver interface; each piece lives where the knowledge it needs lives. HyperFEM (extension on GridapPETSc, short branches off main): S0′, S1, S2, S4-B. HyperFEMKernels: S3, S4-A, S5, exposed as Gridap LinearSolvers with the coarse solver injected (no PETSc dependency). Discarded: everything in HyperFEMKernels; PETSc SNES end to end.
• D3 Hardware (2026-10-23): criteria measured on the M5 Max laptop (128 GB, 18 threads) at 50k/200k/1M DOFs; LU only up to where it fits; cluster (Julia 1.13) only for MPI scaling and GPU.
• D4 Reference problems (2026-10-29, refined 2026-10-30): EM and hMRE; the ×5 criterion is measured on EM (plate P1, box Q2); hMRE enters through its 3D / vacuum cases and the u–φ blocks; hmre_cell (3D periodic cell) is the illustrative case of the gains. Visco secondary; aorta CFD out of the block.
• D5 Convergence criteria (2026-10-29): follow the third-medium note — relative with the consistent tangent, absolute (atol_projected, on the scale of the reactions) in projected mode, curvature criterion relative and only when the linear solve converged.
• D6 h-multigrid (2026-10-31): S4 split into S4-A (matrix-free p-multigrid) and S4-B (PETSc PCMG on a hierarchy obtained by uniform refinement of a coarse Gmsh mesh, for Q1/P1 models). S4-B starts with a test of Gridap's refine on Gmsh tetrahedra and hexahedra (gate). External review in design/solvers-s4b-review.md; its open questions are raised when S4-B starts.
• S2 (DS2-1…DS2-6, 2026-10-31): positive-real form (φ rows negated) + GMRES; one V-cycle per block; default split additive by measurement; block AMGSolver for CI.
• Gate 0 (2026-10-31): the GAMG default is the best tuning found; the iteration growth on the 3D cell comes from the ×100 stiffness contrast, not from AMG on P1 elasticity. GAMG unchanged.
• R4a (DR4-1…DR4-8, 2026-10-31/11-05): Eisenstat–Walker forcing term (choice 2 only, ηmax = η0 = 1e-2, never tighter than the solver's tolerance at the start of the Newton), default off; no automatic fallback to a direct solver; one-pass residual_and_jacobian! deferred; README recommends PETScSolver(U; reuse=3) + forcing=EisenstatWalker() for multi-field production runs; projected-mode tolerance is an open question for S5 (settled by DS5-4).
• S4-A (DS4A-1…DS4A-11, 2026-11-07/08): the caller passes the spaces of every level, HyperFEM provides level_spaces; mechanics (single field) only — per-block p-multigrid for (u, φ) deferred; acceptance on an isotropic unit cube (mech_cube), the slender box reported; coarse solver one GAMG V-cycle with rigid-body modes (LU/MUMPS optional, injected, no PETSc in the engine); GMRES by default, CG valid (symmetric cycle); Chebyshev degree 2, range 20 (by the smoother sweep); the Dirichlet projection reuses the last usable levels; Q3 → Q1 directly; LOR Q1 deferred. Released as 0.2.2 (0.2.1 had gone out with Second derivatives (∇∇) of the unknown, test and trial functions; fallback by triangulation; cheap update #47).
• S4-B (DS4B-1…DS4B-9, 2026-11-08): Gridap's refine(model, 2) for every element type (step 0 gate passed on Gmsh tetrahedra and hexahedra, affine and distorted); Dirichlet tags must be closed sets on the coarse mesh; plain models in the hierarchy, glues kept apart; prolongation free → free by insertion from the refinement rule; Galerkin coarse operators; pure-Julia V-cycle (option J) in HyperFEM, the PETSc PCMG route (P) optional (only needed for MPI); Chebyshev on inverted nodal blocks with a Lanczos estimate in the block inner product; with reuse only the coarse levels lag, the fine smoother follows the current tangent; multi-field through FieldSplitSolver with one V-cycle per block (reported, not acceptance).
• S5 (DS5-1…DS5-5, 2026-11-09/10): forcing term recommended in the README for h-multigrid / field split, defaults unchanged; linear_indefinite(ns) hook reported by the multigrid solvers, FieldSplitSolver, PETSc (KSP_DIVERGED_INDEFINITE_*) and CG in MatrixFreeSolver (Krylov.jl linesearch = true); with HybridTangent an indefinite report switches to the projected tangent at the same state, ConsistentTangent unchanged; projected-mode solves at ηmax floored by atol_projected (adopted: projected Newton iterations did not grow); HybridTangent(; stay_projected = 2) stays projected after the second indefinite report of a solve (2026-11-10); electro-mechanics: HessianRegularization(model::ElectroMechano) projects the total ∂²Ψ/∂F∂F (Maxwell term included), formulation in E, never in D; acceptance on the EM pull-in, the hMRE beam is not an instability case as it stands (iterative solvers fail first by conditioning).
• S3 (DS3-1…DS3-5, 2026-11-05): Krylov.jl; LaggedAssembled(solver; every = 3) (confirmed by measurement); LOR Q1 moved to S4-A; Krylov in Julia, PETSc only inside LaggedAssembled (MatShell in phase 3); Newton ↔ solver hooks stay in HyperFEM, HyperFEMKernelsExt adds the MatrixFreeSolver methods.

Stages

Execution order: S0 → S0′ → R1 → R2 → S1 → R3 → transient models → S2 → R4a → S3 → S4-A → S4-B → S5 (R stages tracked in #23).

Results so far

Laptop M5 Max, 6 Julia threads, PETSc serial, backend :compiled; per Newton iteration or solve! as stated. Full tables in design/solvers.md and the status documents; raw rows in benchmark/results/.

• Since the start of the plan (3D hMRE cell, 104k): HyperFEM v0.0.7 (Gridap + LU) ≈ 1125 s / ~76 GB → :compiled + fieldsplit 78 s / 3.8 GB (×14 time, ×20 peak memory).
• Iterations: P1 mechanics +39 % over 53k → 809k; em_plate_c +37 % over 71k → 1.08M; 3D cell +42 % over 29k → 480k, driven by the ×100 stiffness contrast (flat at contrast 1).
• Same Newton trajectory and solution as LU/MUMPS in every converged run.
• Matrix-free product vs SpMV on Q2 hexahedra: ×5 (57k–974k), ×6 on Q3; on P1 tetrahedra ×10 slower (P1 stays assembled; S4-B is its route).
• p-multigrid Krylov iterations on the cube Q2: 25.5 → 26.6 → 27.8 from 25k to 957k (+9 %); GAMG on the assembled Q2: 60 → 64 → 66.
• h-multigrid CG iterations on P1 tetrahedra: 32 at 414k and at 893k (GAMG 53 → 85); Q1 hexahedra 16.7–16.8 from 101k to 799k (GAMG 28 → 30). Memory does not discriminate on the assembled tangent (both dominated by the matrix).
• The GAMG references are PETSc in serial (no MPI in this build).

Open items carried forward

• Thin-structure elasticity AMG (K_uu cond ≈ 4e7 on em_plate/em_box): MUMPS (or lagged MUMPS matrix-free) until a dedicated task (semi-coarsening, line smoothers, BDDC).
• Monolithic GAMG returns a wrong update with a converged flag on the 3D cell.
• Interface-aware coarsening for high stiffness contrast (or S4-B with the interface on the coarse mesh, or BDDC).
• Multi-field (u, φ) p-multigrid: per-block p-multigrid inside a field split — now possible through FieldSplitSolver (S4-B), which accepts any LinearSolver per block.
• Block field splits are not h-independent on the EM plate (h-multigrid blocks 41 → 56 GMRES iterations from 71k to 555k, GAMG blocks 55 → 98): the coupling the block splits ignore.
• h-multigrid fine products are ~60 % of a CG iteration: symmetric (upper-triangle) storage or the matrix-free fine operator are the next levers; the PETSc PCMG route for MPI.
• p-multigrid cycle is coarse-limited (one GAMG V-cycle on Q1) and its setup is ≈ 23 % of a Newton iteration at 957k: a second coarse cycle and lagging the λmax estimate are the next levers.
• Rebinding a MatrixFreeSolver to an operator of another concrete type (review B8, pre-existing).
• Matrix-free on simplices needs the matrix-product formulation of the assembly kernels.
• Preconditioner application through PETScSolver.solve! converts vectors per call (MatShell/PETSc vectors in phase 3).

Traceability

Every PR of the block (in this repository or in HyperFEM) references this issue and names its stage (S1, S4-A, …). Measurements go in benchmark/ with their origin (machine, threads, Julia version). Decisions are recorded here when taken.

Generated by Claude Code