No abstract

We introduce a family of implementations of low order, additive, geometric multilevel solvers for systems of Helmholtz equations arising from Schrödinger equations. Both grid spacing and arithmetics may comprise complex numbers and we thus can apply complex scaling to the indefinite Helmholtz operator. Our implementations are based upon the notion of a spacetree and work exclusively with a finite number of precomputed local element matrices. They are globally matrix-free.Combining various relaxation factors with two grid transfer operators allows us to switch from additive multigrid over a hierarchical basis method into a Bramble-Pasciak-Xu (BPX)-type solver, with several multiscale smoothing variants within one code base. Pipelining allows us to realise full approximation storage (FAS) within the additive environment where, amortised, each grid vertex carrying degrees of freedom is read/written only once per iteration. The codes realise a single-touch policy. Among the features facilitated by matrix-free FAS is arbitrary dynamic mesh refinement (AMR) for all solver variants. AMR as enabler for full multigrid (FMG) cycling-the grid unfolds throughout the computation-allows us to reduce the cost per unknown per order of accuracy.The present paper primary contributes towards software realisation and design questions. Our experiments show that the consolidation of single-touch FAS, dynamic AMR and vectorisation-friendly, complex scaled, matrix-free FMG cycles delivers a mature implementation blueprint for solvers of Helmholtz equations in general. For this blueprint, we put particular emphasis on a strict implementation formalism as well as some implementation correctness proofs.

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

This paper presents the general purpose framework Peano for the solution of partial differential equations (PDE) on adaptive Cartesian grids. The strict structuredness and inherent multilevel property of these grids allows for very low memory requirements, efficient (in terms of hardware performance) implementations of parallel multigrid solvers on dynamically adaptive grids, and arbitrary spatial dimensions. This combination of advantages distinguishes Peano from other PDE frameworks. We describe shortly the underlying octree-like grid type and its most important properties. The main part of the paper shows the framework concept of Peano and the implementation of a Navier-Stokes solver as one of the main currently implemented application examples. Various results ranging from hardware and numerical performance to concrete application scenarios close the contribution.

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 present a family of spacetree-based multigrid realizations using the tree's multiscale nature to derive coarse grids. They align with matrix-free geometric multigrid solvers as they never assemble the system matrices which is cumbersome for dynamically adaptive grids and full multigrid. The most sophisticated realizations use BoxMG to construct operator-dependent prolongation and restriction in combination with Galerkin/Petrov-Galerkin coarse-grid operators. This yields robust solvers for nontrivial elliptic problems. We embed the algebraic, problem-and grid-dependent multigrid operators as stencils into the grid and evaluate all matrix-vector products in-situ throughout the grid traversals. While such an approach is not literally matrix-free-the grid carries the matrix-we propose to switch to a hierarchical representation of all operators. Only differences of algebraic operators to their geometric counterparts are held. These hierarchical differences can be stored and exchanged with small memory footprint. Our realizations support arbitrary dynamically adaptive grids while they vertically integrate the multilevel operations through spacetree linearization. This yields good memory access characteristics, while standard colouring of mesh entities with domain decomposition allows us to use parallel manycore clusters. All realization ingredients are detailed such that they can be used by other codes.

SUMMARYThis paper presents a new e cient way to implement multigrid algorithms on adaptively reÿned grids. To cope with todays demands in high-performance computing, we cannot do without such highly sophisticated numerical methods. But if we do not implement them very carefully, we lose a lot of e ciency in terms of memory usage: using trees for the storage of hierarchical multilevel data causes a large amount of non-local (in terms of the physical memory space) data accesses, and often requires the storage of pointers to neighbours to allow the evaluation of discrete operators (di erence stencils, restrictions, interpolations, etc.). The importance of this problem becomes clear if we remember that storage and not the CPUs is the bottleneck on modern computers.We established a cache-oblivious and storage-minimizing algorithm based on the concept of spacetree grids combined with a cell-oriented operator evaluation, a linear ordering of grid cells along a space-ÿlling curve, and a sophisticated construction of linearly processed data structures for vertex data. In this context, we could show that the implementation of a dynamically adaptive F-cycle is, ÿrst, very natural and, second, does not cause any overhead in terms of storage usage and access as adaptivity and multilevel data do not disturb the linear processing order of our data structures.

In this paper, we present a parallel multigrid PDE solver working on adaptive hierarchical cartesian grids. The presentation is restricted to the linear elliptic operator of second order, but extensions are possible and have already been realised as prototypes. Within the solver the handling of the vertices and the degrees of freedom associated to them is implemented solely using stacks and iterates of a Peano space-filling curve. Thus, due to the structuredness of the grid, two administrative bits per vertex are sufficient to store both geometry and grid refinement information. The implementation and parallel extension, using a spacefilling curve to obtain a load balanced domain decomposition, will be formalised. In view of the fact that we are using a multigrid solver of linear complexity O(n), it has to be ensured that communication cost and, hence, the parallel algorithm's overall complexity do not exceed this linear behaviour. This work has partially been funded by DFG's research unit FOR493 and the DFG project HA 1517/25-1/2.