利用GPU对OpenFOAM解算器PCG进行加速的方法技术

本发明专利技术公开了一种利用GPU对OpenFOAM解算器PCG进行加速的方法，包括：1)将OpenFOAM网格化所生成的格式为ldu的稀疏矩阵转化为更加适合于GPU运算的ELLB稀疏矩阵存储格式；2)将格式转换后的稀疏矩阵与向量进行乘法，根据不同的存储样式采用不同的相乘方法。3)根据PCG解算器的原理，对涉及到稀疏矩阵向量乘法的运算时引入ELLB稀疏矩阵存储格式和相乘方法。本发明专利技术的优点是：充分利用GPU的硬件特性，通过引入ELLB稀疏矩阵存储格式，提高该格式下稀疏矩阵乘法的效率，从而提高OpenFOAM中里利用解算器PCG求解实际问题的效率。PCG求解实际问题的效率。PCG求解实际问题的效率。

【技术实现步骤摘要】
利用GPU对OpenFOAM解算器PCG进行加速的方法


[0001]本专利技术涉及计算机软件及计算机应用
，特别涉及一种利用GPU对CFD开源软件OpenFOAM解算器PCG进行加速的方法。

技术介绍

[0002]计算流体动力学(CFD)相较于传统的理论流体学以及实验流体学，有着实验环境搭建相对简单、计算耗时较低以及能够验证现实中无法试验的理想或极端条件的优势。目前CFD的计算主要是通过相关的软件框架进行，OpenFOAM是一个常用CFD的一个开源软件，有着源码开放，能自定义解算器等优点。求解CFD问题关键在于求解所构建的微分形式的控制方程，此类方程往往不能求出其解析解，只能通过数值离散的方法，转化为线性方程组Ax＝b的求解。而要获得对应的线性方程组，必须先利用网格划分计算域，网格的数目决定着求解的精度。要得到更高的精度，需要划分更多的网格，其计算量也越来越多，所耗费的时间也更长。为了提高线性方程组的求解速度，利用GPU进行并行加速。
[0003]OpenFOAM中常见的求解线性方程组的方法是预条件共轭梯度算法(PCG)。PCG算法主要由稀疏矩阵向量乘法(SpMV)、向量间的运算以及预条件子的求解等几个运算组成。其中，SpMV运算占是整体PCG算法的很大部分开销。OpenFOAM中的稀疏矩阵存储格式为ldu格式，不利于GPU进行并行运算，因此需要对矩阵进行格式转换。常见的稀疏矩阵存储格式包括：CSR、ELL、HBY等，这些矩阵相较于OpenFOAM中原有的矩阵存储格式更适合于GPU运算，但没有充分利用GPU的硬件特性。因此，通过利用GPU并行优化SpMV从而提高OpenFOAM中PCG算法求解网格化所构成的线性方程组的速度。
[0004]黄敏[黄敏,丁萍,罗海飚.共轭梯度法在GPU及Xeon Phi下的并行优化及比较[J].华南理工大学学报(自然科学版),2015,43(11):35
‑
46+53.]等人对CSR稀疏矩阵存储方式进行改良，采用非零元素个数升序排序的方法，使GPU中邻近wrap的数据处理量尽可能接近，避免计算资源的浪费和解决负载不平衡的问题，进而提高SpMV的效率，从而提高PCG的效率。
[0005]改进后的CSR矩阵在GPU中进行并行计算时，可能由于矩阵每行元素存储数目的不一致而导致没有充分利用GPU的计算资源。同时，该方法所优化的PCG算法仅验证了解决公共数据集的问题，没有进行实际问题求解的效果验证。

技术实现思路

[0006]本专利技术针对现有技术的缺陷，提供了一种利用GPU对OpenFOAM解算器PCG进行加速的方法，提高OpenFOAM求解CFD问题的效率。
[0007]为了实现以上专利技术目的，本专利技术采取的技术方案如下：
[0008]一种利用GPU对OpenFOAM解算器PCG进行加速的方法，包括以下步骤：
[0009]S1：将OpenFOAM所离散后的线性方程组Ax＝b的矩阵A由ldu格式先转化为ELLPACK格式，再转化:ELLB格式。
[0010]S2：对value以及col_idx进行分块处理，分为长为ω，分块后子块的矩阵宽度为σ的矩阵块进行存储。同时，记录每行元素的数目，并进行归并求和，作为e_num数组，用于求解辅助矩阵。
[0011]S3：为了区分每一个块的起始位置，引入tile_ptr向量，该向量代表每一个块的第一个元素所在的行的标号。同时，引入辅助向量bit_flag确定每一行开始的元素。
[0012]S4：根据整个矩阵中的非零元素数目nnz及矩阵的整体行数数目row_num求出每行的元素平均数目avg：
[0013][0014]根据avg的大小，确定分块后子块的宽度大小为σ：
[0015][0016]提高读取数据时的数据利用率，将分块后的子块长度ω设为32。将分块后的矩阵按列进行存储。
[0017]S5：对于块矩阵中的辅助向量bit_flag，采用二进制的方式进行存储，缩减所需要的内存空间，用一个unsigned int元素代表一个矩阵块的行开始数。根据矩阵的每行非零元素数目e_num与矩阵宽度σ的元素的大小，求出每一行非零元素在分块后所占的行数，根据行数求出每一个矩阵块的bit_flag。
[0018]S6：对于整体辅助向量tile_ptr，根据原本向量e_num的值，进行二分查找，进而求出每一个矩阵块所对应的tile_ptr的值。
[0019]S7：基于ELLB的SpMV并行优化；
[0020]对于ELLB稀疏矩阵存储格式，对于每一行元素，分配一个线程进行处理。通过设置矩阵宽度σ，使得每个线程计算所需要的开销是一致的。
[0021]SpMV的步骤如下：
[0022]计算每一行数据的元素和；
[0023]根据bit_flag向量针对每一块的值进行处理；
[0024]根据tile_ptr的值求解最终的y值。
[0025]S8：将并行优化的SpMV算法集成至OpenFOAM的PCG算法中；
[0026]进一步地，S7中在对每一个块的每一行进行求和以后，所得到的每一块的数据根据bit_flag的数据可以分为两种类别：
[0027]一种是块中的bit_flag数组中的每一个位的值全为0或首位为1且其他位为零。一种是块中的bit_flag数组中的包含多个值为1的位，或包含一个值为1的位但位置不是bit_flag的首位。
[0028]对于第一种类别，在对块中的一行每一个元素进行求和后，对整块所求的的和进行归并处理，求出整体和sum。再根据bit_flag的首位元素是否是1以及tile_ptr将整体求和的结果sum加至对应的y元素上。
[0029]对于第二种类别，也是先对块中的每一行元素进行求和，再根据块中所对应的bit_flag元素以及tile_ptr元素，将对应的和加至对应的y上。
的拥有25个非零元素的矩阵为例，从ELLPACK转化为ELLB存储格式如图2所示。
[0051]为了提高该算法GPU并行化的效率，为矩阵分块的σ以及ω提供一种自动化的选择方案，以满足SIMT的计算模式。
[0052]首先，对于矩阵宽度σ，该值与每行的平均非零元素有关。根据整个矩阵中的非零元素数目nnz及整体行数m求出每行的元素平均数目avg：
[0053][0054]根据avg的大小，我们确定每一个块行数的大小σ：
[0055][0056]对于矩阵高度与GPU的硬件结构有关，CUDA中的最小的逻辑单位是thread，而最小硬件执行单位是wrap，一个wrap包含32个thread，GPU在读取数据时是以wrap为单位读取的，因此，为了提高读取数据时的数据利用率，将矩阵高度ω设为32。为了满足SIMT的计算模式，将分块后的矩阵按列进行存储，这样在每一个wrap读取数据时，可以读取相邻内存的数据，减少数据的重复读取，提高运算效率。
[0057]对于块矩阵中的辅助向量bit_本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种利用GPU对OpenFOAM解算器PCG进行加速的方法，其特征在于，包括以下步骤：S1：将OpenFOAM所离散后的线性方程组Ax＝b的矩阵A由ldu格式先转化为ELLPACK格式，再转化:ELLB格式；S2：对value以及col_idx进行分块处理，分为长为ω，分块后子块的矩阵宽度为σ的矩阵块进行存储；同时，记录每行元素的数目，并进行归并求和，作为e_num数组，用于求解辅助矩阵；S3：为了区分每一个块的起始位置，引入tile_ptr向量，该向量代表每一个块的第一个元素所在的行的标号；同时，引入辅助向量bit_flag确定每一行开始的元素；S4：根据整个矩阵中的非零元素数目nnz及矩阵的整体行数数目row_num求出每行的元素平均数目avg：根据avg的大小，确定分块后子块的宽度大小为σ：提高读取数据时的数据利用率，将分块后的子块长度ω设为32；将分块后的矩阵按列进行存储；S5：对于块矩阵中的辅助向量bit_flag，采用二进制的方式进行存储，缩减所需要的内存空间，用一个unsignedint元素代表一个矩阵块的行开始数；根据矩阵的每行非零元素数目e_num与矩阵宽度σ的元素的大小，求出每一行非零元素在分块后所占的行数，根据行数求出每一个矩阵块的bit_flag；S6：对于整体辅助向量tile_ptr，根据原本向量e_num的值，进行二分查找，进而求出每一个矩阵块所对应的tile_ptr的值；S7：基于ELLB的SpMV并行优化；对于ELLB稀疏矩阵存储格式，对于每一行元素，分配一个线程进行处理；通过设置矩阵宽度σ，使得每个线程计算所需要的开销是一致的；SpMV的步骤如下：计算每一行数据的元素和；根据bit_flag向量针对每一块的值进行处理；根据tile_ptr的值求解最终的y值；S8：将并行优化的SpMV算法集...

【专利技术属性】
技术研发人员：王卓薇，陈浩楠，程良伦，
申请(专利权)人：广东能哥知识科技有限公司，
类型：发明
国别省市：