Source-linked AI summary
Fast matrix-free evaluation of discontinuous Galerkin finite element operators
Martin Kronbichler, Katharina Kormann
TL;DR
Matrix-free DG operator evaluation needs efficient treatment of tensor-product quadrature, data access, geometry, and communication on general quadrilateral and hexahedral meshes. The paper develops and compares sum-factorized kernels and vectorized layouts for cells and faces, reaching up to 60% of arithmetic peak locally, while full evaluation is constrained by memory bandwidth and communication.
Problem
A systematic performance analysis of generic, matrix-free DG operator evaluation for cells and faces, including data access and parallelization choices, is needed.
Method
The paper combines sum factorization, vectorization over several cells, optimized data layouts, and MPI strategies for matrix-free DG operator evaluation.
Results
Local kernels reach up to 60% of arithmetic peak on Intel Haswell and Broadwell and up to 50% on Intel Knights Landing, while MPI exchange can consume up to a third of single-node operator-evaluation time.
Takeaways & Limitations
Matrix-free DG evaluation can achieve high local-kernel throughput, but full-operator performance depends strongly on memory bandwidth, geometry handling, vector movement, and MPI exchange.
Takeaways & Limitations
The face-integral data-access strategy incurs additional data transfer, and the considered vectorization scheme increases indirect addressing in some alternatives.
Abstract
from arXiv · showhide
We present an algorithmic framework for matrix-free evaluation of discontinuous Galerkin finite element operators based on sum factorization on quadrilateral and hexahedral meshes. We identify a set of kernels for fast quadrature on cells and faces targeting a wide class of weak forms originating from linear and nonlinear partial differential equations. Different algorithms and data structures for the implementation of operator evaluation are compared in an in-depth performance analysis. The sum factorization kernels are optimized by vectorization over several cells and faces and an even-odd decomposition of the one-dimensional compute kernels. In isolation our implementation then reaches up to 60\% of arithmetic peak on Intel Haswell and Broadwell processors and up to 50\% of arithmetic peak on Intel Knights Landing. The full operator evaluation reaches only about half that throughput due to memory bandwidth limitations from loading the input and output vectors, MPI ghost exchange, as well as handling variable coefficients and the geometry. Our performance analysis shows that the results are often within 10\% of the available memory bandwidth for the proposed implementation, with the exception of the Cartesian mesh case where the cost of gather operations and MPI communication are more substantial.
1 Introduction
This work develops a systematic, performance-focused framework for matrix-free discontinuous Galerkin operator evaluation on quadrilateral and hexahedral meshes. It combines optimized tensor-product kernels, vectorized data layouts, and MPI parallelization for general meshes and variable coefficients.
- Scope: The framework targets quadrilateral and hexahedral meshes with moderate polynomial degrees 2 ≤p ≤10, including general geometries and variable coefficients.The authors also report results beyond this range to check for sudden performance breakdowns.
- Motivation: The paper fills a gap by systematically analyzing DG matrix-free operator evaluation and its performance envelopes independently of user interfaces.It targets generic software interfaces rather than specializations tied to particular equations or parallelization schemes.
- Contributions: Efficient local kernels use shape-value interpolation and derivative computation with vectorization over several elements.The contribution focuses on optimized tensor-product operations for cell and face integrals.
- Contributions: The work analyzes degree-of-freedom storage and efficient MPI parallelization as integral parts of high-performance operator evaluation.These choices address both data movement and communication in the matrix-free implementation.
- Scalability: The resulting algorithms are designed for massively parallel use because operator-evaluation communication is restricted to nearest neighbors.Selected implementations are available in the deal.II finite element library and interoperate with continuous finite element code paths.
2 DG algorithm
The paper evaluates DG operators matrix-free by looping over cells and faces, using numerical quadrature, sum factorization, and tensor-product kernels for efficient interpolation and integration.
- Face operators: Face integration uses traces from neighboring cells, numerical fluxes, face quadrature weights, and opposite normals to form conservative local contributions.Boundary faces instead use exterior states determined by boundary conditions, including mirror-principle data for Dirichlet conditions.
- DG formulations: The framework covers stationary advection and Laplacian DG discretizations, whose linear operators map coefficient values u to tested integrals y.The approach is based on numerical integration rather than final cell matrices, so it extends to nonlinear equations.
- Quadrature formulation: Numerical integration transforms tensor-product basis functions from reference coordinates to physical elements and evaluates cell terms at tensor-product quadrature points.Geometry enters through the element mapping and Jacobian during quadrature-point evaluation.
- Operator evaluation: Matrix-free DG evaluation loops over cells, interior faces, and boundary faces, gathering local values, computing quadrature-point operations, integrating with test functions, and assembling global results.The workflow also includes MPI import and residual compression for distributed execution.
- Sum factorization: Sum factorization applies one-dimensional interpolation and derivative matrices along coordinate directions, replacing naive multidimensional evaluation with small matrix-matrix multiplications.For general bases, the transformation reduces gradient tensor-product calls from d^2 to 2d and lowers the arithmetic cost compared with naive evaluation.
3 Performance-driven design of sum factorization kernels
The kernels combine sum factorization, vectorization across cells, cache-aware loop organization, register optimization, and even-odd decomposition to accelerate DG operator evaluation. Across Intel systems, isolated kernels achieve substantial fractions of arithmetic peak, while cache capacity, memory behavior, compiler choices, and face-integral costs constrain performance.
- Vectorization strategy: Vectorization across several cells outperforms within-cell and automatic vectorization for low and moderate polynomial degrees.Across-cell vectorization avoids cross-lane exchanges and separate gathers, whereas automatic vectorization is slower by factors of 2.1 to 3.7.
- Vectorization strategy: Cache behavior limits preferred across-cell vectorization mainly when temporary data exceed cache capacity, with within-cell vectorization becoming superior beyond k = 21 in 3D.Across-cell vectorization is reported as uncritical for the common range p < 15, while the 3D local vector data exceed L2 cache beyond k = 21.
- Optimization of loop kernels: Register-oriented loop optimization increases throughput by approximately a factor of three for the small matrix-matrix products.The implementation keeps one-dimensional input data in registers and loads interpolation-matrix entries from L1 cache at full vector width.
- Optimization of loop kernels: Even-odd decomposition produces higher and more regular throughput than full-matrix loop-blocking approaches for one-dimensional kernels.The decomposition reduces arithmetic and naturally supports two-way loop unrolling; further blocking may help at p > 20.
- Compute performance on CPUs and Xeon Phi: Up to 60% of arithmetic peak is reached on Haswell and Broadwell, while KNL approaches 50% for low-to-moderate polynomial degrees.Haswell reaches 340–350 GFLOP/s around p = 10; Broadwell reaches 690 GFLOP/s in 2D and 670 GFLOP/s in 3D, while KNL reaches around 1 TFLOP/s for low to moderate degrees.
- Compute performance on CPUs and Xeon Phi: Combined cell and face integrals reach 620 GFLOP/s on Broadwell and 820 GFLOP/s on KNL in the tested periodic same-cell configuration.Face integrals contribute more than two thirds of the combined work in this benchmark.
4 Data access patterns and MPI parallelization
The full DG operator is shaped by vector access, face/cell data movement, and MPI communication, making memory bandwidth and communication central performance limits. Interleaving cell and face work, reducing communication, and avoiding separate face storage improve throughput, but three-dimensional evaluations retain larger gaps from compute-only performance.
- 4.1 Vector access analysis for cell integrals only: Memory bandwidth limits 2D and low-to-medium-degree cell-integral throughput, while Haswell and KNL are mostly compute limited in 3D.The measured envelope is closely described by arithmetic throughput and memory-bandwidth limits, with intermediate gaps attributed partly to cache effects.
- Data access strategy: Interleaving face and cell integrals reuses cached vector data and improves performance when cell integrals are already memory bound.This design is also beneficial for curved geometries and high-bandwidth systems such as KNL because they have lower FLOP/byte ratios.
- 4.5 Index access and ghosting: Slim MPI communication improves throughput by more than 30% for the p = 6 Laplacian and almost 70% for p = 6 advection at 8.9 million degrees of freedom.The improvement is less pronounced on smaller meshes with larger processor-boundary surfaces, while on much larger meshes communication falls below 10% of operator time.
- 4.5 Index access and ghosting: A single-loop cell-and-face implementation is almost twice as fast as face-separate storage at p = 5 and remains 25% faster at p = 11.Face-separate storage becomes superior only above p > 15, despite reaching more than 90 GB/s memory throughput through p = 11.
- 4.6 Performance on Haswell, Broadwell and Knights Landing: In 3D, MPI communication and pack/unpack consume up to one third of compute time, while the Laplacian reaches 1.8–2.1 billion degrees of freedom per second after data exchange is removed.The reported range applies to 4 ≤k ≤9 on 28 Broadwell and 64 KNL cores; remaining gaps include gather/scatter costs and memory stalls.
5 Representation of geometry
The paper compares geometry representations for matrix-free DG operator evaluation, balancing storage, computation, arithmetic intensity, and memory traffic across cells and faces. Precomputed face data generally improves complex-geometry throughput, while collocation derivatives can outperform loaded Jacobians for cell integrals at suitable degrees.
- The four geometry strategies range from storing mapping support points or physical quadrature points to precomputing Jacobians or equation-specific effective coefficients.These choices trade storage and memory access against on-the-fly computation for general geometries and variable coefficients.
- Precomputed face geometry delivers better throughput than computing full cell-and-face geometry, which remains below 350 million degrees of freedom per second.The comparison concerns a 3D Laplacian on a hypershell mesh with 8–57 million degrees of freedom.
- For p > 5, precomputation is preferable for cell geometry because gather-based computation from a continuous finite element field is competitive mainly at low degrees.
- Cell Jacobians computed from precomputed quadrature-point positions can outperform precomputed final coefficients at low degrees, but derivative temporaries spill toward main memory for p > 10.
- Precomputed geometry variants operate near the memory-bandwidth limit, whereas collocation derivatives achieve higher arithmetic intensity with both memory-bound and compute-bound components.At higher degrees, the collocation derivative approach is especially mixed in its performance behavior.
6 Conclusions and future developments
The paper develops optimized sum-factorization kernels and data structures for matrix-free DG cell and face integrals on quadrilateral and hexahedral meshes. The kernels achieve high isolated arithmetic throughput, while full evaluation is constrained by memory traffic and communication, motivating interleaved cell-face processing.
- The implementations support explicit time integration and iterative solution of linear and nonlinear systems whose dominant work is matrix-vector products or residual evaluations.
- The local kernels reach up to 60% of arithmetic peak on Haswell and Broadwell and up to 50% on Knights Landing.The optimizations vectorize over several cells and use even-odd decomposition to reduce one-dimensional interpolation arithmetic.
- Full operator evaluation is mostly memory-bandwidth bound for low and moderate degrees up to ten, so cell and face integrals should be interleaved.
- MPI data exchange can consume up to one-third of operator evaluation time on a single node, reducing throughput from around 2 billion degrees of freedom per second.