（池州学院机电工程学院，安徽 池州 246306）
摘要：近年来，随着通用图形处理器技术的发展，单个GPU的浮点性能比CPU高一个量级以上。采用GPU实现CFD计算为快速仿真分析提供了有力的途径，但目前GPU上的CFD计算方法尚未成熟，仍有大量的相关数值计算问题有待研究和完善。本文采用CUDA技术，针对NVIDIA Tesla K40c GPU，编写了基于有限体积法的二维CFD求解程序。着重研究了网格排序、存储器对齐、数据存储格式对计算效率的影响。
关键词：网格排序算法；行压缩稀疏矩阵格式；广义极小残量法
1．引言
1999年8月，NVIDIA在其原有的TNT显示芯片的基础上发布了全球首款图像处理器（Graphics Process Unit, GPU）——GeForce 256。作为一款图形处理器，其编程接口支持开源的OpenGL标准以及微软DirectX7.0标准，分担了原来由中央处理器（Central Process Unit, CPU）承担的大量图形处理相关的计算任务。自此以后，在游戏工业以及视景仿真行业等相关产业的推动下，GPU得到了迅猛发展[1]。
相比CPU而言，GPU虽然诞生时间较短，但其发展速度已经超过了CPU，按照摩尔定律，CPU的计算能力每18-24个月翻一倍，而GPU的计算能力几乎每12个月就翻一倍。目前主流GPU的单精度浮点数计算能力是当前CPU的10倍，而其内存带宽也已经达到了CPU的5倍[2]。从热力学上而言，通过在单一CPU上集成更多的晶体管实现CPU性能的提高是有极限的，因为在同一芯片上集成过多的晶体管会导致单位面积的热功率过高，使芯片过热失效。
本文基于CUDA架构讨论了利用GPU加速Roe格式有限体积法求解Euler方程的可能性。
2. CFD及其CUDA计算加速技术
2.1 GPU与CUDA简介
如前所述，自1999年NVIDIA发布第一款GPU以来，GPU经过十几年的发展，其计算性能已经超过了CPU，其使用范围也从传统的图形计算延伸到通用计算当中。
本文使用的NVIDIA（英伟达）Tesla k40图形卡基于Kepler架构，是目前最新单颗GPU性能最高的GPU。因为使用单指令多数据（single instruction multiple data）特性并行计算，使其计算能力大大提升。单精度浮点计算能力达到4290GFLOPs，双精度浮点计算能力达到1430GFLOPs，其计算核心达到2880个，片上显存12GB，内存带宽达到288Gbytes/sec，是高性能科学计算的理想选择。
本文中在GPU上运行的程序代码基于C语言的拓展集-CUDA（Compute Unified Device Architecture:计算统一架构），目前最新版本为v7.5。CUDA允许程序员作为一个独立的设备从CPU使用GPU及其随机访问存储器，该设备能够同时执行以块（Block）组织的大量线程（Thread）。多个线程区块组织在一起形成一个线程块网格（Grid）。所有的线程在GPU上执行内核（Kernel）代码。在硬件水平上，每一个流式多处理器单元可以同时执行超过1024个线程或者32个线程束（warp，每32个并行执行的线程集合）。每一个线程块对应一个流式多处理器单元及其16KB共享内存。由16384个（16KB）可读写的32位寄存器组成的共享内存决定并限制了同一流式多处理单元上线程区块的数量（同一流式多处理单元上最多8个线程区块）。此外，在内核执行期间，同一流式多处理单元上的所有区块的线程均可访问同一共享内存、只读常量内存以及纹理内存空间。
2.2．CFD程序的CUDA计算加速原理和方法
CFD程序中，大量存在着的数组计算大部分存在并行加速的可能。比如数组内积的计算、矩阵与数组乘法的计算，由于数组内积的计算过程除了最后的求和部分之外，之前的各项乘法计算是彼此独立，互不相关的，则与上述数组加法的实现过程类似，在GPU上将原有的两个数组的对应项的乘积计算封装在内核中，并行执行，而将最后的求和过程保留到CPU上计算。类似可以实现稀疏矩阵对数组乘法的计算。上述数组与矩阵计算在我们之后使用GMRES算法求解方程组时需要大量运用。
3. GMRES算法理论简介
对于具体的问题其流程如下：
3.2重启动gmres算法
实际问题当中，系数矩阵特征值没有特定的关系，在特征值较大时，一般GMRES算法的收敛性很差，这意味着，即使 已经很大，子空间迭代深度很高，最终的解仍然无法收敛。随着 的增加，GMRES过程对于计算内存与计算时间的需求也会提高，由于这一原因，当 过大时，GMRES算法已经失去了实用意义。可行的解决方案之一是使用重启动GMRES算法，即给定最大子空间深度，如果已经达到最大子空间深度而当前近似解不满足收敛要求，则重新使用当前近似解为初始解开始子空间迭代直至近似解精度达到收敛要求。
3.3 预条件
如前所述，实际问题当中系数矩阵因为具体的实际问题而不同，其特征值并不会具有特定的一般性的关系，而GMRES算法在系数矩阵特征值很大时，收敛性会变差。因此通常情况下，对于实际问题直接使用GMRES算法，并不能起到良好的加速效果，甚至在矩阵规模较小时，还可能会适得其反。使用重启动GMRES算法，可以通过不断取近似程度高的初始解缓解上述问题，但系数矩阵的特征值并未改变。预条件处理就是对当前方程进行一定条件的变换，从而将原有的大特征值系数矩阵问题转换为一个小特征值系数矩阵问题，通过求解变换后的收敛性好的矩阵方程间接得到原方程的解。预条件处理中所做的变换就是通过某一预条件矩阵将系数矩阵变换为某一具有较小特征值矩阵过程。常见的有左预条件处理，右预条件处理，分解预条件处理[25]。
对于左预条件处理，本质是找到一个预条件子矩阵 使得 的特征值较小，则原方程等价于
（1）．diag预条件
Diag预条件认为实际问题中矩阵对角项，对方程的解的影响较大，对矩阵特征值的贡献起主要作用，因而将系数矩阵的对角项作为预条件子M。
4. 预条件gmres算法的GPU实现
4.1 数据结构与计算流程
由于NVIDIA发布的CUDA库可以看作C语言的一个拓展集，因此使用CUDA基于GPU实现的GMRES算法与上述GMRES的CPU实现整体基本相同，区别之处在于过程中部分的矩阵计算操作可以写作使用CUDA的内核函数实现的并行算法，与此同时需要注意变量在CPU与GPU上的存储以及设备之间数据的复制操作。
4.2 结果对比
图4.1改进后的GMRES过程在CPU与GPU上的运行时间对比以及残差收敛情况
在矩阵规模较小时，CPU上矩阵计算占用较大优势，原因在于此时数据复制到GPU的时间将会占总用时的较大比例，使得GPU计算的优势不明显，甚至落后于CPU。但是当数据量增长时，GPU由于有强大的多线程计算能力，其计算用时增长缓慢，而CPU计算用时则随计算维度增加平方级增长。因此在大数据量的前提下，GPU上的GMRES算法其优势是明显的。
结论
本文主要探讨了基于CUDA库编写的GPU并行CFD算法的设计与实现。使用经过网格排序后的网格进行计算，合理选择稀疏矩阵存储格式以及数据组织方式，并且对所使用的GMRES算法进行改进可以提升程序的运行过程。为此最终实现的二维Euler方程求解器中，我们使用了RCM算法对所使用的非结构网格进行排序，使用CSR格式存储计算中的稀疏矩阵，使用对齐结构体数组存储数据，并在GMRES算法中加入了预条件处理和重启过程。
参考文献
[1]Sanders J, Kandrot E,聂雪军.GPU高性能编程CUDA实战[M].机械工业出版社, 2011.
[2]Nvidia C. Programming guide[J].2008.
[3]刘峰. GPU加速的云的生成和动态模拟[D].浙江大学,2005.
[4]高宏俊.基于流体模型和 GPU 的云的实时仿真[D].中国优秀硕士学位论文全文数据库,2007.
[5]周世哲,满家巨.基于多重网格法的实时流体模拟[J].计算机辅助设计与图形学学报,2007,19(7):935-940.
论文作者:张李平
论文发表刊物:《信息技术时代》2018年3期
论文发表时间:2019/1/2
标签：矩阵论文; 算法论文; 特征值论文; 条件论文; 线程论文; 数组论文; 系数论文; 《信息技术时代》2018年3期论文;