ExaHyPE: An Engine for Parallel Dynamically Adaptive Simulations of Wave Problems

TL;DR

ExaHyPE is a versatile, high-performance engine for solving hyperbolic PDEs that supports adaptive mesh refinement and can be customized for various scientific applications, enabling rapid development of simulation codes.

Contribution

This paper introduces ExaHyPE, a novel engine that simplifies and accelerates the development of high-accuracy hyperbolic PDE simulations on supercomputers.

Findings

Supports adaptive mesh refinement for accuracy

Enables simulation of complex wave phenomena

Reduces development time from weeks to hours

Abstract

ExaHyPE ("An Exascale Hyperbolic PDE Engine") is a software engine for solving systems of first-order hyperbolic partial differential equations (PDEs). Hyperbolic PDEs are typically derived from the conservation laws of physics and are useful in a wide range of application areas. Applications powered by ExaHyPE can be run on a student's laptop, but are also able to exploit thousands of processor cores on state-of-the-art supercomputers. The engine is able to dynamically increase the accuracy of the simulation using adaptive mesh refinement where required. Due to the robustness and shock capturing abilities of ExaHyPE's numerical methods, users of the engine can simulate linear and non-linear hyperbolic PDEs with very high accuracy. Users can tailor the engine to their particular PDE by specifying evolved quantities, fluxes, and source terms. A complete simulation code for a new hyperbolic PDE can often be realised within a few hours - a task that, traditionally, can take weeks, months, often years for researchers starting from scratch. In this paper, we showcase ExaHyPE's workflow and capabilities through real-world scenarios from our two main application areas: seismology and astrophysics.

Tables

Table 1. Table 1: An overview of the main kernels implementing the algorithm, all are implemented for 2D and 3D and in linear and non-linear versions.

ADER-DG
spaceTimePredictor	Linear: Cauchy-Kowalevsky,
	Nonlinear: weak element-local space-time DG
Integrals	Implements face, surface, and volume integration
FV
Godunov	Implements the Godunov FV scheme
MUSCL-Hancock	Implements the MUSCL-Hancock FV scheme
slopeLimiter	Various slope limiters , default: minmod
Both (ADER-DG + FV)
solutionUpdate	Updates the solution vector
riemannSolver	Riemann solver, default: Rusanov
stableTimeStepSize	Calculates next stable time step size
Limiter
projectOnFVLimiterSpace	Projects DG solution on FV subcells
projectOnDGSpace	Projects FV solution back onto DG grid
discreteMaximumPrinciple	Limiter criterion, checks discrete maximum principle

Table 2. Table 2: An overview over the optional and required user-implemented functions.

Required Functions
adjustPointSolution	Set values in the state vector.
	Required: initial values. Optional: constraints
eigenvalues	Eigenvalues to calculate time step restriction
boundaryValues	Define boundary conditions
Problem-dependent Functions
flux	flux function
nonConservativeProduct	non-conservative products
source, pointSource	source terms
multiplyMaterialParameterMatrix	material parameters
viscousFlux	Flux terms with access to gradients
viscousEigenvalues	Eigevalues for time step restriction in viscous terms
Optional Functions
init	Initialisation of external libraries
refinementCriterion	Only for applications with AMR
	define additional areas to refine
isPhysicallyAdmissable	Only for applications with Limiter
	Gives physical admissibility criteria

Table 3. Table 3: Material parameters of the LOH.1 benchmark.

$x$	$c_p\left[m/s\right]$	$c_s\left[m/s\right]$	$\rho \left[kg/m^3\right]$
	$4000$	$2000$	$2600$
$\ge 1km$	$6000$	$3464$	$2700$

Table 4. Table 4: Convergence of the ADER-DG scheme for the Michel accretion test at t = 1.0 𝑡 1.0 t=1.0 .

	$N=2$		$N=3$		$N=6$
#cells	Error $L^1$	Order	Error $L^1$	Order	Error $L^1$	Order
$3^3$	2.56e-02		2.31e-02		3.23e-02
$9^3$	4.22e-03	1.65	2.27e-03	2.11	2.28e-05	6.45
${27}^3$	4.66e-04	2.01	3.09e-05	3.91	6.24e-08	5.52
${81}^3$	4.06e-05	2.21	4.47e-07	3.86	1.09e-11	7.88

Table 5. Table 5: Hybrid scaling of the GRMHD application. A N = 6 𝑁 6 N=6 ADER-DG approximation plus Finite Volumes limiting is used. Up to two levels of adaptive mesh refinement are employed.

		regular		$\Delta l=1$		$\Delta l=2$
# cores (nodes)	Threads	time	efficiency	time	efficiency	time	efficiency
56 (2)	14	91.93	1.00	139.43	1.00	395.90	1.00
112 (4)	7	37.68	1.22	57.03	1.29	151.11	1.31
224 (8)	4	30.88	0.77	34.17	0.89	139.42	0.93
448 (16)	2	17.67	0.43	26.66	0.65	57.02	0.85
784 (28)	1	13.88	0.32	22.89	0.44	34.17	0.67

Equations

$$\frac{\partial }{\partial t}Q+\nabla \cdot F(Q,\nabla Q)+B(Q)\cdot \nabla Q=S(Q)+\sum _{i=1}^{n_{\text{ps}}}{\delta }_i,$$

$$\frac{\partial }{\partial t}\underset{=Q}{\underbrace{\left(\begin{array}{c}\sigma \\ \rho v\end{array}\right)}}+\underset{=B(Q)\cdot \nabla Q}{\underbrace{\left(\begin{array}{cc}E(\lambda ,\mu )& 0\\ 0& 0\end{array}\right)\cdot \nabla \left(\begin{array}{c}v\\ \sigma \end{array}\right)}}+\nabla \cdot \underset{=F(Q)}{\underbrace{\left(\begin{array}{c}0\\ \sigma \end{array}\right)}}=0,$$

$$\frac{\partial }{\partial t}\underset{=Q}{\underbrace{\left(\begin{array}{c}h\\ hu\\ hv\\ b\end{array}\right)}}+\nabla \cdot \underset{=F(Q)}{\underbrace{\left(\begin{array}{cc}hu& hv\\ h{u}^2& huv\\ huv& h{v}^2\\ 0& 0\end{array}\right)}}+\underset{=B(Q)\cdot \nabla Q}{\underbrace{\left(\begin{array}{c}0\\ hg{\partial }_x(b+h)\\ hg{\partial }_y(b+h)\\ 0\end{array}\right)}}=0,$$

$$\frac{\partial }{\partial t}\underset{=Q}{\underbrace{\left(\begin{array}{c}\rho \\ \rho v\\ \rho E\end{array}\right)}}+\nabla \cdot \underset{=F(Q,\nabla Q)}{\underbrace{\left(\begin{array}{c}\rho v\\ v\otimes \rho v+Ip+\sigma (Q,\nabla Q)\\ v\cdot \left(I\rho E+Ip+\sigma (Q,\nabla Q)\right)-\kappa \nabla (T)\end{array}\right)}}\underset{=S(Q)}{\underbrace{=\left(\begin{array}{c}0\\ -gk\rho \\ 0\end{array}\right)}},$$

$$\frac{\partial }{\partial t}\underset{Q}{\underbrace{\left(\begin{array}{c}\sqrt{\gamma }D\\ \sqrt{\gamma }{S}_j\\ \sqrt{\gamma }\tau \\ \sqrt{\gamma }{B}^j\\ \varphi \\ {\alpha }_j\\ \beta \\ {\gamma }_m\end{array}\right)}}+\nabla \cdot \underset{=F(Q)}{\underbrace{\left(\begin{array}{c}\alpha v^{i}D-{\beta }^{i}D\\ \alpha T_j^i-{\beta }^i{S}_j\\ \alpha (S^i-v^{i}D)-{\beta }^i\tau \\ (\alpha v^i-{\beta }^i)B^j-(\alpha v^j-{\beta }^j)B^i\\ 0\\ 0\\ 0\\ 0\end{array}\right)}}+\underset{=B(Q)\cdot \nabla Q}{\underbrace{\left(\begin{array}{c}0\\ \sqrt{\gamma }(\tau {\partial }_j\alpha -\frac{1}{2}{T}^{ik}{\partial }_j{\gamma }_{ik}-T_i^j{\partial }_j{\beta }^i)\\ \sqrt{\gamma }(S^j{\partial }_j\alpha -\frac{1}{2}{T}^{ik}{\beta }^j{\partial }_j{\gamma }_{ik}-T_i^j{\partial }_j{\beta }^j)\\ -{\beta }^j{\partial }_i(\sqrt{\gamma }{B}^i)+\alpha \sqrt{\gamma }{\gamma }^{ji}{\partial }_i\varphi \\ \sqrt{\gamma }\alpha c_h^2{\partial }_j(\sqrt{\gamma }{B}^i)-{\beta }^j{\partial }^j\varphi \\ 0\\ 0\\ 0\end{array}\right)}}=0,$$

$$\frac{\partial Q}{\partial t}+\nabla \cdot F(Q)=0\mathrm{on}\mathrm{\Omega }\subset {\mathbb{R}}^{\mathrm{d}},\mathrm{d}=2,3.$$

$$Q(x,0)=Q_0(x),\forall x\in \Omega ,Q(x,t)=Q_B(x,t),\forall x\in \partial \Omega ,\forall t\in {\mathbb{R}}_0^{+}.$$

$${\int }_{T_i}{\int }_{t^n}^{t^{n+1}}{\theta }_h\frac{\partial q_h}{\partial t}\text{d}x\text{d}t+{\int }_{T_i}{\int }_{t^n}^{t^{n+1}}{\theta }_h\nabla \cdot F(q_h)\text{d}x\text{d}t=0.$$

$${\int }_{T_i}{\theta }_h\Delta q_h\text{d}x=-{\int }_{T_i}{\int }_{t^n}^{t^{n+1}}\nabla {\theta }_h\cdot F(q_h)\text{d}x\text{d}t+{\int }_{\partial T_i}{\int }_{t^n}^{t^{n+1}}{\theta }_{h}G(q_h,F(q_h))\text{d}s\text{d}t$$

$$\Delta t\le \frac{{\text{CFL}}_N}{d(2N+1)}\frac{h}{|{\lambda }_{\text{max}}|},$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ j\\ E\end{array}\right)+\nabla \cdot \left(\begin{array}{c}j\\ \frac{1}{\rho }j\otimes j+pI\\ \frac{1}{\rho }j(E+p)\end{array}\right)=0,$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\alpha \rho \\ \alpha j\\ \alpha E\\ \alpha \end{array}\right)+\nabla \cdot \left(\begin{array}{c}\alpha j\\ \frac{\alpha }{\rho }j\otimes j+\alpha pI\\ \frac{\alpha }{\rho }j(E+p)\\ 0\end{array}\right)+\left(\begin{array}{c}0\\ -p\nabla \alpha \\ 0\\ 0\end{array}\right)=0,$$

$$p=(\gamma -1)\left(E-\frac{1}{2\rho }\Vert j{\Vert }_2^2\right),$$

$$j(x,0)=(1,0),\rho (x,0)=1.0,E(x,0)=2.5,$$

$$j(x,0)=0,\rho (x,0)=\{\begin{array}{ll}1& \text{if }\Vert x-x_0{\Vert }^2<r^2,\\ \frac{1}{8}& \text{else,}\end{array},E(x,0)=\{\begin{array}{ll}1& \text{if }\Vert x-x_0{\Vert }^2<r^2,\\ \frac{1}{10}& \text{else.}\end{array}$$

$$\begin{array}{r}\frac{\partial \mathbf{\sigma }}{\partial t}-E(\lambda ,\mu )\cdot \frac{1}{\alpha }\nabla (\alpha \mathbf{v})+\frac{1}{\alpha }E(\lambda ,\mu )\cdot \mathbf{v}\otimes \nabla \alpha =0,\\ \frac{\partial \alpha \mathbf{v}}{\partial t}-\frac{\alpha }{\rho }\nabla \cdot \mathbf{\sigma }-\frac{1}{\rho }\sigma \nabla \alpha =0,\end{array}$$

$$\frac{\partial \alpha }{\partial t}=0,\frac{\partial \lambda }{\partial t}=0,\frac{\partial \mu }{\partial t}=0,\frac{\partial \rho }{\partial t}=0.$$

$$\frac{d\alpha }{dt}+v\nabla \alpha =0.$$

$$\frac{\partial \sigma }{\partial x}=\frac{1}{J}\left(\frac{\partial }{\partial q}\left(J{q}_x\sigma \right)+\frac{\partial }{\partial r}\left(J{r}_x\sigma \right)+\frac{\partial }{\partial s}\left(J{s}_x\sigma \right)\right)$$

$$\frac{\partial v}{\partial x}=q_x\frac{\partial v}{\partial q}+r_x\frac{\partial v}{\partial r}+s_x\frac{\partial v}{\partial s}$$

$$Q=\left(\begin{array}{c}h\\ hu\\ hv\\ b\end{array}\right)=\left(\begin{array}{c}\max \left\{0,\frac{1}{10}\left(x_0\cos (\omega \cdot t)+x_1\sin (\omega \cdot t)+\frac{3}{4}\right)-b\right\}\\ \frac{1}{2}\omega \sin (\omega \cdot t)h\\ \frac{1}{2}\omega \cos (\omega \cdot t)h\\ \frac{1}{10}(x_0^2+x_1^2)\end{array}\right),$$

$$\begin{array}{rl}\rho (x,t)& =1,\\ v(x,t)& =\exp (-1\mu t)\left(\begin{array}{c}\sin (z)+\cos (y)\\ \sin (x)+\cos (z)\\ \sin (y)+\cos (x)\end{array}\right),\\ p(x,t)& =-\exp (-2\mu t)\left(\cos (x)\sin (y)+\sin (x)\cos (z)+\sin (z)\cos (y)\right)+C,\end{array}$$

$$\frac{\partial }{\partial z}p\left(z\right)=-g\rho (z).$$

$${\Theta }^{\prime }=\{\begin{array}{ll}A& r\le a,\\ A\exp \left(-\frac{(r-a{)}^2}{s^2}\right)& r>a,\end{array}$$

$$\displaystyle\begin{aligned} &\text{warm:}\qquad&&A=$0.5\text{\,}$,\quad&&a=$150\text{\,}\mathrm{m}$,\quad&&s=$50\text{\,}\mathrm{m}$,\quad&&x_{c}=$500\text{\,}\mathrm{m}\mathrm{,}$\quad&&z_{c}=$300\text{\,}\mathrm{m}$,\\ &\text{cold:}\qquad&&A=$-0.15\text{\,}$,\quad&&a=$0\text{\,}\mathrm{m}$,\quad&&s=$50\text{\,}\mathrm{m}$,\quad&&x_{c}=$560\text{\,}\mathrm{m}$,\quad&&z_{c}=$640\text{\,}\mathrm{m}$.\end{aligned}$$

Full text

††journal: Computer Physics Communications

**PROGRAM SUMMARY

**

Program Title: ExaHyPE-Engine

Licensing provisions: BSD 3-clause

Programming languages: C++, Python, Fortran

*Nature of Problem:

*The ExaHyPE PDE engine offers robust algorithms to solve linear and non-linear hyperbolic systems of PDEs written in first order form. The systems may contain both conservative and non-conservative terms.

*Solution method:

ExaHyPE employs the discontinuous Galerkin (DG) method combined with explicit one-step ADER (arbitrary high-order derivative) time-stepping. An a-posteriori limiting approach is applied to the ADER-DG solution, whereby spurious solutions are discarded and recomputed with a robust, patch-based finite volume scheme. ExaHyPE uses dynamical adaptive mesh refinement to enhance the accuracy of the solution around shock waves, complex geometries, and interesting features.*

1 Introduction

The study of waves has always been an important subject of research. It is not difficult to see why: Earthquakes and tsunamis have a direct and serious impact on the daily lives of millions of people. A better understanding of electromagnetic waves has enabled ever faster wireless communication. The study of gravitational waves has allowed new insight into the composition and history of our Universe. Those physical phenomena, despite arising in different fields of physics and engineering, can be modelled in a similar way from a mathematical perspective: as a system of hyperbolic partial differential equations (PDE). The consortium behind the ExaHyPE project (“An Exascale Hyperbolic PDE Engine”) translated this structural similarity into a software engine for modelling and simulating a wide range of hyperbolic PDE systems.

ExaHyPE is intended as an engine, i.e. it allows only a limited number of numerical schemes on a fixed mesh infrastructure, but provides high flexibility in terms of the PDE system to be solved. The consortium focuses on two challenging scenarios, long-range seismic risk assessment, see e.g. [1, 2]; and the search for gravitational waves emitted by binary neutron stars, see e.g [3, 4, 5, 6].

A user with a given application only needs to implement the PDE system and problem-specific well-posed initial and boundary conditions. To exploit dynamic mesh refinement the user only needs to implement suitable criteria for mesh refinement and admissibility of solutions. From a user’s perspective, complicated PDE systems can be implemented without considering the complex issues of designing a performance-oriented high-order solver on a parallel compute cluster. Thus, the Engine enables medium-sized interdisciplinary research teams to quickly realize extreme-scale simulations of grand challenges modelled by hyperbolic conservation laws.

ExaHyPE implements a high-order discontinuous Galerkin (DG) approach. DG schemes were first introduced by Reed et al. [7] for the neutron transport equation and subsequently extended to general hyperbolic systems in a series of papers by Cockburn et al. [8, 9, 10, 11, 12]. Within our DG framework, higher-order accuracy in time and space is achieved using the Arbitrary high-order DERivative (ADER) approach first introduced for purely linear constant-coefficient equations in [13] and then extended to non-linear systems in [14]. In the original ADER approach, high-order accuracy in time is achieved by using the Cauchy-Kowalevsky procedure to replace time derivatives with spatial derivatives. This procedure is very efficient for linear problems [15], but becomes cumbersome for non-linear problems. Moreover, this procedure cannot deal with stiff source terms. To tackle these shortcomings, an alternative ADER approach was introduced by Dumbser et al. [16]. In this approach, the Cauchy-Kowalevsky procedure is replaced by an implicit solution of a cell-local space-time weak formulation of the PDE. This removes the problem dependency of the approach and allows the handling of stiff source terms. ExaHyPE provides an implementation of both ADER-DG variants.

In non-linear hyperbolic PDE systems, discontinuities and steep gradients can arise even from smooth initial conditions. High-order methods may then produce spurious oscillation which decrease the approximation quality and may even render the computed solution unphysical. To remedy this issue, a wide range of limiters for DG schemes have been proposed. Notably, these include approaches based on artificial viscosity [17, 18], filtering [19] and WENO- or HWENO-based reconstruction [20, 21]. The approach taken within the ExaHyPE engine is multi-dimensional optimal-order-detection (MOOD). This approach was initially applied to finite volume schemes [22, 23, 24] and has recently been extended to DG schemes [25]. In this approach, the solution is checked a-posteriori for certain admissibility or plausibility criteria and is marked troubled if it does not meet them. Troubled cells are then recalculated with a more robust finite volume scheme. The MOOD approach bypasses limitations of a-priori detection of troubled zones. We identify problematic regions after each time step and roll back to a more robust scheme on demand. This means that we do not need to reliably know all areas with issues a-priori. This approach allows for good resolution of shocks and other discontinuities [26].

The restriction to an adaptive cartesian mesh represents one of ExaHyPE’s fundamental design choices. Problems are discretised on tree-structured fully adaptive Cartesian meshes provided by the Peano framework [27, 28]. Dynamical adaptive mesh refinement (AMR) enhances the shock-capturing abilities of the ADER-DG scheme further [29] and allows good resolution of local features. ExaHyPE allows for complex geometry, either handled by a curvilinear mesh transformation, or by a diffuse interface approach. While the former approach is linear operation it is restricted to geometries that can be mapped smoothly to a cube. The latter approach extends a given PDE system by a parameter representing the volume fraction of material in the cell and thus determines the physical boundary through a diffuse interface instead of boundary-fitted unstructured meshes. Complex geometries can thus be readily represented and extensions to moving geometries are also possible [2, 30, 31].

ExaHyPE comprises the following key features:

High-order ADER discontinuous Galerkin (ADER-DG) with a-posteriori subcell limiting and finite volume (FV) schemes;

• 2.

Dynamic mesh refinement on Cartesian grids in two and three dimensions;

• 3.

A simple API that allows users to quickly realise complex applications;

• 4.

User-provided code can be written in Fortran or C++;

• 5.

Automatically generated architecture- and application-specific optimised ADER-DG routines;

• 6.

Shared memory parallelisation through Intel’s Threading Building Blocks (TBB);

• 7.

Distributed memory parallelisation with MPI.

Furthermore, ExaHyPE offers a wide range of post-processing and plotting facilities such as support for the output formats vtk or tecplot. The software can be compiled with Intel and GNU compilers and provides a switch for choosing different compilation modes: The release mode aggressively optimises the application, while the assertion and debug compilation modes activate the assertions within the code and print additional output. Users and developers can write log filters to filter output relevant to them. For continuous testing a Jenkins server [32] verifies that the code compiles with all parallelisation features in all different compilation modes including its ability to resolve complex geometries.

Our paper starts with an overview of the problem setting using several examples. We then briefly sketch the solution methods used in ExaHyPE, focusing on the ADER-DG algorithm. Next we use a simple example to demonstrate the workflow when using ExaHyPE and describe the engine architecture. We conclude with a sequence of numerical examples from various application areas to demonstrate the capabilities of the engine.

2 Problem Formulation

We consider hyperbolic systems of balance laws that may contain both conservative and non-conservative terms. They have to be given in the following first-order form:

[TABLE]

where $Q:\Omega\subset\mathbb{R}^{d}\mapsto\mathbb{R}^{\nu}$ is the state vector of the $\nu$ conserved variables, $\Omega\subset\mathbb{R}^{d}$ is the computational domain, $F(Q)$ is the flux tensor that may also depend on the gradient of Q in order to model viscous effects, and $B(Q)$ represents its non-conservative part. Finally, $S(Q)$ is the source term and $\delta_{i}$ are the given $n_{\text{ps}}$ point sources.

Hyperbolic systems in the form (1) can be used to model a wide range of applications that involve waves. In the following, we demonstrate the versatility of our formulation by introducing equations from three different application areas. Numerical experiments for these examples are provided in Section 6.

2.1 Waves in Elastic Media

Linear elastodynamics describes waves propagating through elastic heterogeneous media by relating displacements, velocities, stress and strain. The momentum equations of motion are derived from Hooke’s law and the conservation of momentum. Following the form of equation (1) we write them as

[TABLE]

where $\rho$ denotes the mass density, $v$ the velocity and the $\sigma$ stress tensor, which can be written in terms of its six independent components as $\sigma=(\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy},\sigma_{xz},\sigma_{yz})$. In this paper we consider isotropic materials, i.e. the material matrix $E(\lambda,\mu)$ depends only on the two Lamé constants $\lambda$ and $\mu$ of the material. However, the formulation holds also for more general anisotropic materials.

These equations can be used to simulate seismic waves, such as those radiated by earthquakes. In this context the restriction of ExaHyPE to Cartesian meshes seems to be restrictive. However, adaptive Cartesian meshes can be extended to allow the modelling of complex topography. Two methods have been implemented in ExaHyPE to represent complex topographies. The first approach treats ExaHyPE’s adaptive Cartesian mesh as reference domain that is mapped to a complex topography via high-order curvilinear transformations [33, 34]. The second approach, a diffuse interface method, represents the topography as a smooth field [2]. These approaches are described in more detail in Section 6.2.

2.2 Shallow Water Equations

In atmospheric and oceanic modelling of coastal areas, horizontal length scales are typically significantly greater than the vertical length scale. In this case, fluid flow can be modelled with the two-dimensional shallow water equations instead of the more complicated three-dimensional Navier-Stokes equations.

Following the form of equation (1) they can be written as

[TABLE]

where $h$ denotes the height of the water column, $(u,v)$ the horizontal flow velocity, $g$ the gravity and $b$ the bathymetry.

Hyperbolic systems of balance laws have non-trivial equilibrium solutions in which flux and source terms cancel. A well-balanced numerical scheme is capable of maintaining such an equilibrium state. In Section 3 we describe the ADER-DG scheme used in ExaHyPE and in Section 6 we will describe how to keep the scheme well balanced while allowing wetting and drying.

2.3 Compressible Navier-Stokes Equations

ExaHyPE is extensible to non-hyperbolic equations. As an example of such an extension we show the compressible Navier-Stokes equations. They are used to model the dynamics of a viscous fluid, and are given by

[TABLE]

where $\rho$ denotes the density, $\rho v$ the momentum, $\rho E$ the energy density, $T$ the temperature and $p$ the pressure (this term includes gravitational effects). The temperature diffusion is given by $\kappa\nabla T$ with constant $\kappa$, these effects depend on the gradient of $Q$. In the source term, the vector $k$ is the unit vector in $z$-direction and $g$ is the gravitation of Earth. The viscous effects are modelled by the stress tensor $\sigma(Q,\nabla Q)$. Thus, the flux from (1) has been modified to allow $\nabla Q$ as input. More details on the implementation of these equations can be found in [35].

2.4 General Relativistic Magneto-Hydrodynamics

The equations of classical magnetohydrodynamics (MHD) are used to model the dynamics of an electrically ideally conducting fluid with comparable hydrodynamic and electromagnetic forces. When modelling astrophysical objects with strong gravitational fields, e.g. neutron stars, it becomes necessary to model the background space-time as well. We use the standard $3+1$ split to decompose the four dimensional space-time manifold into 3D hyper-surfaces parameterised by a time coordinate $t$. The background space-time is introduced into the equations in the form of a non-conservative product.

Following the form of equation (1) they can be written as

[TABLE]

where $i,j=1,2,3$ and $m=1...6$.

The curved space-time is parameterised by several hyper-surface variables: lapse $\alpha$, spatial metric tensor $\gamma$ shift vector $\beta$ and extrinsic curvature $K$. The spatial metric tensor is given as a vector of its six independent components $\gamma=(\gamma_{11},\gamma_{12},\gamma_{13},\gamma_{22},\gamma_{23},\gamma_{33})$ and has the determinant $\sqrt{\gamma}:=\det{\gamma}$. Further, $D=W\rho$ is the conserved density, which is related to the rest mass density $\rho$ by the Lorentz factor $W$, $v^{i}$ is the fluid velocity, $T$ is the Maxwell 3-energy momentum tensor, $S$ is the conserved momentum, $B$ is the magnetic field and $\tau$ is the conserved energy density. Finally, $\phi$ is an artificial scalar introduced to ensure a divergence-free magnetic field, and $c_{h}$ is the characteristic velocity of the divergence cleaning. For more details on this formulation see e.g. [36, 37, 38].

3 Solver Components

In ExaHyPE’s engine concept the numerical method is given and in general the user does not need to interact with any solver components. A layer of generated glue code separates the user from the kernel calls. The main kernels are given in Table 1, of these only the slope-limiter and Riemann solver can be modified by the user. However, all modifications to the kernels are optional since stable defaults are provided. This section briefly summarises the numerical algorithms used. For brevity, let us use a pared down form of (1), in which only the flux term is nontrivial:

[TABLE]

Assume that (1) is subject to appropriate initial and boundary conditions:

[TABLE]

3.1 Discretisation

ExaHyPE allows for Cartesian grids in two and three dimensions. The computational domain is divided into a grid $\Omega=\bigcup_{i}T_{i}$ using a space-tree construction scheme [27, 39]. Each cell $T_{i}$ is recursively refined giving an adaptive Cartesian grid as shown in Figure 1. Meshes can be defined with an arbitrary numbers of elements in each direction on the coarsest level. Refinement is based on tripartitioning.

The rationale behind a commitment to tripartitioning is two-fold. On the one hand, we rely on the Peano AMR framework [39] which internally orders all cells along a Peano space-filling curve. This yields a high temporal and spatial locality of the cell and face accesses; a property the Peano curve shares with other space-filling curves. However, the Peano curve also allows us to realise all temporary data structures through cache-friendly stacks [27]. On the other hand, tripartitioning is an intrinsic match for patch-based Finite Volume methods which work with odd numbers of volumes per axis: A subdivision of a cell here makes a former volume centre coincide with a volume centre on the next finer level. This simplifies the structure of inter-resolution transfer operators and preserves them over many resolution levels.

3.2 Finite Volumes

Classically systems of hyperbolic PDE have been solved using finite difference or finite volume schemes. In a finite volume method cell averages are calculated. Volume integrals in a PDE that contain a divergence term are converted to surface integrals using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each element. However, to achieve high order accuracy in a finite volume scheme, large stencils and expensive recovery or reconstruction procedures are needed. Examples include essentially non-oscillatory (ENO) or weighted ENO (WENO) schemes, see e.g. [40, 41].

As shown in Table 1 ExaHyPE provides access to classical Godunov type schemes as well as MUSCL schemes [42, 43]. Further, users can use parts of the generic compute routines while implementing other parts on their own. For instance, in both the ADER-DG and Finite Volume schemes, users can overwrite the Riemann Solver with their own implementation.

3.3 ADER-DG

Instead of representing the solution as cell averages, DG methods represent the solution within each cell by a (high-order) polynomial. The ADER-DG method is a one-step predictor-corrector scheme. The integration in time is initially performed only within the element, neglecting the element interfaces and then a single correction step is performed to take the element interfaces into account [44].

The method operates on a weak form of equation (1)

[TABLE]

This formulation has replaced the solution $Q$ with a discrete function $q_{h}$, represented by Lagrange polynomials $\theta_{h}$ in time and space. This solution is refered to as the space-time predictor and intuitively corresponds to a Taylor-expansion in space and time. Hence, it is equal to the time-evolution of the discrete solution within the control element for a smooth solution. The space-time test function is in the space of piecewise polynomials, which is constructed as tensor products of Lagrange polynomials over Gauss-Legendre or Gauss-Lobatto points.

Integrating over each element $T_{i}$ and over the current time interval $[t^{n},t^{n+1}]$ then gives the element-local weak formulation. Note that the formulation does not take into account any information from neighboring elements. The predictor is subsequently corrected using contributions from neighboring cells using a Riemann-solver.

However, first we need to compute $q_{h}$. This can be evaluated directly by the Cauchy-Kovalevsky algorithm for linear problems. For non-linear models, we use the ADER-DG method proposed by Dumbser et al. [45]. It is a fixed-point iteration and can be seen as a discrete counterpoint to the well-known Picard iteration, see [26].

There are three phases per ADER-DG time step:

Per grid cell $T_{i}$ and time interval $[t^{n},t^{n+1}]$, we first implicitly solve (5), predicting the local evolution. The concurrent solves of (5) do not take into account any information from neighbouring elements and thus yield jumps along the cell faces in $q_{h}$ and $F(q_{h})$. 2. 2.

We traverse all faces of the grid and compute a numerical normal flux $G(q_{h},F(q_{h}))$ from both adjacent cells. ExaHyPE uses a Rusanov flux by default; users can replace it with any other Riemann solver. 3. 3.

In the corrector step, we traverse the cells again and solve

[TABLE]

for $\Delta q_{h}=q_{h}(t^{n+1})-q_{h}(t^{n})$, see [16]. The one-step update (6) is derived by multiplying with a spatial testfunction and partially integrating (1). Equation 6 can be easily inverted given that the ansatz and test space typically yield a diagonal mass matrix.

3.4 Time-Step Restrictions

Nonlinear effects and mesh adaptation require adjustments to the time step size during a simulation. This is expressed by the CFL condition, which gives an upper bound on the stable time step size for explicit DG schemes:

[TABLE]

where $h$ and $|\lambda_{\text{max}}|$ are the mesh size and the maximum signal velocity, respectively, and $\text{CFL}_{N}<1$ is a stability factor that depends on the polynomial order [16].

3.5 A-Posteriori Limiting

The unlimited ADER-DG algorithm will suffer from numerical oscillations (Gibbs phenomenon) in the presence of steep gradients or shock waves. Therefore a limiter must be applied. The approach followed in ExaHyPE is based on the a-posteriori MOOD method of Loubère et al. [25]. The solution is checked a-posteriori for certain admissibility or plausibility criteria and is recalculated with a robust FV scheme if it does not meet them. The FV patch size is chosen to have an order of $2N+1$, this is the smallest cell size that does not violate the CFL condition (7).

In contrast to the original approach, ExaHyPE’s approach incorporates the observation that cells usually require a recalculation with FV multiple time steps in a row after the initial check failed. Therefore, ExaHyPE implements the a-posteriori limiting ADER-DG method as a hybrid ADER-DG-FV method (Figure 2).

As a-posteriori detection criteria we use

Physical Admissibility Criteria: Depending on the system of PDEs being studied certain physical constraints can be placed on the solution. For the shallow water equations these are positivity of the water height. These criteria are supplied by the user along with the PDE terms. 2. 2.

Numerical Admissibility Criteria: To identify shocks we use a relaxed discrete maximum principle (DMP) in the sense of polynomials, for details see e.g. [25].

We call DG cells that do not satisfy the above criteria troubled cells. If a cell is flagged as troubled and has not been troubled in the previous time step, the scheme goes back to the old time step and recomputes the solution in all troubled cells (and their direct neighbours) with the FV scheme.

4 Engine Architecture and Programming Workflow

This section briefly walks through the workflow of setting up an application in ExaHyPE using the shallow water equations given in (2) as an example. The architecture of ExaHyPE is illustrated in Figure 3. ExaHyPE is a solver engine, domain-specific code has to be written by the user to obtain simulation code. In Figure 3 a turquoise colour is used to highlight files written by the user. To write an ExaHyPE application users typically start from a specification file. The specification file is passed to the ExaHyPE toolkit, which creates glue code, empty application-specific classes and optionally application and architecture tailored core routines. The application specific classes are filled by the user with the PDE terms. This code can be written in C++ or Fortran. The generated glue code and the initially empty templates make up the ExaHyPE user solver.

The ExaHyPE core relies on the Peano framework (green) for its dynamically adaptive Cartesian meshes. In addition, it provides an efficient mesh traversal loop that ExaHyPE’s algorithms plug into. Peano itself is a third-party component.

The number of dependencies in the ExaHyPE core is minimal. However, the architecture may not fulfil the requirements of all applications, so the user can extend and connect further software fragments to the ExaHyPE core. In red we show an optional dependency, libxsmm [46]. This package provides efficient kernels for small matrix multiplications, which are used by the optimised ADER-DG routines described in Section 5.2. Similarly, further software packages can be added by the user as needed.

In Table 2, a brief summary of the components of the user solver is given. In the following Section we highlight the parts of the specification file that need to be modified for each component and show the implementation of the user functions for the shallow water equations. The solver components can be set up flexibly by modifying the specification file. Users only need to write application-specific code that sets up their PDE system.

4.1 Code Generation and Compilation

After the specification file has been written it is handed over to the ExaHyPE toolkit. A minimal ExaHyPE specification file is shown below.

exahype-project SWE output-directory const = ./SWE

computational-domain dimension const = 2 width = 1.0, 1.0 offset = 0.0, 0.0 end-time = 1.0 end computational-domain

solver ADER-DG MySWESolver variables const = h:1,hu:1,hv:1,b:1 order const = 3 maximum-mesh-size = 0.1 time-stepping = global type const = non-linear terms const = flux,ncp optimisation const = generic, usestack end solver end exahype-project

To prepare this example for the simulation, run

./Toolkit/toolkit.sh Demonstrators/SWE.exahype

The toolkit generates a Makefile, glue code, as well as various helper files. Among them is one C++ class per solver that was specified in the specification file. Within each implementation file, the user can specify initial conditions, mesh refinement control, etc. For example, to set the eigenvalues the following function is generated in the file MySWESolver.cpp.

void SWE::MySWESolver::eigenvalues(const double* const Q, const int direction, double* const lambda) { // @todo Please implement/augment if required lambda[0] = 1.0; lambda[1] = 1.0; lambda[2] = 1.0; lambda[3] = 1.0; }

This function can then be filled with the eigenvalues of the PDE system under consideration. Similarly, the other functions flux, initial conditions and boundary conditions can be defined in the file. For the shallow water equations this would be:

void SWE::MySWESolver::eigenvalues(const double* const Q, const int direction, double* const lambda) { ReadOnlyVariables vars(Q); Variables eigs(lambda);

const double c = std::sqrt(gravity*vars.h()); double u_n = Q[direction + 1] * 1.0/vars.h();

eigs.h() = u_n + c; eigs.hu() = u_n - c; eigs.hv() = u_n; eigs.b() = 0.0; }

In this example we have used named variables to enhance readability. However, ExaHyPE also allows the user to access the vectors directly as can be seen in the generated function above.

The whole build environment is generated. A simple make will create the ExaHyPE executable. ExaHyPE’s specification files always act as both specification and configuration file, i.e when running the code the specification file is passed in again.

./ExaHyPE-SWE ./SWE.exahype

A successful run yields a sequence of vtk files that you can open with Paraview or VisIt. In this example, we plot two quantities. Such an output is shown in Figure 4.

Global metrics such as integrals can also be realised using a plotter with no output variables, i.e. variable const = 0.

5 Parallelisation and Optimisation Features

ExaHyPE relies on the Peano framework for mesh generation. Peano implements cache-efficient, tree-structured, adaptive mesh refinement. To traverse the cells or vertices the mesh traversal automaton of Peano runs along the Peano space-filling-curve (SFC), see Figure 1. The action of a PDE operator is mapped to the SFC traversal by plugging into events triggered by the traversal automaton, these mappings can be generated by a toolkit.

Peano also takes care of the distributed-memory and shared-memory parallelisation. Domain decomposition for distributed-memory parallel simulations is realised by forking off or merging subtrees. The domain decomposition can be influenced via load balancing callbacks. The shared-memory parallelisation relies on identifying regular substructures in the tree and employing parallel-for loops in these areas [48]. Recently, we introduced a runtime tasking system to Peano and ExaHyPE that introduces additional multi-threading concurrency [49] in highly adaptive mesh regions and allows overlapping computations with MPI communication.

ExaHyPE builds on the Peano framework and it inherits a user model from Peano: Our engine removes the responsibility for algorithmic issues from the user. The time step, the mesh traversal and the parallelisation are implemented in a generic way. The user only controls the PDE system being solved, by specifying how many quantities are needed, which terms are required, and ultimately what all terms look like. Our goal is to allow users to focus on the physics only and to hide away as many implementation details as possible.

5.1 High-Level Optimisations

ExaHyPE contains several high level algorithmic optimisations that users can switch on and off at code startup through the specification file. To gain access to these optimisations, we add the following optional section to the minimal specification file shown in Section 4.

global-optimisation fuse-algorithmic-steps = all spawn-predictor-as-background-thread = on spawn-amr-background-threads = on end global-optimisation

A discussion of each individual algorithmic tuning is beyond scope of this paper. Some techniques are:

Step fusion: Each time step of an ExaHyPE solver consists of three phases: computation of local ADER-DG predictor (and projection of the prediction onto the cell faces), solve of the Riemann problems at the cell faces, and update of the predicted solution. We may speed up the code if we fuse these four steps into one grid traversal.

• 2.

Spawn background jobs: TBB has a thread pool in which idle threads wait for tasks, as soon as new tasks are spawned the threads in this pool start this task. We refer to these threads as background threads. Since spawning and scheduling these tasks has a overhead we wait until a certain number of tasks has accumulated and schedule these tasks together. Certain space-time-predictor computations can be spawned as a background job. Costly AMR operations such as the imposition of initial conditions and evaluation of refinement criteria can also be performed as a background jobs.

• 3.

**Modify storage precision: **ExaHyPE internally can store data in less-than-IEEE double precision. This can reduce the memory footprint of the code significantly.

Most users will not modify these options while they prototype. When they start production runs, they can tweak the engine instantiation through them. The rationale behind exposing these control values is simple: It is not clear at the moment which option combinations robustly improve performance. An auto-tuning/machine learning approach to find the optimal parameter combinations could be considered..

ExaHyPE allows running multiple simulations on the same computational mesh. This can be facilitated by adding multiple solver descriptions to the specification file. Each solver uses its own base grid and refinement criterion.

5.2 Optimised ADER-DG Routines

One of ExaHyPE’s key ideas is to use tailored, extremely optimised code routines whenever it comes to the evaluation of the ADER-DG scheme’s steps and other computationally expensive routines. AMR projections and the space-time predictor (5) are typically the dominant steps. If the user decides to use these optimised variants in the specification file, the toolkit calls a specific code generator Python3 module and links to its output in the generated glue code so that the calls to the generic routines are replaced by calls to the generated optimised ones.

SIMD operations, notably introduced with AVX-512 on KNL and Skylake, become increasingly critical to fully exploit the potential of modern CPUs. Therefore, on Intel machines, the optimised routines’ main goal is to either directly use SIMD or enable as much auto-vectorisation from the compiler as possible. To that end the optimised routines use Intel’s libxsmm [46], which is the second third-party building block in the basic ExaHyPE architecture. We map all tensor operations required by the ADER-DG algorithm to general matrix multiplications (gemms) and apply architecture-specific vectorised matrix operations. Furthermore, the code generation allows the introduction of data padding and alignment to get the most out of the compiler auto-vectorisation.

While improving SIMD vectorisation is our current prime use case for the optimised kernels, we provide alternative optimisation dimensions: Instead of a vectorisation of the native ADER-DG loops, the optimised routines can also be configured to work with vectorised PDE formulations. This implies that the PDE terms are computed through SIMD operations instead of scalar ones and requires some work from the user. This has been found to greatly improve performance [50].

Instead of the standard Picard iteration, we provide a continuous extension Runge-Kutta scheme (CERK) to yield an initial guess to the Picard iteration, leading to a much lower number of required iterations.

6 Numerical Results

In this section, we show the ExaHyPE engine in action. The applications in this section are taken from our introductory discussion in Section 2. Only the Euler equations are added as an all-time classic starting point for solvers of hyperbolic systems. All specification and source files used to generate the results in this section are made publicly available, and all of the tests themselves can be run on a standard laptop. The scaling tests shown require the use of a larger cluster.

6.1 Euler Equations

The non-linear compressible Euler equations model the flow of an inviscid fluid with constant density. Solutions of the Euler equations are sometimes used as approximations to real fluids problem, e.g. the lift of a thin airfoil. They are given by

[TABLE]

where $\rho$ denotes the mass density, $j\in\mathbb{R}^{d}$ denotes the momentum density, $E$ denotes the energy density, $p$ denotes the fluid pressure, to be given by an equation of state, and $\otimes$ denotes the outer product.

We extend the model by a color function denoting the volume fraction of material in a cell starting from the Baer-Nunziato model [51]:

[TABLE]

In the above PDE system $\alpha$ denotes the volume fraction of material present. The pressure is given by:

[TABLE]

where $\gamma=1.4$ is the ratio of specific heats .

Note that the addition of $\alpha$ has introduced a non-conservative term into the equation. In the above we have given a reduced Baer-Nunziato model, similar to the approach presented in [30, 31, 52, 2]. This model allows simulation of fluid-structure-interaction (FSI) problems on adaptive Cartesian grids, without requiring boundary-fitted grids.

In Figure 5 we present computational results for flow over a thin airfoil. Initial conditions were:

[TABLE]

As boundary conditions we set an inflow boundary at the left, an outflow boundary at the right and $v\cdot\nabla n=0$ at the top and bottom boundaries. Results are shown for a NACA 4612 airfoil [53]. The results generate the expected bow shock.

To demonstrate the shared-memory scalability of the code we use a variant of the Sod shock tube problem, the circular, spherical explosion problem, also referred to as a “multi-dimensional Sod shock tube”. This test case uses the following initial conditions

[TABLE]

and wall boundary conditions.

In Figure 6 we show the shared-memory scalability in two and three dimensions. This scaling test was run on the $14$ cores of an Intel Xeon E5-v3 processor with $2.2$ GHz core frequency (SuperMUC Phase 2). The tests were run on one processor with 14 cores using Intel’s TBB [54] for parallelisation. We consider a regular base grid and allow a dynamic refinement criterion to add up to two additional levels $\Delta\ell\in\{1,2\}$ of cells around the shock front. In these tests all kernel level optimisations that are available in ExaHyPE have been used. As a result, scalability is more challenging, however, we are interested mainly in reducing the overall run-time.

The scalability improves as the work per element increases, i.e. it improves with higher polynomial degree $N$ or for increasing dimension. In general good scalability can be observed on up to 14 cores at higher orders. All ExaHyPE codes employ a hybrid parallelisation strategy with at least two MPI ranks per node [55]. Results for this hybrid parallelisation strategy are provided in Section 6.6.

6.2 Elastic Wave Equation

The restriction of ExaHyPE to Cartesian meshes seems restrictive in the context of most seismic applications. In this section, we introduce two methods for incorporating complex geometries into such a mesh without modifying the underlying numerical method. We would like to highlight that both of these methods omit manual mesh-generation, which is typically required as a pre-processing step in computational seismology applications and poses a major bottleneck. Hexahedral mesh generation can easily consume weeks to months and is limited for complex geometries of boundary and interface conditions, while form-fitted unstructured tetrahedral meshes allow for automatised meshing and complex geometries [56, 57], however, pose numerical challenges, e.g. in form of misshaped sliver elements [58].

6.2.1 Diffuse Interface Approach

As in the Euler equations we can extend the elastic wave equation by a parameter $\alpha$, which reprents the volume fraction of the solid medium present in a control volume [2]. Diffuse interfaces completely avoid the problem of mesh generation, since all that is needed for the definition of the complex surface topography is to set a scalar colour function to unity inside the regions covered by the solid and to zero outside.

The diffuse interface model is given by:

[TABLE]

where ${E}(\lambda,\mu)$, $\lambda$, $\mu$, $\rho$ and the stress tensor $\bm{\sigma}$ are defined as in the Section 2.1.

The material parameters are assumed to remain constant, i.e.

[TABLE]

Physically, $\alpha$ represents the volume fraction of the solid medium present in a control volume. The equations become non-linear in those regions in which $0<\alpha<1$.

This formulation can be extended to allowing moving materials. Instead of solving $\frac{d\alpha}{dt}=0$, we can solve

[TABLE]

In this way the free surface boundary is allowed to move according the local velocity field [2].

To verify the accuracy of this method we solve the layer over homogeneous halfspace (LOH.1) benchmark problem described by Day et al. [59]. This problem is a well-known reference benchmark for seismic wave propagation in numerical codes. The LOH.1 benchmark considers way propagation in a hexahedral geometry filled with two materials that are stacked on top of each other. The first material is characterised by a lower density and smaller seismic wave speeds. The exact material parameters are defined in Table 3. The parameters $\lambda$ and $\mu$ of the equation can be derived from these values.

A point source is placed $2$ km below the surface at the centre of the domain, such that the resulting wave propagates through the change of material.

In the diffuse interface method we limit the free surface to be able to resolve the local discontinuity of $\alpha$. The left picture in Figure 7 shows the resulting distribution of limited (red) and unlimited areas (blue) and indicates the transition of both (light blue and red). The right picture shows the velocity field in x direction at $t=1.1$. Our implementation of the diffuse interface method is able to successfully resolve both the changing material parameters and the absorbing surface boundary conditions.

6.2.2 Curvilinear meshes

Our second approach models complex topographies by applying a curvilinear transformation to the elements of an adaptive Cartesian mesh. To do so we generate surface quadrature nodes depending on the topography, create a 2D curvilinear interpolation of those quadrature nodes on domain boundaries with topography curves and domain edges as constraints, and finally generate a Jacobian mapping $J$ at each node. Taking this mapping into account, the linear elastic wave equations can be written as:

[TABLE]

where $J$ is the Jacobian matrix of the mapping from mesh to topography. This approach allows us to model meshes with complex topographies including faults and inner branches [33, 34]. The ADER-DG algorithm as described previously, can then be applied directly to this version of the elastic wave equation.

Figure 8 shows a snapshot of a numerical experiment that simulates propagation of seismic waves in the area around the mountain Zugspitze, in Germany. Both the curvilinear approach and the diffuse interface method are able to resolve the complicated topographies in this scenario which is motivated by the AlpArray experiment.

Figure 9 shows the shared-memory scalability of the linear ADER-DG implementation on SuperMUC’s Phase 2 when running the LOH.1 benchmark comparing the scalability of the two meshing approaches. We consider a regular base grid with $27^{3}$ cells and allow a dynamic refinement criterion to add one additional level of cells to the layer over the halfspace. Memory constraints prevented further refinement. As in the non-linear test case the scalability improves as the work per element increases. However, due to the lower amount of total work in the linear setting a higher polynomial degree needs to be reached to attain scalability in this case. This means that the overall scalability of the diffuse interface approach is higher, however, this comes at the price of a higher overall computational cost due to the non-linearity of the formulation.

The question of which method should be used is highly problem dependent. The curvilinear method is purely linear and as such computationally cheap. The non-linear diffuse interface approach, on the other hand, requires additional limiting on the surface. However, the time step size for the DIM is independent of the simulated topography, while in the curvilinear method it is highly dependent on the distortion introduced by the topography. Simulations with a relatively smooth, flat surface are expected to be faster with the curvilinear method, while highly varying topographies are better discretised with the DIM. DIM has the added advantage of allowing moving free surface boundaries, a feature that is difficult to realise with curvilinear meshes.

6.3 Shallow Water Equations

We demonstrate the dynamic mesh refinement capabilities of ExaHyPE via the shallow water equations (2). Even more importantly, we present a non-toy problem which uses ExaHyPE’s 2D facilities. To solve (2), we use the a-posteriori limiting ADER-DG method and equip its FV limiter with a recently developed HLLEM Riemann solver [60]. The solver allows for wetting and drying. This also demonstrates that users can inject their own Riemann solvers in ExaHyPE.

The same scenario is available in ExaHyPE using a Finite Volume scheme instead of an ADER-DG scheme with FV limiter [47, 61]. In both implementations the same Riemann solver can be used. The Riemann solver considers a dry tolerance $\epsilon>0$ below which cells are marked as dry. The constant $\epsilon$ is necessary to avoid negative water heights, which usually lead to unphysical non-linear effects. If one of the cells is flooded the jump in bathymetry is taken into account by the Riemann solver. When reconstructing the DG solution from the FV limiter the water level is reconstructed first and then the bathymetry is subtracted. This ensures that the non-linear reconstruction process does not produce artificial waves. In tsunami simulations, this scheme allows us to simulate areas close to the coast with the FV subcell limiter, while areas in the deep ocean are processed by the high order ADER-DG method.

To test the shallow water equations given in (2) we use the oscillating lake scenario. Here, a water droplet travels in circular motion over a dry basin. The topography of the basin is resolved using the final variable $b$, which remains constant throughout the simulation. The analytical solution, which is used as an initial condition is given by

[TABLE]

where $\omega=\sqrt{0.2g}$ and $g=9.81$ denotes the gravitational constant. The initial conditions are supplemented by wall boundary conditions. In this setup we use polynomial order $N=3$ and a base grid of $7\times 7$ elements. We adaptively refine mesh elements where the water height is small but not zero. We allow up to two levels of adaptive refinement.

This scenario provides a challenging test case for numerical codes due to the continuous wetting and drying of cells. There exists an analytical solution for this benchmark scenario allowing us to verify the well-balancedness of a scheme, as well as in the resolution of drying or inundated cells.

In the top row of Figure 10, we plot the water height at $t=0,1,2,3$. Despite the frequent wetting and drying the scheme performs well. In those areas in which the DG solution is used, we note the accuracy of the higher order scheme on a coarse grid. In those areas in which a FV solution is necessary, the dynamic AMR is useful to retain accuracy and avoid diffusion. The bottom row of Figure 10 shows the locations in which the FV limiter is active and adaptive mesh refinement is used. Furthermore, it shows the error in the water height. We observe that the FV limiter is used only in those cells in which it is required, the limiter locations move with the solution. The error remains below $10^{-3}$ in all cells.

In Figure 11, we show the shared-memory scalability for the oscillating lake. This scaling test was run on one $2.2$ GHz $14$ core Intel Xeon E5-v3 processor of SuperMUC Phase 2. We consider a regular base grid and allow a dynamic refinement criterion to add up to two additional levels of cells around the water droplet. In this area the FV limiter is active. As in the previous test case the scalability improves for higher polynomial degree $N$.

6.4 General Relativistic Magneto-Hydrodynamics

To uncover and validate ExaHyPE’s high convergence order, we require a test case which does not contain discontinuities, as the limiter reduces the convergence order. Furthermore, the analytical solution has to be known. We simulate the spherical accretion onto a stationary black hole. This well-known test case has a known analytical solution [62]. Details of the setup in the context of ADER-DG methods can also be found in [6]. This setup is challenging since GRMHD is a non-linear equation of $19$ variables containing both fluxes and non-trivial non-conservative products. Further complexity arises from the the closeness of the singularity at the critical radius to the computational domain.

We use the analytical solution as the external state vector for the Riemann solver at the boundary of the domain $\partial\Omega$. The test is performed in Kerr-Schild coordinates on the spatial domain $[2.2,12.2]^{3}$, away from the critical radius.

In Figure 12 we show the solution at time $t=1.0$ and the convergence of the error in terms of $L^{1}$-norm for three different polynomial degrees $N=2,3$ and $6$. Table 4 lists $L^{1}$-errors and convergence order for each polynomial degree. This shows that our scheme converges to the expected order of $N+1$.

For the next test, we move to a simulation for which the analytical solution is not known. In this test, a stable non-rotating and non-magnetised neutron star is simulated in three space dimensions by solving the GRMHD equations in the Cowling approximation, i.e. assuming a static background spacetime. The initial state has been obtained as a solution to the Tolman-Oppenheimer-Volkoff (TOV) equations. The corresponding fluid and metric variables are compatible with the Einstein field equations. We set the magnetic field to zero for TOV stars.

The TOV equations constitute a non-linear ODE system in the radial space coordinate, that has been solved numerically by means of a fourth order Runge-Kutta scheme on a very fine grid with step size $dr=0.001$. The parameters of the problem have been chosen to be: $\rho_{c}=1.28\cdot 10^{-3}$, adiabatic exponent $\Gamma=2$ and a constant atmospheric pressure $p_{\text{atm}}=10^{-12}$.

The star sits at the origin of the domain $[-15,15]^{3}$. We apply reflection boundary conditions at the three simulation boundary surfaces $x=0$, $y=0$ and $z=0$. At the surfaces $x=R$, $y=R$ and $z=R$, we apply exact boundary conditions, evaluating the initial data at the boundary.

In Figure 13, we show the shared memory of scalability of the the non-linear hybrid ADER-DG - FV implementation of the GRMHD equations for the TOV star run on SuperMUC phase 2. We consider a regular base grid with $27^{3}$ cells and allow a dynamic refinement criterion to add one additional level of cells along a sphere around the TOV star.

6.5 Compressible Navier-Stokes

The compressible Navier-Stokes equations (3) demonstrate that ExaHyPE is not only able to solve hyperbolic PDEs, but is also capable of handling viscous effects. We use the numerical flux of [63] to integrate this into our ADER-DG framework.

To test the equations, we utilise the Arnold-Beltrami-Childress (ABC) flow

[TABLE]

where the constant $C=100/1.4$ governs the Mach number and $\mu=0.01$ is the viscosity [64]. This equation has an analytical solution in the incompressible limit. We impose this solution at the boundary. In Figure 14 we show a comparison between the analytical solution and a simulation with order $N=2$ on a regular mesh of $27^{3}$ cells. We see a good agreement with the analytical solution a time $t=1.0$.

In addition, we evaluate our scenario for the colliding bubbles scenario of [65]. The initial conditions of this scenario are obtained by an atmosphere with a background state that is in hydrostatic balance, i.e. where

[TABLE]

Initially, the domain has the same potential temperature $\Theta$. This is then perturbed by

[TABLE]

where $r$ is the Euclidean distance from the centre of the bubble $(x_{c},z_{c})$ and the spatial position $(x,z)$. We have two bubbles, with constants

[TABLE]

The simulation runs for $600\text{\,}\mathrm{s}$ and uses a viscosity of $\mu=0.01$ for regularisation. We use order $N=5$ and a mesh with up to two adaptive AMR levels. The results of this simulation are shown in Figure 15 (left). There is an excellent agreement with the reference solution of [65]. For further details on the setup, we refer to [35]. In Figure 15 (right), we show the shared-memory scalability using the ABC-flow on SuperMUC phase 2. In this test case the shared memory scalability is very good for all tested polynomial orders.

6.6 Hybrid Parallelisation Strategy

In the previous sections we showed the shared memory scalability of various test cases in two and three dimensions. In general we observed that the scalability improves as the work per element increases, i.e. it improves with higher polynomial degree $N$, for increasing dimension, or PDE complexity. In this section we demonstrate the effectiveness of a hybrid parallelisation strategy using the GRMHD TOV star setup described in Section 6.4. We consider a regular base grid of $79^{3}$ cells and allow up to two levels of adaptive mesh refinement. As before all kernel level optimisations have been switched on.

This scaling test was run on up to $56$ nodes of SuperMUC Phase 2 consisting of two $2.2$ GHz $14$ core Intel Xeon E5-v3 processors each. The tests were run using both Intel’s TBB and MPI for parallelisation. As a reference value we use a test run on $56$ cores and calculate the corresponding parallel efficiency compared to this baseline. Table 5 contains the results of this strong scaling test. Here, efficiency measures the speedup obtained divided by the expected optimal speedup. The efficiency is over $1.0$ in some cases can be explained by the chosen baseline of $56$ cores.

7 Conclusions and Future Work

This paper introduces a software engine that allows users to write higher-order ADER-DG codes for hyperbolic PDE systems with both conservative and non-conservative terms and a-posteriori limiters. The engine, ExaHyPE, has been designed to work on a wide range of computer systems from laptops to large high-performance compute clusters. Its core vision is to provide an engine rather than a framework: Technical details both on the computational science and numerical side are hidden from the user. A specification file plus very few routines realising the PDE terms are typically the only customisation points modified by user codes. Users can focus on the physics. To achieve this high level of abstraction, writing codes in ExaHyPE requires the use of a pre-specified set of numerical methods. The code thus clearly stands in the tradition of software packages such as Clawpack [66], as opposed to more generic, general-purpose software packages such as AMReX, deal.II or DUNE.

There are three general directions for future work. First, the sketched application areas have to make an impact in their respective domain. This comprises application-specific performance engineering and studies of real-world setups and data. It notably also comprises the coupling of engine applications with other code building blocks. Examples for first work into this direction are [33, 34, 2]. Second, we have to continue to investigate and to invest into the methods under the hood of ExaHyPE. Examples for such new ingredients on our agenda are accelerator support, local time stepping, and more dynamic load balancing and autotuning [67, 68]. Finally, we plan extensions of the core paradigm of the engine. The development of ExaHyPE from scratch has been possible as we restricted ourself to ADER-DG only. In the future, we plan to elaborate to which degree we are able to add particle-based (Particle-in-Cell methods or simple tracers) algorithms or multigrid solvers (for elliptic subterms) on top of ExaHyPE. This will open new user communities to the engine. First feasibility studies for this [69, 70, 71] already exist.

8 Acknowledgements and Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 671698, www.exahype.eu.

The ExaHyPE team acknowledges additional support by the European Union’s Horizon 2020 research and innovation program (ChEESE, grant no. 823844), A.-A.G. specifically thanks the ERC (starting grant number: 852992, TEAR), the German Research Foundation (DFG) (projects no. KA 2281/4-1, GA 2465/2-1, GA 2465/3-1), KONWIHR – the Bavarian Competence Network for Technical and Scientific High Performance Computing (project NewWave), BaCaTec (project no. A4), and KAUST-CRG (GAST, grant no. ORS-2016- CRG5-3027 and FRAGEN, grant no. ORS-2017-CRG6 3389.02) for additional support.

The authors gratefully acknowledge the support by the Leibniz Supercomputing Centre (www.lrz.de), which also provided the computing resources on SuperMUC (grant no. pr48ma and grant no. pr63qo) .

We would especially like to thank the many people who have made contributions to ExaHyPE, in particular our previous team members Benjamin Hazelwood, Angelika Schwarz, Vasco Varduhn and Olindo Zanotti.