Techniques, Tricks, and Algorithms for Efficient GPU-Based Processing of Higher Order Hyperbolic PDEs

Abstract

GPU computing is expected to play an integral part in all modern Exascale supercomputers. It is also expected that higher order Godunov schemes will make up about a significant fraction of the application mix on such supercomputers. It is, therefore, very important to prepare the community of users of higher order schemes for hyperbolic PDEs for this emerging opportunity.

Not every algorithm that is used in the space-time update of the solution of hyperbolic PDEs will take well to GPUs. However, we identify a small core of algorithms that take exceptionally well to GPU computing. Based on an analysis of available options, we have been able to identify weighted essentially non-oscillatory (WENO) algorithms for spatial reconstruction along with arbitrary derivative (ADER) algorithms for time extension followed by a corrector step as the winning three-part algorithmic combination. Even when a winning subset of algorithms has been identified, it is not clear that they will port seamlessly to GPUs. The low data throughput between CPU and GPU, as well as the very small cache sizes on modern GPUs, implies that we have to think through all aspects of the task of porting an application to GPUs. For that reason, this paper identifies the techniques and tricks needed for making a successful port of this very useful class of higher order algorithms to GPUs.

Application codes face a further challenge—the GPU results need to be practically indistinguishable from the CPU results—in order for the legacy knowledge bases embedded in these applications codes to be preserved during the port of GPUs. This requirement often makes a complete code rewrite impossible. For that reason, it is safest to use an approach based on OpenACC directives, so that most of the code remains intact (as long as it was originally well-written). This paper is intended to be a one-stop shop for anyone seeking to make an OpenACC-based port of a higher order Godunov scheme to GPUs.

We focus on three broad and high-impact areas where higher order Godunov schemes are used. The first area is computational fluid dynamics (CFD). The second is computational magnetohydrodynamics (MHD) which has an involution constraint that has to be mimetically preserved. The third is computational electrodynamics (CED) which has involution constraints and also extremely stiff source terms. Together, these three diverse uses of higher order Godunov methodology, cover many of the most important applications areas. In all three cases, we show that the optimal use of algorithms, techniques, and tricks, along with the use of OpenACC, yields superlative speedups on GPUs. As a bonus, we find a most remarkable and desirable result: some higher order schemes, with their larger operations count per zone, show better speedup than lower order schemes on GPUs. In other words, the GPU is an optimal stratagem for overcoming the higher computational complexities of higher order schemes. Several avenues for future improvement have also been identified. A scalability study is presented for a real-world application using GPUs and comparable numbers of high-end multicore CPUs. It is found that GPUs offer a substantial performance benefit over comparable number of CPUs, especially when all the methods designed in this paper are used.

Appendix A

The OpenACC implementation of a HOG scheme in the C/C++ language is described here. The overall structure is identical to that of the Fortran code structure described in Sect. 2.2. The reader who is familiar with C/C++ should be able to pick up all the OpenACC insights from the figures shown in this Appendix. In this Appendix, we also highlight the differences between Fortran and C/C++ syntax and use of the OpenACC directives.

Figure A1 shows the header file containing the preprocessor definitions of the domain sizes “nx, ny, nz” and the function declarations. The data type and array sizes are specified in the function interface as required by the C/C++ syntax, which is different from Fortran. It can be noted that the variable “dt_next” in the function “Update()” is passed by reference as a pointer since this is a scalar output variable. Such special care is not needed in a Fortran subroutine, since all the variables are passed by reference by default. Similar to the Fortran code, the function “U_skinny_to_U()” moves the minimal data from CPU to GPU at the beginning of every timestep.

Fig. A1

Header file “hydro.h”, containing the definitions of the domain sizes “nx, ny, and nz” in the x-, y-, and z-directions, respectively; and also the function declarations, which are essential for a hydrodynamic simulation. This header file is common to all the C/C++ functions and the main function

Figure A2 shows the overall structure of a hydro code in C/C++ , which is analogous to the Fortran structure shown in Fig. 1. The OpenACC directives have a different structure in C/C++ , they start with “#pragma acc” and they do not require an “!$acc end” statement. So, the analogous form for “!$acc end data” in Fig. 1 is not present in Fig. 2.

Fig. A2

Structure of a C/C++ version of the hydro code (in pseudocode format) at second order with OpenACC extensions to allow it to perform well on a GPU. The overall structure is identical to that of the Fortran version. The exceptions are, the arrays dimensions starts from “0” in C and the zone indexing is reversed from “(i, j, k)” of Fortran to “[k][j][i]” of C, to efficiently operate on the row-major array storage format in C. The OpenACC pragmas are identical to the Fortran version, except the C/C +  + version does not need the ending “#pragma acc end” pragma

Figure A3 shows the C/C++ code for the “U_skinny_to_U()” function. This is equivalent to the Fortran code shown in Fig. 2. The OpenACC directives that are outside the triply nested loop are identical in C/C++ and Fortran. The Fortran has the “:” operator for the whole array operation, whereas the C/C++ does not have any special operator for this. Therefore, we need another “for” loop inside the triply nested loop. This inner for-loop will be executed serially. It has been marked with the OpenACC directive “#pragma acc loop seq”.

Fig. A3

Structure of the “U_Skinny_to_U” function. The array “U_skinny” contains the minimum data that should be carried from the host to the device. This subroutine just copies the contents of “U_skinny” to “U (:, :, :, :, 1)”. The OpenACC pragmas are identical to that of the Fortran version

Figure A4 shows the overall structure of a “Limiter()” function in C/C++ . This is analogous to the Fortran subroutine shown in Fig. 3. The OpenACC directives are identical, except for the “!$acc end” directive, which is not needed in C/C++ . Unlike Fortran, the C/C++ does not allow functions to be defined inside another function. So, the “MC_limiter()” function is defined as an inline function, outside the “Limiter()”. The “MC_limiter()” is called inside a triply nested loop serially. So, it is marked with “#pragma acc routine seq”.

Fig. A4

C/C++ version of the application of a simple monotone-centered limiter in the x-direction to all the fluid variables. The inline function “MC_Limiter” is made to accommodate a compression factor which can indeed be changed based on the flow field being limited. The OpenACC directives are identical to that of the Fortran code, except the C version does not need the ending “#pragma acc end” pragma. Also, the inline function must be identified as “acc routine seq” pragma, since the C does not differentiate between a function and a subroutine like Fortran does

Figure A5 shows the overall structure of a “Predictor()” function in C/C++ . This is equivalent to the Fortran subroutine, shown in Fig. 4. The overall structure and the OpenACC directives are identical to the Fortran code. As mentioned before, the C/C++ lacks the “:” operator for array operations. Due to this, “for” loops are used for the array operations, inside the triply nested loop. These are marked with “#pragma acc loop seq” to indicate their serial execution.

Fig. A5

C/C++ version of the predictor step for building the fifth mode, i.e., the mode that contains the time rate of change. The OpenACC directives are identical to that of the Fortran code, except in the C version the inner for-loops need to be marked with the “acc loop seq” pragma indicating their sequential execution. This is not needed in the Fortran version, where the “:” operator is used for such operations

Figure A6 shows the C/C++ code for the “Make_Flux()” function. This is equivalent to the Fortran code shown in Fig. 5. The OpenACC directives are identical to the Fortran code. An additional “#pragma acc loop seq” is used for the inner loop, which performs the operations on each fluid component.

Fig. A6

C/C +  + version of the pseudocode for the subroutine to build the flux along the x-faces. The “U_L” variable stores the value of “U” obtained from the zone which is at the left to the face and the “U_R” variable stores the value of “U” obtained from the zone which is at the right to the face. Thus, the facial values of the conserved variable “U” are stored in “U_L, U_R”. This is then fed into the Riemann solver function to build the flux for the corresponding face. The OpenACC directives are identical to that of the Fortran code, except in the C version the inner for-loops need to be marked with the “acc loop seq” pragma indicating their sequential execution. This is not needed in the Fortran version, where the “:” operator is used for such operations. In addition, it can be noted that the index of small arrays such as “U_L, U_R, flux_x_ptwise” starts with “0” instead of “1”

Figure A7 shows the equivalent C/C++ code for the “Make_dU_dt()” function. This is analogous to the Fortran shown in Fig. 6. The overall structure and the OpenACC directives are identical between the C/C++ and Fortran.

Fig. A7

C/C++ version of the function which builds the “dU_dt” update from the fluxes obtained from the “Make_Flux_X, Make_Flux_Y and Make_Flux_Z” functions. The “:” operator of the Fortran is replaced with a sequential inner for-loop

Figure A8 shows the C/C++ code for the “Update_U_Timestep()” function. Here we perform the time update for the conserved variable “U” and estimate the next timestep, “dt_next”. This is analogous to the Fortran code shown in Fig. 7. In the C/C +  + version, the scalar variable “dt_next” is passed by reference as a pointer. In such a case, the function and the OpenACC directives are unaware of the data structure of the pointer variable. Therefore, to keep it simple, an additional temporary variable “dt1” is used for the reduction operation. The variable “dt1” is copied to GPU at the beginning, using the directive “#pragma acc data copyin(dt1)”. At the end of the triply nested loop, the final data stored in “dt1” is copied back to the CPU, using the OpenACC directive “#pragma acc update host(dt1)”. Then, it is copied to the “dt_next” variable.

Fig. A8

C/C++ version of the function which updates the “U” from the “dU_dt”. In addition to this, here the “U_skinny” is updated from the new zone-centered values of “U”. The OpenACC directives are identical to that of the Fortran code, except in the C version the inner for-loop need to be marked with the “acc loop seq” pragma indicating its sequential execution. Also in this function, the next timestep “dt_next” is estimated, using the CFL number, zone size, and U. The smallest among all the zones is found using the “fmin()” function. To correctly parallelize this reduction operation, the “reduction(min: dt1)” pragma is used