Accelerating High-Order Finite Element Simulations at Extreme Scale with FP64 Tensor Cores

TL;DR

This paper demonstrates that the FP64 tensor cores on NVIDIA GPUs can be used to further accelerate large-scale finite element scientific computing applications, achieving significant speedups in key kernels of MFEM, a scalable open-source finite element library widely used in HPC applications.

Abstract

Finite element simulations play a critical role in a wide range of applications, from automotive design to tsunami modeling and computational electromagnetics. Performing these simulations efficiently at the high resolutions needed for practical applications and scientific insights necessitates the use of high-order methods and large-scale supercomputing. While much progress has been made in porting finite element codes to GPU systems in recent years, additional improvements in the efficiency and computational speed of GPU-accelerated high-order finite element simulations are in constant demand. In this paper, we demonstrate that the FP64 tensor cores on NVIDIA GPUs can be used to further accelerate such simulations, achieving significant speedups in key kernels of MFEM, a scalable open-source finite element library widely used in HPC applications. By integrating FP64 tensor cores with kernel fusion optimizations, we were able to achieve up to 2$\times$ performance gains and up to 83% energy efficiency gains on NVIDIA's Grace Hopper GH200 and Grace Blackwell GB200 architectures. To the best of our knowledge, this is the first time that FP64 tensor cores have been directly programmed to accelerate large-scale finite element scientific computing applications. We demonstrate the performance of the optimized kernels at exascale by showing near-perfect weak scaling efficiency and 90% strong scaling efficiency across nearly 10,000 GPUs on the Alps system. The new algorithms and MFEM enhancements directly benefit complex production codes, including the 2025 Gordon Bell Prize-winning application for real-time tsunami forecasting.

Figures (5)

• Figure 1: Finite element operators, $A$, where $P$ handles parallel scattering and MPI communication across global true degrees of freedom (T-vectors), $G$ manages mesh topology to local subdomains (L-vectors) and elements (E-vectors), $B$ encodes geometry mappings and tensor-product basis functions, and $D$ encapsulates physics at quadrature points (Q-vectors)---with numerical kernels reducing to batched small dense tensor contractions implemented as matrix-matrix multiplications (DGEMMs) for performance portability via vendor-optimized BLAS.
• Figure 2: Shared memory banks for matrix A/B/C of an $m = 25$, $n=5$, $k=4$ matrix-matrix multiplication with the maps shown in (\ref{['eq:fm-fn-fk']}). In each of the phases, the 16 lanes access 16 distinct banks thus there are no shared memory bank conflicts.
• Figure 3: Throughput, in billion degrees of freedom (GDOF) per second, for the finite element kernels corresponding to the off-diagonal blocks in (\ref{['eq:block-mult']}) on a single NVIDIA GH200 Grace Hopper Superchip (left) and a single NVIDIA GB200 Grace Blackwell Superchip (right). Fusing loops into a single kernel and using tensor-core-enabled DMMA instructions ("DMMA Fused PA") together yields a 2$\times$ performance gain compared to the original ("PA") kernel.
• Figure 4: Strong scalability of the MFEM finite element solver on Alps, from 36 nodes (144 NVIDIA GH200 Superchips) to 2,304 nodes (9,216 GH200). Black numbers above each bar indicate the corresponding parallel efficiencies. The three kernel implementations---PA, Fused MF and Fused PA---(left plot) and their respective FP64 tensor-core-accelerated versions (DMMA; right plot) achieve excellent strong scaling (86--91%) over a 64$\times$ increase of nodes.
• Figure 5: Weak scalability of the MFEM finite element solver on Alps, from 36 nodes (144 GH200) to 2,304 nodes (9,216 GH200). Numbers along the graph lines indicate parallel efficiency. Each of the three kernel implementations---PA (top), Fused MF (middle) and Fused PA (bottom)---and their respective FP64 tensor-core-accelerated versions (DMMA) achieve ideal (linear) weak scaling over a 64$\times$ increase of nodes.