CFD efficient GPU computing method without synchronous ILU precondition

Abstract

The invention discloses a CFD (Computational Fluid Dynamics) efficient GPU (Graphics Processing Unit) computing method without a synchronous ILU (Inter-Link Unit) precondition, which belongs to the field of high-performance parallel computing and comprises the following steps of: taking a block size as a block row group, taking the block row group as a parallel granularity, and initializing a dependency relationship for the block row group; determining the index of each thread processing task line, and finding the diagonal position of the line; a busy waiting mechanism is used for judging the change of the dependency relationship between the block row groups, and accumulation operation is carried out on each row; the dependency relationship in the block is solved, and the latest value is transmitted to the thread in the block; and completing calculation of all corresponding solution vectors x of the block row group, setting the dependency relationship of the block row group to be solved, and performing notification. According to the method, the solving efficiency is improved, the synchronization overhead is reduced, the data access performance is improved, and the method can meet the complex calculation requirement of a large-scale unstructured grid.

Description

Technical Field

[0001] The present invention relates to the field of high-performance parallel computing, and more specifically, to an efficient CFD GPU computing method with a synchronization-free ILU preconditioner. Background Art

[0002] With the development of computational fluid dynamics (CFD) technology, the efficiency and accuracy of solving the compressible fluid dynamics equations have become increasingly important in engineering applications. During the CFD solution process, especially when dealing with unstructured grid problems, the solution of large sparse linear systems is often encountered. To improve the solution speed of these large linear systems, preconditioner techniques, especially the incomplete LU (ILU) preconditioner, are widely used to accelerate the convergence of iterative solvers.

[0003] The core idea of the ILU preconditioner is to perform an approximate LU decomposition on the sparse matrix A, decomposing it into a lower triangular matrix L and an upper triangular matrix U to accelerate the convergence of iterative solutions. This decomposition method does not retain all non-zero elements but reduces fill-in according to different strategies (such as ILU(0) only retaining the non-zero structure of the original matrix) to keep the preconditioner sparse. Then, during the iterative solution process, L and U are used for forward and backward substitution, that is, solving the triangular systems of equations Ly = b and Ux = y respectively.

[0004] Although the ILU preconditioner technique has been widely used to improve the efficiency of linear solvers, existing research and optimization have mostly focused on the decomposition stage, with relatively less attention paid to the solution stage. Especially in unstructured CFD applications, traditional ILU solution methods (such as those based on scalar sparse row formats) exhibit problems such as low efficiency and poor parallel scalability in multi-variable coupling problems. This is because in complex unstructured grids, data dependencies and irregular matrix sparse structures limit the parallel performance of traditional solution methods.

[0005] In addition, although many researchers have proposed parallel optimization schemes for solving sparse triangular systems of equations based on GPUs, many existing methods for solving sparse triangular systems of equations are usually developed based on the CSR (Compressed Sparse Row) format. Such methods perform well in some applications of scalar sparse matrices, but in the multi-physical variable coupling scenarios of CFD, matrices are usually stored in the Block Sparse Row (BSR) format, which combines several consecutive elements of the matrix into matrix blocks for storage. Compared with the scalar sparse row format (CSR), the BSR format is particularly suitable for processing systems with multiple variables at each grid point, such as the multiple physical variables (e.g., velocity, pressure, etc.) corresponding to each grid point in three-dimensional compressible fluid problems. By combining the physical quantities at each grid point into matrix blocks, the BSR format can better utilize data locality and reduce the overhead of storage and memory access.

[0006] Currently, for solving sparse triangular systems of equations in the BSR format, the existing cuSPARSE math library provides corresponding implementations that support forward and backward substitution in the BSR format on GPUs. However, although the cuSPARSE library has a certain influence in the field of GPU computing, its performance in the BSR format is still limited, and it fails to fully utilize the latest technologies and optimization strategies. Especially when dealing with large sparse matrices, the computational efficiency is not ideal. Therefore, there is still much room for improvement and application prospects for further optimization of sparse triangular solving in the BSR format, especially in combination with the latest technologies such as the asynchronous parallel algorithm under heterogeneous computing architectures, which is also an effective method to further improve the performance of unstructured CFD. Summary of the Invention

[0007] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a high-efficiency CFD GPU computing method with an asynchronous ILU preconditioner, which improves the solving efficiency, reduces the synchronization overhead, enhances the data access performance, and can adapt to the complex computing requirements of large-scale unstructured grids.

[0008] The purpose of the present invention is achieved through the following solutions: A high-efficiency CFD GPU computing method with an asynchronous ILU preconditioner includes the following steps: S1. Take the rows of the block size as block row groups, with the block row groups as the parallel granularity, and initialize the dependency relationship for the block row groups; S2. Determine the index of the task row processed by each thread and find the diagonal position of the row; S3. Use the busy waiting mechanism to judge the change of the dependency relationship between block row groups and perform an accumulation operation on each row; S4. Resolve the dependencies within the block and pass the latest value to the threads within the block; S5. Complete the calculation of all corresponding solution vectors x for the block row group, set the dependencies of the block row group as resolved, and give a notification.

[0009] Further, in step S1, using the rows of the block size as the block row group, with the block row group as the parallel granularity, and initializing the dependencies for the block row group, specifically including the following sub-steps: Form a row group with consecutive rows of the block size, allocate it to a warp on the GPU for processing, and set a dependency identifier get_value for each block row group, initialized to 0 indicating that the dependency has not been processed.

[0010] Further, in step S2, determining the index of the task row processed by each thread and finding the diagonal position of the row, specifically including the following sub-steps: A thread in the warp is responsible for one row of the block row group, and determines the index of the currently processed row through the thread coordinates. Among them, the lower triangular matrix L is executed in the order from the first row to the last row, and the upper triangular matrix U is executed in the order from the last row to the first row. Find the diagonal position of this row as the end position.

[0011] Further, in step S3, using the busy waiting mechanism to judge the change of dependencies between block row groups and perform an accumulation operation on each row, specifically including the following sub-steps: According to the value of get_value, if it is 0, it means the dependency is not resolved and perform busy waiting. If it is 1, it means the dependency has been resolved. Store the new solution vector x value in the shared memory, multiply all non-zero elements within this block by the corresponding x value, and perform an accumulation operation in a loop.

[0012] Further, in step S4, resolving the dependencies within the block and passing the latest value to the threads within the block, specifically including the following sub-steps: The first row within the block calculates the corresponding x value according to the accumulated sum obtained in step S3, and passes this x value to the subsequent rows within the block through the operation of the __shfl__sync() function of the GPU. The non-zero elements corresponding to this x in the subsequent rows continue to perform the multiply-accumulation operation. When a certain row calculates to the diagonal position, calculate the x value of the corresponding row, and repeat the calculation until all rows within the block are calculated.

[0013] Further, in step S5, completing the calculation of all corresponding solution vectors x for the block row group, setting the dependencies of the block row group as resolved, and giving a notification, specifically including the following sub-steps: After completing the x calculations for all rows within a block, first use the GPU's __threadfence() operation to ensure that all memory operations are complete, and then set the get_value of this block row group to 1, indicating that the dependencies of this block row group have been resolved.

[0014] The beneficial effects of the present invention include: (1) The present invention designs a preconditioner ILU based on GPU, optimizes the sparse triangular solution operations in the forward and backward substitutions of ILU, and adapts to the complex calculation requirements of large-scale unstructured grids.

[0015] (2) The present invention adopts strategies such as reasonable fine-grained task mapping relationships, efficient memory access patterns, and lightweight dependency control mechanisms, improving the solution efficiency, reducing the synchronization overhead, enhancing the data access performance, and being able to fully utilize the parallel computing power of the GPU.

[0016] (3) The present invention is based on the solution of sparse triangular equations in BSR format and performs better than the existing official GPU library cuSPARSE in a large number of matrix tests. As a preconditioner for unstructured CFD, the present invention significantly improves the convergence efficiency of the calculation. Description of the Drawings

[0017] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0018] Figure 1 It is a simple example diagram of the ILU process of the embodiment of the present invention, including the original matrix structure, the decomposed lower triangular matrix L and upper triangular matrix U, and the BSR storage format; Figure 2 It is a diagram of the optimized block-level parallel execution without synchronization in the embodiment of the present invention; Figure 3 It is a diagram of the test results on multiple matrices in the embodiment of the present invention, with the benchmark being the official GPU library cuSPARSE; Figure 4 It is a diagram of the performance comparison of the specific implementation of the embodiment of the present invention in the unstructured CFD solvers (flowstar and SU2) for grids m6-80w and aircraft. Detailed Implementation Manner

[0019] All features disclosed in all embodiments in this specification, or steps in all methods or processes implicitly disclosed, except for mutually exclusive features and / or steps, can be combined and / or extended, replaced in any manner.

[0020] Term Explanation The __shfl__sync() function is a synchronization function in CUDA used for efficient data exchange between threads within the same warp. It controls the participating threads through a mask to achieve low-latency inter-thread communication.

[0021] The __threadfence() function is a synchronization function in CUDA used to ensure that the writes of the current thread to global or shared memory are visible to all other threads, and is used to achieve memory consistency within the device scope.

[0022] The specific implementation process of the present invention is as follows: In view of the problems in the background, the present invention provides an efficient GPU computing method for CFD with a synchronization-free ILU preconditioner. In the inventive concept, aiming at the problems that there are few optimization algorithms for the current GPU-based block sparse row (BSR) format sparse triangular solution, and the existing implementation schemes (such as cuSPARSE) fail to fully utilize the latest parallel computing technologies, which limit their efficiency in large-scale unstructured grid applications, by optimizing the sparse triangular solution operations in the preconditioner ILU forward and backward substitutions, using the CPU-GPU heterogeneous architecture and the synchronization-free parallel strategy, the purpose of improving the solution efficiency, reducing the synchronization overhead, enhancing the data access performance, and adapting to the complex computing requirements of large-scale unstructured grids is achieved.

[0023] In a preferred embodiment, the present invention adopts a strategy of reasonable fine-grained task mapping relationship, efficient memory access mode, and lightweight dependency control mechanism to accelerate the solution of the preconditioner ILU sparse triangular equation based on GPU, and integrates the optimization method into the unstructured CFD application, specifically including the following steps: Step S1: Take the rows of the block size (block row group) as the parallel granularity, and initialize the dependency relationship for the block row group; form a row group with consecutive rows of the block size, allocate it to a warp on the GPU for processing, and set a dependency relationship identifier get_value for each block row group, initialized to 0 indicating that the dependency has not been processed.

[0024] Step S2: Determine the index of the task row processed by each thread, and find the diagonal position of this row; one thread in the warp is responsible for one row of the block row group, and determines the index of the row being processed through the thread coordinates. The lower triangular matrix L is executed in the order from the first row to the last row, and the upper triangular matrix U is executed in the order from the last row to the first row. Find the diagonal position of this row according to the row index as the end position.

[0025] Step S3: The busy-waiting mechanism determines the change in the dependency relationship between block row groups and performs an accumulation operation on each row; according to the value of get_value, if it is 0, it means the dependency is not resolved and busy-waiting is performed. If it is 1, it means the dependency has been resolved. The new x value is stored in the shared memory, and all non-zero elements within the block are multiplied by the corresponding x value, and the accumulation operation is performed in a loop.

[0026] Step S4: Resolve the dependencies within the block and pass the latest value to other threads within the block; the first row within the block calculates the corresponding x value based on the previously obtained accumulation sum, and through the operation of the GPU's __shfl__sync() function, this x value is passed to the subsequent rows within the block. The non-zero elements corresponding to this x in the subsequent rows continue to perform the multiply-accumulation operation. When a certain row calculates to the diagonal position, the x value corresponding to that row is calculated, and the above operations are repeated until all rows within the block are calculated.

[0027] Step S5: Complete the calculation of all corresponding x values for the block row group, set the dependency relationship of the block row group to resolved, and notify other threads; after completing the calculation of all x values corresponding to the rows within the block, first use the GPU's __threadfence() operation to ensure that all memory operations are completed, and then set the get_value of this block row group to 1, indicating that the dependency of this block row group has been resolved.

[0028] As Figure 1 shown, it is a simple example of the ILU process of the embodiment of the present invention, including the original matrix structure, the decomposed lower triangular matrix L and upper triangular matrix U, and the BSR storage format; As Figure 2 shown, it is an example of block-level parallel execution without synchronization after optimization in the embodiment of the present invention; As Figure 3 shown, it is the test result of the embodiment of the present invention on multiple matrices, with the benchmark being the GPU official library cuSPARSE; As Figure 4 shown, it is the performance comparison of the embodiment of the present invention's solution specifically implemented in unstructured CFD solvers (flowstar and SU2) for grids m6 - 80w and aircraft.

[0029] The units involved in the embodiments of the present invention can be implemented in software or in hardware, and the described units can also be set in a processor. Among them, the names of these units do not constitute a limitation to the unit itself in some cases.

[0030] According to one aspect of an embodiment of the present invention, there is provided a computer program product or a computer program, the computer program product or the computer program including computer instructions stored in a computer-readable storage medium. A processor of a computer device reads the computer instructions from the computer-readable storage medium, and the processor executes the computer instructions, so that the computer device executes the methods provided in the above various alternative implementations.

[0031] As another aspect, an embodiment of the present invention further provides a computer-readable medium, which may be included in the electronic device described in the above embodiment; or may exist alone without being assembled into the electronic device. The above computer-readable medium carries one or more programs, and when the above one or more programs are executed by an electronic device, the electronic device implements the methods described in the above embodiments.

Claims

1. An efficient CFD GPU computing method without a synchronous ILU preconditioner, characterized in that It includes the following steps: S1. Use the rows of the block size as a block row group, with the block row group as the parallel granularity, and initialize the dependency relationship for the block row group; S2. Determine the index of the task row processed by each thread and find the diagonal position of the row; S3. Use the busy waiting mechanism to judge the change of the dependency relationship between block row groups and perform an accumulation operation on each row; S4. Solve the dependency relationship within the block and pass the latest value to the threads within the block; S5. Complete the calculation of all corresponding solution vectors x of the block row group, set the dependency relationship of the block row group to be solved, and give a notification.

2. The CFD efficient GPU computing method with a non-synchronous ILU preconditioner according to claim 1, wherein In step S1, the step of using the rows of the block size as a block row group, with the block row group as the parallel granularity, and initializing the dependency relationship for the block row group specifically includes the following sub-steps: Form a row group with consecutive rows of the block size, allocate it to be processed by a warp on the GPU, and set a dependency flag get_value for each block row group, initialized to 0 indicating that the dependency has not been processed.

3. The CFD efficient GPU computing method with a non-synchronous ILU preconditioner according to claim 2, wherein In step S2, the step of determining the index of the task row processed by each thread and finding the diagonal position of the row specifically includes the following sub-steps: A thread in the warp is responsible for one row of the block row group, and determines the index of the currently processed row through the thread coordinates. Among them, the lower triangular matrix L is executed in the order from the first row to the last row, and the upper triangular matrix U is executed in the order from the last row to the first row. Find the diagonal position of this row as the end position.

4. The CFD efficient GPU computing method without a synchronous ILU preconditioner according to claim 3, characterized in that In step S3, the step of using the busy waiting mechanism to judge the change of the dependency relationship between block row groups and perform an accumulation operation on each row specifically includes the following sub-steps: According to the value of the dependency flag get_value, if it is 0, it means the dependency is not solved and perform busy waiting. If it is 1, it means the dependency has been solved. Store the new solution vector x value into the shared memory, and multiply all non-zero elements within this block by the corresponding x value, and perform an accumulation operation in a loop.

5. The CFD efficient GPU computing method with an asynchronous ILU preconditioner according to claim 4, characterized in that In step S4, the step of solving the dependency relationship within the block and passing the latest value to the threads within the block specifically includes the following sub-steps: The first row within the block calculates the corresponding x value according to the accumulated sum obtained in step S3, and through the operation of the __shfl__sync() function of the GPU, passes this x value to the subsequent rows within the block. The non-zero elements corresponding to this x in the subsequent rows continue to perform multiplication and accumulation operations. When a certain row is calculated to the diagonal position, calculate the x value of the corresponding row, and repeat the calculation until all rows within the block are calculated.

6. The CFD efficient GPU computing method without a synchronous ILU preconditioner according to claim 5, characterized in that In step S5, the step of completing the calculation of all corresponding solution vectors x of the block row group, setting the dependency relationship of the block row group to be solved, and giving a notification specifically includes the following sub-steps: After completing the calculation of x corresponding to all rows within the block, first use the __threadfence() operation of the GPU to ensure that all memory operations are completed, and then set the dependency flag get_value of this block row group to 1, indicating that the dependency of this block row group has been solved.

Images