Simulation of a turbulent flow subjected to favorable and adverse pressure gradients

Abstract

This paper reports the results from a direct numerical simulation of an initially turbulent boundary layer passing over a wall-mounted “speed bump” geometry. The speed bump, represented in the form of a Gaussian distribution profile, generates a favorable pressure gradient region over the upstream half of the geometry, followed by an adverse pressure gradient over the downstream half. The boundary layer approaching the bump undergoes strong acceleration in the favorable pressure gradient region before experiencing incipient or very weak separation within the adverse pressure gradient region. These types of flows have proved to be particularly challenging to predict using lower-fidelity simulation tools based on various turbulence modeling approaches and warrant the use of the highest fidelity simulation techniques. The present direct numerical simulation is performed using a flow solver developed exclusively for graphics processing units. Simulation results are utilized to examine the key phenomena present in the flowfield, such as relaminarization/stabilization in the strong acceleration region succeeded by retransition to turbulence near the onset of adverse pressure gradient, incipient/weak separation and development of internal layers, where the sense of streamwise pressure gradient changes at the foot, apex and tail of the bump.

Appendix

1.1 GPU code scaling performance

The flow solver implements strategies that overlap communication with computation, which are crucial for achieving good performance on GPUs. This overlap is achieved as follows. The GPUs first compute the near-boundary information needed by their neighbors and copy the data to their host central processing units (CPUs) for MPI communication. The GPUs then compute the interior points while the host CPUs handle the MPI communication in parallel. The host CPUs run several OpenMP threads (one dedicated thread per face of a grid block) to perform the MPI data exchange. Assuming a large enough GPU workload, by the time the GPUs have finished computing interior points, the MPI communication among the host CPUs has been completed and the host CPUs have copied the exchanged MPI data back to their corresponding GPUs. The GPUs can then update the points near their block interface using the data received from their neighbors and can proceed further. The data copies between the GPU and the host CPU take place asynchronously, meaning that the GPU can do the computing (as long as there is sufficient work) and the data exchange with the host CPU simultaneously. Moreover, the GPU can run several independent tasks in parallel, as long as the resources needed by the computational kernels are available. The GPU operations are synchronized at critical points.

During the course of the developmental work, it was realized that the chosen explicit algorithms run extremely fast on the GPU since they are essentially made up of many independent multiply–add operations, at which the GPU excels. However, the host-assisted MPI communication among the GPUs was found to be a major bottleneck, as already known from others’ similar experience in GPU code development. We observed that, unless we assign a large workload to the GPU and overlap computation with communication, MPI communication over the network cannot keep up with GPU computation. When the workload per GPU is too small, an inefficiency arises because the GPU finishes its work way too quickly and then sits idle while waiting for the communication to catch up. This idle time is nothing but a wasted opportunity to have the GPU do useful work. Hence, in order to have the GPU run at full efficiency and make the best use of limited computational resources, the present strategy is to assign the maximum possible workload to the GPU and overlap computation with communication as much as possible. To achieve the maximum possible workload, we adjust the grid block size per GPU in order to nearly or fully max out the available GPU global memory, which is 16 GB per GPU on the Summit system. This corresponds to about 24–25 million points per GPU on Summit. The GPUs at the NASA Advanced Supercomputing Division (NAS) facility have 32 GB of global memory, which allows nearly 50 million points per GPU on that system.

In this work, the host-assisted MPI communication approach is chosen over the so-called “GPU-aware MPI” option, which accomplishes direct communication among the GPUs without the complication of sending the data through the host CPUs. This choice was made because the performance of the GPU-aware MPI implementation available at the time of code development was found to be rather poor and unacceptable. Although the implementation of host-assisted MPI communication is more involved, assigning a large workload to the GPUs and having the host CPUs handle the MPI communication overlaps computation with communication effectively and hides the communication cost, as noted above. It is therefore well worth the additional complication since it improves the overall code performance.

Because of the slow communication issue relative to computation, we do not expect to see good strong scaling performance on the current systems. Communication simply would not be able to keep up with computation if the GPUs were assigned too little work. This issue can only be resolved by a much faster communication network, which presently does not exist. Hence, the only way to make the best use of limited resources at present is to run the code in its maximum efficiency mode. Although the maximum efficiency mode may not minimize the wall-clock run time, it will certainly minimize the total number of node-hours needed for a given problem size, which is the more relevant performance metric given the fact that the available resources are limited.

A weak scaling study has been performed on the Summit system, in which the total workload per GPU is kept fixed while increasing the total number of GPUs. Each node on Summit contains 6 NVIDIA Tesla V100 GPUs, with 16 GB of global memory per GPU. Our test case for the scaling study involves the flow over a flat plate. For the weak scaling study runs, the total number of GPUs is varied from 150 (25 nodes) up to 6000 (1000 nodes). Note that 1000 nodes make up $21.7\%$ of the entire system. For each case, we assign a grid block of $256 \times 320 \times 288$ points to each GPU. Test runs for each node-count case were performed to measure the total time taken to run 1000 time steps and compute the corresponding average time per step. The third-order, three-stage explicit Runge–Kutta scheme is used for time advancement. Figure 31 shows how the average time taken per step, in milliseconds, varies as the number of Summit nodes are increased. Note that the horizontal axis is plotted in logarithmic scale. Depending upon how the nodes assigned to a particular job are distributed within the system and the network route among the nodes, as well as the network traffic due to other jobs running on the system, code execution speeds may vary from one run to another. We therefore see a slight variation in the average time per step as the number of nodes changes. The average time per step is nearly constant and remains between 149 and 150 ms. This is the payoff for utilizing a strategy that assigns a large workload to the GPU and overlaps GPU computation with MPI communication as much as possible. The “mean” value of the time taken per step (i.e., the value averaged over all node-count runs) is about 149.5 ms for the given grid block size per GPU. This corresponds to roughly one nanosecond per grid point per Summit node per time step.

Fig. 31

Weak scaling performance on the Summit system

The overall performance of the GPU code was measured to be $10\%$ of the peak theoretical double-precision performance of the V100 GPU, which is rated at 7.8 Tera-floating-point-operations (TFLOPs) per second. The most compute-intensive individual kernels (or subroutines) were found to achieve as high as $17\%$ of the peak. To determine the performance limiting factor, the code was profiled using the available performance analysis tools. For a grid block of about 24 million points per GPU, it was found that, at each time step, the code moves about 71 GB worth of total data from the global memory to the registers of the processing units, and about another 39 GB worth of total computed data from the registers back to the global memory. Thus, there is a data movement of about 110 GB at each time step, between the global memory and the registers. Now, in order to illustrate what limits the code performance, suppose the code only moves this much data between the main memory and registers, but does not do any computing at all. How long would this data movement alone take? To answer this question, we first note that the V100 GPU architecture provides 900 GB/s of peak memory bandwidth. For simplicity, let us assume that the data movement takes place at the peak bandwidth. That would mean that at least $110/900 \approx 0.122$ s would be needed to move that amount of data. Our profiling measurements show that the most compute-intensive kernels generally achieve about 80 to $90\%$ of the peak memory bandwidth, so the actual data movement takes a bit longer than this estimate. With an overall average rate of $85\%$ of the peak bandwidth, the data movement would take about 0.144 s. We also know that for the given grid block size per GPU, the code takes about 0.15 s to perform all operations and advance the simulation for one time step. The profiling measurements show that the actual compute time is around $19\%$ of the total elapsed time. Even though the actual computations overlap with the memory operations, the data movement between the memory and registers still constitutes a significant chunk of the time taken per computational time step according to the above timings. These observations lead us to the conclusion that our code performance is bound by the available memory bandwidth. In other words, the memory bandwidth is not sufficient to transfer data into and out of the registers at the rate demanded by the processing units. Thus, the well-known adage of the computer science world, “The FLOPs are free, you are paying for the memory bandwidth!” is still very much valid in our case.

As noted above, with about 80 to $90\%$ of the peak memory bandwidth achieved by the most compute-intensive kernels, our memory-bound GPU code is not far off from its maximum possible performance; hence, reaching a much greater percentage of the peak theoretical FLOPs per second performance of the V100 architecture is not feasible. We anticipate that potential memory bandwidth improvements in future-generation GPU architectures should enable our code to achieve higher FLOP counts per second on those systems.

Fig. 32

Instantaneous snapshots of turbulent channel flow

Fig. 33

Turbulent channel flow mean streamwise velocity and Reynolds stress component profiles in wall units

1.2 Performance comparison to CPU flow solver

We now provide the performance comparison between this new GPU code and our previous CPU code, which was most recently used in the simulation of flow separation problems [23, 24]. The speedup factors depend on how the comparisons are made. For example, the “node-to-node” speedup factor, which is derived from the performance comparison based on one Summit node with 6 GPUs versus one dual-socket Intel Skylake CPU node with 40 cores, comes out to about $75\times$. This comparison is for the version of the CPU code based on the third-order, three-stage explicit Runge–Kutta scheme, same as that employed in the GPU code. Note that the CPU code uses high-order compact finite-difference [1] and filtering schemes [7, 26], whose implementations are well-optimized for the CPU. We estimate that the implementation of explicit finite-difference and filtering schemes, employed in the GPU code, into the CPU solver would result in a performance improvement in only about $25\%$. Even in such a case, the node-to-node speedup factor would still be about $60\times$. Basing the comparison on the number of GPUs versus the number of CPU cores would lead to an equivalent speedup factor of $60 \times 40/6 = 400\times$, meaning that one GPU is worth 400 CPU cores. We also note that one dual-socket Intel Skylake CPU, which contains 40 cores total, is comparable in price and power consumption to a V100 GPU. If we were to base the comparison on one V100 GPU versus one dual-socket Intel Skylake CPU, the corresponding speedup factor would be $400/40=10\times$.

The CPU code also has a second-order implicit time integration version. It is about 2.5 times more costly per time step than the explicit integration scheme. The explicit time integration scheme is normally run at a CFL number of around 0.8. Because of time accuracy concerns, the second-order implicit time integration scheme should not normally be run at a CFL number greater than 5 or so. This corresponds to a factor of about 6.25 increase in the time step with the implicit scheme. The increased computational cost per time step of the implicit scheme gives a speedup factor of about $6.25/2.5 = 2.5$ over the explicit scheme. Thus, comparing the GPU code to the implicit version of CPU code as is, we obtain a node-to-node speedup factor of $75/2.5 = 30\times$. The corresponding comparison based on one V100 GPU versus one dual-socket Intel Skylake CPU would yield a speedup factor of $30/6=5\times$.

As noted earlier, it is possible to further accelerate the GPU code by switching to a second-order, single-stage explicit scheme developed by Verstappen and Veldman [25]. The second-order scheme runs at half the time step of the third-order Runge–Kutta scheme but requires only one right-hand-side computation per time step as opposed to three right-hand-side evaluations needed for the Runge–Kutta scheme. It therefore provides an acceleration factor of $3/2=1.5\times$ over the Runge–Kutta scheme. With the second-order explicit time advancement scheme implemented in the GPU code, the performance comparison to the implicit version of CPU code (which is also second-order accurate in time), based on one V100 GPU versus one dual-socket Intel Skylake CPU, would provide a speedup factor of $7.5\times$.

We should note here that a significant effort was also put into optimizing the implicit time integration version of the CPU code discussed above. The explicit time integration version of the CPU code is essentially the conversion from the optimized implicit version. Hence, the above performance comparisons are between the fastest running versions of the GPU and CPU codes, and are as fair as currently possible.

We also note that the connection speed between the GPU and the host CPU has some effect on the overall code performance. Because of the slower GPU to CPU connection on the NVIDIA Tesla V100 GPU nodes installed at NAS (which is PCI-Express, while Summit nodes have the faster NVLink⁷), the GPU flow solver runs about $5\%$ slower on the NAS GPU nodes relative to Summit. The block interface data to be exchanged among the GPUs during the computation are first copied to their respective host CPUs, which then perform the actual data exchange via MPI communication and copy the exchanged data back to the GPUs. Although the code overlaps communication with computation as much as possible, the GPU to CPU connection speed and specific MPI implementation details still play a role on the overall performance.

1.3 GPU code validation

For code validation, a turbulent channel flow problem was considered. The Reynolds number of the fully developed turbulent channel flow is $\mathrm{Re}_{\tau } = \rho _\text {bulk} u_{\tau } h / \mu _\text {wall} = 590$, where $\rho _\text {bulk}$ is the bulk density, $u_{\tau }$ is the wall friction velocity, h is the channel half-height and $\mu _\text {wall}$ is the viscosity on the wall. The domain size is $2\pi h$ in the streamwise direction, x, 2h in the wall-normal direction, y, and $\pi h$ in the spanwise direction, z. The flow is periodic both in the streamwise and spanwise directions and is bounded by solid walls at $y=0$ and 2h. Because of the imposed streamwise periodicity, a source term is added to the streamwise momentum and energy equations to drive the flow at a constant mass flow rate. The Mach number based on bulk velocity and sound speed on the wall is set to 0.2 in order to facilitate a comparison with incompressible DNS results available in the literature.

The grid used in the DNS contains $768 \times 512 \times 768$ points along x, y and z directions, respectively. The total number of points is about 302 million. The grid resolution in wall units is $\varDelta x^+ \approx 4.8$ in the streamwise direction and $\varDelta z^+ \approx 2.4$ in the spanwise direction. In the vertical direction, $\varDelta y^+ \approx 0.23$ on the wall and $\varDelta y^+ \approx 5.4$ at the channel centerline. For the computation, 18 NVIDIA Tesla V100 GPUs are used on the Summit system. Each GPU solves a grid block of $256^3$ points. The second-order time integration scheme [25] is used with a CFL number of 0.4. Explicit filtering is applied at every 5 time steps with a filtering parameter of $\sigma = 0.09$. Note that the finer overall grid resolution used here, relative to that in the speed bump case, and the simpler configuration keeps the simulation stable with less frequent filtering. To ensure full convergence of the time-averaged results, the flow statistics are averaged over $1009 h/u_\text {bulk}$.

Figure 32 depicts the instantaneous snapshots of the turbulent channel flow. As seen here, the flow structures appear very smooth and provide evidence that the minimal filtering applied in the present case is able to provide a solution free of any numerical wiggles. Figure 33 compares our mean streamwise velocity and Reynolds stress component profiles with the data from Moser et al. [11] as well as Vreman and Kuerten [27]. Both of these groups used spectral methods to solve the incompressible flow equations for the turbulent channel flow at the same $\mathrm{Re}_{\tau }$. As seen in the comparisons, our results show an excellent agreement with their results.