Generating high performance matrix kernels for earthquake simulations with viscoelastic attenuation

SeisSol is an earthquake simulator. It uses a discontinuous Galerkin method on unstructured grids with ADER time discretisation and achieves arbitrary high convergence order in space and time. In previous work, SeisSol was tuned for elastic rheological models on supercomputers with modern commodity hardware and shows high sustained peak performance and excellent parallel efficiency. We extend the high performance implementation to viscoelastic rheological models. These are commonly used to model the intrinsic attenuation of seismic waves, which is an important physical effect in the presence of high-frequency seismic waves. While the equation system gets three times larger, we show that we only have an overhead of about 75% compared to elastic rheological models. A strong scaling study shows that we achieve about 20% of peak performance. The numerical scheme leads to a series of small matrix chain products that have to be computed for each element. SeisSol uses a code generator for highly efficient matrix multiplications, however, several manual steps were necessary for tuning each matrix chain product. We present a framework that optimises code generation for matrix chain products in an almost fully automatic way.