CN109726441B - 体和面混合gpu并行的计算电磁学dgtd方法 - Google Patents

Abstract

本发明属于计算电磁学、高性能计算领域，具体为一种基于体和面混合GPU并行的计算电磁学时域间断伽辽金方法。本发明基于图形处理器GPU来进行并行计算，将电场和磁场计算过程分别划分为三个核函数：体积分核函数、面积分核函数、时间更新核函数。对于DGTD时域电磁学问题半离散格式，体积核函数用于求其右端的体右端项，面积核函数用于求其右端的求和项中的面右端项，时间更新核函数将其所有右端项求和，并使用蛙跳时间格式求出关于时间t的微分方程进行时间推进。最终，本发明提供了一种新的DGTD方法求解时域电磁学问题路径，且效果优于现有技术。

Description

体和面混合GPU并行的计算电磁学DGTD方法

技术领域

本发明属于计算电磁学、高性能计算领域，涉及一种时域间断伽辽金方法(DGTD)的并行加速技术，具体为一种基于体和面混合GPU并行的计算电磁学时域间断伽辽金方法。

背景技术

自从上世纪中叶以来，电磁波在军事、科研和日常生活中都得到了广泛的应用，人们对电磁学的研究也在一步步深入。最近几十年，随着技术的发展，电磁学问题已经涉及医学、材料、通信、环境等学科，电磁环境越来越复杂，单纯的理论研究无法得出一个解析解，因此研究分析困难，甚至无法研究，单纯的实验则面临着成本高、周期长等等问题，于是计算电磁学应运而生。

作为一种“数字实验”，计算电磁学使用计算机，利用电磁学的理论进行数值模拟，具有低成本、短周期的特点，能够给涉及电磁学的设计和研究提供非常多的帮助。而在计算电磁学领域中，时域间断伽辽金法(DGTD)是一种求解时域电磁学问题的新兴方法。

DGTD求解无源时域电磁学问题时，可以得到半离散格式为

公式(1)中的

分别为第i单元磁场和电场的质量矩阵，

为第i单元的刚度矩阵，

分别为第i单元的按照基函数展开的电场和磁场矢量的展开系数，

为第i单元的第f个面的通量矩阵，n_if为第i单元的第f个面的单位外法向量，

和

分别为此面的数值通量，且有

公式(2)中{A}_il＝A_i+A_l，

和

分别为第i单元的阻抗和导纳，下标l为与第i单元相邻的单元的编号，α∈[0,1]是一个可变参数，当α＝0时为中心通量，当α＝1时为迎风通量。

该方法源于时域有限体积法和有限元法，且相比于这两种方法，DGTD有更好的灵活性和精度，更重要的是，该方法具有天然的良好并行性，可以较为容易地与并行计算技术相结合。但DGTD相比于有限体积法的一个缺点是每个单元内的计算量更大，因此传统的串行实现的DGTD的计算时间非常长；多核CPU实现的DGTD，一方面加速效果严重受限于CPU的核数，另一方面设备成本很高，因此往往难以实现大规模的并行运算。

近十年来，图形处理器(GPU)并行计算技术获得了高性能计算领域的广泛关注，这在于GPU在有良好的并行性能的同时，还有相对于大规模CPU并行比较低廉的价格和更低的能耗，此外从传统C代码到CUDA C代码的移植比较容易，相对较简单的移植就能够看到不错的加速比。由于GPU所具有的单指令多线程(SIMT)架构能够实现大量重复工作的并行，而这一点和DGTD非常契合。因此，将DGTD与GPU并行相结合，将能够大幅提高DGTD的性能，减少计算的时间成本。在此之前基于GPU并行的DGTD方法中，在整个计算公式(1)以及(2)的流程中均采用以体网格单元为并行计算的基本单元。而由于计算数值通量时不同于计算体积分，需要的是本单元和相邻单元的场数据，线程会对内存进行不可预期的访问，从而影响并行的加速效果。

发明内容

针对上述存在问题或不足，为解决现有现有DGTD求解时域电磁学问题路径单一效果不佳的问题，本发明提供了体和面混合GPU并行的计算电磁学时域间断伽辽金(DGTD)方法。

该方法的具体步骤如下：

S1、在主机端为读取的三角形面网格数组按照边界条件类型进行排序，同一边界条件的面的编号相连，内部面置于最后，之后计算质量矩阵的逆的系数矩阵，大小为d_i*d_i，计算面通量矩阵的系数矩阵，大小为d_f*d_f，d_i为每单元的自由度(DOF)，d_f为每个面的DOF，并对电场数组:E_x、E_y、E_z和磁场数组H_x、H_y、H_z进行初始化，每个场的数组大小为d_i*K，K为单元总数；

S2、为GPU端的电磁场数组，所有单元的体右端项，所有面的面右端项以及电磁场数组分配好GPU全局内存空间，为质量矩阵的逆的系数矩阵、面通量矩阵的系数矩阵分配好GPU常量内存，内存分配具体数量为：体右端项三个分量的数组rhsx、rhsy、rhsz的大小分别为为d_i*K，面右端项又分为一个面在其左侧单元的右端项的三个分量rhsx_l、rhsy_l、rhsz_l和在其右侧单元的右端项的三个分量rhsx_r、rhsy_r、rhsz_r，每个数组的大小为d_i*K_f；质量矩阵的逆矩阵和通量矩阵的系数矩阵大小分别为d_i*d_i，d_f*d_f；K_f为面的总数，d_f为每个面的DOF，并将在主机端初始化完毕的电场、磁场数组从主机内存复制到GPU显存中；

S3、在主机端判断时间步n是否达到仿真步数上限，若是，结束计算，释放所有已分配的GPU显存空间，否则继续；

S4、在GPU端，先计算每个单元的刚度矩阵，每个线程需要计算刚度矩阵的一行，大小为d_i，再计算每个单元的电场体右端项

且线程的分配方式如下：

依据每个线程计算一个DOF的原则分配，每个线程块负责处理N_T个单元，使得一个线程块内的线程数量取：

其中上式中16为CUDA线程束大小的一半，符号

代表A向下取整。本步骤核函数计算结果存储在对应的体右端项三个坐标分量的数组rhsx、rhsy、rhsz中。

S5、在GPU端，先由通量矩阵的系数矩阵计算出每个面的通量矩阵，再计算每个面对其左侧单元、右侧单元的电场面右端项：

和

且线程的分配方式如下：

依据每个线程计算一个DOF的原则分配，每个线程块负责处理N_F个面的左、右右端项，使得一个线程块内的线程数量取：

对于不同的边界条件类型，由于面在S1中排过序，可以通过面的编号判断边界类型并且执行相应计算，需要注意的是，边界面只存在一个相邻的单元，且我们统一将这个单元当作左侧单元，因此边界面只存在左面右端项，本步骤核函数的计算结果根据是对左侧还是右侧单元的贡献而分别存储在左、右面对应的右端项三个坐标分量的数组rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S6、在GPU端，先由质量矩阵的逆的系数矩阵计算每个单元质量矩阵的逆，再组合每个单元电场的体右端项和面右端项为总右端项，并用质量矩阵的逆乘以总右端项得到

最后计算

得到新时刻的电场；本步骤核函数的线程分配与S4一致，计算完毕后将GPU端的电场数组使用cudaMemcpyAsync()函数异步拷贝到主机端进行输出或者后处理；

S7、在GPU端，先计算每个单元的刚度矩阵，每个线程需要计算刚度矩阵的一行，再计算每个单元的磁场体右端项

且线程的分配方式同S4，本步骤核函数计算结果存储在对应的体右端项三个坐标分量的数组rhsx、rhsy、rhsz中；

S8、在GPU端，先由通量矩阵的系数矩阵计算每个面的通量矩阵，再计算每个面对其左侧单元、右侧单元的磁场面右端项：

和

且线程的分配方式同S5，边界条件面的处理同S5，本步骤核函数的计算结果根据是对左侧还是右侧单元的贡献而分别存储在对应的左、右面的右端项三个坐标分量的数组rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S9、在GPU端，先由质量矩阵的逆矩阵的系数矩阵计算每个单元质量矩阵的逆，再组合每个单元磁场的体右端项和面右端项为总右端项，并用刚度矩阵的逆乘总右端项得到

最后计算

得到新时刻的磁场；本步骤核函数的线程分配方式同S4，计算完毕后将GPU端的磁场数组使用cudaMemcpyAsync()异步拷贝到主机端进行输出或者后处理，一个时间步长计算完成，转S3。

在上述步骤中，步骤S4和S7在计算体右端项时，每个线程从全局内存读取一个DOF的场。

在上述步骤中，步骤S6和S9中使用异步拷贝函数cudaMemcpyAsync()能够使得当前不在更新的场数据的拷贝和另一个场的计算同时进行，从而隐藏数据拷贝时间。

在上述步骤中，步骤S5和S8的线程分配是基于面网格的，每个面对其左右两个单元的右端项都有不同的贡献，如果用下标L和R来分别标识一个面的左、右两个单元的物理量，n_LR表示当前计算的面上由左单元指向右单元的单位法向量，而n_RL＝-n_LR表示当前面上由右单元指向左单元的单位法向量，则一个面关于

的迎风通量可以写为

进而面右端项可以写为

将n_RL＝-n_LR带入公式(4)可以得到

其中S_LL,S_RR,S_LR,S_RL分别代表由大小为d_f*d_f扩展为d_i*d_i的面通量矩阵，扩展方式为根据该面的在左右单元中的局部编号，将原通量矩阵的每行每列放到扩展矩阵中该面的自由度在两个单元中所对应的自由度的行列中，下标代表矩阵是相对于左右哪个单元扩展的，第一个下标是对行的扩展，第二个下标是对列的扩展。在这其中，每个线程从全局内存中读取一个面通量矩阵的一行，每个线程先读取左单元中一个DOF的电磁场，在计算完公式(5)中需要使用左单元场的四项之后，读取右单元中一个DOF的电磁场并计算其余四项。这是由于公式(5)中左、右单元的面右端项中的每一个项都被使用了两遍，通过安排计算顺序，可以在使用最小的本地内存的情况下，避免重复读取全局内存中的场数据，同时避免一些重复计算。关于

的面右端项的计算与之同理。

步骤S5和S8中，计算每个面只有等于面DOF数量的线程，而每个面对应的右端项的数量等于体DOF数量，分配等于面DOF数量的线程是因为只有这些DOF才会对单元有贡献，其余DOF的贡献为零，而分配每个面右端项数组的数量等于体DOF数量是为了在S6和S9步骤中更方便地组合右端项。

由于体DOF数大于面DOF数，需要面与DOF对应关系的数组FACE_MAP，来根据一个面在单元内的局部编号，来确定哪些DOF需要分配到各个线程来计算。具体而言，由于面的DOF一一对应于体的DOF，所以当一个线程所计算的面的DOF对其所在体的DOF有贡献时进行计算，没有贡献时置零。比如对于一阶情况而言，0号面的三个DOF对应于所在体中编号1,2,3的DOF，所以1,2,3号DOF需要被S5和S8中每个面对应的3个线程计算，而0号DOF不用计算直接置零。

步骤S6和S9中，计算一个单元的总右端项时，需要用到四面体的拓扑关系，即围成该四面体的四个面的编号以及此单元是一个面的左侧单元还是右侧单元。每个线程根据该线程所计算的单元编号和DOF，从全局内存获取一个体右端项

再根据四个面编号和DOF，从全局内存中获取一个

和一个

同时用标识符F[0]＝1来标明此单元是当前计算的面的左单元，否则F[0]＝0，F[1]＝1来标明此单元是当前计算的面的右单元，否则F[1]＝0，再计算

由于F[0]和F[1]中有且仅有一个为1一个为0，所以上式可以将正确的面右端项整合到总右端项上，如此重复4次，每个面右端项都组合到了总右端项中。

步骤S4到S9中，需要用到单元或者面的矩阵，在本发明的线程分配方式下，每个单元或者面只需要处理一个DOF，因此也只需要矩阵的一行，所以在计算时，不使用共享内存存储矩阵，每个线程只计算并存储本单元DOF对应的一行矩阵，这样可以减少不必要的计算量，并且完全避免使用共享内存时可能存在的存储体冲突。比如某一线程的局部编号是2，则计算第2个DOF的值，那么该线程只需要矩阵第2行的值。

本发明基于图形处理器(GPU)来进行并行计算，并且采用英伟达统一设备计算架构CUDA实现，将电场和磁场计算过程分别划分为三个核函数：体积分核函数、面积分核函数、时间更新核函数。对于DGTD时域电磁学问题半离散格式，体积核函数用于其右端的第一项(体右端项)，面积核函数用于求其右端的求和项中的子项(面右端项)，时间更新核函数将其所有右端项求和，并使用蛙跳时间格式求出关于时间t的微分方程进行时间推进。

综上所述，本发明提供了一种新的DGTD方法求解时域电磁学问题路径，且效果优于现有技术。

附图说明

图1为本发明算法整体流程图；

图2为面积分核函数依据面编号来区分边界类型流程图；

图3为按照体单元并行的线程块、线程分布以及对应的数据在全局内存中的布置方式、线程访问内存示意图；

图4为按照面并行的线程块、线程分布以及对应的数据在全局内存中的布置方式、线程访问内存示意图；

图5为实施例计算500步时的电场幅值云图。

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

实施例中所计算的问题为5×2.5×3.75mm³的矩形金属空腔，划分的有限元网格包含51874个四面体单元，10797个结点，106452个面，只有理想导体(PEC)边界一种边界条件，5408个PEC边界面。基函数采用一阶基函数，即每个单元的DOF为4，每个面的DOF为3，使用TE101模式的场进行初始化。

对于本实施例，需要先从网格获取单元、结点和面的数组，再将面数组按照边界信息排序。此外如附图3，需要给3个体右端项数组和6个场数组各分配51874*4＝207496个大小的空间，如附图4，需要给6个面右端项数组各分配106452*4＝425808个大小的空间。然后为质量矩阵的逆矩阵分配GPU常量内存，大小为4*4＝16，为通量矩阵分配GPU常量内存，大小为3*3＝9，为面与DOF对应的数组FACE_MAP，大小为4*3＝12，常量内存需定义为全局定义，需要加CUDA的__constant__修饰符。此外，需要定义两个CUDA流stream1和stream2，用于实现拷贝和计算的并发执行，可以使用cudaStreamCreate(&p_stream)来创建流对象。

在开始计算前，使用cudaMemcpy()函数将单元、面和结点信息、刚度矩阵、初始化的电场、各种系数矩阵拷贝到GPU端。之后开始时间推进。

对于体积分计算，取每个线程块计算4个单元最合适，理论上应划分为51874/4＝12969个线程块，实际为了通用性，线程块划分为两个维度，第一维固定为128，第二维计算得102，线程划分则取一维共4*4＝16个线程/线程块。在每个线程中，只需要获取刚度矩阵的一行(4个)即可，而电磁场则存在共享内存中，整个线程块共用。

同理，对于面积分计算，每个面DOF为3，每个线程块计算16个面合适，理论上应划分为106452/16＝6654个线程块，同样取二维，第一维固定为128，第二维计算得52个，线程划分为3*16＝48个线程/线程块。在每个线程中，只需要获取通量矩阵的一行(3个)即可，每个线程也只需要一个电磁场的值，因此这里的电磁场不放在共享内存中，由于数值通量会被多个线程用到，因此将计算好的数值通量存在共享内存中。

在面积分核函数中涉及到FACE_MAP矩阵，用于对应面与DOF，其第0行为0号顶点对面上的自由度：1,2,3；第1行为1号顶点对面的自由度：0,2,3；第2行为2号顶点对面的自由度：0,1,3；第3行为3号顶点对面的自由度：0,1,2。矩阵中每一行存在的编号就是需要计算的自由度的编号，不存在的编号就是直接置零的编号。

面积分核函数还涉及边界条件的处理，对于本例，只有PEC一个边界类型，因此面积分只有两种情况需要判断。处理方式是，每个线程在核函数一开始就获取该线程处理的面的编号，如果是在PEC编号范围内，执行PEC边界的计算，否则执行内部面的计算，由于面网格依照边界类型排过序，在这样的方式下，最多有一个线程束可能出现线程分化，比如本例中有5408个边界面，每个面有3个线程，因此共计5408*3＝16224个线程，前507个线程束正好计算完所有的PEC边界面，因此本例中未出现线程分化。

时间更新核函数的线程分配方式同体积分。在计算上，每个线程只需要计算单元质量矩阵的一行(4个)即可，右端项需要存储在共享内存中。在组合右端项时，先获取该单元的体贡献，再对该单元的四个面进行循环，将四个面的面贡献累加到体贡献上来得到总贡献，最后使用蛙跳格式更新当前计算的场。

在上述的计算过程中，使用cudaMemcpyAsync()函数在每一步的电磁场计算完毕后，发起电磁场的异步拷贝，以重叠计算和拷贝时间，减少一定的时间消耗。

图5为实施例计算500步时的电场幅值云图，可以确定计算结果正确无误。经过测试，本例在CPU为Intel Xeon ES 2637v4，GPU为Tesla K40c的计算平台上，采用VS10环境、双精度运算，相比于同等的单CPU串行代码有约35倍加速比。而对于使用二阶基函数、不同网格数量的谐振腔问题或者平面波问题的算例，使用本发明的方法测试中获得了相对于串行代码的48到159倍不等加速比，具有明显的加速效果。相比Stylianos Dosopoulos等人报道的使用双精度运算、单GPU相对单CPU的19倍加速比有一定优势。

Claims (1)

1.体和面混合GPU并行的计算电磁学DGTD方法，具体步骤如下：

S1、在主机端为读取的三角形面网格数组按照边界条件类型进行排序，同一边界条件的面的编号相连，内部面置于最后，之后计算质量矩阵的逆的系数矩阵，大小为d_i*d_i，计算面通量矩阵的系数矩阵，大小为d_f*d_f，d_i为每单元的DOF，d_f为每个面的DOF；并对电场数组:E_x、E_y、E_z和磁场数组H_x、H_y、H_z进行初始化，每个场的数组大小为d_i*K，K为单元总数；其中DOF为自由度；

S2、为GPU端的电磁场数组，所有单元的体右端项，所有面的面右端项以及电磁场数组分配好GPU全局内存空间，为质量矩阵的逆的系数矩阵、面通量矩阵的系数矩阵分配好GPU常量内存，内存分配具体数量为：

体右端项三个分量的数组rhsx、rhsy、rhsz的大小分别为d_i*K，面右端项又分为一个面在其左侧单元的右端项的三个分量rhsx_l、rhsy_l、rhsz_l和在其右侧单元的右端项的三个分量rhsx_r、rhsy_r、rhsz_r，每个数组的大小为d_i*K_f；质量矩阵的逆矩阵和通量矩阵的系数矩阵大小分别为d_i*d_i，d_f*d_f；K_f为面的总数，d_f为每个面的DOF，并将在主机端初始化完毕的电场、磁场数组从主机内存复制到GPU显存中；

S3、在主机端判断时间步n是否达到仿真步数上限，若是，结束计算，释放所有已分配的GPU显存空间，否则继续；

S4、在GPU端，先计算每个单元的刚度矩阵，每个线程需要计算刚度矩阵的一行，大小为d_i，再计算每个单元的电场体右端项

为第i单元的刚度矩阵，

为第i单元按照基函数展开的磁场矢量的展开系数，且线程的分配方式如下：

依据每个线程计算一个DOF的原则分配，每个线程块负责处理N_T个单元，使得一个线程块内的线程数量取：

上式中16为CUDA线程束大小的一半，符号

代表A向下取整，本步骤计算结果存储在对应的体右端项三个坐标分量的数组rhsx、rhsy、rhsz中；

S5、在GPU端，先由通量矩阵的系数矩阵计算出每个面的通量矩阵，再计算每个面对其左侧单元、右侧单元的电场面右端项：

和

v∈{x,y,z}，

为第i单元的通量矩阵，其中L表示面的左侧单元，R表示面的右侧单元，n为单位外法向量，

为数值通量；且线程的分配方式如下：

依据每个线程计算一个DOF的原则分配，每个线程块负责处理N_F个面的左、右右端项，使得一个线程块内的线程数量取：

本步骤的计算结果根据是对左侧还是右侧单元的贡献而分别存储在左、右面对应的右端项三个坐标分量的数组rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S6、在GPU端，先由质量矩阵的逆的系数矩阵计算每个单元质量矩阵的逆，再组合每个单元电场的体右端项和面右端项为总右端项，并用质量矩阵的逆乘以总右端项得到

为第i单元电场的质量矩阵，

为第i单元的第f个面的通量矩阵，n_if为第i单元的第f个面的单位外法向量，

为此面的数值通量；最后计算

得到新时刻的电场；

本步骤核函数的线程分配与S4一致，计算完毕后将GPU端的电场数组使用cudaMemcpyAsync()函数异步拷贝到主机端进行输出或者后处理；

S7、在GPU端，先计算每个单元的刚度矩阵，每个线程需要计算刚度矩阵的一行，再计算每个单元的磁场体右端项

为第i单元按照基函数展开的电场矢量的展开系数，且线程的分配方式同S4，本步骤核函数计算结果存储在对应的体右端项三个坐标分量的数组rhsx、rhsy、rhsz中；

S8、在GPU端，先由通量矩阵的系数矩阵计算每个面的通量矩阵，再计算每个面对其左侧单元、右侧单元的磁场面右端项：

和

v∈{x,y,z}，且线程的分配方式同S5，边界条件面的处理同S5，本步骤核函数的计算结果根据是对左侧还是右侧单元的贡献而分别存储在对应的左、右面的右端项三个坐标分量的数组rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S9、在GPU端，先由质量矩阵的逆矩阵的系数矩阵计算每个单元质量矩阵的逆，再组合每个单元磁场的体右端项和面右端项为总右端项，并用刚度矩阵的逆乘总右端项得到

为第i单元磁场的质量矩阵，

为此面的数值通量，最后计算

得到新时刻的磁场；本步骤核函数的线程分配方式同S4，计算完毕后将GPU端的磁场数组使用cudaMemcpyAsync()异步拷贝到主机端进行输出或者后处理，一个时间步长计算完成，转S3；

所述步骤S4和S7在计算体右端项时，每个线程从全局内存读取一个DOF的场；

所述步骤S5和S8的线程分配是基于面网格的，每个面对其左右两个单元的右端项都有不同的贡献，如果用下标L和R来分别标识一个面的左、右两个单元的物理量，n_LR表示当前计算的面上由左单元指向右单元的单位法向量，而n_RL＝-n_LR表示当前面上由右单元指向左单元的单位法向量，则一个面关于

的迎风通量可以写为

进而面右端项可以写为

将n_RL＝-n_LR带入公式(4)可以得到

其中Z为阻抗，S_LL,S_RR,S_LR,S_RL分别代表由大小为d_f*d_f扩展为d_i*d_i的面通量矩阵，扩展方式为根据该面的在左右单元中的局部编号，将原通量矩阵的每行每列放到扩展矩阵中该面的自由度在两个单元中所对应的自由度的行列中，下标代表矩阵是相对于左右哪个单元扩展的，第一个下标是对行的扩展，第二个下标是对列的扩展；

在这其中，每个线程从全局内存中读取一个面通量矩阵的一行，每个线程先读取左单元中一个DOF的电磁场，在计算完公式(5)中需要使用左单元场的四项之后，读取右单元中一个DOF的电磁场并计算其余四项。

Images