LCS.jl: A High-Performance, Multi-Platform Computational Model in Julia for Turbulent Particle-Laden Flows

Abstract: Multiphase turbulent flow phenomena are observed not only in industrial devices but also in environmental flows, and direct numerical simulation (DNS) plays a key role in their investigation. Many numerical models have been developed; nevertheless, few models are highly optimized for GPU platforms, which represent the current mainstream in high-performance computing (HPC). In this study, we developed LCS.jl (Lagrangian Cloud Simulator in Julia), a single-source and multi-platform multiphase turbulence simulation model implemented in Julia language and KernelAbstractions.jl. Validation results confirmed that the present fluid and particle statistics agree well with those obtained in prior studies. A GPU-native particle communication algorithm based on prefix-scan reduced the particle communication cost from approximately 78% (CPU-delegated) to 10% of total execution time. LCS.jl achieved computational performance equivalent to the Fortran implementation in many-processes computations. For GPUs, strong scaling efficiency was maintained above 85% (up to 256 GPUs) and weak scaling efficiency above 90% (up to 216 GPUs) on TSUBAME4.0 (a GPU supercomputer at the Institute of Science Tokyo). LCS.jl achieved a maximum speedup of 18.0x on GPUs over CPUs. A trial heterogeneous execution achieved a 72% reduction in execution time compared to the CPU-only configuration even in configurations where the GPU was not the primary compute device. These results demonstrate that LCS.jl is a multiphase turbulence simulation platform that achieves both portability and scalability across a variety of computational resource configurations.

• The paper introduces a novel Julia-based DNS model for turbulent particle-laden flows that leverages GPU-native optimizations to eliminate particle communication bottlenecks.
• It employs fourth-order finite differences and an RK2 scheme to accurately resolve fluid-particle dynamics, validated against established turbulence statistics.
• LCS.jl achieves over 85% strong scaling efficiency on GPUs and up to 18× speedup over CPU execution, demonstrating its multi-platform performance and adaptability.

LCS.jl: High-Performance, Multi-Platform Simulation for Turbulent Particle-Laden Flows

Introduction

The development of LCS.jl, a high-performance direct numerical simulation (DNS) model for turbulent particle-laden flows fully implemented in Julia, addresses core challenges in simulating environmental and industrial multiphase turbulence at high Reynolds numbers. Traditionally, Euler–Lagrangian DNS frameworks for such problems have been restricted by CPU-centric codebases and suboptimal scaling on modern high-performance computing (HPC) architectures dominated by GPUs. LCS.jl introduces a single-source, multi-platform design facilitated by Julia and KernelAbstractions.jl, enabling efficient deployment and strong scaling across both CPU and GPU clusters. The implementation strategically targets both algorithmic efficiency—particularly in particle communication—and performance portability, supporting current and future heterogeneous HPC environments.

Model Formulation and Numerical Methods

LCS.jl employs incompressible Navier–Stokes equations for the carrier fluid, discretized on a uniform Cartesian mesh using conservative fourth-order finite differences for convective and viscous terms, and second-order differences for pressure. Time integration leverages a two-stage, second-order Runge–Kutta (RK2) scheme, with a pressure–velocity coupling resolved via a highly simplified marker-and-cell (HSMAC) algorithm. Turbulent forcing is achieved using the Reduced Communication Forcing (RCF) method, ensuring energy injection at large scales with reduced inter-process communication.

The particle phase evolves according to a nonlinear Stokes drag model, where inertial particles are tracked with individualized ODEs for position and velocity. The adaptive time stepping enforces convective and viscous CFL constraints. Fluid-particle coupling is handled by trilinear interpolation. Validation against reference statistics—root-mean-square velocity, Reynolds number, and higher-order moments—demonstrates numerical fidelity of the fluid and particle algorithm.

Platform Abstraction and Communication Optimizations

Achieving platform-agnostic code without performance regression is realized using KernelAbstractions.jl for kernel parallelism and explicit memory movement. LCS.jl introduces a 'foraxes' abstraction for loop traversal, generating platform-specific code suited for CPUs, CUDA, AMD, Metal, and other backends, while unifying the user kernel logic (cf. Listing 1 in the source).

Critical to multi-GPU DNS is efficient data exchange. LCS.jl applies advanced communication optimizations in two aspects:

• HALO Exchange: The communication–computation overlap hides boundary data transfer latency by interleaving communication and computational kernels, while time-blocking aggregates multiple iteration steps to reduce communication frequency.
• Particle Communication: The GPU-native prefix-scan algorithm enables dynamic, parallel determination of send buffers for particles migrating across subdomain boundaries. This effectively eliminates the global sequential dependency present in legacy CPU-delegated schemes and dramatically reduces wall-clock time allocated to particle communication.

  Figure 2: Speedup of wall-clock time per time-integration step as a function of GPU count, demonstrating gains from communication–computation overlap, time-blocking, and their combination.

Validation and Physical Results

Physical validation was performed for fluid and particle statistics across a range of grid resolutions ($N^3 = 128^3$ through $2048^3$) and Taylor-microscale Reynolds numbers. The model recovers expected turbulence scaling laws: kinetic energy spectra normalized by Kolmogorov scales show extended inertial ranges aligning closely with Kolmogorov's $k^{-5/3}$ scaling.

Figure 3: Energy spectra $E(k)$ normalized by Kolmogorov scaling across multiple resolutions, confirming resolution of inertial and dissipative ranges.

For the particle phase, the validated radial distribution function at contact, $g(r=R)$, demonstrates Stokes number and Reynolds number dependence consistent with previous DNS studies. Clustering statistics crucial for cloud microphysics modeling are quantitatively reproduced.

Figure 4: Radial distribution function $g(r=R)$ at contact as a function of Reynolds number, showing agreement with reference simulations and correct reproduction of clustering dynamics.

GPU and CPU Performance: Scaling and Speedup

LCS.jl's single-source design achieves performance parity with hand-tuned Fortran codes at scale. On the TSUBAME4.0 supercomputer, strong scaling efficiency on GPUs exceeded 85% up to 256 GPUs and weak scaling efficiency surpassed 90% up to 216 GPUs.

Figure 1: Strong scaling—wall-clock time per step versus device number for fixed-sized DNS on GPUs and CPUs compared against ideal scaling.

Figure 7: Weak scaling—wall-clock time per step for fixed device-local sizes, demonstrating near-ideal scaling on both GPU and CPU clusters.

For comparable workloads, GPU execution of LCS.jl achieves up to 18$\times$ speedup versus CPU execution, after adjusting for optimal thread/process binding and memory bandwidth utilization.

Notably, the GPU-native prefix-scan implementation reduces particle communication from 78% (when CPU delegation is used) to 10% of total execution time per time-integration step, effectively removing the primary scaling bottleneck for large-scale particle-laden turbulent DNS.

Heterogeneous Execution and Resource Flexibility

LCS.jl offers out-of-the-box support for heterogeneous execution wherein computation and analysis modules are dynamically dispatched to different compute resources. A configuration running fluid and particle time stepping on CPUs and statistics computation on a GPU achieved a 72% reduction in total wall-clock time relative to CPU-only execution, even when the GPU was not the primary simulation device. This establishes LCS.jl as a resource-adaptive platform, suitable for gradual or partial GPU cluster upgrades.

Communication Patterns and Topological Effects

The analysis of communication reveals topology-induced performance variations. Power-of-two (e.g., $4\times4\times4$) and non-power-of-two (e.g., $3\times3\times3$) GPU process grids exhibit distinct inter/intra-node communication patterns. Power-of-two topologies retain at least one intra-node direction per GPU, enhancing performance and maintaining expected scaling, whereas non-power-of-two grids can impose full inter-node communication for many ranks, leading to measurable overhead.

Figure 5: Communication patterns demonstrating the benefits of power-of-two domain decompositions for minimizing inter-node exchanges in GPU clusters.

Implications and Future Outlook

The results of LCS.jl have significant implications for DNS of particle-laden turbulence in scientific and engineering domains including atmospheric microphysics, cloud physics, and industrial multiphase reactors:

• Performance Portability: The architecture-agnostic design ensures adaptability to evolving HPC hardware, minimizing refactor cost and facilitating reproducibility and long-term sustainability.
• Algorithmic Scalability: By eliminating the traditional particle communication bottleneck, the platform unlocks DNS at higher Reynolds numbers and larger particle counts, previously inaccessible due to communication overhead.
• Flexible Resource Utilization: Heterogeneous execution strategies can be exploited in mixed-infrastructure environments, integrating into existing workflows and extending the lifecycle of legacy CPU clusters via incremental GPU adoption.
• Open Science: With the public release of LCS.jl, further development, benchmarking, and application across disciplines become possible.

Potential future developments include integration of more complex microphysical models, adaptive mesh refinement, or coupling with weather and climate models. Extensions leveraging Julia's evolving HPC ecosystem, such as unified transpilation frameworks (Nicusan et al., 22 Jul 2025), may further enhance performance and scalability.

Conclusion

LCS.jl represents a significant advance in the simulation of turbulent particle-laden flows, unifying high physical fidelity, scalability, and performance portability across computing platforms. Through algorithmic advancements such as GPU-native prefix-scan communication and rigorous validation, LCS.jl demonstrates that performance and portability need not be mutually exclusive. The modular, open-source Julia implementation positions LCS.jl as a robust platform for research in current and future heterogeneous HPC environments, facilitating exploration of multiphase turbulence regimes at unprecedented resolution and scale.

Reference: "LCS.jl: A High-Performance, Multi-Platform Computational Model in Julia for Turbulent Particle-Laden Flows" (2604.11008).