Multiple-GPU accelerated high-order gas-kinetic scheme for direct numerical simulation of compressible turbulence

The CPU is regarded as host, and GPU is treated as device. Data-parallel, compute-intensive operations running on the host are transferred to device by using kernels, and kernels are executed on the device by many different threads. For CUDA, these threads are organized into thread blocks, and thread blocks constitute a grid. Such computational structures build connection with Nvidia GPU hardware architecture. The Nvidia GPU consists of multiple streaming multiprocessors (SMs), and each SM contains streaming processors (SPs). When invoking a kernel, the blocks of grid are distributed to SMs, and the threads of each block are executed by SPs. In summary, the correspondence between software CUDA and hardware Nvidia GPU is simply shown in Fig.1. In the following, the single-GPU accelerated HGKS will be introduced briefly.

Device implementation is handled automatically by the CUDA, while the kernels and grids need to be specified for HGKS. One updating step of single-GPU accelerated HGKS code can be divided into several kernels, namely, initialization, calculation of time step, WENO reconstruction, flux computation at cell interface in $x$, $y$, $z$, and the update of flow variables. The main parts of the single-GPU accelerated HGKS code are labeled in blue as Algorithm.1. Due to the identical processes of calculation for each kernel, it is natural to use one thread for one computational cell when setting grids. Due to the limited number of threads can be obtained in one block (i.e., the maximum $1024$ threads in one block), the whole computational domain is required to be divided into several parts. Assume that the total number of cells for computational mesh is $N_{x}\times N_{y}\times N_{z}$, which is illustrated in Fig.2. The computational domain is divided into $N_{x}/\text{block}_{x}$ parts in $x$-direction, where the choice of $\text{block}_{x}$ for each kernel is an integer defined according to the requirement and experience. If $N_{x}$ is not divisible by $\text{block}_{x}$, an extra block is needed. When setting grid, the variables ”dimGrid” and ”dimBlock” are defined as

As shown in Fig.2, one thread is assigned to complete the computations of a cell $\Omega_{ijk}$, and the one-to-one correspondence of thread index threadidx, block index blockidx and cell index $(i,j,k)$ is given by

Thus, the single-GPU accelerated HGKS code can be implemented after specifying kernels and grids.

Due to the limited computational resources of single-GPU, it is natural to develop the multiple-GPU accelerated HGKS for large-scale DNS of turbulence. The larger computational scale can be achieved with the increased available device memories. It is not straightforward for the programming with multiple-GPUs, because the device memories are not shared across GPUs and the tasks need to be coordinated appropriately. The computation with multiple GPUs is implemented using multiple CPUs and multiple GPUs, where the GPU-GPU communication is coordinated via the MPI. CUDA-aware MPI library is chosen [34], where GPU data can be directly passed by MPI function and transferred in the efficient way automatically.

For the parallel computation with MPI, the domain decomposition are required for both CPU and GPU, where both one-dimensional (slab) and two-dimensional (pencil) decomposition are allowed. For CPU computation with MPI, hundreds or thousands of CPU cores are usually used in the computation, and the pencil decomposition can be used to improve the efficiency of data communications [23]. In terms of GPU computation with MPI, the slab decomposition is used, which reduces the frequency of GPU-GPU communications and improves performances. The domain decomposition and GPU-GPU communications are presented in Fig.3, where the computational domain is divided into $N$ parts in $x$-direction and $N$ GPUs are used. For simplicity, $P_{i}$ denotes the $i$-th MPI process, which deals with the tasks of $i$-th decomposed computational domain in the $i$-th GPU. To implement the WENO reconstruction, the processor $P_{i}$ exchanges the data with $P_{i-1}$ and $P_{i+1}$ using MPI$\_$RECV and MPI$\_$SEND. The wall boundary condition of $P_{0}$ and $P_{N-1}$ can be implemented directly. For the periodic boundary, the data is exchanged between $P_{0}$ and $P_{N-1}$ using MPI$\_$RECV and MPI$\_$SEND as well. In summary, the detailed data communication for GPU computation with MPI are shown in Algorithm.2. Fig.4 shows the HGKS code frame with multiple-GPUs. The initial data will be divided into $N$ parts according to domain decomposition and sent to corresponding process $P_{1}$ to $P_{N-1}$ using MPI$\_$SEND. As for input and output, the MPI process $P_{0}$ is responsible for distributing data to other processors and collecting data from other processors, using MPI$\_$SEND and MPI$\_$RECV. For each decomposed domain, the computation is executed with the prescribed GPU. The identical code is running for each GPU, where kernels and grids are similarly set as single-GPU code [29].

In following studies, the CPU codes run on the desktop with Intel Core i7-9700 CPU using OpenMP directives and TianHe-II supercomputer system with MPI and Intel mpiifort compiler, The HGKS codes with CPU are all compiled with FP64 precision. For the GPU computations, the Nvidia TITAN RTX GPU and Nvidia Tesla V100 GPU are used with Nvidia CUDA and Nvidia HPC SDK. The detailed parameters of CPU and GPU are given in Table.2 and Table.2, respectively. For TITAN RTX GPU, the GPU-GPU communication is achieved by connection traversing PCIe, and there are $2$ TITAN TRTX GPUs in one node. For Tesla V100 GPU, there are $8$ GPUs inside one GPU node, and more nodes are needed for more than $8$ GPUs. The GPU-GPU communication in one GPU node is achieved by Nvidia NVLink. The communication across GPU nodes can be achieved by Remote Direct Memory Access (RDMA) technique, including iWARP, RDMA over Converged Ethernet (RoCE) and InfiniBand. In this paper, RoCE is used for communication across GPU nodes.

As shown in Table.2, Tesla V100 is equipped with more FP64 precision cores and much stronger FP64 precision performance than TITAN RTX. For FP32 precision performance, two GPUs are comparable. Nvidia TITAN RTX is powered by Turing architecture, while Nvidia Tesla V100 is built to HPC by advanced Ampere architecture. Additionally, Tesla V100 outweighs TITAN RTX in GPU memory size and memory bandwidth. Thence, much excellent performance is excepted in FP64 precision simulation with Tesla V100. To validate the performance with FP32 and FP64 precision, numerical examples will be presented.

Taylor-Green vortex (TGV) is a classical problem in fluid dynamics developed to study vortex dynamics, turbulent decay and energy dissipation process [35, 36, 37]. In this case, the detailed efficiency comparisons of HGKS running on the CPU and GPU will be given. The flow is computed within a periodic square box defined as $-\pi L\leq x,y,z\leq\pi L$. With a uniform temperature, the initial velocity and pressure are given by

In the computation, $L=1,V_{0}=1,\rho_{0}=1$, and the Mach number takes $Ma=V_{0}/c_{0}=0.1$, where $c_{0}$ is the sound speed. The fluid is a perfect gas with $\gamma=1.4$, Prandtl number $Pr=1.0$, and Reynolds number $Re=1600$. The characteristic convective time $t_{c}=L/V_{0}$.

In this section, GPU codes are all compiled with FP64 precision. To test the performance of single-GPU code, this case is run with Nvidia TITAN RTX and Nvidia Tesla V100 GPUs using Nvidia CUDA. As comparison, the CPU code running on the Intel Core i7-9700 is also tested. The execution times with different meshes are shown in Table.3, where the total execution time for CPU and GPUs are given in terms of hours and $20$ characteristic convective time are simulated. Due to the limitation of single GPU memory size, the uniform meshes with $64^{3}$, $128^{3}$ and $256^{3}$ cells are used. For the CPU computation, $8$ cores are utilized with OpenMP parallel computation, and the corresponding execution time $T_{8,CPU}^{total}$ is used as the base for following comparisons. The speedups $S_{gc}$ are also given in Table.3, which is defined as

where $T_{8,CPU}^{total}$ and $T_{1,GPU}^{total}$ denote the total execution time of $8$ CPU cores and $1$ GPU, respectively. Compared with the CPU code, 7x speedup is achieved by single TITAN RTX GPU and 16x speedup is achieved by single Tesla V100 GPU. Even though the Tesla V100 is approximately $15$ times faster than the TITAN RTX in FP64 precision computation ability, Tesla V100 only performs $2$ times faster than TITAN RTX. It indicates the memory bandwidth is the bottleneck for current memory-intensive simulations, and the detailed technique for GPU programming still required to be investigated.

For single TITAN TRX with $24$ GB memory size, the maximum computation scale is $256^{3}$ cells. To enlarge the computational scale, the multiple-GPU accelerated HGKS are designed. Meanwhile, the parallel CPU code has been run on the TianHe-II supercomputer [23], and the execution times with respect to the number of CPU cores are presented in Table.4 for comparison. For the case with $256^{3}$ cells, the computational time of single Tesla V100 GPU is comparable with the MPI code with supercomputer using approximately $300$ Intel Xeon E5-2692 cores. For the case with $512^{3}$ cells, the computational time of CPU code with supercomputer using 1024 Intel Xeon E5-2692 cores is approximately $3$ times longer than that of GPU code using $8$ Tesla V100 GPUs. It can be inferred that the efficiency of GPU code with $8$ Tesla V100 GPUs approximately equals to that of MPI code with $3000$ supercomputer cores, which agree well with the multiple-GPU accelerated finite difference code for multiphase flows [26]. These comparisons shows the excellent performance of multiple-GPU accelerated HGKS for large-scale turbulence simulation. To further show the performance of GPU computation, the scalability is defined as

Fig.6 shows the log-log plot for $n$ and $S_{n}$, where $2$, $4$ and $8$ GPUs are used for the case with $256^{3}$ cells, while $8$, $12$ and $16$ GPUs are used for the case with $512^{3}$ cells. The ideal scalability of parallel computations would be equal to $n$. With the log-log plot for $n$ and $T_{n,GPU}^{total}$, an ideal scalability would follow $-1$ slope. However, such idea scalability is not possible due to the communication delay among the computational cores and the idle time of computational nodes associated with load balancing. As expected, the explicit formulation of HGKS scales properly with the increasing number of GPU. Conceptually, the total computation amount increases with a factor of $16$, when the number of cells doubles in every direction. Taking the communication delay into account, the execution time of $256^{3}$ cells with 1 GPU approximately equals to that of $512^{3}$ cells with 16 GPUs, which also indicates the scalability of GPU code. When GPU code using more than $8$ GPUs, the communication across GPU nodes with RoCE is required, which accounts for the worse scalability using $12$ GPUs and $16$ GPUs. Table.5 shows the execution time in flow solver $T_{n,GPU}^{flow}$ and the communication delay between CPU and GPU $T_{n,GPU}^{com}$. Specifically, $T_{n,GPU}^{flow}$ is consist of the time of WENO reconstruction, flux calculation at cell interface and update of flow variables, while $T_{n,GPU}^{com}$ concludes time for MPI$\_$RECV and MPI$\_$SEND for initialization and boundary exchange, and MPI$\_$REDUCE for global time step. The histogram of $T_{n,GPU}^{com}$ and $T_{n,GPU}^{flow}$ with different number of GPU is shown in Fig.6 with $256^{3}$ and $512^{3}$ cells. With the increase of GPU, more time for communication is consumed and the parallel efficiency is reduced accounting for the practical scalability in Fig.6. Especially, the communication across GPU nodes with RoCE is needed for the GPU code using more than $8$ GPUs, which consumes longer time than the communication in single GPU node with NVLink. The performance of communication across GPU nodes with InfiniBand will be tested, which is designed for HPC center to achieve larger throughout among GPU nodes.

The time history of kinetic energy $E_{k}$, dissipation rate $\varepsilon(E_{k})$ and enstrophy dissipation rate $\varepsilon(\zeta)$ of the above nearly incompressible case from GPU code is identical with the solution from previous CPU code, and more details can be found in [23]. With the efficient multiple-GPU accelerated platform, the compressible Taylor-Green vortexes with fixed Reynolds number $Re=1600$ Mach number $Ma=0.25$, $0.5$, $0.75$ and $1.0$ are also tested. Top view of Q-criterion with Mach number $Ma=0.25$, $0.5$, $0.75$ and $1.0$ at $t=10$ using $256^{3}$ cells are presented in Fig.7. The volume-averaged kinetic energy is given by

where $\Omega$ is the volume of the computational domain. The total viscous dissipation is defined as

where $\mu$ is the coefficient of viscosity and $\bm{\omega}=\nabla\times\bm{U}$. The time histories of kinetic energy and total viscous dissipation are given in Fig.8. Compared with the low Mach number cases, the transonic and supersonic cases with $Ma=0.75$ and $Ma=1.0$ have regions of increasing kinetic energy rather than a monotone decline. This is consistent with the compressible TGV simulations [38], where increases in the kinetic energy were reported during the time interval $2\leq t\leq 4$. A large spread of total viscous dissipation curves is observed in Fig.8, with the main trends being a delay and a flattening of the peak dissipation rate for increasing Mach numbers. The peak viscus dissipation rate for the $Ma=0.25$ curve matches closely to that of $Ma=0.1$. Meanwhile, for the higher Mach number cases, the viscous dissipation rate is significantly higher in the late evolution region (i.e., $t\geq 11$) than that from near incompressible case.

For most GPUs, the performance of FP32 precision is stronger than FP64 precision. In this case, to address the effect of FP32 and FP64 precision in terms of efficiency, memory costs and simulation accuracy, both the near incompressible and compressible turbulent channel flows are tested. Incompressible and compressible turbulent channel flows have been studied to understand the mechanism of wall-bounded turbulent flows [2, 4, 8]. In the computation, the physical domain is $(x,y,z)\in[0,2\pi H]\times[-H,H]\times[0,\pi H]$, and the computational mesh is given by the coordinate transformation

where the computational domain takes $(\xi,\eta,\zeta)\in[0,2\pi H]\times[0,3\pi H]\times[0,\pi H]$ and $b_{g}=2$. The fluid is initiated with density $\rho=1$ and the initial streamwise velocity $U(y)$ profile is given by the perturbed Poiseuille flow profile, where the white noise is added with $10\%$ amplitude of local streamwise velocity. The spanwise and wall-normal velocity is initiated with white noise. The periodic boundary conditions are used in streamwise $x$-direction and spanwise $z$-directions, and the non-slip and isothermal boundary conditions are used in vertical $y$-direction. The constant moment flux in the streamwise direction is used to determine the external force [23]. In current study, the nearly incompressible turbulent channel flow with friction Reynolds number $Re_{\tau}=395$ and Mach number $Ma=0.1$, and the compressible turbulent channel flow with bulk Reynolds number $Re=3000$ and bulk Mach number $Ma=0.8$ are simulated with the same set-up as previous studies [23, 39]. For compressible case, the viscosity $\mu$ is determined by the power law as $\mu=\mu_{w}(T/T_{w})^{0.7}$, and the Prandtl number $Pr=0.7$ and $\gamma=1.4$ are used. For the nearly incompressible turbulent channel flow, two cases $G_{1}$ and $G_{2}$ are tested, where $256\times 128\times 128$ and $256\times 256\times 128$ cells are distributed uniformly in the computational space as Table.7. For the compressible flow, two cases $H_{1}$ and $H_{2}$ are tested as well, where $128\times 128\times 128$ and $128\times 256\times 128$ cells are distributed uniformly in the computational space as Table.7. In the computation, the cases $G_{1}$ and $H_{1}$ are implemented with single GPU, while the cases $G_{2}$ and $H_{2}$ are implemented with double GPUs. Specifically, $\Delta y^{+}_{min}$ and $\Delta y^{+}_{max}$ are the minimum and maximum mesh space in the $y$-direction. The equidistant meshes are used in $x$ and $z$ directions, and $\Delta x^{+}$ and $\Delta z^{+}$ are the equivalent plus unit, respectively. The fine meshes arrangement in case $G_{2}$ and $H_{2}$ can meet the mesh resolution for DNS in turbulent channel flows. The $Q$-criterion iso-surfaces for $G_{2}$ and $H_{2}$ are given Fig.9 in which both cases can well resolve the vortex structures. With the mesh refinement, the high-accuracy scheme corresponds to the resolution of abundant turbulent structures.

Because of the reduction in device memory and improvement of arithmetic capabilities on GPUs, the benefits can be achieved by using FP32 precision compared to FP64 precision. In view of these strength, FP32-based and mixed-precision-based high-performance computing start to be explored [30, 31]. However, whether FP32 precision is sufficient for the DNS of turbulent flows still needs to be discussed. It is necessary to evaluate the effect of precision on DNS of compressible turbulence. Thus, the GPU-accelerated HGKS is compiled with FP32 precision and FP64 precision, and both nearly incompressible and compressible turbulent channel flows are used to test the performance in different precision. For simplicity, only the memory cost and computational time for the nearly incompressible cases are provided in Table.8, where the execution times $T_{2,GPU}^{total}$ are given in terms of hours and the memory cost is in GB. In the computation, two GPUs are used and $10H/U_{c}$ statistical period is simulated. In Table.8, $S_{fp}$ represents the speed up of $T_{2,GPU}^{total}$ with FP32 to $T_{2,GPU}^{total}$ with FP64, and $R_{fp}$ denotes the ratio of memory cost with FP32 to memory cost with FP64. As expected, the memory of FP32 precision is about half of that of FP64 precision for both TITAN RTX and Tesla V100 GPUs. Compared with FP64-based simulation, 2.7x speedup is achieved for TITAN RTX GPU and 1.7x speedup is achieved for Tesla V100 GPU. As shown in Table.8, the comparable FP32 precision performance between TITAN RTX and Tesla V100 also provides comparable $T_{2,GPU}^{total}$ for these two GPUs. We conclude that Turing architecture and Ampere architecture perform comparably for memory-intensive computing tasks in FP32 precision.

To validate the accuracy of HGKS and address the effect with different precision, the statistical turbulent quantities are provided for the nearly incompressible and compressible turbulent channel flows. For FP32-FP64 comparison, all cases restarted from the same flow fields, and $200H/U_{c}$ statistical periods are used to average the computational results for all cases. For the nearly incompressible channel flow, the logarithmic formulation is given by

where the von Karman constant $\kappa=0.40$ and $B=5.5$ is given for the low Reynolds number turbulent channel flow [2]. The mean streamwise velocity profiles with a log-linear plot are given in Fig.11, where the HGKS result is in reasonable agreement with the spectral results [4]. In order to account for the mean property of variations caused by compressibility, the Van Driest (VD) transformation [40] for the density-weighted velocity is considered

where $\langle\cdot\rangle$ represents the mean average over statistical periods and the X- and Z-directions. For the compressible flow, the transformed velocity is expected to satisfy the incompressible log law Eq.(10). The mean streamwise velocity profiles $\langle U\rangle^{+}$ and $\langle U\rangle_{VD}^{+}$ with VD transformation as Eq.(11) are given in Fig.11 in log-linear plot. HGKS result is in reasonable agreement with the reference DNS solutions with the mesh refinements [9]. The mean velocity with FP32 and FP64 precision are also given in Fig.11 and Fig.11. It can be observed that the qualitative differences between different precision are negligible for statistical turbulent quantities with long time averaging.

For the nearly incompressible turbulent channel flow, the time averaged normalized root-mean-square fluctuation velocity profile $U_{rms}^{+}$, $V_{rms}^{+}$, $W_{rms}^{+}$ and normalized Reynolds stress profiles $-<U^{\prime}V^{\prime}>$ are given in Fig.12 for the cases $G_{1}$ and $G_{2}$. The root mean square is defined as $\phi_{rms}^{+}=\sqrt{(\phi-\langle\phi\rangle)^{2}}$ and $\phi^{\prime}=\phi-\langle\phi\rangle$, where $\phi$ represents the flow variables. With the refinement of mesh, the HGKS results are in reasonable agreement with the reference date, which given by the spectral method with $256\times 193\times 192$ cells [4]. It should also be noted that spectral method is for the exact incompressible flow, which provides accurate incompressible solution. The HGKS is more general for both nearly incompressible flows and compressible flows. For the compressible flow, the time averaged normalized root-mean-square fluctuation velocity profiles $U_{rms}^{+}$, $V_{rms}^{+}$, $W_{rms}^{+}$, Reynolds stress $\langle-\rho U^{{}^{\prime}}V^{{}^{\prime}}\rangle/\langle\tau_{w}\rangle$, the root-mean-square of Mach number $M_{rms}^{+}$ and turbulent Mach number $M_{t}$ are presented in Fig.13 for the cases $H_{1}$ and $H_{2}$. The turbulent Mach number is defined as $M_{t}=q/\left\langle c\right\rangle$, where $q^{2}=\langle(U^{{}^{\prime}})^{2}+(V^{{}^{\prime}})^{2}+(W^{{}^{\prime}})^{2}\rangle$ and $c$ is the local sound speed. The results with $H_{2}$ agree well with the refereed DNS solutions, confirming the high-accuracy of HGKS for DNS in compressible wall-bounded turbulent flows. The statistical fluctuation quantities profiles with FP32 precision and FP64 precision are also given in Fig.12 and Fig.13. Despite the deviation of $U_{rms}^{+}$ for different precision, the qualitative differences in accuracy between FP32 and FP64 precision are acceptable for statistical turbulent quantities with such a long time average.

The instantaneous profiles of $U^{+}_{rms}$, $V^{+}_{rms}$, $W^{+}_{rms}$ at $t=20H/U_{c}$ and $200H/U_{c}$ are shown in Fig.14 for the cases $G_{2}$ and $H_{2}$ with FP32 and FP64 precision, where cases restarted with the same initial flow field for each case. Compared with the time averaged profiles, the obvious deviations are observed for the instantaneous profiles. The deviations are caused by the round-off error of FP32 precision, which is approximately equals to or larger than the errors of numerical scheme. Since detailed effect of the round-off error is not easy to be analyzed, FP32 precision may not be safe for DNS in turbulence especially for time-evolutionary turbulent flows. In terms of the problems without very strict requirements in accuracy, such as large eddy simulation (LES) and Reynolds-averaged Navier–Stokes (RANS) simulation in turbulence, FP32 precision may be used due to its improvement of efficiency and reduction of memory. It also strongly suggests that the FP64 precision performance of GPU still requires to be improved to accommodate the increasing requirements of GPU-based HPC.