英语原文共 5 页
基于多GPU的三维直流电阻率模型有限元加速法
Hairil Anwar，Achmad Imam Kistijantoro
（万隆理工学院 电气工程与信息学院，万隆 印度尼西亚）
摘要：本文采用有限元法解决使用多图像处理器（Graphics Processing Unit,GPU）加速三维直流电阻率建模。运用GPU减少有限元法计算过程中最耗时的大型线性方程组的求解。由于共轭梯度解算器常被用于求解大型线性方程组，所以我们利用数据矢量化开发了共轭梯度解算器内核，并以PTX汇编形式编写。我们在GTX 750Ti GPU和特斯拉C2050 GPU上进行了测试，发现在GTX 750 Ti GPU上我们开发的内核比CUSPARSE库具有更好的稀疏矩阵向量性能，但是在特斯拉 C2050GPU上其性能不如CUSPARSE库，这两者与库相比具体的加速比分别为1.4倍和0.7倍，通过使用两个完全相同的GPU，我们发现多GPU实现的加速比是单个GPU的1.9倍，与串行的CPU相比，使用两个GPU可以实现大约10倍的加速比。
关键词：共轭梯度，CUDA，直流电阻率建模，有限元法，多GPU
- 引言
在三维直流电阻率模型中，必须解决由二阶偏微分方程和一组边界条件给出的边值问题有限元法（finite element method，FEM）是一种广泛用于求解偏微分方程的边界值问题的数值方法，同时FEM在电磁仿真中得到了广泛的应用。与有限差分法相比，有限元法对数值解的精度更高。有限元法还提供了域离散化的灵活性。这两种方法都被用来解决大型线性方程组的问题，对于规模较大的线性方程组FEM运行缓慢，因此，减少计算时间是非常重要的，尤其是对于模拟元素数量可以达到10^6阶的大型问题。
GPU是为加速图像计算而设计的硬件。随着目前的发展使执行GPU的通用计算成为可能。由于GPU有着众多核，使其能够并行地执行大型浮点运算，在相同的价格下，GPU提供了比CPU更强的浮点运算能力，与在CPU上实现FEM相比，在GPU上实现FEM的性能有着显著地提高。本文已经在GPU上建立了线性方程组的有限元并行计算方法，同时，对单GPU进行了三维电阻率建模，并且显示了在CPU上实现的加速。
执行一个庞大而复杂的三维直流电阻率模型，需要大内存，而一个GPU的内存有限，通过使用多个GPU，则可以进行大规模的问题建模，提高计算性能。但是在多GPU系统中，由于数据依赖于不同的GPU会导致计算性能的出现问题。
- 直流电阻率模型
在直流电阻率模型中，电阻率介质中电流源的总电势以数值形式确定，用二次势能公式来分析一个点源的直流电阻率模型是由如下式（1）的边界值问题给出的：
（1）
式中，是介质电导率（1/Ohm.m）,u是二次电势（伏特），是初始电势（伏特），为电导率的变化，计算公式如式（2）所示：
（2）
这个问题的边界条件由如下式（3）和式（4）给出：
=0 （3）
（4）
首先，次级势能（u）由具有表面区域的边界条件式（3）和非表面区域边界式（4）的式（1）解得。总电势（V）则通过下式（5）计算：
其中，初电势（）由下式（6）计算得到
本文采用二次电位公式对三维直流电阻率进行建模，有限元法求解了边界值问题，采用四面体元离散化模域，首先使用块元素对模域进行离散化，然后将每个块元素离散为5个四面体单元，而不向块元素中添加任何节点，再利用伽辽金加权余量法得到四面体单元的线性方程组，并且以压缩稀疏行（compress sparse row,CSR）格式的矩阵存储线性方程组，这样便于高效存储和计算。
图一 共轭梯度算法
在有限元法中求解线性方程组是最耗时的计算，共轭梯度法（conjugate gradient，CG）常被用于求解有限元法生成的线性方程组。CG求解器是一种迭代方法，每次迭代需要1个稀疏矩阵-向量乘法（sparse matrix-vector multiplication，SPMV）,2个点积，2个AXPY运算和1个AYPX运算。
- 多GPU方法的实现
通过将CG解算器计算移到GPU来加速三维直流电阻率模型，为了在多GPU系统中实现良好的CG性能，采用了一些实现细节，这些实现的具体方法如下：
- 计算和数据传输之间重叠
- 矩阵存储方式改进
- 使用CUDA事件和非阻塞主机方案的高性能同步
图2 矩阵和向量划分方案
图3 CSR存储方法的改进
矩阵和向量正好适用于要求按行划分和GPU分布的计算，在SPMV操作中，每次CG迭代中每个GPU都需要一个完整的向量p，主机内存从每块GPU中获得被划分的向量p的集合，然后将其复制给每块GPU使其有一个完整的向量p。SPMV操作分两个阶段执行，使得计算和向量p的传输能够重叠。第一阶段SPMV是在A对角线块（A diagonal，Ad）之间进行操作，对每块GPU已经拥有的向量p进行划分。第二阶段SPMV是A的非对角块（A non-diagonal，And）之间的操作，对向量p进行划分，数据传输后归GPU所有，在p的数据传输和SPMV的第一阶段执行重叠。
第一阶段和第二阶段存储不同格式的矩阵A，是为了实现更好的SPMV性能。Ad是以标准CSR格式存储，而And是以修改后的CSR格式存储，它是一种没有被有效地存储在CSR格式中的不完整行的矩阵块，这可以通过添加CSR格式的行信息来排除矩阵中的所有空行，从而有效地存储不完整的行矩阵。
CUDA事件为多GPU系统提供了详细的同步机制，CUDA事件的同步不会阻止主机将工作调度到GPU流中，为了使用这种同步方式，所有的工作都必须在GPU中完成，而主机只负责将工作调度到GPU中。
我们也开发了CG求解器（PTX-v4）所需要的内核，我们的内核在从全局内存中读取数据并存储到GPU全局内存时使用数据向量化。写在PTX组件中的那些内核能够进行数据矢量化，以128位的数据读取和存储在一次处理中相当于4个单精度数据或2个双精度数据。为了使用数据矢量化，必须使用零来填充矩阵和向量，在我们的实现中使用了4的倍数个填充。
4．性能测试
为了验证我们的多GPU实现的性能，我们使用了不同模型大小的三维直流电阻率模型问题，一个用于测试目的的简单双层地球模型问题，在两个不同的系统和配置不同的GPU上进行，第一个系统配备两个GTX 750Ti GPU,第二个系统配备两个特斯拉C2050 GPU。我们开发的内核（PTX-v4）在单精度和双精度方面与CUDA库（CUBLAS和CUSPARSE）性能比较，这是用1个GPU和2个GPU来实施测量对串行CPU实现的加速。
图4 简单双层地球模型
表一 测试模型/问题
表二 测试系统
表三 PTX-v4内核性能与CUBLAS/CUSPARSE 内核的性能
表3显示了在单精度和双精度方面CUDA库性能与PTX-v4内核性能的对比。SPMV内核性能对CG求解器性能影响最大，在GTX 750Ti GPU中PTX-v4 SPMV内核的性能优于CUDA库，单精度可达1.49左右，双精度大约可达1.40，但是这个SPMV内核的性能低于特斯拉C2050 GPU中的CUDA库，在CUDA库中不提供AYPX内核，在使用CUDA库实现时，AYPX内核由CUDA库提供的缩放内核和AXPY内核执行实现。我们可以通过使用更直接的内核来获得更好的性能，PTX-v4 AXPY和DOT内核的性能与CUDA库相当。
图5 A系统的单精度性能比较
图6 A系统的双精度性能比较
图7 B系统的单精度性能比较
图8 B系统的双精度性能比较
图5到图8 显示了使用1个GPU和2个GPU实现CG求解器的CUDA库实现和PTX-v4实现的性能对比，相对于串行CPU实现的性能测量，从图中可以看出，模型的规模对性能的影响较大，规模越大，性能越好。在最大的规模下，2个GPU实现的速度比1个GPU实现的速度快1.9倍，与串行CPU实现相比，使用2个GPU可实现约10倍的加速比。
5. 结论
我们成功地在多GPU系统上实现了有限元法的共轭梯度求解，并且采用有限元法对三维直流电阻率进行了建模。对于规模最大的测试问题，我们采用的两块GPU的多GPU系统实现与使用一块GPU系统实现相比达到了1.9倍的加速比。另外，在单精度和双精度下我们开发的内核PTX-v4和CUDA库（CUBLAS和CUSPARSE）都可以达到这样的加速，而且在GTX 750Ti的GPU上，我们的内核有更好的SPMV性能相对于CUDA库而言，但是在特斯拉C2050的GPU上比CUDA库性能低，这两种不同GPU系统情况下，PTX-v4与CUDA库性能比较分别为大约1.4倍和0.7倍左右，与串行CPU实现相比，使用2个GPU可实现大约10倍的加速。
致谢
第一作者感谢印度尼西亚教育捐赠基金（LPDP）的支持。
参考文献：
[1] Telford, W.M., Geldart, L.P., dan Sheriff R.E. 1990. Applied Geophysics. Cambridge: USA.
[2] C. Polycarpu. Anastasis. 2006. An Introduction to the Finite Element Method in Electromagnetism. Morgan amp; Claypool Publisher : USA.
[3] Wu, Xiaoping, A 3-D finite-element algorithm for DC resistivity modelling using the shifted incomplete Cholesky conjugate gradient method, 2003. Geophys. J. Int. (2003) 154, 947–956.
[4] Mohammad, W. Srigutomo, D. Sutarno and P. Sumintadireja, 'The Modeling of 2D Controlled Source Audio Magnetotelluric (Csamt) Responses Using Finite Element Method,' Journal of Electromagnetic Analysis and Applications, Vol. 4 No. 7, 2012, pp. 293-304.
[5] Kiss, Imre; Gyimothy, S.; Badics, Z.; Pavo, J., 'Parallel Realization of the Element-by-Element FEM Technique by CUDA,' in Magnetics, IEEE Transactions on , vol.48, no.2, pp.507-510, Feb. 2012
[6] CUDA C Programming guide, March 2015.
[7] Dziekonski, A.; Lamecki, A.; Mrozowski, M., 'Jacobi and Gauss-Seidel preconditioned complex conjugate gradient method with GPU acceleration for finite element method,' Microwave Conference (EuMC), 2010 European , vol., no., pp.1305,1308, 28-30 Sept. 2010.
[8] Zimon, J.; Zoworka, M., 'Implementation of selected numerical algorithms for solving sparse matrixes using CUDA technology,' Electrodynamic and Mechatronic System (SELM), 2013 International Symposium on , vol., no., pp.77,78, 15-18 May 2013.
[9] J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics. Antennas, Microwave Circuits and Scattering Applications, ser. IEEE Series on Electromagnetic Wave Theory. Piscataway, NJ: IEEE Press, 1998.
[10] Liu Kun, 'Graphics Processor Unit (GPU) acceleration of Time-Domain Finite Element
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
