Coalesced Matrix-Free Geometric Multigrid on Persistent Cell-Wise Storage

Abstract: We present a geometric multigrid preconditioner for high-order continuous finite elements that operates entirely on redundant, cell-wise stored vectors: the assembled global vector is never formed on any level of the hierarchy. In this storage paradigm the machinery that classically complicates adaptive multigrid dissolves. Hanging-node constraints are never assembled: we prove that the plain tensor-product transfer operators, applied to the \emph{unassembled} residual, algebraically reproduce the classical constrained restriction, including the action of the transposed constraint matrix, and the edge operators of local smoothing reduce to a pointwise masking of the residual, with no splitting of the level operator into interior and edge blocks. As a consequence, the single inter-cell primitive of the whole V-cycle can use a topologically oblivious structured kernel even on adaptively refined meshes. We prove that the resulting cell-wise V-cycle is equivalent, iterate by iterate, to the classical local multigrid method, and therefore inherits its convergence theory. Numerical experiments for the Laplace operator confirm grid-independent convergence that is essentially unaffected by local refinement; on a single GPU, using nothing more than a masked point-Jacobi smoother, the solver sustains up to $1.1$\, GDoF/s per V-cycle in double precision and reaches end-to-end solve throughput on par with patch-smoother-based solvers.

• The paper introduces a matrix-free GMG preconditioner that eliminates global assembled vectors through persistent cell-wise storage.
• It employs cell-wise prolongation and restriction to achieve grid-independent convergence and bandwidth-limited performance on GPUs.
• The method removes traditional hanging-node constraint handling, streamlining adaptive refinement in high-order finite element solvers.

Coalesced Matrix-Free Geometric Multigrid on Persistent Cell-Wise Storage

Introduction and Motivation

The paper "Coalesced Matrix-Free Geometric Multigrid on Persistent Cell-Wise Storage" (2607.03413) innovates within matrix-free high-order finite element solvers by presenting a geometric multigrid (GMG) preconditioner that leverages redundant, persistent cell-wise storage throughout the entire hierarchy. The key architectural advance is the elimination of all explicit global assembled vectors and the classical machinery for handling hanging-node constraints on adaptive, non-conforming meshes. The method operates directly on unassembled, cell-local data at all stages—operator evaluation, inter-grid transfer, smoothing, and Krylov acceleration—such that only a single inter-cell primitive (direct stiffness summation, DSS) mediates all required communication.

This paradigm genuinely alters the complexity and implementation of high-performance adaptive multigrid on modern hardware, particularly GPUs, where memory bandwidth and locality dominate the cost model. The paper establishes both theoretical equivalence with classical local multigrid preconditioners and demonstrates grid-independent iteration counts and near-peak bandwidth-limited throughputs in large-scale tests.

Cell-Wise Storage and Algebraic Structure

The persistent cell-wise storage strategy retains all field representations as redundant vectors with one entry per local degree of freedom (DoF) per cell. No assembled global vector is ever formed. Crucially, this storage induces a primal-dual structure:

• Primal objects (e.g., iterates, corrections) correspond to $C^0$-continuous fields, whose multiple per-cell copies are synchronized.
• Dual objects (e.g., element-local residuals) are nonsynchronized local integral sums, physically meaningful without global assembly.

Formally, the gather operator $\mathcal{G}$ maps an assembled vector to its redundant cell-local representation; its transpose $\mathcal{G}^T$ performs additive assembly. Operator applications, vector algebra, and Krylov dot products are defined entirely via local operations or (when required) via the pairing identity $$\langle \tilde{f}, u \rangle = \langle \mathcal{G}^T \tilde{f}, u_\text{global} \rangle$$.

This structure supports full matrix-free implementation with no need for scatter/gather lists, atomic operations, or indirect addressing—all bottlenecks in high-performance GPU execution.

Multigrid Hierarchy and Redundant Inter-Grid Transfers

The method constructs a multigrid hierarchy through nested mesh and function spaces $(\mathbb{V}_\ell)$ aligned with adaptive refinement, without ever assembling global system matrices. Algorithmically, the V-cycle proceeds using:

• Cell-wise prolongation via tensor-product embeddings, respecting local cell orientation.
• Cell-wise restriction as the adjoint, which when applied to unassembled, cell-local duals, algebraically recovers the effect of classical restriction composed with constraint application—crucial for non-conforming interfaces.

Critically, the following operator identities hold:

$$\mathcal{P}_{\text{cell}, \ell-1}^\ell \mathcal{G}_{\ell-1} = \mathcal{G}_\ell \mathcal{P}_{CG, \ell-1}^\ell, \qquad \mathcal{G}_{\ell-1}^T \mathcal{R}_{\text{cell}, \ell}^{\ell-1} = \mathcal{R}_{CG, \ell}^{\ell-1} \mathcal{G}_\ell^T,$$

ensuring that the local, structure-oblivious transfer operations commute with the abstract global ones.

Treatment of Hanging Nodes and Edge Constraints

The major conceptual advance is the abolition of constraint handling for hanging nodes. Local smoothing on adaptively refined meshes is realized by exact masking on the refinement edge $E_\ell$:

• The smoother acts only on interior DoFs (Dirichlet constraint on $E_\ell$).
• The edge residual, required for coarse grid correction, is evaluated fully via unassembled operations.
• No constraint matrix or explicit interface resolution is required.

The approach is justified by analysis showing that the resulting method is algebraically equivalent to the classical local-multigrid with explicit constraint enforcement, inheriting its convergence guarantees.

Smoother Implementation and DSS Realization

The core primitive for enforcing continuity post-smoothing is direct stiffness summation (DSS):

• DSS assembles all partial sums for shared DoFs by locality-aware, memory-coalesced communication, completely avoiding atomics and connectivity traversals.
• On structured macroblocks, DSS is realized as sequential dimensionally-split passes: each pass assembles one direction, leveraging bandwidth-limited streaming.
• For inter-block boundaries or highly adaptive regions, a complementary set of kernels perform minimal indirect communication, with an impact proportional only to mesh surface area complexity.

The actual smoother is damped Jacobi, but the architectural framework admits any (potentially $p$-robust) overlapping Schwarz or patch smoothers, provided their corrections can be represented as cell-wise conforming fields and assembled by DSS.

Numerical Results: Convergence, Throughput, and Scalability

Several significant performance and convergence benchmarks are reported:

• Grid-independent convergence: V-cycle contraction factors are nearly unchanged both on uniform and adaptively refined (with deep hanging-node interfaces) meshes, with contraction rates ranging from 0.36 ($p=2$) to 0.50 ($\mathcal{G}$0) per cycle.
• Throughput: On A100 GPUs in double precision, sustained rates reach up to $\mathcal{G}$1 GDoF/s per V-cycle and over $\mathcal{G}$2 MDoF/s end-to-end solve rate (including Krylov acceleration) at $\mathcal{G}$3, with negligible loss on adaptive/hanging-node grids.
• Bandwidth-bound execution: Roofline plots and kernel analysis confirm pure bandwidth-boundedness, with all major kernels (evaluation, restriction, prolongation, smoother) saturating the achievable memory bandwidth. Further increasing arithmetic intensity yields no practical throughput improvement.

Figure 1: The tensor-product mesh and cell partitioning underpinning coalesced cell-wise memory and block-structured DSS; adaptive refinement is handled without special case code.

Theoretical and Practical Implications

The work fundamentally demonstrates that:

• All algebraic machinery for mesh conformity/hanging nodes can be subsumed by the redundancy and properties of tensor-product transfer operators and DSS, provided cell-wise storage is used throughout.
• This supports implementation approaches that are completely topologically oblivious, with all adaptive, unstructured, and highly refined regions handled through a uniform pipeline.
• On hardware with high memory bandwidth and thread-level parallelism (GPUs, exascale CPUs), the method is optimal in terms of machine utilization, shifting performance considerations to data movement rather than complex kernel fusion or threading strategies.

The ability to plug in advanced smoothers without changing the storage or communication strategy suggests straightforward extension to $\mathcal{G}$4-robust and robustly preconditioned high-order discretizations, and to more complex operator families (e.g., $\mathcal{G}$5, $\mathcal{G}$6), as well as to hybrid/mixed and unfitted methods.

Future Directions

Key avenues for further investigation and development include:

• Hybridizing with locality-optimized smoothers to reduce geometry streaming and further collapse the per-DoF bandwidth cost, particularly essential on high-order curved meshes.
• Reducing working precision (e.g., fp32 or fp16) now that correctness is uncoupled from arithmetic precision, pending stability analysis.
• Deploying on large-scale, distributed architectures, leveraging the inherent communication locality and kernel fusion this method supports.
• Extending constraint- and valence-free approach to nonconforming methods (DG, CutFEM, SB methods) by adapting the cell-wise inter-grid transfer abstractions.

Conclusion

This paper establishes a rigorous, efficient, and highly practical matrix-free multigrid methodology for high-order FEM with adaptive refinement, demonstrating that persistent cell-wise storage obviates all classical constraint handling and supports peak hardware utilization. The confluence of algebraic equivalence to classical methods, architectural minimalism, and empirical performance situates this method as a new baseline for GPU-accelerated, high-fidelity, adaptive finite element simulations.

---

Reference:

M. Wichrowski, "Coalesced Matrix-Free Geometric Multigrid on Persistent Cell-Wise Storage," (2607.03413).