Speaker
Description
The objective of my talk is the efficient implementation of solvers for linear systems arising from the discretization of partial differential equations (PDEs) with higher-order finite element methods. The classical workflow is to first construct a sparse matrix and a right-hand side through the evaluation of element integrals in one part of the code. The matrix and vector are then passed to a specialized solver software, which usually employs an iterative method of the Krylov family combined with some sophisticated preconditioner. While this abstraction has been highly successful for many problems, research on performance optimizations of this parts reveals that all steps are characterized by a low arithmetic intensity, i.e., a low ratio of arithmetic operations compared to memory access of loading all matrix entries repeatedly into the compute units. Due to the continuous hardware evolution of the last decades, computations have actually become very cheap in comparison to memory access. In fact, modern CPUs and GPUs can sustain more than 100 arithmetic operations for each data element loaded from main (RAM) memory. It is therefore natural to seek for algorithms that avoid storing the big sparse matrices altogether, and rather apply the operator action on the fly within iterative solvers. This leads to matrix-free solvers.
The talk will present the main design characteristics of fast matrix-free operator evaluation for PDEs based on computing the cell and face integrals on the fly with the deal.II finite element library. I will consider both memory access and efficiency of arithmetic work, such as sum-factorization techniques, with respect to modern hardware of the exascale era, with the goal to minimize the time-to-solution. In the talk, applications in computational fluid dynamics will be shown, comparing both classical H1 and L2 conforming methods against H(div) conforming Raviart-Thomas operators. While these methods show large speedups for matrix-vector products in isolation, the overall solver strategy crucially depends on efficient preconditioners that balance low iteration counts with fast operator evaluation, for which either matrix-free and matrix-based ingredients can be attractive. In our research, we consider both multigrid techniques with polynomial or geometric mesh coarsening, block-Jacobi methods with local solvers based on the fast diagonalization method and approximate incomplete matrix factorizations. Finally, the talk will look into the node-level performance optimizations of these algorithms and scalability to large supercomputers.
