Shifting the Sweet Spot: High-Performance Matrix-Free Method for High-Order Elasticity

TL;DR

This work tackles the memory bottleneck in high-order elasticity simulations by developing a deeply optimized matrix-free PA operator integrated with Geometric Multigrid (GMG). It introduces a tensor-product, sum-factorization based $O(p^4)$ algorithm, augmented with Voigt symmetry, macro-kernel fusion, and slice-wise loop optimizations, and validates it within MFEM. Across x86 and ARM platforms, it achieves kernel speedups up to $83\times$ and end-to-end improvements up to $16.8\times$, while enabling large-scale solves with tens of millions of degrees of freedom where FA fails due to memory constraints. The results demonstrate a practical, scalable path for high-order elasticity on mainstream CPUs, bridging theory and practice and pointing to future directions such as mixed-precision computing and extension to more complex physics.

Abstract

In high-order finite element analysis for elasticity, matrix-free (PA) methods are a key technology for overcoming the memory bottleneck of traditional Full Assembly (FA). However, existing implementations fail to fully exploit the special structure of modern CPU architectures and tensor-product elements, causing their performance "sweet spot" to anomalously remain at the low order of $p \approx 2$, which severely limits the potential of high-order methods. To address this challenge, we design and implement a highly optimized PA operator within the MFEM framework, deeply integrated with a Geometric Multigrid (GMG) preconditioner. Our multi-level optimization strategy includes replacing the original $O(p^6)$ generic algorithm with an efficient $O(p^4)$ one based on tensor factorization, exploiting Voigt symmetry to reduce redundant computations for the elasticity problem, and employing macro-kernel fusion to enhance data locality and break the memory bandwidth bottleneck. Extensive experiments on mainstream x86 and ARM architectures demonstrate that our method successfully shifts the performance "sweet spot" to the higher-order region of $p \ge 6$. Compared to the MFEM baseline, the optimized core operator (kernel) achieves speedups of 7x to 83x, which translates to a 3.6x to 16.8x end-to-end performance improvement in the complete solution process. This paper provides a validated and efficient practical path for conducting large-scale, high-order elasticity simulations on mainstream CPU hardware.

Figures (6)

• Figure 1: Conceptual decomposition of the FEM operator into a chain of discrete operators, as implemented in libraries like MFEM.
• Figure 2: Architecture of the hybrid Geometric Multigrid (GMG) preconditioner. On the fine grid levels (right), a matrix-free Chebyshev polynomial smoother is used. On the coarsest level (left), the system is explicitly assembled and solved robustly with an Algebraic Multigrid (AMG) preconditioner. FA, PA, and PAop all operate within this framework, differing only in operator type and whether the integrator is optimized.
• Figure 3: Dataflow patterns for sum-factorization: (a) Naive fiber-wise pattern with large 3D intermediate tensors causing memory bottlenecks. (b) Optimized slice-wise pattern where computation is localized to cache-resident 2D slices, maximizing data locality.
• Figure 4: Comparison of total solve time (left) and peak memory usage (right) for Full Assembly (FA), baseline Partial Assembly (PA), and our optimized operator (PAop). The experiment is conducted on a problem with approximately 6.5 million Degrees of Freedom (DoFs) while varying the polynomial degree $p$. Our PAop method demonstrates superior performance, achieving speedups of up to 16.8x over PA and 2722x over FA. In terms of memory, the FA method's consumption grows exponentially, whereas PA and PAop maintain a nearly constant, low memory footprint.
• Figure 5: Computational throughput (in MDof/s) and corresponding speedup of the baseline PA and our optimized PAop kernels on Kunpeng (left) and AMD (right) platforms. The performance of the baseline PA kernel degrades sharply for $p>3$, a clear sign of cache thrashing as its working set exceeds the cache capacity. In contrast, our PAop kernel sustains high throughput across all polynomial degrees by improving memory locality and reducing the working set size. This results in maximum speedups of 54x on the Kunpeng platform and 83x on the AMD platform, demonstrating the effectiveness and portability of our optimization strategies.