2025/09/12 by Dominik Still, Still, Dominik, Niklas Fehn +5
Computer Science · Engineering · Physics and Astronomy · #65M60 #68W10 #76M10 #Advanced Numerical Methods in Computational Mathematics #FOS: Computer and information sciences #FOS: Mathematics #Mathematical Software (cs.MS) #Matrix Theory and Algorithms #Model Reduction and Neural Networks #Numerical Analysis (math.NA) #Performance (cs.PF)
paper · pdf · doi:10.48550/arxiv.2509.10226
openalex publication_date 2025/09/12 · openalex created_date 2025/10/10 · openalex updated_date 2026/07/28
This study presents novel strategies for improving the node-level performance of matrix-free evaluation of continuous and discontinuous Galerkin spatial discretizations on unstructured tetrahedral grids. In our approach the underlying integrals of a generic finite-element operator are computed cell-by-cell through numerical quadrature using tabulated dense local matrices of shape functions, achieving high throughput for low to moderate-order polynomial degrees. By employing dense matrix-matrix products instead of matrix-vector products for the cell-wise interpolation, the method reaches over 60% of peak performance. The optimization strategies exploit explicit data parallelism to enhance computational efficiency, complemented by a hierarchical mesh reordering algorithm that improves data locality. The matrix-free implementation achieves up to a 6× speedup compared to a global sparse matrix-based approach at a polynomial degree of three. The effectiveness of the method is demonstrated through numerical experiments on the Poisson and Navier--Stokes equations. The Poisson operator is preconditioned by a hybrid multigrid scheme that combines auxiliary continuous finite-element spaces, polynomial and geometric coarsening where possible while employing algebraic multigrid on the coarse mesh. Within the preconditioner, the implementation transitions between the matrix-free and matrix-based strategies for optimal efficiency. Finally, we analyze the strong scaling behavior of the Poisson and Helmholtz operators, demonstrating the method's potential to solve large real-world problems.