Scalability of high-performance PDE solvers

Abstract

Performance tests and analyses are critical to effective high-performance computing software development and are central components in the design and implementation of computational algorithms for achieving faster simulations on existing and future computing architectures for large-scale application problems. In this article, we explore performance and space-time trade-offs for important compute-intensive kernels of large-scale numerical solvers for partial differential equations (PDEs) that govern a wide range of physical applications. We consider a sequence of PDE-motivated bake-off problems designed to establish best practices for efficient high-order simulations across a variety of codes and platforms. We measure peak performance (degrees of freedom per second) on a fixed number of nodes and identify effective code optimization strategies for each architecture. In addition to peak performance, we identify the minimum time to solution at 80% parallel efficiency. The performance analysis is based on spectral and p-type finite elements but is equally applicable to a broad spectrum of numerical PDE discretizations, including finite difference, finite volume, and h-type finite elements.

ExaDG: High-Order Discontinuous Galerkin for the Exa-Scale

“…The benchmark involves a continuous finite element discretization of the Laplacian (3), using matrix-free operator evaluation within a conjugate gradient solver preconditioned by the matrix diagonal. In this study, we consider the case BP5 [29], see also https://ceed.exascaleproject.org/bps/, which integrates the weak form (4) of polynomial degree k using a Gauss-Lobatto quadrature formula with n 1D q = k + 1 quadrature points on a cube with deformed elements. While this integration is not exact, it is the typical spectral element setup with an identity interpolation matrix S i = I in Eq.…”

“…Since the metric terms do not depend on the shape function indices i K and j , and the sum over j does not depend on i K , the summations in the equation can be broken up into (1) an dn q × n dof,ele matrix operation to evaluate the reference element derivative of u (K) at the quadrature points, (2) the application of metric terms as well as other physics terms at n q quadrature points, and (3) an n dof,ele × dn q matrix operation to test by all n dof,ele test functions and perform the summation over the quadrature points. The separation of point-wise physics evaluation at quadrature points is a common abstraction in integration-based matrix-free methods [29,41,49,50].…”

“…• support for both continuous [49] and discontinuous finite elements on uniform and adaptively refined meshes with hanging nodes and deformed elements, • support for arbitrary polynomial expansions on quadrilateral and hexahedral element shapes as well as tensor product quadrature rules, • minimization of arithmetic operations by using available symmetries, such as the even-odd decomposition [69] and a switch between the collocation derivative (5) for n 1D q ≈ k + 1 quadrature points or an alternative variant based on derivatives of the original polynomials as used in [49] and discussed in [29], • flexible implementation of operations at quadrature points, • vectorization across several elements to optimally use SIMD units (AVX, AVX-512, AltiVec) of modern processors, • applicability to modern multi-core CPUs as well as GPUs [51,57], • data access optimizations such as element-based loops for DG elements [50,56], • and MPI implementation with tight data exchange as well as MPI-only and shared-memory models [43,48,54].…”

Exascale applications: skin in the game

“…To achieve this efficiency, high-order methods use mesh elements that are mapped from canonical reference elements (hexes, wedges, pyramids and tetrahedra) and exploit, where possible, the tensor-product structure of the canonical mesh elements and finite-element spaces. Through matrix-free partial assembly, the use of canonical reference elements enables substantial cache efficiency and minimizes extraneous data movement in comparison to traditional low-order approaches [5].…”

A matrix‐free approach for finite‐strain hyperelastic problems using geometric multigrid

“…The finite element method (FEM) with linear or quadratic elements (that reduce the potential for locking), or mixed or enhanced formulations (that avoid locking) are often adopted, whereas higher order elements are rarely employed. This can potentially reduce the advantage of matrix‐free methods using sum‐factorization techniques, which are generally more competitive for higher order elements . However, we note that, even with linear FEM, large‐scale computations with as many as 10 12 unknowns are possible with a matrix‐free approach but infeasible with sparse matrices as there is simply not enough memory to store the sparse tangent matrix even on large supercomputers .…”