Efficient GPU Computing Method for CFD without Synchronous ILU Preconditioner

Abstract

The present invention discloses a CFD high - efficiency GPU computing method without a synchronous ILU pre - conditioner, belonging to the field of high - performance parallel computing, including the steps of: taking the rows of a 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; determining the index of each thread to process the task row and finding the diagonal position of the row; using a busy - waiting mechanism to judge the change of the dependency relationship between block - row groups and performing an accumulation operation on each row; solving the dependency relationship within the block and passing the latest value to the threads within the block; completing the calculation of all corresponding solution vectors x of the block - row group, setting the dependency relationship of the block - row group as resolved, and giving a notification. The present invention improves the solution efficiency, reduces the synchronization overhead, enhances the data access performance, and can adapt to the complex computing requirements of large - scale unstructured grids.

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 synchronous-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 technology, especially the incomplete LU (ILU) preconditioner, is 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 technology has been widely used to improve the efficiency of linear solvers, most existing research and optimizations focus 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 synchronization-free 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 CFD efficient GPU computing method with a synchronization-free ILU preconditioner, which improves the solving efficiency, reduces the synchronization overhead, enhances the data access performance, and can adapt to the complex computational requirements of large-scale unstructured grids.

[0008] The purpose of the present invention is achieved through the following solutions:

[0009] A CFD efficient GPU computing method with a synchronization-free ILU preconditioner, comprising the following steps:

[0010] S1. Use the rows with the block size as block row groups, the block row groups as the parallel granularity, and initialize the dependency relationships for the block row groups;

[0011] S2. Determine the index of the task row processed by each thread and find the diagonal position of the row;

[0012] S3. Use the busy waiting mechanism to determine the change in the dependency relationship between block row groups, and perform an accumulation operation on each row;

[0013] S4. Solve the dependency relationship within the block and pass the latest value to the threads within the block;

[0014] S5. Complete the calculation of all corresponding solution vectors x for the block row group, set the dependency relationship of the block row group to solved, and give a notification.

[0015] Further, in step S1, taking the rows of the block size as block row groups, using the block row groups as the parallel granularity, and initializing the dependency relationship for the block row groups specifically includes the following sub-steps:

[0016] 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, which is initialized to 0 indicating that the dependency has not been processed.

[0017] Further, in step S2, 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:

[0018] 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 according to the row index as the end position.

[0019] Further, in step S3, using the busy waiting mechanism to determine the change in the dependency relationship between block row groups and performing an accumulation operation on each row specifically includes the following sub-steps:

[0020] 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 into the shared memory, and multiply all non-zero elements within the block by the corresponding x value, and perform the accumulation operation in a loop.

[0021] Further, in step S4, solving the dependency relationship within the block and passing the latest value to the threads within the block specifically includes the following sub-steps:

[0022] 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.

[0023] Furthermore, in step S5, after all the corresponding solution vectors x of the block row group are calculated, the dependency relationship of the block row group is set as resolved and a notification is made, which specifically includes the following sub-steps:

[0024] After the calculation of x corresponding to all rows within the block is completed, first use the __threadfence() operation of the GPU 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.

[0025] The beneficial effects of the present invention include:

[0026] (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.

[0027] (2) The present invention adopts strategies such as a reasonable fine-grained task mapping relationship, an efficient memory access mode, and a lightweight dependency control mechanism, improves the solution efficiency, reduces the synchronization overhead, enhances the data access performance, and can fully exert the parallel computing power of the GPU.

[0028] (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

[0029] 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, other drawings can be obtained based on these drawings without creative efforts.

[0030] 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;

[0031] Figure 2 It is an example diagram of block-level parallel execution without synchronization after optimization in the embodiment of the present invention;

[0032] Figure 3 It is a test result diagram of the embodiment of the present invention on multiple matrices, with the benchmark being the official GPU library cuSPARSE;

[0033] Figure 4This is a performance comparison chart of the specific implementation of the solution of the embodiment of the present invention in the unstructured CFD solvers (Flowstar and SU2) for the meshes m6-80w and aircraft. Detailed implementation manners

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

[0035] Term explanations

[0036] The __shfl__sync() function is a synchronization function in CUDA for efficiently exchanging data between threads within the same warp. By controlling the participating threads through a mask, it enables low-latency inter-thread communication.

[0037] The __threadfence() function is a synchronization function in CUDA for ensuring 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.

[0038] The specific implementation process of the present invention is as follows:

[0039] In view of the problems in the background, the present invention provides an efficient CFD GPU computing method with a synchronization-free ILU preconditioner. In the inventive concept, aiming at the problems that there are few optimization algorithms for the sparse triangular solution of the current GPU-based block sparse row (BSR) format, and the existing implementation solutions (such as cuSPARSE) fail to fully utilize the latest parallel computing technologies, restricting their efficiency in large-scale unstructured grid applications, by optimizing the sparse triangular solution operations in the forward and backward substitutions of the preconditioner ILU, and 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.

[0040] In a preferred embodiment, the present invention adopts a strategy of a reasonable fine-grained task mapping relationship, an efficient memory access mode, and a lightweight dependency control mechanism to accelerate the solution of the preconditioner ILU sparse triangular equation based on the GPU, and integrates the optimization method into the unstructured CFD application, which specifically includes the following steps:

[0041] 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, which is initialized to 0 indicating that the dependency has not been processed.

[0042] Step S2: Determine the index of the task line processed by each thread and find the diagonal position of that line; one thread in a warp is responsible for one line of the block row group, and determines the index of the currently processed line 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 that line according to the row index as the end position.

[0043] Step S3: The busy waiting mechanism judges the change of the dependency relationship between block row groups and performs an accumulation operation on each line; 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. Store the new x value in the shared memory, and multiply all non-zero elements in the block by the corresponding x value and perform an accumulation operation in a loop.

[0044] Step S4: Solve the dependencies within the block and pass the latest value to other threads within the block; the first line within the block calculates the corresponding x value according to the previously obtained accumulated sum, and uses the __shfl__sync() function of the GPU to pass this x value to the subsequent lines within the block. The non-zero elements corresponding to this x in the subsequent lines continue to perform the multiply-accumulation operation. When a certain line calculates to the diagonal position, calculate the x value of the corresponding line, and repeat the above operations until all lines within the block are calculated.

[0045] Step S5: Complete the calculation of all corresponding x values of 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 lines within the block, first use the __threadfence() operation of the GPU 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.

[0046] 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 the upper triangular matrix U, and the BSR storage format;

[0047] As Figure 2 shown, it is an example of block-level parallel execution without synchronization after optimization of the embodiment of the present invention;

[0048] As Figure 3 shown, it is the test result of the embodiment of the present invention on multiple matrices, with the GPU official library cuSPARSE as the benchmark;

[0049] As Figure 4 shown, it is the performance comparison of the specific implementation of the embodiment of the present invention in unstructured CFD solvers (Flowstar and SU2) on grids m6-80w and aircraft.

[0050] 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 provided in a processor. Among them, the names of these units do not constitute a limitation to the unit itself in some cases.

[0051] According to one aspect of the embodiments 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, the computer instructions being 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 implementation manners.

[0052] As another aspect, the embodiments of the present invention further provide a computer-readable medium, which may be included in the electronic device described in the above embodiments; or may exist alone without being assembled into the electronic device. The above computer-readable medium carries one or more programs, and when the 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, The steps include the following: S1: Taking the rows of the block size as a block row group, using the block row group as the parallel granularity, and initializing the dependency relationship for the block row group; S2: Determining the index of the task row processed by each thread and finding the diagonal position of the row; S3: Using the busy waiting mechanism to judge the change of the dependency relationship between block row groups and performing an accumulation operation on each row; S4: Solving the dependency relationship within the block and passing the latest value to the threads within the block; S5: Completing the calculation of all corresponding solution vectors x of the block row group, setting the dependency relationship of the block row group as resolved, and giving a notification.

2. The CFD efficient GPU computing method without a synchronous ILU preconditioner according to claim 1, characterized in that, In step S1, the operation of taking the rows of the block size as a block row group, using the block row group as the parallel granularity, and initializing the dependency relationship for the block row group specifically includes the following sub-steps: Forming a row group with consecutive rows of the block size, allocating it to a warp on the GPU for processing, and setting a dependency relationship identifier get_value for each block row group, initializing it to 0 to indicate that the dependency has not been processed.

3. The CFD efficient GPU computing method with a non-synchronous ILU preconditioner according to claim 2, characterized in that In step S2, the operation 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. The diagonal position of this row is found according to the row index and used as the end position.

4. The CFD efficient GPU computing method with a non-synchronous ILU preconditioner according to claim 3, wherein In step S3, the operation of using the busy waiting mechanism to judge the change of the dependency relationship between block row groups and performing an accumulation operation on each row specifically includes the following sub-steps: According to the value of the dependency relationship identifier 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 solution vector x value is stored in the shared memory, and all non-zero elements within this block are multiplied by the corresponding x value, and the accumulation operation is performed in a loop.

5. The CFD efficient GPU computing method with an asynchronous ILU preconditioner according to claim 4, wherein In step S4, the operation 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, 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 this row is calculated, and the calculation is repeated until all rows within the block are calculated.

6. The CFD efficient GPU computing method without a synchronous ILU preconditioner according to claim 5, wherein In step S5, the operation of completing the calculation of all corresponding solution vectors x of the block row group, setting the dependency relationship of the block row group as resolved, 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 relationship identifier get_value of this block row group to 1, indicating that the dependency of this block row group has been resolved.

Images