Multigrid p-Robustness at Jacobi Speeds: Efficient Matrix-Free Implementation of Local p-Multigrid Solvers

TL;DR

This paper introduces a high-performance, matrix-free local p-multigrid solver that achieves robust convergence at Jacobi speeds, combining efficiency with effectiveness in high-order finite element methods.

Contribution

It presents a novel matrix-free implementation of vertex-patch smoothers using sum-factorization and SIMD, enabling efficient, scalable multigrid smoothing.

Findings

Achieves p-robust convergence comparable to patch-based smoothers

Maintains optimal memory scaling even on distorted meshes

Outperforms traditional methods in computational efficiency

Abstract

Vertex-patch smoothers are essential for the robust convergence of geometric multigrid methods in high-order finite element applications, yet their adoption is traditionally hindered by the prohibitive cost of solving local patch problems. This paper presents a high-performance, matrix-free implementation of a p-multigrid local solver that dismantles the trade-off between smoothing effectiveness and computational efficiency. We focus on the practical realization of this iterative approach, leveraging sum-factorization and explicit SIMD vectorization to minimize memory footprint and maximize arithmetic throughput. The performance analysis demonstrates that the solver effectively hides data-fetching latencies and maintains optimal $\mathcal{O}(p^d)$ memory scaling, even when dominated by geometric data on distorted meshes. The result is a robust smoother that rivals the execution speed of simple pointwise smoothers while preserving the convergence benefits of patch-based methods.