Towards heterogeneous parallelism for SPHinXsys

1. Problem Statement

The paper addresses two distinct but interconnected challenges in the field of computational fluid dynamics (CFD) using Smoothed Particle Hydrodynamics (SPH):

• Physical Modeling Challenge: Existing SPH methods struggle to accurately simulate turbulent, wall-bounded flows (with or without separation) when coupled with the Reynolds-Averaged Navier-Stokes (RANS) model. The core issue is the inconsistency between the Lagrangian nature of SPH and the RANS model, caused by intense shearing and near-wall discontinuities. This leads to issues such as dissipation incompatibility, turbulent kinetic energy (TKE) disturbances, and over-prediction of results.
• Computational Challenge: While SPH is powerful, it suffers from high computational costs due to intensive particle interactions. Existing GPU-accelerated SPH libraries often rely on low-level, vendor-specific directives (e.g., CUDA) or OpenCL, which separate compute kernels from the main codebase. This separation makes the software difficult to program, maintain, and extend, particularly for engineers who are not experts in parallel computing architectures. There is a need for a heterogeneous parallelism approach that leverages both multi-core CPUs and many-core GPUs (like NVIDIA GPUs) without sacrificing code maintainability or open-source development benefits.

2. Methodology

The authors propose a dual-pronged methodology: improving the physical accuracy of the SPH-RANS coupling and re-architecting the SPHinXsys library for heterogeneous computing.

A. Physical Modeling Improvements (SPH-RANS Coupling)

To resolve the Lagrangian-RANS inconsistency, the authors introduce specific treatments for the mainstream and near-wall regions:

• Mainstream Treatments:
  • Adaptive Riemann-eddy Dissipation (ARD): Addresses dissipation incompatibility.
  • De-noised Transport Velocity Formulation: Mitigates turbulent kinetic energy disturbances and prevents over-prediction.
• Near-Wall Treatments:
  • Particle-based Wall Model Realization: Improves the accuracy of the wall model.
  • Weighted Near-Wall Compensation Scheme & Constant Wall-Normal Spacing: Enhances stability and accuracy near boundaries.
  • Wall Dummy Particles: Retained to facilitate future coupling with solid dynamics.
• Verification Tool: A Level-Set-based Boundary-Offset Technique (BOT) is developed to ensure consistent wall-normal distances across different resolutions, enabling rigorous convergence tests.

B. Software Architecture (Heterogeneous Parallelism via SYCL)

The paper details the implementation of SYCL (Data Parallel C++ / DPC++) within the open-source C++ library SPHinXsys.

• Unified Codebase: Unlike previous implementations that require separate kernel definitions, this approach uses a single-source ISO C++ standard. The same source code can be executed sequentially on a CPU, in parallel on a CPU, or in parallel on a GPU (device).
• Three-Layer Kernel Pattern: To handle the constraints of device memory and the inability to use standard C++ features (like virtual functions or dynamic memory allocation) on GPUs, the authors developed an extended programming pattern:
  1. Outer Layer: Handles physical problem definition, global variables, and preprocessing (e.g., particle generation).
  2. Kernel-Shell Layer: Defines the SPH method interface. It uses discrete_ and singular_variables to transfer data from the outer layer to the kernel layer, acting as a bridge.
  3. Kernel Layer: Contains the actual computing logic. Objects here are instantiated and dispatched by the shell layer. Data is passed as raw pointers (arrays) compatible with both host and device.
• Execution Policies: The library defines sequenced_policy, parallel_policy (for CPU), and a new parallel_device_policy (for GPU).
• Algorithmic Adaptations:
  • Neighbor Search: Replaced concurrent_vector (unsupported on devices) with atomic operations for cell linked-list construction.
  • Direct Search: Instead of pre-computing neighbor lists, kernel values are computed directly for particles in nearest cells to save GPU memory.
  • Sorting: Replaced quick_sort with radix_sort for particle sorting, as GPUs perform poorly on recursive algorithms.

3. Key Contributions

1. First Satisfactory Convergence in SPH-RANS: The method achieves, for the first time (to the authors' knowledge), satisfactory convergence for both velocity and turbulent kinetic energy in SPH-RANS simulations of wall-bounded turbulent flows.
2. Heterogeneous Parallelism Framework: Successfully ported SPHinXsys to SYCL, allowing the same codebase to run efficiently on CPUs and GPUs without requiring developers to write low-level, hardware-specific kernels.
3. Open-Source Maintainability: The architecture preserves the benefits of open-source development (unit testing, CI/CD, modularity) while enabling high-performance computing. Developers do not need to understand the execution details of the numerical method to add new features.
4. Benchmarking: Extensive validation against standard benchmarks (dam-break, fish-pass, various channel geometries) and a direct comparison with the established GPU solver DualSPHysics.

4. Results

The paper presents results from several benchmark cases, including straight, curved, and Half Converging-Diverging (HCD) channels, as well as a fish-pass simulation.

• Accuracy: The method produces smoothed and reasonably accurate results for complex turbulent flows.
• Performance on GPU:
  • Using single-precision floating-point on an NVIDIA RTX 2080 Ti, the SYCL implementation achieved up to 27x speedup compared to the CPU execution.
  • Direct Search on GPU with single precision proved to be the most efficient configuration.
• Comparison with DualSPHysics (CUDA-based):
  • Runtime: SPHinXsys was twice as fast as DualSPHysics. For a simulation with 56 million particles, DualSPHysics took 32 hours to simulate 2 seconds, whereas SPHinXsys completed the same in 18 hours.
  • Throughput (GPIPS): SPHinXsys achieved twice the Giga Particle Interactions Per Second (GPIPS) compared to DualSPHysics, confirming that the speedup is due to algorithmic efficiency and not just hardware differences.
  • Validation: The simulation snapshots and data agreed well with existing literature benchmarks.

5. Significance

This work represents a significant step forward in both the physics and software engineering of particle-based simulations:

• Bridging the Gap: It effectively bridges particle-based (SPH) and mesh-based (RANS) modeling, making SPH a viable candidate for industrial-scale turbulent fluid-structure interaction (FSI) simulations.
• Developer Accessibility: By abstracting the hardware complexity through SYCL and a unified codebase, the library lowers the barrier to entry for mechanical engineers and researchers. They can focus on numerical discretization and physics rather than parallelization strategies.
• Future-Proofing: The approach allows for development and testing in environments without GPUs (using standard compilers like GNU on Linux), ensuring that code remains portable and robust. Once verified, the code seamlessly operates on DPC++/GPU environments.
• Scalability: The demonstrated ability to handle 56 million particles efficiently suggests that SPHinXsys is ready for large-scale industrial applications where high-fidelity turbulence modeling is required.