ENVIRONMENT AND ENERGY

Introduction

SeisSol [1] is an HPC software for the simulation of earthquake source mechanisms and resulting ground motion. It employs the Discontinuous Galerkin method with ADER time-stepping (ADER-DG) to solve the elastic wave equation in various types of elastic media, including isotropic and anisotropic and viscoelastic materials. It also features coupling to acoustic materials for the simulation of earthquake-tsunami scenarios. For earthquake sourcing, SeisSol simulates the dynamic coseismic breaking of pre-existing faults due to frictional failure as a physics-based model (“dynamic rupture simulation”) to excite seismic waves.

For the GCS project pr83no, we enhanced SeisSol's capabilities to simulate small induced earthquakes, as they may occur during the preparation and operation of Geothermal energy production plants. For the simulation of such scenarios, the interaction of the solid phase (rock) and the injected fluid (water phase) has to be accounted for. We therefore implemented and optimized the simulation of poroelastic materials [2] in SeisSol, we refactored and extended SeisSol’s routines for dynamic rupture simulations, and we performed first benchmark simulations combining these new features.

Results and Methods

Poroelasticity

Poroelastic media consist of an elastic porous skeleton, which is filled with a viscous fluid. The deformation of the skeleton drives fluid flow and vice versa. The interaction of both phases renders the governing equation stiff and requires a more complicated solution procedure for the ADER-DG scheme – using a locally implicit space-time-predictor. Here, a lot of medium-sized linear systems (with ~1,000 unknowns) have to be solved. Since all these linear systems share the same sparsity pattern, a specialized solution algorithm has been developed. This algorithm takes the block-wise sparsity pattern into account and reduces the required floating point operations by up to a factor of 25 compared to straightforward LU decomposition. The code generator YATeTo [3] is able to map the solution algorithm to optimized implementations of matrix-matrix multiplication (GEMM) and employs hardware-aware backends for optimal performance. Figure 1 shows a strong scaling study for a wave-propagation scenario in poroelastic media. We compare global (GTS) and local time stepping (GTS) for different convergence orders. Scenarios computed with higher convergence order achieve a better performance due to the higher computational intensity. Also, global time stepping achieves a higher performance compared to local time stepping due to the simpler update scheme. Despite the reduced performance, local time stepping reduces the time to solution significantly. The parallel efficiency plot highlights the good scaling behaviour of SeisSol. With a node-level performance of up to 1.4 TFLOP/s per node, SeisSol reaches approximately 40% of the available peak performance of SuperMUC-NG, when using polynomials up to degree 6 as basis functions.

Dynamic Rupture Simulation

The dynamic rupture source model combines frictional failure along fault surfaces with the radiation of seismic waves. Both aspects are linked and drive each other. This source mechanism allows researchers to gain insight into the hidden nucleation processes of earthquakes.

SeisSol’s early roots date back almost 20 years. Over these years, it has almost entirely been re-implemented from Fortran to C++  – with the dynamic rupture implementation being one of the last “legacy” parts. In particular to enable GPU support for upcoming machines, like SuperMUC-NG Phase 2, the respective routines have now been redesigned, as well.

SeisSol uses OpenMP (on CPU machines) or SYCL/CUDA (on GPU architectures) for intra-node parallelism. SYCL has been chosen as it allows a performance-portable implementation without code duplication for AMD, Intel and NVidia GPU architectures. On CPUs, padding of the vector of degrees of freedom enables the efficient use of AVX-512 vector instructions. The refactoring thus yields a better performance in comparison to the legacy implementation also on CPU-based supercomputers.

As part of the redesign, we also re-examined the numerical treatment of dynamic rupture simulations. The latter requires the evaluation of a friction law on the prescribed fault surface, which is discretized by a triangular mesh. The friction law is evaluated at quadrature points following a quadrature rule on the triangular surface cells. In addition to the originally used Stroud rule, we implemented the Dunavant rule, which achieves the same accuracy with less integration points, thus reducing the computational effort. The figure shows a strong scaling study from 800 to 6,336 nodes of SuperMUC-NG comparing the legacy Fortran version of SeisSol and the refactored version with both quadrature rules. While on smaller scales, the improvements are more clearly visible and the C++ version outperforms the legacy Fortran version, the picture is more complicated at large scale. We observe large variations in performance between different runs of the same configuration, which we attribute to effects such as power capping or dynamic changes of the clock frequency.