We present a highly scalable 3D fully-coupled Earth & ocean model of earthquake rupture and tsunami generation and perform the first fully coupled simulation of an actual earthquake-tsunami event and a 3D benchmark problem of tsunami generation by a megathrust dynamic earthquake rupture. Multi-petascale simulations, with excellent performance demonstrated on three different platforms, allow high-resolution forward modeling. Our largest mesh has ≈261 billion degrees of freedom, resolving at least 15 Hz of the acoustic wave field. We self-consistently model seismic, acoustic and surface gravity wave propagation in elastic (Earth) and acoustic (ocean) materials sourced by physics-based non-linear earthquake dynamic rupture, thereby gaining insight into the tsunami generation process without relying on approximations that have previously been applied to permit solution of this challenging problem. Complicated geometries, including high-resolution bathymetry, coastlines and segmented earthquake faults are discretized by adaptive unstructured tetrahedral meshes. This inevitably leads to large differences in element sizes and wave speeds which can be mitigated by ADER local time-stepping and a Discontinuous Galerkin discretization yielding high-order accuracy in time and space.

The numerical solution of partial differential equations is at the heart of many grand challenges in supercomputing. Solvers based on high-order discontinuous Galerkin (DG) discretisation have been shown to scale on large supercomputers with excellent performance and efficiency if the implementation exploits all levels of parallelism and is tailored to the specific architecture. However, every year new supercomputers emerge and the list of hardware-specific considerations grows simultaneously with the list of desired features in a DG code. Thus, we believe that a sustainable DG code needs an abstraction layer to implement the numerical scheme in a suitable language. We explore the possibility to abstract the numerical scheme as small tensor operations, describe them in a domain-specific language (DSL) resembling the Einstein notation, and to map them to small General Matrix-Matrix Multiplication routines. The compiler for our DSL implements classic optimisations that are used for large tensor contractions, and we present novel optimisation techniques such as equivalent sparsity patterns and optimal index permutations for temporary tensors. Our application examples, which include the earthquake simulation software SeisSol, show that the generated kernels achieve over 50% peak performance of a recent 48-core Skylake system while the DSL considerably simplifies the implementation.

We present a high-resolution simulation of the 2004 Sumatra-Andaman earthquake, including non-linear frictional failure on a megathrustsplay fault system. Our method exploits unstructured meshes capturing the complicated geometries in subduction zones that are crucial to understand large earthquakes and tsunami generation. These up-to-date largest and longest dynamic rupture simulations enable analysis of dynamic source effects on the seafloor displacements.