This paper introduces a code generator designed for node-level optimized, extreme-scalable, matrix-free finite element operators on hybrid tetrahedral grids. It optimizes the local evaluation of bilinear forms through various techniques including tabulation, relocation of loop invariants, and inter-element vectorization - implemented as transformations of an abstract syntax tree. A key contribution is the development, analysis, and generation of efficient loop patterns that leverage the local structure of the underlying tetrahedral grid. These significantly enhance cache locality and arithmetic intensity, mitigating bandwidth-pressure associated with compute-sparse, low-order operators. The paper demonstrates the generator's capabilities through a comprehensive educational cycle of performance analysis, bottleneck identification, and emission of dedicated optimizations. For three differential operators ($-\Delta$, $-\nabla \cdot (k(\mathbf{x})\, \nabla\,)$, $\alpha(\mathbf{x})\, \mathbf{curl}\ \mathbf{curl} + \beta(\mathbf{x}) $), we determine the set of most effective optimizations. Applied by the generator, they result in speed-ups of up to 58$\times$ compared to reference implementations. Detailed node-level performance analysis yields matrix-free operators with a throughput of 1.3 to 2.1 GDoF/s, achieving up to 62% peak performance on a 36-core Intel Ice Lake socket. Finally, the solution of the curl-curl problem with more than a trillion ($ 10^{12}$) degrees of freedom on 21504 processes in less than 50 seconds demonstrates the generated operators' performance and extreme-scalability as part of a full multigrid solver.

Block Floating Point (BFP) arithmetic is currently seeing a resurgence in interest because it requires less power, less chip area, and is less complicated to implement in hardware than standard floating point arithmetic. This paper explores the application of BFP to mixed- and progressive-precision multigrid methods, enabling the solution of linear elliptic partial differential equations (PDEs) in energy- and hardware-efficient integer arithmetic. While most existing applications of BFP arithmetic tend to use small block sizes, the block size here is chosen to be maximal such that matrices and vectors share a single exponent for all entries. This is sometimes also referred to as a scaled fixed-point format. We provide algorithms for BLAS-like routines for BFP arithmetic that ensure exact vector-vector and matrix-vector computations up to a specified precision. Using these algorithms, we study the asymptotic precision requirements to achieve discretization-error-accuracy. We demonstrate that some computations can be performed using as little as 4-bit integers, while the number of bits required to attain a certain target accuracy is similar to that of standard floating point arithmetic. Finally, we present a heuristic for full multigrid in BFP arithmetic based on saturation and truncation that still achieves discretization-error-accuracy without the need for expensive normalization steps of intermediate results.

Accurate retrodictions of the past evolution of convection in the Earth's mantle are crucial to obtain a qualitative understanding of this central mechanism behind impactful geological events on our planet. They require highly resolved simulations and therefore extremely scalable numerical methods. This paper applies the massively parallel matrix-free finite element framework HyTeG to approximate solutions to stationary Stokes systems and time-dependent, coupled convection problems. It summarizes the underlying mathematical model and verifies the implementation through semi-analytical setups and community benchmarks. The numerical results agree with the expected outcomes from the literature.

Simulating the Earth’s mantle convection at full convective vigor on planetary scales is a fundamental challenge in Geodynamics even for state of the art high- performance computing (HPC) systems. Realistic Earth mantle convection simulations can contribute a decisive link between uncertain input parameters, such as rheology, and testable preconditions, such as dynamic topography. The vertical deflections predicted by such models may then be tested against the geological record. Considering realistic Earth-like Rayleigh numbers (∼ 108) a resolution of the thermal boundary layer on the order of ∼ 10 km is necessary considering the volume of the Earth’s mantle. Simulating Earth’s mantle convection at this level of accuracy requires solving sparse indefinite systems with more than a trillion degrees of freedom, computational feasible on exascale HPC systems. This can only be achieved by mantle convection codes providing high degrees of parallelism and scalability. Earlier approaches from applying a prototype framework using hierarchical hybrid grids (HHG) as solvers for such systems demonstrated the scalability of the underlying concept for future generations of exascale computing systems. In consideration of the TerraNeo project, here we report on the progress of utilizing the improved framework HyTeG (Hybrid Tetrahedral Grids) based on matrix-free multigrid solvers in combination with highly efficient parallelization and scalability. This will allow to solve systems with more than a trillion degrees of freedom on present and future generations of exascale computing systems. We also report on the advances in developing the scalable mantle convection code TerraNeo using the HyTeG framework to realise extreme-scale mantle convection simulations with a resolution on the order of ∼ 1 km.

This paper presents efficient data structures for the implementation of matrix-free finite element methods on block-structured, hybrid tetrahedral grids. It provides a complete categorization of all geometric sub-objects that emerge from the regular refinement of the unstructured, tetrahedral coarse grid and describes efficient iteration patterns and analytical linearization functions for the mapping of coefficients to memory addresses. This foundation enables the implementation of fast, extreme-scalable, matrix-free, iterative solvers, and in particular geometric multigrid methods by design. Their application to the variable-coefficient Stokes system subject to an enriched Galerkin discretization and to the curl-curl problem discretized with N\'ed\'elec edge elements showcases the flexibility of the implementation. Eventually, the solution of a curl-curl problem with $1.6 \cdot 10^{11}$ (more than one hundred billion) unknowns on more than $32000$ processes with a matrix-free full multigrid solver demonstrates its extreme-scalability.

We employ textbook multigrid efficiency (TME), as introduced by Achi Brandt, to construct an asymptotically optimal monolithic multigrid solver for the Stokes system. The geometric multigrid solver builds upon the concept of hierarchical hybrid grids (HHG), which is extended to higher-order finite-element discretizations, and a corresponding matrix-free implementation. The computational cost of the full multigrid (FMG) iteration is quantified, and the solver is applied to multiple benchmark problems. Through a parameter study, we suggest configurations that achieve TME for both, stabilized equal-order, and Taylor-Hood discretizations. The excellent node-level performance of the relevant compute kernels is presented via a roofline analysis. Finally, we demonstrate the weak and strong scalability to up to $147,456$ parallel processes and solve Stokes systems with more than $3.6 \times 10^{12}$ (trillion) unknowns.

Motivated by challenges in the Earth's mantle convection, we present a massively parallel implementation of an Eulerian-Lagrangian method for the advection-diffusion equation in the advection-dominated regime. The advection term is treated by a particle-based characteristics method coupled to a block-structured finite element framework. Its numerical and computational performance is evaluated in multiple two- and three-dimensional benchmarks, including curved geometries, discontinuous solutions, and pure advection, and it is applied to a coupled nonlinear system modeling buoyancy-driven convection in Stokes flow. We demonstrate the parallel performance in a strong and weak scaling experiment, with scalability to up to 147,456 parallel processes, solving for more than 5.2 x 10(10) (52 billion) degrees of freedom per time-step.

This work is based on the seminar titled ‘Resiliency in Numerical Algorithm Design for Extreme Scale Simulations’ held March 1–6, 2020, at Schloss Dagstuhl, that was attended by all the authors. Advanced supercomputing is characterized by very high computation speeds at the cost of involving an enormous amount of resources and costs. A typical large-scale computation running for 48 h on a system consuming 20 MW, as predicted for exascale systems, would consume a million kWh, corresponding to about 100k Euro in energy cost for executing 10 23 floating-point operations. It is clearly unacceptable to lose the whole computation if any of the several million parallel processes fails during the execution. Moreover, if a single operation suffers from a bit-flip error, should the whole computation be declared invalid? What about the notion of reproducibility itself: should this core paradigm of science be revised and refined for results that are obtained by large-scale simulation? Naive versions of conventional resilience techniques will not scale to the exascale regime: with a main memory footprint of tens of Petabytes, synchronously writing checkpoint data all the way to background storage at frequent intervals will create intolerable overheads in runtime and energy consumption. Forecasts show that the mean time between failures could be lower than the time to recover from such a checkpoint, so that large calculations at scale might not make any progress if robust alternatives are not investigated. More advanced resilience techniques must be devised. The key may lie in exploiting both advanced system features as well as specific application knowledge. Research will face two essential questions: (1) what are the reliability requirements for a particular computation and (2) how do we best design the algorithms and software to meet these requirements? While the analysis of use cases can help understand the particular reliability requirements, the construction of remedies is currently wide open. One avenue would be to refine and improve on system- or application-level checkpointing and rollback strategies in the case an error is detected. Developers might use fault notification interfaces and flexible runtime systems to respond to node failures in an application-dependent fashion. Novel numerical algorithms or more stochastic computational approaches may be required to meet accuracy requirements in the face of undetectable soft errors. These ideas constituted an essential topic of the seminar. The goal of this Dagstuhl Seminar was to bring together a diverse group of scientists with expertise in exascale computing to discuss novel ways to make applications resilient against detected and undetected faults. In particular, participants explored the role that algorithms and applications play in the holistic approach needed to tackle this challenge. This article gathers a broad range of perspectives on the role of algorithms, applications and systems in achieving resilience for extreme scale simulations. The ultimate goal is to spark novel ideas and encourage the development of concrete solutions for achieving such resilience holistically.

Programming current supercomputers efficiently is a challenging task. Multiple levels of parallelism on the core, on the compute node, and between nodes need to be exploited to make full use of the system. Heterogeneous hardware architectures with accelerators further complicate the development process. waLBerla addresses these challenges by providing the user with highly efficient building blocks for developing simulations on block-structured grids. The block-structured domain partitioning is flexible enough to handle complex geometries, while the structured grid within each block allows for highly efficient implementations of stencil-based algorithms. We present several example applications realized with waLBerla, ranging from lattice Boltzmann methods to rigid particle simulations. Most importantly, these methods can be coupled together, enabling multiphysics simulations. The framework uses meta-programming techniques to generate highly efficient code for CPUs and GPUs from a symbolic method formulation.

Software development for applications in computational science and engineering has become complex in recent years. This is mainly due to the increasing parallelism and heterogeneity in modern computer architectures and to the more realistic physical and mathematical models that have to be processed. One idea to address this issue is to use code generation techniques. In contrast to manual implementations in a general-purpose computing language, they allow to integrate automatic code transforms to produce efficient code for different models and platforms. As an example the numerical solution of an elliptic partial differential equation via generated geometric multigrid solvers is considered. We present three code generation approaches for it and discuss their advantages and disadvantages with respect to performance, portability, and productivity.

Simulation of mantle convection on planetary scales is considered a grand-challenge application even in the exascale era. The reason being the enormous spatial and temporal scales that must be resolved in the computation as well as the complexities of realistic models and the large parameter uncertainties that need to be handled by advanced numerical methods. This contribution reports on the TerraNeo project which delivered novel matrix-free geometric multigrid solvers for the Stokes system that forms the core of mantle convection models. In TerraNeo the hierarchical hybrid grids paradigm was employed to demonstrate that scalability can be achieved when solving the Stokes system with more than ten trillion (1.1 ⋅ 1013) degrees of freedom even on present-day peta-scale supercomputers. Novel concepts were developed to ensure resilience of algorithms even in case of hard faults and new scheduling algorithms proposed for ensemble runs arising in Multilevel Monte Carlo algorithms for uncertainty quantification. The prototype framework was used to investigate geodynamic questions such as high velocity asthenospheric channels and dynamic topography and to perform adjoint inversions. We also describe the redesign of our software to support more advanced discretizations, adaptivity, and highly asynchronous execution while ensuring sustainability and flexibility for future extensions.

Realistic simulations in engineering or in the materials sciences can consume enormous computing resources and thus require the use of massively parallel supercomputers. The probability of a failure increases both with the runtime and with the number of system components. For future exascale systems, it is therefore considered critical that strategies are developed to make software resilient against failures. In this article, we present a scalable, distributed, diskless, and resilient checkpointing scheme that can create and recover snapshots of a partitioned simulation domain. We demonstrate the efficiency and scalability of the checkpoint strategy for simulations with up to 40 billion computational cells executing on more than 400 billion floating point values. A checkpoint creation is shown to require only a few seconds and the new checkpointing scheme scales almost perfectly up to more than 260, 000 (218) processes. To recover from a diskless checkpoint during runtime, we realize the recovery algorithms using ULFM MPI. The checkpointing mechanism is fully integrated in a state-of-the-art high-performance multi-physics simulation framework. We demonstrate the efficiency and robustness of the method with a realistic phase-field simulation originating in the material sciences and with a lattice Boltzmann method implementation.

In this article, a new generic higher-order finite-element framework for massively parallel simulations is presented. The modular software architecture is carefully designed to exploit the resources of modern and future supercomputers. Combining an unstructured topology with structured grid refinement facilitates high geometric adaptability and matrix-free multigrid implementations with excellent performance. Different abstraction levels and fully distributed data structures additionally ensure high flexibility, extensibility, and scalability. The software concepts support sophisticated load balancing and flexibly combining finite element spaces. Example scenarios with coupled systems of PDEs show the applicability of the concepts to performing geophysical simulations.

ABSTRACT In this article, a new generic higher-order finite-element framework for massively parallel simulations is presented. The modular software architecture is carefully designed to exploit the resources of modern and future supercomputers. Combining an unstructured topology with structured grid refinement facilitates high geometric adaptability and matrix-free multigrid implementations with excellent performance. Different abstraction levels and fully distributed data structures additionally ensure high flexibility, extensibility, and scalability. The software concepts support sophisticated load balancing and flexibly combining finite-element spaces. Example scenarios with coupled systems of partial differential equations show the applicability of the concepts to performing geophysical simulations. GRAPHICAL ABSTRACT