Towards Matrix-Free Patch Smoothers for the Stokes Problem: Evaluating Local p-Multigrid Solvers

TL;DR

This work targets robust, matrix-free smoothing for Stokes equations using vertex-patch domain decompositions. It introduces a local $p$-multigrid framework to solve patch problems and compares Braess-Sarazin, block GMRES, and Schur-based strategies, showing that a single local Braess-Sarazin iteration suffices to maintain global convergence in the presence of mesh distortion and large viscosity contrasts. The results demonstrate that exact local solves are not necessary for robustness, enabling efficient, scalable solvers suitable for high-order discretizations. The approach holds promise for extending to $H( ext{div})$-conforming discretizations and complex multiphysics applications, where matrix-free performance is critical.

Abstract

Vertex-patch smoothers offer an effective strategy for achieving robust geometric multigrid convergence for the Stokes equations, particularly in the context of high-order finite elements. However, their practical efficiency is often limited by the computational cost of solving the local saddle-point problems, especially when explicit matrix factorizations are not feasible. We explore a fully iterative, matrix-free-compatible approach to the local patch solve using $p$-multigrid techniques. We evaluate different local solver configurations: Braess-Sarazin and block-triangular preconditioners. Our numerical experiments suggest that the Braess-Sarazin approach is particularly resilient. We find that a single iteration of the local solver yields global convergence rates comparable to those obtained with exact local solvers, even on distorted meshes and in the presence of large viscosity jumps.

Figures (8)

• Figure 1: Workflow in a generic smoother application of a patch smoother for a single patch $j$. The three panels represent the main steps: gathering data from the global solution and right-hand side, computing the local residual, and applying the local correction and scattering it back to the global solution. Blue dots indicate interior DoFs of the patch, while red dots indicate all DoFs (including boundary) associated with each cell. The dashed rectangle highlights the patch interior. Figure taken from wichrowski2025pMGimplementaion.
• Figure 2: Patch-smoother data flow for a local p-multigrid solve (illustration for $p=3$ on the fine local level). The coarse level shown is $p=1$ and contains a single interior node. The diagram shows gather, evaluate (local residual), pre-smoothing, restriction to the coarse level, an exact coarse correction, prolongation-and-add, and scatter-add. A dashed arrow indicates proceeding to the next local multigrid iteration. For brevity we illustrate only pre-smoothing; post-smoothing is omitted. Figure taken from wichrowski2025pMGimplementaion.
• Figure 3: Illustration of the vertex-based patches used in our 2D experiments. The top row shows structured (Cartesian) patches, while the bottom row depicts unstructured (simplicial) patches. The distortion level $\delta$ represents the maximum displacement of each interior vertex relative to the local mesh size $h$.
• Figure 4: Number of FGMRES iterations (logarithmic y-axis) vs. mesh distortion for local block solver on structured patches in 2D (left) and 3D (right). Solid lines: velocity block $A$ is approximated via a single $p$-multigrid V-cycle; dashed lines: exact velocity inverse; dotted lines: reference.
• Figure 5: Number of FGMRES iterations (logarithmic y-axis) vs. mesh distortion for Braess--Sarazin smoother on structured patches. Both variants use one smoothing step per level. Solid lines: structured patch, dashed lines: simplicial patch, dotted lines: reference Poisson problem iteration counts (available for $p=3$ and $p=7$). Left subplot: 2D; right subplot: 3D. For $p=7$ the solver did not converge in some runs with $\delta=0.3$ on simplicial patches. Dotted lines: reference Poisson problem iteration counts (available for $p=3$ and $p=7$ in 2D)