Storms of high-energy particles (EduHPC'18 Peachy Assignment)
In this problem, we are provided with a surface (layer) on which cycles occur, consisting of two phases:
- Storm
- Relaxation
During the storm phase, particles fall onto the layer at specific positions. Energy is released into the cell where the particle lands, as well as into its adjacent cells.
This energy is then stabilized during the relaxation phase. In this second phase, the value of each cell becomes the average of the cell itself and its two adjacent cells (if they exist).
After each relaxation phase, the value and position of the layer's current relative maximum must be identified.
At the end of the execution, the program will return two arrays with a size equal to the number of storm cycles:
- The first array will contain the maximum value of storm i at position i.
- The second array will contain the position of that maximum.
You can find more information about the problem in handout.pdf.
Details on our implementations can be found in report.pdf (here).
Key highlights are:
- Domain decomposition: uses a 1D decomposition to distribute workload across MPI ranks.
-
Latency hiding: inter-process halo exchanges are performed using non-blocking primitives (
MPI_Isend,MPI_Irecv) to overlap communication with the computation of internal cells. -
Thread synchronization: minimizes overhead by encapsulating the simulation loop in a single
#pragma omp parallelregion. -
False sharing prevention: uses a padded
ThreadMaxstructure (64 bytes) to ensure thread-local data resides on unique cache lines during reduction. -
Memory efficiency: Implements
$O(1)$ "zero-copy" pointer swapping for double-buffering and a NUMA first-touch policy for optimal memory locality. - Cache tiling: employs a cache blocking (tiling) strategy with a tile size of 2048 floats to ensure impact logic operates at L1 speeds.
- Dual-path energy update: optimizes memory bandwidth by routing data through constant memory for smaller storms or the read-only data cache (via
__ldg()) for larger ones. - Warp shuffling: Leverages
__shfl_down_syncprimitives for high-speed intra-warp tree reductions at register speed. - Optimized occupancy: uses a "sweet spot" of 256 threads per block for balanced performance across different domain sizes.
Our tests on a 128-core cluster showed:
- Scalability: - pure MPI reached a peak speedup of 117.86x on 128 cores.
- Efficiency: - the CUDA implementation outperformed the 128-core CPU configuration by over 2x on balanced workloads (1,000,000 elements).
More information on our testing results can also be found in report.pdf (here).
For the MPI / OpenMP implementation, the project is divided into two main files:
energy_storms_mpi_omp.cenergy_storms_mpi_omp_core.c
For the CUDA implementation:
energy_storms_cuda.cppenergy_storms_cuda_core.cu
In both cases, we recommend using make to run the tests. You can check the configuration in the Makefile and modify the compilation and execution flags in student.mk.
You can run the commands defined in the Makefile, such as make run_mpi or make run_cuda, to execute Test #2 with MPI / OpenMP or CUDA, respectively.
At some point we needed more flexibility for testing, so we created a test.sh script for MPI / OpenMP and testcuda.sh for CUDA, which allows us to run the executable and specify exactly which tests to launch each time.
You can change some variables for testing:
- For MPI / OpenMP, you can change the number of processes with the variable
DEFAULT_PROCSintest.sh. - For CUDA, you can change the block size with the constant
BLOCK_SIZEinenergy_storms_cuda_core.cu.
Then, you can use the following commands:
# MPI / OpenMP
./test.sh [test numbers]
# CUDA
./testcuda.sh [test numbers]For example, running ./test.sh 6 7 will execute Test #6 and Test #7 with the number of processes specified in the script.
(to run the CUDA code, you will obviously need an Nvidia GPU (or you can try some tools like hipify on AMD) - otherwise, you can test it on platforms like Google Colab).