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

MFEM is a widely used finite-element library, but its native linear-elasticity Partial Assembly (PA) path still applies an $O((p+1)^6)$ contraction in the element operator, leaving the CPU operator-throughput sweet spot near $p\approx 2$ in our baseline measurements. This work closes this implementation gap for MFEM linear elasticity on affine tensor-product hexahedral meshes by integrating four well-established tensor-product PA optimizations (sum factorization, Voigt notation, macro-kernel fusion, and slice-wise loop reorganization) into MFEM's native linear-elasticity PA path. The resulting operator is evaluated in high-order GMG-PCG solves using MFEM's geometric multigrid (GMG) components. On AMD EPYC 7713, the optimized operator achieves $7\text{--}83\times$ kernel speedup and $3.6\text{--}16.8\times$ end-to-end speedup across $p\in\{1,2,4,8\}$. At fixed problem size, the kernel-time operator throughput peaks around $p=6$ and remains high at $p=8$, shifting the operator-throughput sweet spot to $p\ge 6$. The same trend is reproduced on Huawei~Kunpeng~920 (ARMv8.2). These results are accompanied by per-stage ablation and hardware-counter characterization; the implementation will be released on GitHub.