In most classical approaches of computational geophysics for seismic wave propagation problems, complex surface topography is either accounted for by boundary-fitted unstructured meshes, or, where possible, by mapping the complex computational domain from physical space to a topologically simple domain in a reference coordinate system. However, all these conventional approaches face problems if the geometry of the problem becomes sufficiently complex. They either need a mesh generator to create unstructured boundary-fitted grids, which can become quite difficult and may require a lot of manual user interactions in order to obtain a high quality mesh, or they need the explicit computation of an appropriate mapping function from physical to reference coordinates. For sufficiently complex geometries such mappings may either not exist or their Jacobian could become close to singular. Furthermore, in both conventional approaches low quality grids will always lead to very small time steps due to the Courant-Friedrichs-Lewy (CFL) condition for explicit schemes. In this paper, we propose a completely different strategy that follows the ideas of the successful family of high resolution shock-capturing schemes, where discontinuities can actually be resolved anywhere on the grid, without having to fit them exactly. We address the problem of geometrically complex free surface boundary conditions for seismic wave propagation problems with a novel diffuse interface method (DIM) on adaptive Cartesian meshes (AMR) that consists in the introduction of a characteristic function 0 ≤ α ≤ 1 which identifies the location of the solid medium and the surrounding air (or vacuum) and thus implicitly defines the location of the free surface boundary. Physically, α represents the volume fraction of the solid medium present in a control volume. Our new approach completely avoids the problem of mesh generation, since all that is needed for the definition of the complex surface topography is to set a scalar color function to unity inside the regions covered by the solid and to zero outside. The governing equations are derived from ideas typically used in the mathematical description of compressible multiphase flows. An analysis of the eigenvalues of the PDE system shows that the complexity of the geometry has no influence on the admissible time step size due to the CFL condition. The model reduces to the classical linear elasticity equations inside the solid medium where the gradients of α are zero, while in the diffuse interface zone at the free surface boundary the governing PDE system becomes nonlinear. We can prove that the solution of the Riemann problem with arbitrary data and a jump in α from unity to zero yields a Godunov-state at the interface that satisfies the free-surface boundary condition exactly, i.e. the normal stress components vanish. In the general case of an interface that is not aligned with the grid and which is not infinitely thin, a systematic study on the distribution of the volume fraction function inside the interfac...

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.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.

We study the performance behaviour of a seismic simulation using the ExaHyPE engine with a specific focus on memory characteristics and energy needs. ExaHyPE combines dynamically adaptive mesh refinement (AMR) with ADER-DG. It is parallelized using tasks, and it is cache efficient. AMR plus ADER-DG yields a task graph which is highly dynamic in nature and comprises both arithmetically expensive tasks and tasks which challenge the memory's latency. The expensive tasks and thus the whole code benefit from AVX vectorization, though we suffer from memory access bursts. A frequency reduction of the chip improves the code's energy-to-solution. Yet, it does not mitigate burst effects. The bursts' latency penalty becomes worse once we add Intel R Optane TM technology, increase the core count significantly, or make individual, computationally heavy tasks fall out of close caches. Thread overbooking to hide away these latency penalties becomes contra-productive with non-inclusive caches as it destroys the cache and vectorization character. In cases where memory-intense and computationally expensive tasks overlap, ExaHyPE's cache-oblivious implementation nevertheless can exploit deep, non-inclusive, heterogeneous memory effectively, as main memory misses arise infrequently and slow down only few cores. We thus propose that upcoming supercomputing simulation codes with dynamic, inhomogeneous task graphs are actively supported by thread runtimes in intermixing tasks of different compute character, and we propose that future hardware actively allows codes to downclock the cores running particular task types.

Tobias (2020) 'Enclave tasking for discontinuous Galerkin methods on dynamically adaptive meshes.', SIAM Journal on Scientic Computing (SISC)., 42 (3). C69-C96.

With the advent of manycore systems, shared memory parallelisation has gained importance in high performance computing. Once a code is decomposed into tasks or parallel regions, it becomes crucial to identify reasonable grain sizes, i.e. minimum problem sizes per task that make the algorithm expose a high concurrency at low overhead. Many papers do not detail what reasonable task sizes are, and consider their findings craftsmanship not worth discussion. We have implemented an autotuning algorithm, a machine learning approach, for a project developing a hyperbolic equation system solver. Autotuning here is important as the grid and task workload are multifaceted and change frequently during runtime. In this paper, we summarise our lessons learned. We infer tweaks and idioms for general autotuning algorithms and we clarify that such a approach does not free users completely from grain size awareness.