Search

ExaGRyPE describes a suite of solvers for numerical relativity, based upon ExaHyPE 2, the second generation of our Exascale Hyperbolic PDE Engine. The presented generation of ExaGRyPE solves the Einstein equations in the CCZ4 formulation under a 3+1 foliation and focuses on black hole spacetimes. It employs a block-structured Cartesian grid carrying a higher-order order Finite Difference scheme with adaptive mesh refinement, it facilitates massive parallelism combining message passing, domain decomposition and task parallelism, and it supports the injection of particles into the grid as data probes or tracers. We introduce the ExaGRyPE-specific building blocks within ExaHyPE 2, and discuss its software architecture and compute-n-feel: For this, we formalize the creation of any specific simulation with ExaGRyPE as a sequence of lowering operations, where abstract logical tasks are successively broken into finer tasks until we obtain an abstraction level that can be mapped onto a C++ executable. The overall program logic is fully specified via a domain-specific Python interface, we map this logic onto a more detailed set of numerical tasks, subsequently lower this representation onto technical tasks that the underlying ExaHyPE engine uses to parallelize the application, before eventually the technical tasks are mapped onto task graphs including the actual PDE term evaluations, initial conditions, boundary conditions, and so forth. We end up with a rigorous separation of concerns which shields ExaGRyPE users from technical details and hence simplifies the development of novel physical models. We present the simulations and data for the gauge wave, static single black holes and rotating binary black hole systems, demonstrating that the code base is mature and usable. However, we also uncover domain-specific numerical challenges that need further study by the community in future work.

The C++ programming language and its cousins lean towards a memory-inefficient storage of structs: The compiler inserts helper bits into the struct such that individual attributes align with bytes, and it adds additional bytes aligning attributes with cache lines, while it is not able to exploit knowledge about the range of integers, enums or bitsets to bring the memory footprint down. Furthermore, the language provides neither support for data exchange via MPI nor for arbitrary floating-point precision formats. If developers need to have a low memory footprint and MPI datatypes over structs which exchange only minimal data, they have to manipulate the data and to write MPI datatypes manually. We propose a C++ language extension based upon C++ attributes through which developers can guide the compiler what memory arrangements would be beneficial: Can multiple booleans be squeezed into one bit field, do floats hold fewer significant bits than in the IEEE standard, or does the code require a user-defined MPI datatype for certain subsets of attributes? The extension offers the opportunity to fall back to normal alignment and padding rules via plain C++ assignments, no dependencies upon external libraries are introduced, and the resulting code remains standard C++. Our work implements the language annotations within LLVM and demonstrates their potential impact, both upon the runtime and the memory footprint, through smoothed particle hydrodynamics (SPH) benchmarks. They uncover the potential gains in terms of performance and development productivity.

The OpenMP API offers both task-based and data-parallel concepts to scientific computing. While it provides descriptive and prescriptive annotations, it is in many places deliberately unspecific how to implement its annotations. As the predominant OpenMP implementations share design rationales, they introduce "quasi-standards how certain annotations behave. By means of a task-based astrophysical simulation code, we highlight situations where this "quasi-standard" reference behaviour introduces performance flaws. Therefore, we propose prescriptive clauses to constrain the OpenMP implementations. Simulated task traces uncover the clauses' potential, while a discussion of their realization highlights that they would manifest in rather incremental changes to any OpenMP runtime supporting task priorities.

The C programming language and its cousins such as C++ stipulate the static storage of sets of structured data: Developers have to commit to one, invariant data model -- typically a structure-of-arrays (SoA) or an array-of-structs (AoS) -- unless they manually rearrange, i.e.~convert it throughout the computation. Whether AoS or SoA is favourable depends on the execution context and algorithm step. We propose a language extension based upon C++ attributes through which developers can guide the compiler what memory arrangements are to be used. The compiler can then automatically convert (parts of) the data into the format of choice prior to a calculation and convert results back afterwards. As all conversions are merely annotations, it is straightforward for the developer to experiment with different storage formats and to pick subsets of data that are subject to memory rearrangements. Our work implements the annotations within Clang and demonstrates their potential impact through a smoothed particle hydrodynamics (SPH) code.

Discrete Element Methods (DEM), i.e.~the simulation of many rigid particles, suffer from very stiff differential equations plus multiscale challenges in space and time. The particles move smoothly through space until they interact almost instantaneously due to collisions. Dense particle packings hence require tiny time step sizes, while free particles can advance with large time steps. Admissible time step sizes can span multiple orders of magnitudes. We propose an adaptive local time stepping algorithm which identifies clusters of particles that can be updated independently, advances them optimistically and independently in time, determines collision time stamps in space-time such that we maximise the time step sizes used, and resolves the momentum exchange implicitly. It is combined with various acceleration techniques which exploit multiscale geometry representations and multiscale behaviour in time. The collision time stamp detection in space-time in combination with the implicit solve of the actual collision equations avoids that particles get locked into tiny time step sizes, the clustering yields a high concurrency level, and the acceleration techniques plus local time stepping avoid unnecessary computations. This brings a scaling, adaptive time stepping for DEM for real-world challenges into reach.

We discuss three SYCL realisations of a simple Finite Volume scheme over multiple Cartesian patches. The realisation flavours differ in the way how they map the compute steps onto loops and tasks: We compare an implementation that is exclusively using a sequence of for-loops to a version that uses nested parallelism, and finally benchmark these against a version modelling the calculations as task graph. Our work proposes realisation idioms to realise these flavours within SYCL. The results suggest that a mixture of classic task and data parallelism performs if we map this hybrid onto a solely data-parallel SYCL implementation, taking into account SYCL specifics and the problem size.

The ISO C++17 standard introduces \emph{parallel algorithms}, a parallel programming model promising portability across a wide variety of parallel hardware including multi-core CPUs, GPUs, and FPGAs. Since 2019, the NVIDIA HPC SDK compiler suite supports this programming model for multi-core CPUs and GPUs. ExaHyPE is a solver engine for hyperbolic partial differential equations for complex wave phenomena. It supports multiple numerical methods including Finite Volumes and ADER-DG, and employs adaptive mesh refinement with dynamic load balancing via space-filling curves as well as task-based parallelism and offloading to GPUs. This study ports ExaHyPE's tasks over blocks of Finite Volumes to the ISO C++ parallel algorithms programming model, and compares its performance and usability against an OpenMP implementation with offloading via OpenMP target directives. It shows that ISO C++ is a feasible programming model for non-trivial applications like our task-based AMR code. The realisation is bare of vendor-specific or non-C++ extensions. It however is slower than its OpenMP counterpart. \vspace{-1cm}

The spherical collapse scenario has great importance in cosmology since it captures several crucial aspects of structure formation. The presence of self-similar solutions in the Einstein-de Sitter (EdS) model greatly simplifies its analysis, making it a powerful tool to gain valuable insights into the real and more complicated physical processes involved in galaxy formation. While there has been a large body of research to incorporate various additional physical processes into spherical collapse, the effect of modified gravity (MG) models, which are popular alternatives to the $\Lambda CDM$ paradigm to explain the cosmic acceleration, is still not well understood in this scenario. In this paper, we study the spherical accretion of collisional gas in a particular MG model, which is a rare case that also admits self-similar solutions. The model displays interesting behaviours caused by the enhanced gravity and a screening mechanism. Despite the strong effects of MG, we find that its self-similar solution agrees well with that of the EdS model. These results are used to assess a new cosmological hydrodynamical code for spherical collapse simulations introduced here, which is based on the hyperbolic partial differential equation engine ExaHyPE 2. Its good agreement with the theoretical predictions confirms the reliability of this code in modelling astrophysical processes in spherical collapse. We will use this code to study the evolution of gas in more realistic MG models in future work., Comment: 20 pages, 13 figures

Resilient algorithms in high-performance computing are subject to rigorous non-functional constraints. Resiliency must not increase the runtime, memory footprint or I/O demands too significantly. We propose a task-based soft error detection scheme that relies on error criteria per task outcome. They formalise how ``dubious'' an outcome is, i.e. how likely it contains an error. Our whole simulation is replicated once, forming two teams of MPI ranks that share their task results. Thus, ideally each team handles only around half of the workload. If a task yields large error criteria values, i.e.~is dubious, we compute the task redundantly and compare the outcomes. Whenever they disagree, the task result with a lower error likeliness is accepted. We obtain a self-healing, resilient algorithm which can compensate silent floating-point errors without a significant performance, I/O or memory footprint penalty. Case studies however suggest that a careful, domain-specific tailoring of the error criteria remains essential.

Space-filling curves (SFCs) are used in high performance computing to distribute a computational domain or its mesh, respectively, amongst different compute units, i.e.~cores or nodes or accelerators. The part of the domain allocated to each compute unit is called a partition. Besides the balancing of the work, the communication cost to exchange data between units determines the quality of a chosen partition. This cost can be approximated by the surface-to-volume ratio of partitions: the volume represents the amount of local work, while the surface represents the amount of data to be transmitted. Empirical evidence suggests that space-filling curves yield advantageous surface-to-volume ratios. Formal proofs are available only for regular grids. We investigate the surface-to-volume ratio of space-filling curve partitions for adaptive grids and derive the maximum surface-to-volume ratio as a function of the number of cells in the partition. In order to prove our main theorem, we construct a new framework for the study of adaptive grids, notably introducing the concepts of a shape and of classified partitions. The new methodological framework yields insight about the SFC-induced partition character even if the grids refine rather aggressively in localised areas: it quantifies the obtained surface-to-volume ratio. This framework thus has the potential to guide the design of better load balancing algorithms on the long term.

Simulations of many rigid bodies colliding with each other sometimes yield particularly interesting results if the colliding objects differ significantly in size and are non-spherical. The most expensive part within such a simulation code is the collision detection. We propose a family of novel multiscale collision detection algorithms that can be applied to triangulated objects within explicit and implicit time stepping methods. They are well-suited to handle objects that cannot be represented by analytical shapes or assemblies of analytical objects. Inspired by multigrid methods and adaptive mesh refinement, we determine collision points iteratively over a resolution hierarchy, and combine a functional minimisation plus penalty parameters with the actual comparision-based geometric distance calculation. Coarse surrogate geometry representations identify "no collision" scenarios early on and otherwise yield an educated guess which triangle subsets of the next finer level potentially yield collisions. They prune the search tree, and furthermore feed conservative contact force estimates into the iterative solve behind an implicit time stepping. Implicit time stepping and non-analytical shapes often yield prohibitive high compute cost for rigid body simulations. Our approach reduces these cost algorithmically by one to two orders of magnitude. It also exhibits high vectorisation efficiency due to its iterative nature.

The orchestration of complex algorithms demands high levels of automation to use modern hardware efficiently. Task-based programming with OpenMP 5.0 is a prominent candidate to accomplish this goal. We study OpenMP 5.0's tasking in the context of a wave equation solver (ExaHyPE) using three different architectures and runtimes. We describe several task-scheduling flaws present in currently available runtimes, demonstrate how they impact performance and show how to work around them. Finally, we propose extensions to the OpenMP standard.

In an era where we can not afford to checkpoint frequently, replication is a generic way forward to construct numerical simulations that can continue to run even if hardware parts fail. Yet, replication often is not employed on larger scales, as na\"ively mirroring a computation once effectively halves the machine size, and as keeping replicated simulations consistent with each other is not trivial. We demonstrate for the ExaHyPE engine -- a task-based solver for hyperbolic equation systems -- that it is possible to realise resiliency without major code changes on the user side, while we introduce a novel algorithmic idea where replication reduces the time-to-solution. The redundant CPU cycles are not burned "for nothing". Our work employs a weakly consistent data model where replicas run independently yet inform each other through heartbeat messages whether they are still up and running. Our key performance idea is to let the tasks of the replicated simulations share some of their outcomes, while we shuffle the actual task execution order per replica. This way, replicated ranks can skip some local computations and automatically start to synchronise with each other. Our experiments with a production-level seismic wave-equation solver provide evidence that this novel concept has the potential to make replication affordable for large-scale simulations in high-performance computing.

The accurate assembly of the system matrix is an important step in any code that solves partial differential equations on a mesh. We either explicitly set up a matrix, or we work in a matrix-free environment where we have to be able to quickly return matrix entries upon demand. Either way, the construction can become costly due to non-trivial material parameters entering the equations, multigrid codes requiring cascades of matrices that depend upon each other, or dynamic adaptive mesh refinement that necessitates the recomputation of matrix entries or the whole equation system throughout the solve. We propose that these constructions can be performed concurrently with the multigrid cycles. Initial geometric matrices and low accuracy integrations kickstart the multigrid, while improved assembly data is fed to the solver as and when it becomes available. The time to solution is improved as we eliminate an expensive preparation phase traditionally delaying the actual computation. We eliminate algorithmic latency. Furthermore, we desynchronise the assembly from the solution process. This anarchic increase of the concurrency level improves the scalability. Assembly routines are notoriously memory- and bandwidth-demanding. As we work with iteratively improving operator accuracies, we finally propose the use of a hierarchical, lossy compression scheme such that the memory footprint is brought down aggressively where the system matrix entries carry little information or are not yet available with high accuracy.

Today's hardware's explosion of concurrency plus the explosion of data we build upon in both machine learning and scientific simulations have multifaceted impact on how we write our codes. They have changed our notion of performance and, hence, of what a good code is: Good code has, first of all, to be able to exploit the unprecedented levels of parallelism. To do so, it has to manage to move the compute data into the compute facilities on time. As communication and memory bandwidth cannot keep pace with the growth in compute capabilities and as latency increases---at least relative to what the hardware could do---communication-avoiding techniques gain importance. We characterise and classify the field of communication-avoiding algorithms. A review of some examples of communication-avoiding programming by means of our new terminology shows that we are well-advised to broaden our notion of "communication-avoiding" and to look beyond numerical linear algebra. An abstraction, generalisation and weakening of the term enriches our toolset of how to tackle the data movement challenges. Through this, we eventually gain access to a richer set of tools that we can use to deliver proper code for current and upcoming hardware generations.

Balancing the workload of sophisticated simulations is inherently difficult, since we have to balance both computational workload and memory footprint over meshes that can change any time or yield unpredictable cost per mesh entity, while modern supercomputers and their interconnects start to exhibit fluctuating performance. We propose a novel lightweight balancing technique for MPI+X to accompany traditional, prediction-based load balancing. It is a reactive diffusion approach that uses online measurements of MPI idle time to migrate tasks temporarily from overloaded to underemployed ranks. Tasks are deployed to ranks which otherwise would wait, processed with high priority, and made available to the overloaded ranks again. This migration is non-persistent. Our approach hijacks idle time to do meaningful work and is totally non-blocking, asynchronous and distributed without a global data view. Tests with a seismic simulation code developed in the ExaHyPE engine uncover the method's potential. We found speed-ups of up to 2-3 for ill-balanced scenarios without logical modifications of the code base and show that the strategy is capable to react quickly to temporarily changing workload or node performance.

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.

Multigrid solvers face multiple challenges on parallel computers. Two fundamental ones read as follows: Multiplicative solvers issue coarse grid solves which exhibit low concurrency and many multigrid implementations suffer from an expensive coarse grid identification phase plus adaptive mesh refinement overhead. We propose a new additive multigrid variant for spacetrees, i.e. meshes as they are constructed from octrees and quadtrees: It is an additive scheme, i.e. all multigrid resolution levels are updated concurrently. This ensures a high concurrency level, while the transfer operators between the mesh levels can still be constructed algebraically. The novel flavour of the additive scheme is an augmentation of the solver with an additive, auxiliary damping parameter per grid level per vertex that is in turn constructed through the next coarser level---an idea which utilises smoothed aggregation principles or the motivation behind AFACx: Per level, we solve an additional equation whose purpose is to damp too aggressive solution updates per vertex which would otherwise, in combination with all the other levels, yield an overcorrection and, eventually, oscillations. This additional equation is constructed additively as well, i.e. is once more solved concurrently to all other equations. This yields improved stability, closer to what is seen with multiplicative schemes, while pipelining techniques help us to write down the additive solver with single-touch semantics for dynamically adaptive meshes.

Simple stencil codes are and remain an important building block in scientific computing. On shared memory nodes, they are traditionally parallelised through colouring or (recursive) tiling. New OpenMP versions alternatively allow users to specify data dependencies explicitly and to outsource the decision how to distribute the work to the runtime system. We evaluate traditional multithreading strategies on both Broadwell and KNL, study the arising assignment of tasks to threads and, from there, derive two efficient ways to parallelise stencil codes on regular Cartesian grids that fuse colouring and task-based approaches., Comment: 10 pages

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