Stencil Computations Research Articles

High-performance xPU Stencil Computations in Julia

High-performance xPU Stencil Computations in Julia

SOMA: Observability, monitoring, and in situ analytics for exascale applications

SummaryWith the rise of exascale systems and large, data‐centric workflows, the need to observe and analyze high performance computing (HPC) applications during their execution is becoming increasingly important. HPC applications are typically not designed with online monitoring in mind, therefore, the observability challenge lies in being able to access and analyze interesting events with low overhead while seamlessly integrating such capabilities into existing and new applications. We explore how our service‐based observation, monitoring, and analytics (SOMA) approach to collecting and aggregating both application‐specific diagnostic data and performance data addresses these needs. We present our SOMA framework and demonstrate its viability with LULESH, a hydrodynamics proxy application. Then we focus on Astaroth, a multi‐GPU library for stencil computations, highlighting the integration of the TAU and APEX performance tools and SOMA for application and performance data monitoring.

CPU–GPU heterogeneous code acceleration of a finite volume Computational Fluid Dynamics solver

This research focuses on accelerating the finite-volume Computational Fluid Dynamics (CFD) solver, SENSEI, through concurrent CPU–GPU heterogeneous computing, leveraging multiple CPUs and GPUs. An overview of SENSEI, its discretization, and the heterogeneous computing workflow utilizing MPI and OpenACC are provided. A performance model for CPU–GPU heterogeneous computing, incorporating ghost cell exchange usually applied in stencil computations, is introduced to estimate the performance. This work explores the scaling performance of CPU–GPU heterogeneous computing in comparison to a pure multi-CPU/GPU setup using a supersonic inlet test case. The results emphasize the advantages of leveraging both CPU and GPU computational power concurrently. Highlighting the collaborative synergy of CPUs and GPUs as co-workers, this work shows improved performance over exclusive use of pure CPUs or GPUs, with the proposed performance model ensuring a fair estimation of these advantages. This work concludes by offering suggestions for application users interested in accelerating scientific computing simulations and providing feedback for hardware architects designing improved CPU–GPU heterogeneous systems for heterogeneous computing.

A Fast Algorithm for Aperiodic Linear Stencil Computation using Fast Fourier Transforms

Stencil computations are widely used to simulate the change of state of physical systems across a multidimensional grid over multiple timesteps. The state-of-the-art techniques in this area fall into three groups: cache-aware tiled looping algorithms, cache-oblivious divide-and-conquer trapezoidal algorithms, and Krylov subspace methods. In this article, we present two efficient parallel algorithms for performing linear stencil computations. Current direct solvers in this domain are computationally inefficient, and Krylov methods require manual labor and mathematical training. We solve these problems for linear stencils by using discrete Fourier transforms preconditioning on a Krylov method to achieve a direct solver that is both fast and general. Indeed, while all currently available algorithms for solving general linear stencils perform Θ( NT ) work, where N is the size of the spatial grid and T is the number of timesteps, our algorithms perform o ( NT ) work. To the best of our knowledge, we give the first algorithms that use fast Fourier transforms to compute final grid data by evolving the initial data for many timesteps at once. Our algorithms handle both periodic and aperiodic boundary conditions and achieve polynomially better performance bounds (i.e., computational complexity and parallel runtime) than all other existing solutions. Initial experimental results show that implementations of our algorithms that evolve grids of roughly 10 7 cells for around 10 5 timesteps run orders of magnitude faster than state-of-the-art implementations for periodic stencil problems, and 1.3× to 8.5× faster for aperiodic stencil problems. Code Repository: https://github.com/TEAlab/FFTStencils

TurboStencil: You only compute once for stencil computation

Stencil computation is an important kind of computational mode that widely used in numeric. It iteratively updates the values of the spatial grid points over multiple time steps according to a given pattern. Existing techniques suffer from high complexity of massive iterative computations. Convolution and fast Fourier transform (FFT) provide the possibility to avoid massive iterations and reduce time complexity of stencil computation. However, current convolution-based fast stencil algorithms cannot effectively solve the problems with aperiodic boundary conditions that are common in practical applications. In this paper, we present a novel algorithm, TurboStencil, for linear stencil computations with aperiodic boundary condition. TurboStencil provides a padding method, eliminating the effects of boundary conditions, to enable convolution for all grid points. For symmetric stencil, TurboStencil only computes once by applying FFT, and thus achieves the time complexity of O(NlogN), where N is the grid data size. For asymmetric stencil, TurboStencil also only computes several times by employing a divide-and-conquer method and FFT, and exhibits a lower complexity than existing stencil algorithms. Experimental results demonstrate that TurboStencil outperforms the state-of-the-art convolution-based fast stencil algorithm by up to 777.1×, 43.2×, and 9.4× for symmetric stencil, and 12.9×, 2.0×, and 1.3× for asymmetric stencil, respectively, on 1D, 2D, and 3D benchmarks.

A sharp interface immersed boundary method for thin-walled geometries in viscous compressible flows

Compressible flows at low Mach numbers are prevalent in various engineering applications, such as combustion, aeroacoustics, heat transfer, and more. In this study, we present a sharp interface immersed boundary method (IBM) tailored for low Mach number viscous compressible flows. The IBM is developed using a robust interpolation scheme, ensuring stability for complex geometries with zero-thickness walls. To satisfy the boundary condition, the interpolation needs to use points from the same side, thus preventing interaction across the zero-thickness wall. Additionally, a simple modification is proposed to fully eliminate interaction between two opposite sides when performing stencil calculations for high order schemes. This is achieved by gradually reducing the order when the stencil is close to the wall. Furthermore, the IBM efficiently addresses the issue of moving objects in compressible flows, specifically the treatment of fresh cells, without necessitating any special treatment. Results showcase the method’s applicability to aerodynamic simulations at low Mach numbers, simulations of compression and expansion processes, and direct numerical simulations of aeroacoustics.

SASA: A Scalable and Automatic Stencil Acceleration Framework for Optimized Hybrid Spatial and Temporal Parallelism on HBM-based FPGAs

Stencil computation is one of the fundamental computing patterns in many application domains such as scientific computing and image processing. While there are promising studies that accelerate stencils on FPGAs, there lacks an automated acceleration framework to systematically explore both spatial and temporal parallelisms for iterative stencils that could be either computation-bound or memory-bound. In this article, we present SASA, a scalable and automatic stencil acceleration framework on modern HBM-based FPGAs. SASA takes the high-level stencil DSL and FPGA platform as inputs, automatically exploits the best spatial and temporal parallelism configuration based on our accurate analytical model, and generates the optimized FPGA design with the best parallelism configuration in TAPA high-level synthesis C++ as well as its corresponding host code. Compared to state-of-the-art automatic stencil acceleration framework SODA that only exploits temporal parallelism, SASA achieves an average speedup of 3.41× and up to 15.73× speedup on the HBM-based Xilinx Alveo U280 FPGA board for a wide range of stencil kernels.

GPU performance analysis for viscoacoustic wave equations using fast stencil computation from the symbolic specification

Seismic forward modeling is a computationally and data-intensive stage in the seismic processing workflow. By profiling the kernels of seismic forward modeling algorithms, it was observed that they need to access a wide variety of memory locations, in addition to the computational cost of performing floating-point operations for the numerical solution of wave equations. In this context, the Roofline model was used to analyze six representative computing kernels in seismic modeling on GPU environment to indicate bottlenecks in the performance and suggest improvements of these wave equation propagators. Based on this, six viscoacoustic equations were implemented using the Devito tool. Experimental data have shown that optimizations in increasing data reuse and decreasing off-chip memory traffic can significantly improve performance.

A compression-based memory-efficient optimization for out-of-core GPU stencil computation

A compression-based memory-efficient optimization for out-of-core GPU stencil computation

Tuning a CPU-Based Stencil Computation in a DPC++ Multi-Device Environment

Reverse Time Migration (RTM) uses the finite-difference (FD) method to compute numerical approximations for the acoustic wave equation. It is a computational bottleneck for RTM applications and therefore needs to be optimized to guarantee timely results and efficiency when allocating resources for hydrocarbon exploration. This article describes our experience reengineering a migrated CUDA-based RTM code to SYCL into a multidevice RTM.

Enhancing DEVITO GPU Allocator Using Unified Memory by NVIDIA

DEVITO is a framework whose objective is to implement optimized stencil computing. Its execution can be carried out both in the CPU and in GPU. For this reason, the data must be manipulated correctly so that, in case of executions in the GPU, they are present in the memory of the GPU at the time of the execution. Natively, DEVITO transfers data every time the operator is executed from OpenACC pragmas. This approach results in performance degradation when the operator is executed repeatedly. To prevent redundant copies and alleviate this bottleneck, an allocator based on unified memory was implemented, which makes manual data transfer between CPU and GPU unnecessary, significantly reducing data transfer time in GPU applications.

An AMR-Based Liquid Film Simulation with Surfactant Transport Using PLIC-HF Method

In this study, an AMR-PLIC-HF method is proposed and implemented by GPU parallel computing based on CUDA programming language and NVIDIA GPU. The present method improves the computation efficiency without compromising the accuracy and conservation of the volume. To satisfy the requirements of stencil points of the PLIC-HF method, an extended stencil computation method based on the tree-based AMR method is proposed and implemented. The Weakly Compressible Scheme (WCS) is used in the present work as a fluid solver. An evolving pressure projection method is adopted to suppress the oscillation induced by the reflection of acoustic waves. The Langmuir model is introduced into the solver to calculate surfactant transport and the Marangoni effect caused by the gradient of the interface concentration of the surfactant. The single vortex flow results verify the accuracy of the AMR-PLIC method. A single bubble rising problem with two different physical property settings is simulated. The results show good agreement with the results given by incompressible solvers. This verifies the accuracy of the two-phase flow solver including the AMR-PLIC-HF method and the WCS. The generation and rupture of liquid film by a single bubble freely rising to an interface is simulated by the present solver with a 1024×2048 AMR grid as the finest resolution. This simulation successfully calculates surfactant transport and the Marangoni effect.

EPSILOD: efficient parallel skeleton for generic iterative stencil computations in distributed GPUs

Iterative stencil computations are widely used in numerical simulations. They present a high degree of parallelism, high locality and mostly-coalesced memory access patterns. Therefore, GPUs are good candidates to speed up their computation. However, the development of stencil programs that can work with huge grids in distributed systems with multiple GPUs is not straightforward, since it requires solving problems related to the partition of the grid across nodes and devices, and the synchronization and data movement across remote GPUs. In this work, we present EPSILOD, a high-productivity parallel programming skeleton for iterative stencil computations on distributed multi-GPUs, of the same or different vendors that supports any type of n-dimensional geometric stencils of any order. It uses an abstract specification of the stencil pattern (neighbors and weights) to internally derive the data partition, synchronizations and communications. Computation is split to better overlap with communications. This paper describes the underlying architecture of EPSILOD, its main components, and presents an experimental evaluation to show the benefits of our approach, including a comparison with another state-of-the-art solution. The experimental results show that EPSILOD is faster and shows good strong and weak scalability for platforms with both homogeneous and heterogeneous types of GPU.

Casper: Accelerating Stencil Computations Using Near-Cache Processing

Stencil computations are commonly used in a wide variety of scientific applications, ranging from large-scale weather prediction to solving partial differential equations. Stencil computations are characterized by three properties: (1) low arithmetic intensity, (2) limited temporal data reuse, and (3) regular and predictable data access pattern. As a result, stencil computations are typically bandwidth-bound workloads, which only experience limited benefits from the deep cache hierarchy of modern CPUs. In this work, we propose Casper, a near-cache accelerator consisting of specialized stencil computation units connected to the last-level cache (LLC) of a traditional CPU. Casper is based on two key ideas: (1) avoiding the cost of moving rarely reused data throughout the cache hierarchy, and (2) exploiting the regularity of the data accesses and the inherent parallelism of stencil computations to increase overall performance.With minimal changes in LLC address decoding logic and data placement, Casper performs stencil computations at the peak LLC bandwidth. We show that by tightly coupling lightweight stencil computation units near LLC, Casper improves performance of stencil kernels by 1.65× on average (up to 4.16×) compared to a commercial high-performance multi-core processor, while reducing system energy consumption by 35% on average (up to 65%). Casper provides 37× (up to 190×) improvement in performance-per-area compared to a state-of-the-art GPU.

Exploiting temporal data reuse and asynchrony in the reverse time migration

Reverse Time Migration (RTM) is a state-of-the-art algorithm used in seismic depth imaging in complex geological environments for the oil and gas exploration industry. It calculates high-resolution images by solving the three-dimensional acoustic wave equation using seismic datasets recorded at various receiver locations. Reverse Time Migration’s computational phases are predominantly composed of stencil computational kernels for the finite-difference time-domain scheme, applying the absorbing boundary conditions, and I/O operations needed for the imaging condition. In this paper, we integrate the asynchronous Multicore Wavefront Diamond (MWD) tiling approach into the full RTM workflow. Multicore Wavefront Diamond permits to further increase data reuse by leveraging spatial with Temporal Blocking (TB) during the stencil computations. This integration engenders new challenges with a snowball effect on the legacy synchronous RTM workflow as it requires rethinking of how the absorbing boundary conditions, the I/O operations, and the imaging condition operate. These disruptive changes are necessary to maintain the performance superiority of asynchronous stencil execution throughout the time integration, while ensuring the quality of the subsurface image does not deteriorate. We assess the overall performance of the new MWD-based RTM and compare against traditional Spatial Blocking (SB)-based RTM on various shared-memory systems using the SEG Salt3D model. The MWD-based RTM achieves up to 70% performance speedup compared to SB-based RTM. To our knowledge, this paper highlights for the first time the applicability of asynchronous executions with temporal blocking throughout the whole RTM. This may eventually create new research opportunities in improving hydrocarbon extraction for the petroleum industry.

HEGrid: A high efficient multi-channel radio astronomical data gridding framework in heterogeneous computing environments

The challenge to fully exploit the potential of existing and upcoming scientific instruments like large single-dish radio telescopes is to process the collected massive data effectively and efficiently. As a ”quasi 2D stencil computation” with the ”Moore neighborhood pattern”, gridding is the most computationally intensive step in data reduction pipeline for radio astronomy studies, enabling astronomers to create correct sky images for further analysis. However, the existing gridding frameworks can either only run on multi-core CPU architecture or do not support high-concurrency, multi-channel data gridding. Their performance is then limited, and there are emerging needs for innovative gridding frameworks to process data from large single-dish radio telescopes like the Five-hundred-meter Aperture Spherical Telescope (FAST). To address those challenges, we developed a High Efficient Gridding framework, HEGrid, by overcoming the above limitations. HEGrid is the first effort to solve the gridding of multi-channel data from the large single-dish radio telescope by multi-pipeline concurrency in the CPU–GPU heterogeneous environment. Specifically, we propose and construct the gridding pipeline in heterogeneous computing environments and achieve multi-pipeline concurrency for high performance multi-channel processing. Furthermore, we propose pipeline-based co-optimization to alleviate the potential negative performance impact of possible intra- and inter-pipeline low computation and I/O utilization, including component share-based redundancy elimination, thread-level data reuse and overlapping I/O and computation. Our experiments are based on both simulated datasets and actual FAST observational datasets. The results show that HEGrid outperforms other state-of-the-art gridding frameworks by up to 5.5x and has robust hardware portability, including AMD Radeon Instinct GPU and NVIDIA GPU.

A method for efficient radio astronomical data gridding on multi-core vector processor

Gridding is the performance-critical step in the data reduction pipeline for radio astronomy research, allowing astronomers to create the correct sky images for further analysis. Like the 2D stencil computation, gridding iteratively updates the output cells by convolution, where the value at each output cell in the space is computed as a weighted sum of neighboring point values. Existing state-of-the-art works have achieved performance improvement of gridding by using multi-core CPUs and GPUs in real-world applications, and their study proved that gridding is a type of scientific computation with high-density computing characteristics. However, low computational performance or high power consumption becomes the main limitation for their processing of large-scale astronomical data. The high-density computing feature of gridding provides opportunities to accelerate it on the multi-core vector processor with vector-SIMD architectures. However, existing works’ (such as those implemented on CPUs or GPUs) task parallelization and data transfer strategies are inefficient to perform gridding directly on the vector processor without any dedicated mapping algorithm.M-DSP is a multi-core vector processor with vector-SIMD architectures designed for the next-generation exascale supercomputer, delivering high performance with ultra-low power consumption. In this paper, we present, for the first time, a novel method to achieve efficient gridding on the M-DSP. Specifically, we propose a gridding workflow designed for the vector-SIMD architectures and present a vectorized version of the gridding convolution algorithm to fully exploit the computational power of the M-DSP. In addition, centering on the processor architectures, we propose task-based parallelization strategies for block and line computing as well as different data loading strategies to achieve high parallel performance and high data transfer efficiency. Experimental results show that our work on M-DSP exhibits very competitive performance compared to other methods running on CPUs or GPUs. This demonstrates the efficiency of our method and the fact that the vector-SIMD architecture is beneficial for scientific computing with ”high density” characteristics, which can exploit its wide vector core and achieve higher performance than its competitors.

Stencil Calculations with Algorithmic Skeletons for Heterogeneous Computing Environments

The development of parallel applications is a difficult and error-prone task, especially for inexperienced programmers. Stencil operations are exceptionally complex for parallelization as synchronization and communication between the individual processes and threads are necessary. It gets even more difficult to efficiently distribute the computations and efficiently implement communication when heterogeneous computing environments are used. For using multiple nodes, each having multiple cores and accelerators such as GPUs, skills in combining frameworks such as MPI, OpenMP, and CUDA are required. The complexity of parallelizing the stencil operation increases the need for abstracting from the platform-specific details and simplify parallel programming. One way to abstract from details of parallel programming is to use algorithmic skeletons. This work introduces an implementation of the MapStencil skeleton that is able to generate parallel code for distributed memory environments, using multiple nodes with multicore CPUs and GPUs. Examples of practical applications of the MapStencil skeleton are the Jacobi Solver or the Canny Edge Detector. The main contribution of this paper is a discussion of the difficulties when implementing a universal Skeleton for MapStencil for heterogeneous computing environments and an outline of the identified best practices for communication intense skeletons.

Scalable communication for high-order stencil computations using CUDA-aware MPI

Modern compute nodes in high-performance computing provide a tremendous level of parallelism and processing power. However, as arithmetic performance has been observed to increase at a faster rate relative to memory and network bandwidths, optimizing data movement has become critical for achieving strong scaling in many communication-heavy applications. This performance gap has been further accentuated with the introduction of graphics processing units, which can provide by multiple factors higher throughput in data-parallel tasks than central processing units. In this work, we explore the computational aspects of iterative stencil loops and implement a generic communication scheme using CUDA-aware MPI, which we use to accelerate magnetohydrodynamics simulations based on high-order finite differences and third-order Runge–Kutta integration. We put particular focus on improving intra-node locality of workloads. Our GPU implementation scales strongly from one to 64 devices at 50%–87% of the expected efficiency based on a theoretical performance model. Compared with a multi-core CPU solver, our implementation exhibits 20–60× speedup and 9–12× improved energy efficiency in compute-bound benchmarks on 16 nodes.

Memory-Aware Functional IR for Higher-Level Synthesis of Accelerators

Specialized accelerators deliver orders of a magnitude of higher performance than general-purpose processors. The ever-changing nature of modern workloads is pushing the adoption of Field Programmable Gate Arrays (FPGAs) as the substrate of choice. However, FPGAs are hard to program directly using Hardware Description Languages (HDLs). Even modern high-level HDLs, e.g., Spatial and Chisel, still require hardware expertise. This article adopts functional programming concepts to provide a hardware-agnostic higher-level programming abstraction. During synthesis, these abstractions are mechanically lowered into a functional Intermediate Representation (IR) that defines a specific hardware design point. This novel IR expresses different forms of parallelism and standard memory features such as asynchronous off-chip memories or synchronous on-chip buffers. Exposing such features at the IR level is essential for achieving high performance. The viability of this approach is demonstrated on two stencil computations and by exploring the optimization space of matrix-matrix multiplication. Starting from a high-level representation for these algorithms, our compiler produces low-level VHSIC Hardware Description Language (VHDL) code automatically. Several design points are evaluated on an Intel Arria 10 FPGA, demonstrating the ability of the IR to exploit different hardware features. This article also shows that the designs produced are competitive with highly tuned OpenCL implementations and outperform hardware-agnostic OpenCL code.