Method of accelerating OpenFOAM solver PCG using GPU

Abstract

The present invention discloses a method for accelerating the OpenFOAM solver PCG by using a GPU, comprising: 1) converting a sparse matrix in the format of ldu generated by OpenFOAM gridding into an ELLB sparse matrix storage format more suitable for GPU operation; 2) multiplying the sparse matrix after the format conversion with a vector, and adopting different multiplication methods according to different storage styles. 3) According to the principle of the PCG solver, the ELLB sparse matrix storage format and multiplication method are introduced for operations involving sparse matrix-vector multiplication. The advantages of the present invention are: making full use of the hardware characteristics of the GPU, by introducing the ELLB sparse matrix storage format, the efficiency of sparse matrix multiplication under the format is improved, thereby improving the efficiency of solving practical problems using the solver PCG in OpenFOAM.

Description

Technical Field

[0001] The present invention relates to the field of computer software and computer application technology, and in particular to a method for accelerating a CFD open source software OpenFOAM solver PCG by using a GPU. Background Art

[0002] Compared with traditional theoretical fluid dynamics and experimental fluid dynamics, computational fluid dynamics (CFD) has the advantages of relatively simple experimental environment construction, low computational time, and the ability to verify ideal or extreme conditions that cannot be tested in reality. At present, CFD calculations are mainly performed through related software frameworks. OpenFOAM is a commonly used open source software for CFD, which has the advantages of open source code and custom solvers. The key to solving CFD problems is to solve the differential form of the constructed control equations. Such equations often cannot be solved analytically, and can only be converted into the solution of the linear equation system Ax=b through numerical discretization methods. In order to obtain the corresponding linear equation system, the computational domain must be divided by grids first, and the number of grids determines the accuracy of the solution. To obtain higher accuracy, more grids need to be divided, and the amount of calculation is increasing, and the time consumed is also longer. In order to improve the speed of solving the linear equation system, GPU is used for parallel acceleration.

[0003] The common method for solving linear equations in OpenFOAM is the preconditioned conjugate gradient algorithm (PCG). The PCG algorithm mainly consists of several operations such as sparse matrix-vector multiplication (SpMV), operations between vectors, and solving preconditioners. Among them, the SpMV operation accounts for a large part of the overall PCG algorithm overhead. The sparse matrix storage format in OpenFOAM is ldu format, which is not conducive to parallel computing of GPU, so the matrix needs to be converted. Common sparse matrix storage formats include: CSR, ELL, HBY, etc. These matrices are more suitable for GPU computing than the original matrix storage format in OpenFOAM, but they do not fully utilize the hardware characteristics of the GPU. Therefore, by using GPU parallel optimization SpMV, the speed of the PCG algorithm in OpenFOAM to solve the linear equations composed of gridding can be improved.

[0004] Huang Min [Huang Min, Ding Ping, Luo Haibiao. Parallel optimization and comparison of conjugate gradient method on GPU and Xeon Phi [J]. Journal of South China University of Technology (Natural Science Edition), 2015, 43(11): 35-46+53.] et al. improved the CSR sparse matrix storage method and adopted the method of ascending sorting of the number of non-zero elements to make the data processing volume of adjacent wraps in the GPU as close as possible, avoiding the waste of computing resources and solving the problem of load imbalance, thereby improving the efficiency of SpMV and thus improving the efficiency of PCG.

[0005] When the improved CSR matrix is ​​parallelized in the GPU, the GPU computing resources may not be fully utilized due to the inconsistency in the number of elements stored in each row of the matrix. At the same time, the PCG algorithm optimized by this method has only been verified for solving the problem of public data sets, and has not been verified for the effect of solving actual problems. Summary of the invention

[0006] In view of the defects of the prior art, the present invention provides a method for accelerating the OpenFOAM solver PCG by using a GPU, thereby improving the efficiency of OpenFOAM in solving CFD problems.

[0007] In order to achieve the above invention object, the technical solution adopted by the present invention is as follows:

[0008] A method for accelerating the OpenFOAM solver PCG using a GPU comprises the following steps:

[0009] S1: Convert the matrix A of the linear equation system Ax=b discretized by OpenFOAM from ldu format to ELLPACK format, and then to ELLB format.

[0010] S2: Divide value and col_idx into blocks of length ω and width σ for storage. At the same time, record the number of elements in each row and merge and sum them as the e_num array for solving the auxiliary matrix.

[0011] S3: In order to distinguish the starting position of each tile, the tile_ptr vector is introduced, which represents the label of the row where the first element of each tile is located. At the same time, the auxiliary vector bit_flag is introduced to determine the element at the beginning of each row.

[0012] S4: Calculate the average number of elements in each row avg based on the number of non-zero elements nnz in the entire matrix and the total number of rows in the matrix row_num:

[0013]

[0014] According to the size of avg, the width of the sub-block after segmentation is determined as σ:

[0015]

[0016] To improve data utilization when reading data, the sub-block length ω after block division is set to 32. The matrix after block division is stored by column.

[0017] S5: For the auxiliary vector bit_flag in the block matrix, it is stored in binary format to reduce the required memory space, and an unsigned int element is used to represent the starting number of a row of a matrix block. According to the number of non-zero elements in each row of the matrix e_num and the size of the elements of the matrix width σ, the number of rows occupied by non-zero elements in each row after block division is calculated, and the bit_flag of each matrix block is calculated according to the number of rows.

[0018] S6: For the overall auxiliary vector tile_ptr, perform a binary search according to the value of the original vector e_num, and then find the value of tile_ptr corresponding to each matrix block.

[0019] S7: Parallel optimization of SpMV based on ELLB;

[0020] For the ELLB sparse matrix storage format, a thread is assigned to process each row of elements. By setting the matrix width σ, the overhead required for calculation by each thread is consistent.

[0021] The steps of SpMV are as follows:

[0022] Calculate the element sum of each row of data;

[0023] Process the value of each block according to the bit_flag vector;

[0024] Solve the final y value based on the value of tile_ptr.

[0025] S8: Integrate the parallel optimized SpMV algorithm into the PCG algorithm of OpenFOAM;

[0026] Furthermore, after summing each row of each block in S7, the data of each block can be divided into two categories according to the data of bit_flag:

[0027] One is that the value of each bit in the bit_flag array in the block is all 0 or the first bit is 1 and the other bits are 0. The other is that the bit_flag array in the block contains multiple bits with a value of 1, or contains a bit with a value of 1 but the position is not the first bit of the bit_flag.

[0028] For the first category, after summing each element in a row of a block, the sums of the entire block are merged to find the overall sum. Then, based on whether the first element of bit_flag is 1 and tile_ptr, the overall sum is added to the corresponding y element.

[0029] For the second category, each row of elements in the block is summed first, and then the corresponding sum is added to the corresponding y according to the corresponding bit_flag element and tile_ptr element in the block.

[0030] Furthermore, the specific steps of S8 are:

[0031] S81: Configure the function interface. The SpMV algorithm inherits the lduMatrix::solver class.

[0032] S82: Generate a dynamic function library, build a cross-framework compilation model, use NVCC in CUDA to compile the EPCG algorithm to generate a dynamic library, and then use WMAKE in OpenFOAM to call the dynamic library to generate a new dynamic library.

[0033] S84: The example calls gPCG. If the optimized PCG needs to be used for solving the example, the controlDict file and the fvSolution file under the system file need to be modified.

[0034] As a preference, based on the CSR sparse matrix storage format, the CSR row_ptr vector is converted into a tile_ptr vector, the original col and val matrices are divided into blocks, and bit_flag is introduced.

[0035] Compared with the prior art, the advantages of the present invention are:

[0036] (1) A matrix sparse storage format ELLB is proposed to effectively improve the efficiency of solving SpMV.

[0037] (2) The parallel accelerated SpMV algorithm is integrated into the PCG solver of OpenFOAM to improve the speed of solving CFD problems.

[0038] (3) Make full use of the hardware characteristics of GPU to parallel optimize SpMV and thus improve the running speed of the PCG algorithm in OpenFOAM. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 ELLPACK sparse matrix storage format diagram of an embodiment of the present invention;

[0040] Figure 2 It is a diagram of an ELLB sparse matrix storage format according to an embodiment of the present invention;

[0041] Figure 3 It is a PCG algorithm execution relationship diagram of an embodiment of the present invention;

[0042] Figure 4 This is a schematic diagram of the dynamic library compilation relationship of an embodiment of the present invention;

[0043] Figure 5 Schematic diagram of heat conduction calculation example of an embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0045] A method for accelerating the OpenFOAM solver PCG using a GPU, comprising:

[0046] 1. The parallel optimization of PCG is mainly to parallelize sparse matrix-vector multiplication. The optimization can be carried out from the following aspects:

[0047] Matrix format conversion module

[0048] The matrix in OpenFOAM is stored in LDU format. For this storage method, the operation of reading a row of non-zero data with k numbers and the vector inner product in the matrix requires 5k random address accesses. In order to improve the efficiency of data transmission, the present invention converts the matrix A of the linear equation system Ax=b discretized by OpenFOAM from ldu format to ELLPACK format ( Figure 1 ), and then converted into the ELLB format which is more suitable for GPU computing.

[0049] Based on the ELLPACK format, this format sets the number of elements in each row in ELLPACK to a multiple of σ. At the same time, value and col_idx are divided into blocks and stored in matrix blocks with a length of ω and a width of σ to improve the efficiency of sparse matrix storage. At the same time, the number of elements in each row is recorded and merged and summed as the e_num array for solving the auxiliary matrix.

[0050] In order to distinguish the starting position of each tile, the tile_ptr vector is introduced, which represents the row number of the first element of each tile. At the same time, the auxiliary vector bit_flag is introduced to determine the element at the beginning of each row. Take an 8×8 matrix with 25 non-zero elements as an example, the conversion from ELLPACK to ELLB storage format is as follows Figure 2 shown.

[0051] In order to improve the efficiency of GPU parallelization of the algorithm, an automatic selection scheme is provided for the matrix blocks σ and ω to meet the SIMT computing mode.

[0052] First, for the matrix width σ, this value is related to the average number of non-zero elements in each row. Based on the number of non-zero elements nnz in the entire matrix and the total number of rows m, the average number of elements in each row avg is calculated:

[0053]

[0054] According to the size of avg, we determine the size of each block row number σ:

[0055]

[0056] The height of the matrix is ​​related to the hardware structure of the GPU. The smallest logical unit in CUDA is thread, and the smallest hardware execution unit is wrap. A wrap contains 32 threads. The GPU reads data in units of wrap. Therefore, in order to improve the data utilization when reading data, the matrix height ω is set to 32. In order to meet the SIMT calculation mode, the matrix after block is stored by column, so that when each wrap reads data, the data of the adjacent memory can be read, reducing the repeated reading of data and improving the calculation efficiency.

[0057] For the auxiliary vector bit_flag in the block matrix, it is stored in binary format to reduce the required memory space. Since the size of the thread bundle in CUDA is 32, an unsigned int element can be used to represent the starting number of a row of a matrix block. According to the number of non-zero elements in each row of the matrix e_num and the size of the elements of the matrix width σ, the number of rows occupied by non-zero elements in each row after block division is calculated, and the bit_flag of each matrix block is calculated based on this number. For the overall auxiliary vector tile_ptr, a binary search is performed based on the value of the original vector e_num to calculate the value of tile_ptr corresponding to each matrix block.

[0058] Parallel optimization of SpMV based on ELLB

[0059] For the ELLB sparse matrix storage format, since the elements of each row are independent of each other, a thread is assigned to process each row of elements. By setting the width σ of each block matrix, the overhead required for calculation by each thread is consistent. The basic steps of SpMV are as follows:

[0060] (1) Calculate the sum of the elements of each row of data;

[0061] (2) Process the value of each block according to the bit_flag vector;

[0062] (3) Solve the final y value based on the value of tile_ptr.

[0063] After summing each row of each block, the data of each block can be divided into two categories according to the data of bit_flag:

[0064] One is that the value of each bit in the bit_flag array in the block is all 0 or the first bit is 1 and the other bits are 0. The other is that the bit_flag array in the block contains multiple bits with a value of 1, or contains a bit with a value of 1 but the position is not the first bit of the bit_flag.

[0065] For the first case, after summing each element in a row of the block, the sums of the entire block are merged to find the overall sum. Then, based on whether the first element of bit_flag is 1 and tile_ptr, the overall sum is added to the corresponding y element.

[0066] For the second case, each row of elements in the block is summed first, and then the corresponding sum is added to the corresponding y according to the corresponding bit_flag element and tile_ptr element in the block.

[0067] 2. Integrate the parallel optimized SpMV algorithm into the PCG algorithm of OpenFOAM.

[0068] Configuration function interface

[0069] According to the original code of the PCG algorithm in OpenFOAM, the optimized algorithm needs to inherit the lduMatrix::solver class. For the specific structure, see Figure 3 shown.

[0070] Generate dynamic function library

[0071] Since the parallel optimization of PCG uses GPU for expansion acceleration, a cross-framework compilation model is built. The EPCG algorithm is compiled using NVCC in CUDA to generate a dynamic library. The dynamic library is then called using WMAKE in OpenFOAM to generate a new dynamic library. The specific structure is shown in Figure 4 shown.

[0072] Example calling gPCG

[0073] If the optimized PCG needs to be used for solving the example, the controlDict file and the fvSolution file under the system file need to be modified.

[0074] Add the library file import to the controlDict file:

[0075] libs: (

[0077] "libexternalgsolvers.so" / / Generated by wmake

[0078] "libgslovers.so" / / Generated by nvcc )

[0080] Modify the specified solver in the fvSolution file. For example, if the solution of variable P uses the PCG algorithm, modify it to:

[0081] solvers

[0082] {

[0083] P

[0084] {

[0085] solver gPCG; / / Use GPU to optimize the PCG algorithm

[0086] …

[0087] }

[0088] }

[0089] This embodiment uses the parallelized PCG solver to solve the heat conduction problem. The heat flux on the lower surface does not change with time, while the heat flux on the upper surface changes with time. Its geometric model is as follows: Figure 5 As shown;

[0090] The model has a length of 1 mm, a width of 1 mm, and a height of 0.025 mm. There is a time-varying heat flux in the upper surface ring. The LaplacianFOAM in OpenFOAM is used for simulation calculation, and the solver used in LaplacianFOAM is PCG.

[0091] Experiments show that compared with the serial PCG solver and the parallel PCG solver using CSR as the matrix storage format, the performance of this solver is improved by an average of 10.1 times and 2.5 times, and the finer the mesh division, the more obvious the performance improvement. Figure 5 The performance of the geometric model in the example is improved by 12.4 times and 3.2 times after the mesh is divided into 100*100*100.

[0092] The method according to the present invention described above can be implemented in hardware, firmware, or as software or computer code that can be stored in a recording medium (such as a CD ROM, RAM, floppy disk, hard disk or magneto-optical disk), or as a computer code that is originally stored in a remote recording medium or a non-temporary machine-readable medium downloaded through a network and will be stored in a local recording medium, so that the method described herein can be stored in such software processing on a recording medium using a general-purpose computer, a special-purpose processor or programmable or special-purpose hardware (such as an ASIC or FPGA). It can be understood that a computer, a processor, a microprocessor controller or programmable hardware includes a storage component (e.g., RAM, ROM, flash memory, etc.) that can store or receive software or computer code, and when the software or computer code is accessed and executed by a computer, a processor or hardware, the processing method described herein is implemented. In addition, when a general-purpose computer accesses the code for implementing the processing shown here, the execution of the code converts the general-purpose computer into a special-purpose computer for executing the processing shown here.

[0093] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the implementation methods of the present invention, and should be understood that the protection scope of the present invention is not limited to such special statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A method for accelerating the OpenFOAM solver PCG using a GPU, characterized in that: The following steps are involved: S1: Convert the matrix A of the linear equation system Ax=b discretized by OpenFOAM from ldu format to ELLPACK format, and then to ELLB format; S2: Divide value and col_idx into blocks, and store them in matrix blocks with a length of ω and a matrix width of σ. At the same time, record the number of elements in each row and merge and sum them as the e_num array for solving the auxiliary matrix. S3: In order to distinguish the starting position of each block, the tile_ptr vector is introduced, which represents the label of the row where the first element of each block is located; at the same time, the auxiliary vector bit_flag is introduced to determine the element at the beginning of each row; S4: Calculate the average number of elements in each row avg based on the number of non-zero elements nnz in the entire matrix and the total number of rows in the matrix row_num: According to the size of avg, the width of the sub-block after segmentation is determined as σ: To improve the data utilization rate when reading data, the sub-block length ω after block division is set to 32; the matrix after block division is stored by column; S5: For the auxiliary vector bit_flag in the block matrix, it is stored in binary mode to reduce the required memory space, and an unsigned int element is used to represent the starting number of a row of a matrix block; according to the number of non-zero elements e_num in each row of the matrix and the size of the element of the matrix width σ, the number of rows occupied by the non-zero elements in each row after block division is calculated, and the bit_flag of each matrix block is calculated according to the number of rows; S6: For the overall auxiliary vector tile_ptr, perform a binary search based on the value of the original vector e_num, and then find the value of tile_ptr corresponding to each matrix block; S7: Parallel optimization of SpMV based on ELLB; For the ELLB sparse matrix storage format, a thread is assigned to process each row of elements; by setting the matrix width σ, the overhead required for calculation by each thread is consistent; The steps of SpMV are as follows: Calculate the element sum of each row of data; Process the value of each block according to the bit_flag vector; Solve the final y value according to the value of tile_ptr; S8: Integrate the parallel optimized SpMV algorithm into the PCG algorithm of OpenFOAM.

2. The method of accelerating the OpenFOAM solver PCG using a GPU according to claim 1, characterized in that: After summing each row of each block in S7, the data of each block can be divided into two categories according to the data of bit_flag: One is that the value of each bit in the bit_flag array in the block is all 0 or the first bit is 1 and the other bits are zero; the other is that the bit_flag array in the block contains multiple bits with a value of 1, or contains a bit with a value of 1 but the position is not the first bit of the bit_flag; For the first category, after summing each element in a row of the block, merge the sums of the entire block to find the overall sum; then add the overall sum result sum to the corresponding y element according to whether the first element of bit_flag is 1 and tile_ptr; For the second category, each row of elements in the block is summed first, and then the corresponding sum is added to the corresponding y according to the corresponding bit_flag element and tile_ptr element in the block.

3. The method of accelerating the OpenFOAM solver PCG using a GPU according to claim 1, characterized in that: The specific steps of S8 are: S81: Configuration function interface, SpMV algorithm inherits lduMatrix::solver class; S82: Generate a dynamic function library, build a cross-framework compilation model, use NVCC in CUDA to compile the EPCG algorithm to generate a dynamic library, and then use WMAKE in OpenFOAM to call the dynamic library to generate a new dynamic library; S84: The example calls gPCG. If the optimized PCG needs to be used for solving the example, the controlDict file and the fvSolution file under the system file need to be modified.

4. The method of accelerating the OpenFOAM solver PCG using a GPU according to claim 1, characterized in that: Based on the CSR sparse matrix storage format, the CSR row_ptr vector is converted into a tile_ptr vector, the original col and val matrices are divided into blocks, and bit_flag is introduced.

Images