High-performance Implementation of Matrix-free High-order Discontinuous Galerkin Methods

Achieving a substantial part of peak performance on todays and future high-performance computing systems is a major challenge for simulation codes. In this paper we address this question in the context of the numerical solution of partial differential equations with finite element methods, in particular the discontinuous Galerkin method applied to a convection-diffusion-reaction model problem. Assuming tensor product structure of basis functions and quadrature on cuboid meshes in a matrix-free approach a substantial reduction in computational complexity can be achieved for operator application compared to a matrix-based implementation while at the same time enabling SIMD vectorization and the use of fused-multiply-add. Close to 60\% of peak performance are obtained for a full operator evaluation on a Xeon Haswell CPU with 16 cores and speedups of several hundred (with respect to matrix-based computation) are achieved for polynomial degree seven. Excellent weak scalability on a single node as well as the roofline model demonstrate that the algorithm is fully compute-bound with a high flop per byte ratio. Excellent scalability is also demonstrated on up to 6144 cores using message passing.