libCEED: Fast algebra for high-order element-based discretizations

implementations


Summary and statement of need
Finite element methods are widely used to solve partial differential equations (PDE) in science and engineering, but their standard implementation (Arndt et al., 2020;Kirk et al., 2006;Logg et al., 2012) relies on assembling sparse matrices. Sparse matrix multiplication and triangular operations perform a scalar multiply and add for each nonzero entry, just 2 floating point operations (flops) per scalar that must be loaded from memory (Williams et al., 2009). Modern hardware is capable of nearly 100 flops per scalar streamed from memory (Rupp, 2020) so sparse matrix operations cannot achieve more than about 2% utilization of arithmetic units. Matrix assembly becomes even more problematic when the polynomial degree p of the basis functions is increased, resulting in O(p d ) storage and O(p 2d ) compute per degree of freedom (DoF) in d dimensions. Methods pioneered by the spectral element community (Deville et al., 2002;Orszag, 1980) exploit problem structure to reduce costs to O(1) storage and O(p) compute per DoF, with very high utilization of modern CPUs and GPUs. Unfortunately, highquality implementations have been relegated to applications and intrusive frameworks that are often difficult to extend to new problems or incorporate into legacy applications, especially when strong preconditioners are required.
libCEED, the Code for Efficient Extensible Discretization , is a lightweight library that provides a purely algebraic interface for linear and nonlinear operators and preconditioners with element-based discretizations. libCEED provides portable performance via run-time selection of implementations optimized for CPUs and GPUs, including support for just-in-time (JIT) compilation. It is designed for convenient use in new and legacy software, and offers interfaces in C99 (International Standards Organisation, 1999), Fortran77 (ANSI, 1978, Python (Python, 2021), Julia (Bezanson et al., 2017), and Rust (Rust, 2021). Users and library developers can integrate libCEED at a low level into existing applications in place of existing matrix-vector products without significant refactoring of their own discretization infrastructure. Alternatively, users can utilize integrated libCEED support in MFEM (Anderson et al., 2020;MFEM, 2021).
In addition to supporting applications and discretization libraries, libCEED provides a platform for performance engineering and co-design, as well as an algebraic interface for solvers research like adaptive p-multigrid, much like how sparse matrix libraries enable development and deployment of algebraic multigrid solvers.

Concepts and interface
Consider finite element discretization of a problem based on a weak form with one weak derivative: find u such that where the functions f 0 and f 1 define the physics and possible stabilization of the problem (Brown, 2010) and the functions u and v live in a suitable space. Integrals in the weak form are evaluated by summing over elements e, where E e restricts to element e, B e evaluates solution values and derivatives to quadrature points, f acts independently at quadrature points, and W e is a (diagonal) weighting at quadrature points. By grouping the operations W e and f into a point-block diagonal D and stacking the restrictions E e and basis actions B e for each element, we can express the global residual in operator notation (Figure 1), where P is an optional external operator, such as the parallel restriction in MPI-based (Gropp et al., 2014) solvers. Inhomogeneous Neumann, Robin, and nonlinear boundary conditions can be added in a similar fashion by adding terms integrated over boundary faces while Dirichlet boundary conditions can be added by setting the target values prior to applying the operator representing the weak form. Similar face integral terms can also be used to represent discontinuous Galerkin formulations. libCEED's native C interface is object-oriented, providing data types for each logical object in the decomposition.

D CeedQFunction
User-defined action at quadrature points B CeedBasis Basis evaluation to quadrature (dense/structured) E CeedElemRestriction Restriction to each element (sparse/boolean) A CeedOperator Linear or nonlinear operator acting on L-vectors libCEED implementations ("backends") are free to reorder and fuse computational steps (including eliding memory to store intermediate representations) so long as the mathematical properties of the operator A are preserved. A CeedOperator is composed of one or more operators defined as in Figure 1, and acts on a CeedVector, which typically encapsulates zerocopy access to host or device memory provided by the caller. The element restriction E requires mesh topology and a numbering of DoFs, and may be a no-op when data is already composed by element (such as with discontinuous Galerkin methods). The discrete basis B is the purely algebraic expression of a finite element basis (shape functions) and quadrature; it often possesses structure that is exploited to speed up its action. Some constructors are provided for arbitrary polynomial degree H 1 Lagrange bases with a tensor-product representation due to the computational efficiency of computing solution values and derivatives at quadrature points via tensor contractions. However, the user can define a CeedBasis for arbitrary element topology including tetrahedra, prisms, and other realizations of abstract polytopes, by providing quadrature weights and the matrices used to compute solution values and derivatives at quadrature points from the DoFs on the element.
The physics (weak form) is expressed through CeedQFunction, which can either be defined by the user or selected from a gallery distributed with libCEED. These pointwise functions do not depend on element resolution, topology, or basis degree (see Figure 2), in contrast to systems like FEniCS where UFL forms specify basis degree at compile time. This isolation is valuable for hp-refinement and adaptivity (where h commonly denotes the average element size and p the polynomial degree of the basis functions; see Babuška & Suri (1994)) and pmultigrid solvers; mixed-degree, mixed-topology, and h-nonconforming finite element methods are readily expressed by composition. Additionally, a single source implementation (in vanilla C or C++) for the CeedQFunctions can be used on CPUs or GPUs (transparently using the NVRTC (2021), HIPRTC, or OCCA (OCCA Development Site, 2021) run-time compilation features).
libCEED provides computation of the true operator diagonal for preconditioning with Jacobi and Chebyshev as well as direct assembly of sparse matrices (e.g., for coarse operators in multigrid) and construction of p-multigrid prolongation and restriction operators. Preconditioning matrix-free operators is an active area of research; support for domain decomposition methods and inexact subdomain solvers based on the fast diagonalization method (Lottes & Fischer, 2005) are in active development. Figure 2: A schematic of element restriction and basis applicator operators for elements with different topology. This sketch shows the independence of Q-functions (in this case representing a Laplacian) on element resolution, topology, and basis degree.

High-level languages
libCEED provides high-level interfaces in Python, Julia, and Rust, each of which is maintained and tested as part of the main repository, but distributed through each language's respective package manager.
The Python interface uses CFFI, the C Foreign Function Interface (C Foreign Function Interface for Python, 2021). CFFI allows reuse of most C declarations and requires only a minimal adaptation of some of them. The C and Python APIs are mapped in a nearly 1:1 correspondence. For instance, a CeedVector object is exposed as libceed.Vector in Python, and supports no-copy host and GPU device interperability with Python arrays from the NumPy (Harris et al., 2020) or Numba (Lam et al., 2015) packages. The interested reader can find more details on libCEED's Python interface in Barra et al. (2020).
The Julia interface, referred to as LibCEED.jl, provides both a low-level interface, which is generated automatically from libCEED's C header files, and a high-level interface. The high-level interface takes advantage of Julia's metaprogramming and just-in-time compilation capabilities to enable concise definition of Q-functions that work on both CPUs and GPUs, along with their composition into operators as in Figure 1.
The Rust interface also wraps automatically-generated bindings from the libCEED C header files, offering increased safety due to Rust ownership and borrow checking, and more convenient definition of Q-functions (e.g., via closures). . libCEED provides a dynamic interface such that users only need to write a single source (no need for templates/generics) and can select the desired specialized implementation at run time. Moreover, each process or thread can instantiate an arbitrary number of backends on an arbitrary number of devices.

Performance benchmarks
The Exascale Computing Project (ECP) co-design Center for Efficient Exascale Discretization (CEED, 2021) has defined a suite of Benchmark Problems (BPs) to test and compare the performance of high-order finite element implementations (Fischer et al., 2020;Kolev et al., 2021). Figure 4 compares the performance of libCEED solving BP3 (CG iteration on a 3D Poisson problem) or CPU and GPU systems of similar (purchase/operating and energy) cost. These tests use PETSc  for unstructured mesh management and parallel solvers with GPU-aware communication ; a similar implementation with comparable performance is available through MFEM. : Performance for BP3 using the xsmm/blocked backend on a 2-socket AMD EPYC 7452 (32-core, 2.35GHz) and the cuda/gen backend on LLNL's Lassen system with NVIDIA V100 GPUs. Each curve represents fixing the basis degree p and varying the number of elements. The CPU enables faster solution of smaller problem sizes (as in strong scaling) while the GPU is more efficient for applications that can afford to wait for larger sizes. Note that the CPU exhibits a performance drop when the working set becomes too large for L3 cache (128 MB/socket) while no such drop exists for the GPU. (This experiment was run with release candidates of PETSc 3.14 and libCEED 0.7 using gcc-10 on EPYC and clang-10/CUDA-10 on Lassen.)

Demo applications and integration
To highlight the ease of library reuse for solver composition and leverage libCEED's full capability for real-world applications, libCEED comes with a suite of application examples, including problems of interest to the fluid dynamics and solid mechanics communities. The fluid dynamics example solves the 2D and 3D compressible Navier-Stokes equations using SU/SUPG stabilization and implicit, explicit, or IMEX time integration; Figure 5 shows vortices arising in the "density current" (Straka et al., 1993) when a cold bubble of air reaches the ground. The solid mechanics example solves static linear elasticity and hyperelasticity with load continuation and Newton-Krylov solvers with p-multigrid preconditioners; Figure 6 shows a twisted Neo-Hookean beam. Both of these examples have been developed using PETSc, where libCEED provides the matrix-free operator and preconditioner ingredient evaluation and PETSc provides the unstructured mesh management and parallel solvers.  libCEED also includes additional examples with PETSc, MFEM, and Nek5000 (Nek5000, 2021).
If MFEM is built with libCEED support, existing MFEM users can pass -d ceed-cuda: /gpu/cuda/gen to use a libCEED CUDA backend, and similarly for other backends. The libCEED implementations, accessed in this way, currently provide MFEM users with the fastest operator action on CPUs and GPUs (CUDA and HIP/ROCm) without writing any libCEED Q-functions.