From 5902ba3110c4c9eaa61f8afb347b7259f93379a4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Apr 2026 14:17:46 +1000 Subject: [PATCH] Default Stokes velocity subsolve to FGMRES (#147) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the FGMRES fieldsplit configuration that gthyagi diagnosed and validated at scale on Gadi (issue #147). Credit for finding the underlying CG/FCG / standard-GMRES failure mode in the spherical Kramer free-slip Nitsche benchmark, isolating the trigger conditions, and validating FGMRES as the working configuration goes to gthyagi; this commit only promotes that configuration to the UW3 default. Changes the default fieldsplit velocity KSP from CG to FGMRES across both the SNES_Stokes_SaddlePt constructor and the strategy setter (used when solver.strategy is reassigned). Adds an explicit ksp_max_it=200 safety bound to both pressure and velocity subsolves. Why --- The Schur fieldsplit preconditioner uses inner Krylov solvers (GAMG with mg_levels_ksp_max_it=3 and _converged_maxits=true) whose application is non-linear. CG/FCG and standard GMRES build their residual recurrence assuming a linear preconditioner, so the recurrence drifts from the true residual under non-stationary preconditioning. Combined with the weakly-indefinite coarse operators that GAMG produces on partition-stressed problems with sharp internal sources (Kramer case1's internal-boundary delta forcing), the velocity block trips PETSc's "DIVERGED_PC_FAILED / indefinite matrix" check or the GMRES residual-recursion mismatch. FGMRES handles non-stationary preconditioning by storing the preconditioned vectors explicitly rather than relying on the m-term recurrence; it is the right default for any solver that composes nested Krylov methods, independent of bilinear-form symmetry. gthyagi's at-scale validation (#147) showed FGMRES completing at np=144, cellsize=1/32 on Gadi where CG/FCG and standard GMRES both failed. Locally re-validated: pytest -m "level_1 and tier_a" on amr-dev passes 56/3/0. Compatibility ------------- - pc_side is intentionally not set: FGMRES is right-preconditioned by construction, and not setting pc_side keeps the defaults compatible with users who override the velocity KSP to CG/FCG (which require left preconditioning). - The 8 existing tests that explicitly set fieldsplit_velocity_ksp_type=fcg after solver construction continue to work — user overrides cleanly replace the default. - solver.strategy = "default" re-applies the new FGMRES defaults via the same code path. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 31 +++++++++++++++---- 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index e1a965ffc..e5c483e3b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3990,8 +3990,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): p_name = "pressure" # pressureField.clean_name v_name = "velocity" # velocityField.clean_name - # Works / mostly quick + # Pressure subsolve — flexible GMRES + GASM. (FGMRES is right- + # preconditioned by construction; no need to set pc_side explicitly.) self.petsc_options[f"fieldsplit_{p_name}_ksp_type"] = "fgmres" + self.petsc_options[f"fieldsplit_{p_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{p_name}_ksp_rtol"] = self._tolerance self.petsc_options[f"fieldsplit_{p_name}_pc_type"] = "gasm" # self.petsc_options[f"fieldsplit_{p_name}_pc_gasm_type"] = "basic" @@ -4003,8 +4005,20 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # self.petsc_options[f"fieldsplit_{p_name}_pc_gamg_type"] = "agg" # self.petsc_options[f"fieldsplit_{p_name}_pc_gamg_repartition"] = True - # Great set of options for gamg - self.petsc_options[f"fieldsplit_{v_name}_ksp_type"] = "cg" + # Velocity subsolve — flexible GMRES + GAMG. + # + # FGMRES (not CG/FCG) is the default to remain robust under non-stationary + # preconditioning and weakly-indefinite coarse operators. The mg_levels + # KSP is bounded but variable-iteration (mg_levels_ksp_converged_maxits) + # so the GAMG application is non-linear; CG/FCG's residual recurrence + # cannot accommodate this, and at scale (large parallel partitions, sharp + # internal sources, free-slip Nitsche) PETSc reports DIVERGED_PC_FAILED + # / indefinite matrix or GMRES residual-recursion mismatch. See issue + # #147 for the spherical-Kramer benchmark on Gadi that drove this change; + # gthyagi's validated FGMRES configuration completed at np=144, + # cellsize=1/32 where CG/FCG and standard GMRES both failed. + self.petsc_options[f"fieldsplit_{v_name}_ksp_type"] = "fgmres" + self.petsc_options[f"fieldsplit_{v_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{v_name}_ksp_rtol"] = self._tolerance * 0.1 self.petsc_options[f"fieldsplit_{v_name}_pc_type"] = "gamg" self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_type"] = "agg" @@ -4396,8 +4410,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): p_name = "pressure" # pressureField.clean_name v_name = "velocity" # velocityField.clean_name - # Works / mostly quick + # Pressure subsolve — flexible GMRES + GASM. (FGMRES is right- + # preconditioned by construction; no need to set pc_side explicitly.) self.petsc_options[f"fieldsplit_{p_name}_ksp_type"] = "fgmres" + self.petsc_options[f"fieldsplit_{p_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{p_name}_ksp_rtol"] = self._tolerance self.petsc_options[f"fieldsplit_{p_name}_pc_type"] = "gasm" self.petsc_options[f"fieldsplit_{p_name}_pc_gasm_type"] = "basic" @@ -4409,8 +4425,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # self.petsc_options[f"fieldsplit_{p_name}_pc_gamg_type"] = "agg" # self.petsc_options[f"fieldsplit_{p_name}_pc_gamg_repartition"] = True - - self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "cg" + # Velocity subsolve — see corresponding block in __init__ for the + # rationale (FGMRES default for robustness to non-stationary GAMG + # preconditioning and weakly-indefinite coarse operators; issue #147). + self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fgmres" + self.petsc_options[f"fieldsplit_velocity_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_velocity_pc_type"] = "gamg" self.petsc_options[f"fieldsplit_velocity_pc_gamg_type"] = "agg" self.petsc_options[f"fieldsplit_velocity_pc_gamg_repartition"] = True