Search

A variety of methods is available to quantify uncertainties arising within the modeling of flow and transport in carbon dioxide storage, but there is a lack of thorough comparisons. Usually, raw data from such storage sites can hardly be described by theoretical statistical distributions since only very limited data is available. Hence, exact information on distribution shapes for all uncertain parameters is very rare in realistic applications. We discuss and compare four different methods tested for data-driven uncertainty quantification based on a benchmark scenario of carbon dioxide storage. In the benchmark, for which we provide data and code, carbon dioxide is injected into a saline aquifer modeled by the nonlinear capillarity-free fractional flow formulation for two incompressible fluid phases, namely carbon dioxide and brine. To cover different aspects of uncertainty quantification, we incorporate various sources of uncertainty such as uncertainty of boundary conditions, of parameters in constitutive relations, and of material properties. We consider recent versions of the following non-intrusive and intrusive uncertainty quantification methods: arbitrary polynomial chaos, spatially adaptive sparse grids, kernel-based greedy interpolation, and hybrid stochastic Galerkin. The performance of each approach is demonstrated assessing expectation value and standard deviation of the carbon dioxide saturation against a reference statistic based on Monte Carlo sampling. We compare the convergence of all methods reporting on accuracy with respect to the number of model runs and resolution. Finally, we offer suggestions about the methods’ advantages and disadvantages that can guide the modeler for uncertainty quantification in carbon dioxide storage and beyond.

In the age of data collection, machine learning algorithms have to be able to efficiently cope with vast data sets. This requires scalable algorithms and efficient implementations that can cope with heterogeneous hardware. We propose a new, performance-portable implementation of a well-known, robust, and versatile multi-class classification method that supports multiple Graphics Processing Units (GPUs) from different vendors. It is based on a performance-portable implementation of the approximate k-nearest neighbors (k-NN) algorithm in SYCL. The k-NN assigns a class to a data point based on a majority vote of its neighborhood. The naive approach compares a data point x to all other data points in the training data to identify the k nearest ones. However, this has quadratic runtime and is infeasible for large data sets. Therefore, approximate variants have been developed. Such an algorithm is the Locality-Sensitive Hashing (LSH) algorithm, which uses hash tables together with locality-sensitive hash functions to reduce the data points that have to be examined to compute the k-NN. To the best of our knowledge, there is no distributed LSH version supporting multiple GPUs from different vendors available so far despite the fact that k-NNs are frequently employed. Therefore, we have developed the library. It provides the first hardware-independent, yet efficient and distributed implementation of the LSH algorithm that is suited for modern supercomputers. The implementation uses C++17 together with SYCL 1.2.1, which is an abstraction layer for OpenCL that allows targeting different hardware with a single implementation. To support large data sets, we utilize multiple GPUs using the Message Passing Interface (MPI) to enable the usage of both shared and distributed memory systems. We have tested different parameter combinations for two locality-sensitive hash function implementations, which we compare. Our results show that our library can easily scale on multiple GPUs using both hash function types, achieving a nearly optimal parallel speedup of up to 7.6 on 8 GPUs. Furthermore, we demonstrate that the library supports different SYCL implementations—ComputeCpp, hipSYCL, and DPC++—to target different hardware architectures without significant performance differences.

Robust prediction of the behavior of complex physical and engineering systems relies on approximating solutions in terms of physical and stochastic domains. For higher resolution and accuracy, simulation models must increase the number of deterministic and stochastic variables and therefore further increase the dimensionality of the problem. Sparse grids are an established technique to tackle higher-dimensional problems. Their efficient tensor product structure allows the creation of accurate surrogates from few model evaluations. Classical approaches use hat functions, resulting in non-differentiable surrogates, or global basis functions, resulting in potential instabilities. Therefore, we propose using modified not-a-knot B-splines to overcome both problems. Additionally, we use established spatially adaptive refinement criteria to reduce the number of model evaluations even further. We compare our technique to other data-driven uncertainty quantification methods in a real-world benchmark for probabilistic risk assessment for carbon dioxide storage in geological formations.

We present a novel approach for the visualization of fracture processes in peridynamics simulations. In peridynamics simulation, materials are represented by material points linked with bonds, providing complex fracture behavior. Our approach first extracts the cracks from each time step by means of height ridge extraction. To avoid deterioration of the structures, we propose an approach to extract ridges from these data without resampling. The extracted crack geometries are then combined into a spatiotemporal structure, with special focus on temporal coherence and robustness. We then show how this structure can be used for various visualization approaches to reveal fracture dynamics, with a focus on physical mechanisms. We evaluate our approach and demonstrate its utility by means of different data sets.

The present paper uses Langevin dynamics (LD) to investigate the aggregation of soot nano-particles in turbulent flows. Interparticle forces are included, and the computation of the individual particles by LD is retained even after aggregate formation such that collision events and locations can be based on center-to-center particle distances without invoking any modelling assumptions of aggregate shape and/or collision frequency. We focus on the interactions between the specific hydrodynamic conditions and the particle properties and their effect on the resulting agglomerates’ morphologies. The morphology is characterized by the fractal dimension, D f . Computations of particle aggregation in homogeneous isotropic turbulence and in shear flows dominated by counter-rotating vortices with a wide range of turbulence intensities and particle sizes indicate that the evolution of the agglomerates’ shapes can be adequately parameterized by the size of the agglomerates and the Knudsen and Peclet numbers, the latter being based on the smallest turbulence scales. The computations further suggest that the shapes of agglomerates of certain sizes are relatively independent of time and relatively insensitive to larger turbulence structures. The fractal dimensions are modelled as functions of radius of gyration, Kn and Pe. The fitted expressions show good agreement with the LD simulations and represent the entire growth process of the agglomerates. A direct comparison of selected aggregates with experimental data shows very good qualitative agreement. A thorough quantitative validation of the evolution of the computed aggregate characteristics is, however, presently hindered by the challenges for and therefore lack of suitable experiments under appropriately controlled conditions.

We present a highly scalable demonstration of a portable asynchronous many-task programming model and runtime system applied to a grid-based adaptive mesh refinement hydrodynamic simulation of a double white dwarf merger with 14 levels of refinement that spans 17 orders of magnitude in astrophysical densities. The code uses the portable C++ parallel programming model that is embodied in the HPX library and being incorporated into the ISO C++ standard. The model represents a significant shift from existing bulk synchronous parallel programming models under consideration for exascale systems. Through the use of the Futurization technique, seemingly sequential code is transformed into wait-free asynchronous tasks. We demonstrate the potential of our model by showing results from strong scaling runs on National Energy Research Scientific Computing Center’s Cori system (658,784 Intel Knight’s Landing cores) that achieve a parallel efficiency of 96.8% using billions of asynchronous tasks.

This paper shows a new approach to estimate the critical traction for Mode I crack opening before crack growth by numerical simulation. For quasi-static loading, Linear Elastic Fracture Mechanics predicts the critical traction before crack growth. To simulate the crack growth, we used bond-based peridynamics, a non-local generalization of continuum mechanics. We discretize the peridynamics equation of motion with a collocation by space approach, the so-called EMU nodal discretization. As the constitutive law, we employ the improved prototype micro brittle material model. This bond-based material model is verified by the Young’s modulus from classical theory for a homogeneous deformation for different quadrature rules. For the EMU-ND we studied the behavior for different ratios of the horizon and nodal spacing to gain a robust value for a large variety of materials. To access this wide range of materials, we applied sparse grids, a technique to build high-dimensional surrogate models. Sparse grids significantly reduce the number of simulation runs compared to a full grid approach and keep up a similar approximation accuracy. For the validation of the quasi-static loading process, we show that the critical traction is independent of the material density for most material parameters. The bond-based IPMB model with EMU nodal discretization seems very robust for the ratio $${\delta }/{\varDelta X}=3$$ for a wide range of materials, if an error of 5 % is acceptable.

We study the global communication of the numerical solution of high-dimensional PDEs.We design two optimal communication schemes for the sparse grid combination technique.We present a new communication model based on the system's latency and bandwidth.The communication model predicts the performance of the communication schemes.Experimental results on several current supercomputers confirm the predictions. The numerical treatment of high-dimensional partial differential equations is among the most compute-hungry problems and in urgent need for current and future high-performance computing (HPC) systems. It is thus also facing the grand challenges of exascale computing such as the requirement to reduce global communication. To cope with high dimensionalities we employ a hierarchical discretization scheme, the sparse grid combination technique. Based on an extrapolation scheme, the combination technique additionally mitigates the need for global communication: multiple and much smaller problems can be computed independently for each time step, and the global communication shrinks to a reduce/broadcast step in between. Here, we focus on this remaining synchronization step of the combination technique and present two communication schemes designed to either minimize the number of communication rounds or the total communication volume. Experiments on two different supercomputers show that either of the schemes outperforms the other depending on the size of the problem. Furthermore, we present a communication model based on the system's latency and bandwidth and validate the model with the experiments. The model can be used to predict the runtime of the reduce/broadcast step for dimensionalities that are yet out of scope on current supercomputers.

The optimization of performance of complex simulation codes with high computational demands, such as Octo-Tiger, is an ongoing challenge. Octo-Tiger is an astrophysics code simulating the evolution of star systems based on the fast multipole method using adaptive octrees as the central data structure. Octo-Tiger was implemented using high-level C++ libraries, specifically HPX and Vc, which allows its use on different hardware platforms. Recently, we have demonstrated excellent scalability in a distributed setting.In this paper, we study the node-level performance of Octo-Tiger on an Intel Knights Landing platform. We focus on Octo-Tiger's fast multipole method, as it is the computationally most demanding component. By using HPX and a futurization approach, we can efficiently traverse the adaptive octrees in parallel. On the core-level, threads process sub-grids using multiple 1074-element stencils.In numerical experiments, simulating the time evolution of a rotating star on an Intel Xeon Phi 7250 Knights Landing processor, Octo-Tiger shows good parallel efficiency. For the fast multipole algorithm, we achieved up to 408 GFLOPS, resulting in a speedup of 2x compared to a 24-core Skylake-SP platform, using the same high-level abstractions.

With future exascale computers expected to have millions of compute units distributed among thousands of nodes, system faults are predicted to become more frequent. Fault tolerance will thus play a key role in HPC at this scale. In this paper we focus on solving the 5-dimensional gyrokinetic Vlasov-Maxwell equations using the application code GENE as it represents a high-dimensional and resource-intensive problem which is a natural candidate for exascale computing. We discuss the Fault-Tolerant Combination Technique, a resilient version of the Combination Technique, a method to increase the discretization resolution of existing PDE solvers. For the first time, we present an efficient, scalable and fault-tolerant implementation of this algorithm for plasma physics simulations based on a manager-worker model and test it under very realistic and pessimistic environments with simulated faults. We show that the Fault-Tolerant Combination Technique - an algorithm-based forward recovery method - can tolerate a large number of faults with a low overhead and at an acceptable loss in accuracy. Our parallel experiments with up to 32k cores show good scalability at a relative parallel efficiency of 93.61%. We conclude that algorithm-based solutions to fault tolerance are attractive for this type of problems.

Comparing different accelerated cluster architectures by a single application is a tough piece of work because this application has to be optimized with respect to platform-dependent features. In this work, we demonstrate such an optimization for a data mining algorithm which solves regression and classification problems on vast data sets. Our technique is based on least squares regression, and its major component is the iterative matrix-free solution of a linear system of equations. By processing data sets ranging from several hundreds of thousands instances to multi-million data points in strong-scaling and weak-scaling settings, we are able to estimate the amount of parallelism needed to unleash the performance of classic CPU-based machines and clusters employing Intel Xeon Phi coprocessors and NVIDIA Kepler GPUs. Only in strong-scaling experiments, GPUs and coprocessors suffer from their tremendous amount of needed parallelism and get outperformed by dual socket Intel Sandy Bridge nodes at large scale more than 64 nodes/accelerators. However, in weak-scaling scenarios, a speed-up larger than 2X over an entire CPU node can be achieved by a single accelerator. Copyright © 2015 John Wiley & Sons, Ltd.

The analysis of simulation results and the verification against experimental data is essential to develop and interpret simulation models for impact damage. We present two visualization techniques to post-process particle-based simulation data, and we highlight new aspects for the quantitative comparison with experimental data. As the underlying simulation model we consider the particle method Peridynamics, a non-local generalization of continuum mechanics. The first analysis technique is an extended component labeling algorithm to extract the fragment size and the corresponding histograms. The distribution of the fragment size can be obtained by real-world experiments as demonstrated in Schram and Meyer (Simulating the formation and evolution of behind armor debris fields. ARL-RP 109, U.S. Army Research Laboratory, 2005), Vogler et al. (Int J Impact Eng 29:735–746, 2003). The second approach focuses on the visualization of the stress after an impact. Here, the particle-based data is re-sampled and rendered with standard volume rendering techniques to address the interference pattern of the stress wave after reflection at the boundary. For the extraction and visual analysis, we used the widely-used Stanford bunny as a complex geometry. For a quantitative study with a simple geometry, the edge-on impact experiment (Schradin, Scripts German Acad Aeronaut Res 40:21–68, 1939; Strassburger, Int J Appl Ceram Technol 1:1:235–242, 2004; Kawai et al., Procedia Eng 103:287–293, 2015) can be applied. With these new visualization approaches, new insights for the quantitative comparison of fragmentation and wave propagation become intuitively accessible.

The solution of higher-dimensional problems, such as the simulation of plasma turbulence in a fusion device as described by the five-dimensional gyrokinetic equations, is a grand challenge for current and future high-performance computing. The sparse grid combination technique is a promising approach to the solution of these problems on large-scale distributed memory systems. The combination technique numerically decomposes a single large problem into multiple moderately-sized partial problems that can be computed in parallel, independently and asynchronously of each other. The ability to efficiently combine the individual partial solutions to a common sparse grid solution is a key to the overall performance of such large-scale computations. In this work, we present new algorithms for the recombination of distributed component grids and demonstrate their scalability to 180, 225 cores on the supercomputer Hazel Hen.

We present a novel data-driven approach to propagate uncertainty. It consists of a highly efficient integrated adaptive sparse grid approach. We remove the gap between the subjective assumptions of the input’s uncertainty and the unknown real distribution by applying sparse grid density estimation on given measurements. We link the estimation to the adaptive sparse grid collocation method for the propagation of uncertainty. This integrated approach gives us two main advantages: First, the linkage of the density estimation and the stochastic collocation method is straightforward as they use the same fundamental principles. Second, we can efficiently estimate moments for the quantity of interest without any additional approximation errors. This includes the challenging task of solving higher-dimensional integrals. We applied this new approach to a complex subsurface flow problem and showed that it can compete with state-of-the-art methods. Our sparse grid approach excels by efficiency, accuracy and flexibility and thus can be applied in many fields from financial to environmental sciences.

As data has become easier to collect and precise sensors have become ubiquitous, data mining with large data sets has become an important problem. Because sparse grid data mining scales only linearly in the number of data points, large data mining problems have been successfully addressed with this method. Still, highly efficient algorithms are required to process very large problems within a reasonable amount of time.

Substantial modifications of both the choice of the grids, the combination coefficients, the parallel data structures and the algorithms used for the combination technique lead to numerical methods which are scalable. This is demonstrated by the provision of error and complexity bounds and in performance studies based on a state of the art code for the solution of the gyrokinetic equations of plasma physics. The key ideas for a new fault-tolerant combination technique are mentioned. New algorithms for both initial- and eigenvalue problems have been developed and are shown to have good performance.

Investigating the interplay between muscular activity and motion is the basis to improve our understanding of healthy or diseased musculoskeletal systems. To be able to analyze the musculoskeletal systems, computational models are used. Albeit some severe modeling assumptions, almost all existing musculoskeletal system simulations appeal to multibody simulation frameworks. Although continuum-mechanical musculoskeletal system models can compensate for some of these limitations, they are essentially not considered because of their computational complexity and cost. The proposed framework is the first activation-driven musculoskeletal system model, in which the exerted skeletal muscle forces are computed using 3-dimensional, continuum-mechanical skeletal muscle models and in which muscle activations are determined based on a constraint optimization problem. Numerical feasibility is achieved by computing sparse grid surrogates with hierarchical B-splines, and adaptive sparse grid refinement further reduces the computational effort. The choice of B-splines allows the use of all existing gradient-based optimization techniques without further numerical approximation. This paper demonstrates that the resulting surrogates have low relative errors (less than 0.76%) and can be used within forward simulations that are subject to constraint optimization. To demonstrate this, we set up several different test scenarios in which an upper limb model consisting of the elbow joint, the biceps and triceps brachii, and an external load is subjected to different optimization criteria. Even though this novel method has only been demonstrated for a 2-muscle system, it can easily be extended to musculoskeletal systems with 3 or more muscles.