CN107562690A - 一种基于gpu的aim‑po混合方法 - Google Patents
一种基于gpu的aim‑po混合方法 Download PDFInfo
- Publication number
- CN107562690A CN107562690A CN201710638068.6A CN201710638068A CN107562690A CN 107562690 A CN107562690 A CN 107562690A CN 201710638068 A CN201710638068 A CN 201710638068A CN 107562690 A CN107562690 A CN 107562690A
- Authority
- CN
- China
- Prior art keywords
- node
- matrix
- gpu
- aim
- area
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Withdrawn
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于GPU的AIM‑PO混合算法，包括以下步骤：1.对自适应积分法中的扩展系数矩阵以及近区阻抗矩阵填充进行并行化求解；2.对自适应积分法中的矩阵方程进行并行化求解；3.引入基于GPU的无堆栈KD‑Tree对PO区域进行遮挡面判断；4.对AIM‑PO混合方法中的PO区域电流系数进行并行化求解；5.对附加激励向量进行并行化求解。
Description
技术领域
本发明涉及一种基于GPU的AIM-PO混合方法，属于计算电磁学技术领域。
背景技术
在频率域电磁计算方法中，自适应积分法是对传统矩量法的一种改进，有效的降低了算法的空间复杂度和时间复杂度。在求解超大电尺寸物体的电磁问题时，单一的全波方法(矩量法，自适应积分法等)由于具有较高的复杂度已经难以适用，因此发展出了全波方法和高频近似混合方法：迭代AIM-PO混合算法。
迭代AIM-PO混合算法将天线平台划分为AIM区域和PO区域。一般将天线及其附近区域划分为AIM区域，其他部分为PO区域。在对区域进行划分时，边界上的RWG基函数的两个三角形可能分别属于AIM区域和PO区域，此时，应该将该未知量归于AIM区域来处理，剩余的部分使用PO来近似。对计算目标进行合理的区域划分后，AIM区域由于具有精细的天线结构和较强的表面电流分布，故采用自适应积分法进行计算。
通过物理光学法的相关理论可知，PO区域将目标散射场近似成表面的感应电磁流的积分来计算。在利用迭代AIM-PO算法求解电大尺寸平台的天线问题时，则是将AIM区域的表面电流视为激励源，PO区域被照射后便可以根据PO的近似表达式求出PO区域的表面电流：
通过一系列的数学推导，PO区域的表面电流系数可以表示成：
IPO＝A*IAIM
式中，A*是一个n×k的转换矩阵，IAIM是利用自适应积分法解得的AIM区域的表面电流系数。
根据AIM对PO区域的照射作用得到了PO区域的表面电流IPO，但是PO区域的表面电流对AIM区域同样有不可忽略的扰动作用，这表明整个系统此时并没有达到一个稳定的状态。为了获得一个稳定的状态，需要考虑PO区域的作用，将PO区域的表面电流对AIM区域的扰动作用等效为一个附加激励向量，该附加激励向量可以表示成：
ΔV＝Z*IPO
式中Z*是一个n×k的转换矩阵，IPO是解得的PO区域表面电流系数。
附加激励向量ΔV求解出来之后，便可以将该激励附加到AIM矩阵方程的右端，则新的激励向量可以表示为：
Vn＝Vn-1+ΔV
式中，n≥1，Vn表示第n次迭代时的激励向量。新的激励向量求解出来之后，便可以重新代入到AIM区域的矩阵方程中求解，计算出AIM区域受到PO区域表面电流扰动后的新的表面电流分布。这便完成了一次迭代，但是AIM区域的表面电流的变化又会对PO区域照射后产生新的表面电流。因此，为了使整个系统达到一个相对稳定的状态，需要进行多次迭代，不断地对AIM区域和PO区域的表面电流进行修正，直到两个区域的电流在连续的两次迭代过程中的误差达到一个可以接受的范围内，即满足以下条件：
在上式中，Ii,AIM和Ii-1,AIM分别是第i次和第i-1次迭代求解得到的AIM区域的表面电流系数。当第i次和第i-1次求解得到的电流误差在可接受的迭代精度ε之内时，便认为整个系统达到了一个稳定的状态，此时便可以停止迭代过程，完成AIM区域和PO区域表面电流系数的计算。
在现有的AIM-PO混合方法中，存在着以下问题：
(1)自适应积分法求解速度较慢，尤其在阻抗矩阵的填充以及矩阵方程的迭代过程中消耗较多的求解时间。
(2)PO区域三角面元可见性探查效率低下。
(3)PO区域电流系数以及附加激励向量求解的过程中由于转换矩阵A*和Z*的规模较大，求解时间较长。
这些问题导致AIM-PO算法在求解电大尺寸的电磁问题时计算时间长，效率低下，因此针对该方法发展快速方法是非常必要的。
发明内容
本发明为解决现有技术中存在的AIM-PO算法在求解电大尺寸的电磁问题时计算时间长，效率低下等问题，提出一种利用GPU加速求解迭代的自适应积分法-物理光学发混合方法。
为了解决上述技术问题，本发明提供了基于GPU的AIM-PO混合算法，包括以下步骤：
1.对自适应积分法中的扩展系数矩阵以及近区阻抗矩阵填充进行并行化求解；
2.对自适应积分法中的矩阵方程进行并行化求解；
3.引入基于GPU的无堆栈KD-Tree对PO区域进行遮挡面判断；
4.对AIM-PO混合方法中的PO区域电流系数进行并行化求解；
5.对附加激励向量进行并行化求解。
作为本发明的进一步优化，本发明所述的步骤1)中，利用GPU并行求解扩展系数矩阵，对于每个RWG基函数，需要计算扩展系数用来得到辅助基函数，由于每个RWG基函数的扩展系数向量的计算相互独立，因此在GPU中，为每一个RWG基函数分配一个GPU线程进行多极子展开，用来计算扩展系数矩阵；
近区阻抗矩阵填充包括以下步骤：
1.1、判断两个RWG基函数是否满足近区条件，如果满足则记录该元素在矩阵中的位置索引；
1.2、将步骤1.1中记录的近区阻抗矩阵元素的位置索引传入到GPU端，为每一个近区阻抗元素分配一个GPU线程用来计算其阻抗值。
作为本发明的进一步优化，本发明所述的判断两个RWG基函数满足近区条件为：对近区阻抗矩阵元素Zmn，首先需要判断第m个RWG基函数和第n个RWG基函数构成的基函数对之间的距离，该距离小于近区耦合距离。
作为本发明的进一步优化，本发明所述的步骤2)中对矩阵方程求解，采用双共轭梯度法来进行自适应积分法中的矩阵方程求解，利用CUDA平台中的CUSPARSE库完成了矩阵与向量乘积计算。
作为本发明的进一步优化，本发明所述的步骤3)中对PO区域进行遮挡面判断，具体包括以下步骤：
3.1、在PO区域构建一棵KD-Tree，将其进行线索化；
3.2、根据步骤3.1中构建好的KD-Tree，建立节点数组和面元编号数组，并将节点数组和面元编号数组从主机端拷贝到设备端，进行射线与三角面元的可见性探查测试；
其中，步骤3.1中对KD-Tree进行线索化，包括为树中每个叶子节点添加六个线索指针，每个线索指针保存有与该叶子节点包围盒的一个面相邻的节点地址，在对KD-Tree进行遍历测试射线和面元是否相交的过程中，射线通过在节点包围盒的出射点位置以及节点中的线索指针找到下一个需要访问的节点；
其中，步骤3.2中建立节点数组，将KD-Tree的节点进行编号，将指针的地址信息转化为节点的编号信息，使用一个节点数组保存所有的节点，节点内的所有指针信息均通过节点数组的下标索引来表示，所述指针信息包括左、右节点指针和每个叶子节点内的六个线索指针；
步骤3.2中建立面元编号数组：采用一个新的一维面元编号数组来存储节点内的三角面元编号信息，在节点内保存有三角面元个数和节点内第一个三角形面元在面元编号数组中的索引，每个节点内所包含的三角形面元编号信息在面元编号数组中连续存储，通过三角面元个数和节点内第一个三角形面元在面元编号数组中的索引便可以访问任何一个节点内部三角面元的编号；
步骤3.2中进行射线与三角面元的可见性探查测试：根据叶子节点中保存的线索指针来搜索遍历，具体步骤如下：
5.1、根据AIM区域的RWG基函数和PO区域三角形面元的位置坐标确定射线ray；
5.2、将射线ray和KD-Tree根节点进行相交测试，得到入射位置和出射位置；
5.3、当入射位置距离射线源点比出射位置距离源点近时，计算得到射线ray和节点包围盒的交点，如果当前交点所在的节点不是叶子节点，则根据交点的位置寻找叶子节点，如果交点与左或右子节点相交，则测试射线ray和左子节点或右子节点包围盒，更新入射位置、出射位置以及交点位置，并根据更新后的交点位置信息向下一层寻找子节点，直到找到叶子节点为止；
5.4、找到叶子节点后，需要将射线ray和该叶子节点内的所有三角面元进行相交测试，如果相交，计算相交位置，并更新最近位置λmin以及最近位置对应的三角面元编号Tmin，如果未相交，则不需要更新最近位置λmin；
5.5、相交测试完毕后，利用出射位置判断出射面，并根据叶子节点的线索指针定位下一个需要搜索的节点，如果计算出的下一个节点的索引合法，则将射线ray和下一个节点包围盒进行相交测试，重新计算入射位置和出射位置，返回步骤5.3；如果下一个节点的索引不是一个合法的索引，则表明遍历结束；将得到的最近位置对应的三角面元编号Tmin和PO区域当前进行可见面判断的三角面元的编号进行比较，如果相等，则表明该三角面元对于射线ray可见，位于亮区，否则表明未与该三角面元相交或者被其他的三角面元遮挡，位于暗区。
作为本发明的进一步优化，对于AIM区域的每个RWG基函数，在GPU上构建数量与PO区域的三角面元数量一致的一维线程，所述线程可以并行进行射线与三角面元的可见性探查。
作为本发明的进一步优化，本发明所述的步骤4中，PO区域电流系数的求解需要计算矩阵A*和AIM区域电流系数IAIM的乘积，采用矩阵分块并行求解A*和IAIM的乘积，所述A*为一个K×N的矩阵，其中，N为AIM区域RWG基函数数目，K为PO区域三角面元数量，具体包括以下步骤：
4.1、将矩阵A*按行进行划分得：
上式中，A*被按行划分成了m+1个子矩阵Ai，该子矩阵的规模根据GPU全局存储区剩余空间来确定；
4.2、对A*的子矩阵Ai，设其维度为p×K，对Ai中各个元素的填充过程相互独立，在GPU中分配线程网格空间(p/t+1,K/t+1,1)，其中块内的线程为线程空间(t,t,1)，t为线程块的维度,每个线程计算一个矩阵元素，计算完毕后，便可计算Ai和IAIM的乘积，最终将得到的结果组合得PO区域表面电流系数IPO：
所述步骤5中，附加激励向量ΔV的求解需要计算矩阵Z*和PO区域电流系数IPO的乘积，矩阵Z*是一个K×N的矩阵，其规模和矩阵A*相同，具体包括以下步骤：
a、采取分块计算：将矩阵Z*进行分块：
式中，矩阵Z*被按行划分成了m+1个子矩阵Zi；
b、对Zi并行填充后，可得Zi和IPO的乘积，附加激励向量ΔV被表示为：
本发明的有益效果是：本发明与现有技术相比：本发明可大大减少计算时间，加快电磁计算效率。
附图说明
图1为本发明的利用GPU改进自适应积分法流程图；
图2为本发明的利用GPU并行求解扩展系数矩阵示意图；
图3为本发明的线索化KD-Tree示意图；
图4为本实施例中载有两个蝶形天线的飞机平台示意图；
图5为本实施例中基于GPU的迭代AIM-PO算法和FEKO的MoM-PO计算的远场辐射方向图对比；其中，(a)为XOZ平面；(b)为YOZ平面；(c)为XOY平面。
具体实施方式
以下结合附图和实施例做进一步说明。
步骤一：利用GPU并行化求解扩展系数矩阵和填充近区阻抗元素：
图1为利用GPU并行化改进混合方法中的自适应积分法的流程，远区阻抗矩阵元素由扩展系数矩阵Λx，Λy，Λz，Λd和G等多个矩阵的乘积构成，每个RWG基函数的扩展系数向量都需要进行一次矩阵求逆的计算以及多次矩阵向量的乘积运算，每个RWG基函数的扩展系数矩阵的求解过程相互独立，因此该部分的计算具有良好的并行性，适合使用GPU进行并行化处理。在CUDA编程模型中，若RWG基函数个数为N，可以配置包含有m个GPU线程的一维block块，相对应的，block块的数量应该为个，为更加有效的利用线程，每个block块的的线程个数m应为32的倍数，整个计算流程如图2所示。
在AIM中，近区阻抗元素采用传统矩量法来求解，采用RWG基函数和权函数，并利用M点高斯积分将积分离散化得到阻抗矩阵的元素：
其中，M为高斯积分点的数目，lm为第m个RWG基函数公共边的长度，ln为第n个RWG基函数公共边的长度，从三角形面元内部指向该面元的自由顶点，从该面元的自由顶点指向三角形面元内部，wp和wq为对应的高斯积分点的权值。在式(1)中，实际包含了四个部分的积分，分别是 式中的正负号用来标明T+和T-两个三角面元。当阻抗矩阵中源点和场点相同或者两点的距离特别近时，此时Rpq为零或趋近于零，则求解阻抗矩阵的积分会出现奇异性，不能直接利用高斯积分公式求解。对该奇异性元素的求解有多种方法，如加减奇异项、Duffy坐标变换算法等。由于加减奇异项方法适用范围广，计算量相对较小，因此本实施例中采用加减奇异项的方式来解决奇异性。
近区阻抗矩阵Znear是一个N×N的稀疏矩阵，其中N是RWG基函数的个数，对近区阻抗矩阵元素Zmn，首先需要判断第m个RWG基函数和第n个RWG基函数构成的基函数对之间的距离，如果该距离小于近区耦合距离，才需要求解其精确的阻抗值，如果按照普通的矩阵填充规则，每一个CUDA线程用来判断近区条件并计算阻抗元素，则会导致以下问题：
(1)一般采用自适应积分法来计算的物体都具有电大尺寸，此时描述物体也需要更多的三角面元，对应的RWG基函数对的数量也会很大，这将会导致程序中分配过多的线程。
(2)在设计基于GPU的并行算法时，由于GPU的SIMT特性，应该尽量减少程序中的分支结构。由于在CUDA架构中，GPU的基本执行单元是warp，每个warp中通常有32个线程，warp中的每个线程都会执行相同的指令。如果同一个warp中的不同的线程执行不同的分支，则在同一时刻，除了正在执行分支的线程运行之外，其余的分支的线程都将处于阻塞状态。因为在求解阻抗元素前需要判断两个RWG基函数的距离是否满足近区条件，因此在同一个warp中将会出现大量的分支结构，这会极大的降低GPU中线程的利用率。
为了解决上面两个问题，本实施例将近区阻抗矩阵的填充分为两个步骤，第一步是判断两个RWG基函数是否满足近区条件，如果满足则记录该元素在矩阵中的位置索引。第二步则将第一步记录的近区阻抗矩阵元素索引传入到GPU端，可以为每一个近区阻抗元素分配一个GPU线程用来计算阻抗值。通过分为两个步骤计算，可以有效的消除GPU中条件分支，提高线程的利用效率。
步骤二：对自适应积分法中的矩阵方程进行并行化求解：
在双共轭梯度法BICGSTAB迭代求解矩阵方程的过程中，每次迭代过程需要计算两次矩阵与向量的乘积和多次向量乘积。在自适应积分法中，矩阵与向量的乘积需要计算ZnearI和ZfarI，其中Znear为近区阻抗矩阵，Zfar为远区阻抗矩阵，I为电流系数向量。
近区阻抗矩阵Znear和电流系数向量I的乘积可以借助于CUSPARSE并行运算库来实现，首先将近区阻抗矩阵Znear压缩成CSR(Compressed Sparse Row，压缩稀疏行)格式，然后调用库中相关函数来完成稀疏矩阵和向量相乘运算。CUSPARSE库中提供了优化良好的稀疏矩阵相关操作的高性能并行实现，可以最大限度的发挥GPU的强大的并行性能。CUSPARSE库提供了三个等级的与稀疏矩阵以及系数向量相关的代数操作函数：第一级函数主要用于进行近区阻抗矩阵和稠密向量之间的计算；第二级函数主要用于进行近区阻抗矩阵和稠密向量之间的计算；第三级函数主要用于进行近区阻抗矩阵和多个稠密向量之间的计算；这三级函数均支持单精度浮点数、双精度浮点数以及单精度复数和双精度复数。本实施例中主要利用该库中的第二级函数来实现近区阻抗矩阵Znear和电流系数向量I的乘积。
远区阻抗矩阵Zfar和电流系数向量I相乘的过程包括多次的近区阻抗矩阵与电流系数向量I的运算以及一维快速傅里叶变换及其逆变换。由于扩展系数矩阵为近区阻抗矩阵，因此扩展系数矩阵和向量的乘积运算可以利用CUSPARSE库来高效的实现。计算过程中的快速傅里叶变换及其逆变换可以借助于CUDA中的CUFFT库来快速实现，CUFFT库中提供了快速傅里叶变换及其逆变换的GPU并行接口，有着非常理想的加速效果。
步骤三：利用基于GPU的KD-Tree加速PO区域面元可见性探查；
为了对PO区域进行可见面探查，首先在PO区域构建一棵KD-Tree，并将其进行线索化，然后根据构建好的KD-Tree，建立节点数组和面元编号数组，并将这两个数组从主机端拷贝到设备端，进行射线与三角面元的测试。
线索化KD-Tree的思想是为树中每个叶子节点添加六个线索指针，每个线索指针保存有与该叶子节点包围盒的一个面相邻的节点地址，具体如图3所示。在对KD-Tree进行遍历测试射线和三角面元是否相交的过程中，射线可以通过在节点包围盒的出射点位置以及节点中的线索指针信息找到下一个需要访问的节点，线索KD-Tree的节点的遍历操作不再依赖于堆栈，因此可以在GPU中被高效的实现，并且这种遍历操作也减少了非必要的内部节点遍历，加快了判断的速度。
KD-Tree节点内部会存储该节点包围盒内部的三角面元的编号信息，因为每个节点包围盒内部的三角面元数目不一定相同，因此在节点内会存在动态内存的分配，这意味着节点所占用的内存大小不同，导致内存中的KD-Tree节点信息难以拷贝到GPU的全局内存中。在KD-Tree的每个节点中，保存了指向其左右子节点的指针信息，当将内存中的KD-Tree节点信息拷贝到GPU全局内存中后，节点在GPU全局内存中的地址会和主机端内存的地址不同，导致节点内的指针信息失效，而且由于KD-Tree各个节点之间是通过指针相连，其在主机端内存中的内存空间并不连续，因此也不利于将节点信息从主机端拷贝到GPU全局内存中。
为了解决节点存储空间不连续以及内部指针信息失效的问题，可以将KD-Tree的节点进行编号，将指针的地址信息转化为节点的编号信息。实现该方法需要使用一个节点数组保存所有的节点，节点内的所有指针信息包括左、右节点指针和每个叶子节点内的六个线索指针，均利用节点数组的下标索引来表示，通过将所有的节点保存到一维数组中，解决了节点内存空间不连续的问题，通过利用节点数组的下标索引模拟节点指针，解决了将节点信息从主机端内存拷贝到GPU全局内存中指针信息失效的问题。
在KD-Tree的节点中，由于节点内所存储的三角面元编号的个数不同，导致每个节点所占的内存大小不同，为了实现将节点信息拷贝到GPU全局内存中，则需要对KD-Tree原有的数据结构进行修改，取消节点内动态内存的分配。本实施例采用一个新的一维面元编号数组来存储节点内的三角面元编号信息，在节点内部保存三角面元个数和节点内第一个三角形面元在面元编号数组中的索引，每个节点内部所包含的三角形编号信息在数组中连续存储。通过三角面元个数和节点内第一个三角形面元在面元编号数组中的索引便可以访问任何一个节点内部三角面元的编号。一维面元编号数组也很容易的在主机端和设备端进行拷贝操作。
在利用GPU进行射线与三角面元测试时，不需要在递归的进行节点搜索，而是根据节点中保存的线索指针来搜索遍历，具体算法如下：
(1)根据AIM区域基函数和PO区域三角形面元的位置坐标确定射线ray。
(2)将射线ray和KD-Tree根节点进行相交测试，得到入射位置λentry和出射位置λexit。
(3)当λentry<λexit，即入射位置距离射线源点比出射位置距离源点近时，计算射线ray和节点包围盒的交点Pentry，如果当前交点Pentry所在的节点不是叶子节点，则根据交点Pentry的位置寻找叶子节点，如果Pentry与左子节点相交，则测试射线ray和左子节点包围盒，更新λentry、λexit以及Pentry，并根据Pentry向下一层寻找子节点；若Pentry与右子节点相交，则进行相同的操作，直到找到叶子节点为止。
(4)找到叶子节点后，需要将射线ray和该叶子节点内的所有三角面元进行相交测试，如果相交，计算相交位置，并更新最近位置λmin以及最近位置对应的三角面元编号Tmin,否则执行步骤(5)。
(5)测试完毕后，利用出射位置λexit判断出射面，并根据节点的线索指针信息定位下一个需要搜索的节点，如果计算出的下一个节点的索引合法，则将射线ray和下一个节点包围盒进行相交测试，重新计算入射位置λentry和出射位置λexit，并返回步骤(3)。如果下一个节点的索引不是一个合法的索引，则表明遍历结束。此时将得到的Tmin和PO区域当前进行可见面判断的三角面元的编号进行比较，如果相等，则表明该三角面元对于射线ray可见，位于亮区，否则表明未与该三角面元相交或者被其他的三角面元遮挡，位于暗区。
若AIM区域的RWG基函数的数量为N，PO区域的三角面元数量为K，在AIM-PO算法中，则需要进行N×K次射线与三角面元遮挡测试，即每个RWG基函数都需要和PO区域的每个三角面元进行可见性探查测试。一般来说，PO区域的三角面元数量K的规模都比较大，此时对于AIM区域的每个RWG基函数，可以在GPU上构建与PO区域三角面元数量一致的K个一维线程，这K个一维线程可以并行的进行光线和三角面元的可见性探查，可以大大的加快PO区域的遮挡面判断过程。
步骤四：利用GPU并行求解PO区域表面电流系数：
PO区域电流的求解需要计算矩阵A*和AIM区域电流系数IAIM的乘积。A*是一个K×N的矩阵，其中N为AIM区域RWG基函数数目，K为PO区域三角面元数量。对于具有较多未知量的电大尺寸平台而言，矩阵A*的规模往往很大，无法完整的保存在GPU的全局内存或者主机端内存中，因此本实施例中采用矩阵分块的思想并行求解A*和IAIM的乘积。
由于A*是一个K×N的矩阵，因此可以将矩阵A*按行进行划分得到：
上式中，A*被按行划分成了m+1个子矩阵，子矩阵的规模根据GPU全局存储区剩余空间来确定。可得PO区域表面电流系数：
由于受到GPU全局存储空间的限制，在GPU中，每次计算A*中的一个子矩阵的全部元素。对A*的子矩阵Ai，设其维度为p×K，由于对Ai中各个元素的填充过程相互独立，因此可以在GPU中分配线程网格空间(p/t+1,K/t+1,1)，其中块内的线程为线程空间(t,t,1)，利用每个线程计算一个矩阵元素，计算完毕后，便可以计算Ai和IAIM的乘积，最终将得到的结果组合得到PO区域的表面电流IPO。
步骤五：利用GPU并行求解附加激励：
附加激励ΔV的求解需要计算矩阵Z*和PO区域电流系数IPO的乘积，Z*是一个K×N的矩阵，其规模和A*相同，由于受到内存空间大小的限制，在计算Z*和IPO的乘积时，同样采取分块计算的策略。将Z*进行分块：
式中，Z*被按行划分成了m+1个子矩阵，附加激励可以被表示为：
在利用GPU计算子矩阵Zi的过程中，其GPU线程分配策略可以参照PO区域表面电流IPO求解过程中对A*分块求解的线程分配策略。Zi并行填充完毕后，便可以计算Zi和IPO的乘积，将所有的子矩阵和IPO的乘积计算完毕后便可以组合得到附加激励向量ΔV。
为了验证本实施例的方法的效果，进行以下实验：
1、实验环境
本实验环境如表1所示
表1实验环境
| 环境 | 型号或规格 |
| CPU | Intel Xeon CPU E5-2609v3@1.90GHz×2 |
| 内存 | 32.0GB |
| GPU | NVIDIA Tesla K20c |
| 操作系统 | Windows10 |
2、实验算例
针对一个安装在飞机上方的两个蝶形天线计算其辐射特性，该飞机模型及双蝶形天线模型如图4所示，飞机平台尺寸为18.74m×13.52m×3.97m，两个蝶形天线的工作频率均为500MHz，工作波长为0.60m。在本实施例中，将两个蝶形天线划分为全波区域，使用GPU-AIM进行求解，其余区域划分为PO进行计算。
3、实验结果及分析
图5展示了利用GPU-AIM-PO计算得到的远场辐射方向图和FEKO的MoM-PO的结果实验对比，从图中可以看出实验结果吻合良好。
从表2中可以看出，利用GPU进行并行改进后，算法耗时大大缩短，从1267.5s下降到了111.6s，计算时间减少了约11.4倍，体现出了GPU并行化算法的强大优势。与FEKO的MoM-PO方法相比，优势更加明显。
表2算法内存及耗时对比
| 算法 | 内存需求(GB) | 耗时(s) |
| GPU-AIM-PO | 1.43 | 111.6 |
| CPU-AIM-PO | 1.43 | 1267.5 |
| FEKO-MoM-PO | 0.56 | 4549.47 |
Claims (7)
1.基于GPU的AIM-PO混合算法，其特征在于：包括以下步骤：
1)对自适应积分法中的扩展系数矩阵以及近区阻抗矩阵填充进行并行化求解；
2)对自适应积分法中的矩阵方程进行并行化求解；
3)引入基于GPU的无堆栈KD-Tree对PO区域进行遮挡面判断；
4)对AIM-PO混合方法中的PO区域电流系数进行并行化求解；
5)对附加激励向量进行并行化求解。
2.根据权利要求1所述的基于GPU的AIM-PO混合算法，其特征在于：所述步骤1)中，利用GPU并行求解扩展系数矩阵：在GPU中，为每一个RWG基函数分配一个GPU线程进行多极子展开，用来计算扩展系数矩阵；
近区阻抗矩阵填充包括以下步骤：
1.1、判断两个RWG基函数是否满足近区条件，如果满足则记录该元素在矩阵中的位置索引；
1.2、将步骤1.1中记录的近区阻抗矩阵元素的位置索引传入到GPU端，为每一个近区阻抗元素分配一个GPU线程用来计算其阻抗值。
3.根据权利要求2所述的基于GPU的AIM-PO混合算法，其特征在于：所述步骤1.1中判断两个RWG基函数满足近区条件为：对近区阻抗矩阵元素Zmn，首先需要判断第m个RWG基函数和第n个RWG基函数构成的基函数对之间的距离，该距离小于近区耦合距离。
4.根据权利要求1所述的基于GPU的AIM-PO混合算法，其特征在于：所述步骤2)中对矩阵方程求解：采用双共轭梯度法来进行自适应积分法中的矩阵方程求解，利用CUDA平台中的CUSPARSE库完成矩阵与向量乘积计算。
5.根据权利要求1所述的基于GPU的AIM-PO混合算法，其特征在于：所述步骤3)中对PO区域进行遮挡面判断，具体包括以下步骤：
3.1、在PO区域构建一棵KD-Tree，将其进行线索化；
3.2、根据步骤3.1中构建好的KD-Tree，建立节点数组和面元编号数组，并将节点数组和面元编号数组从主机端拷贝到设备端，进行射线与三角面元的可见性探查测试；
其中，步骤3.1中对KD-Tree进行线索化，包括为树中每个叶子节点添加六个线索指针，每个线索指针保存有与该叶子节点包围盒的一个面相邻的节点地址，在对KD-Tree进行遍历测试射线和面元是否相交的过程中，射线通过在节点包围盒的出射点位置以及节点中的线索指针找到下一个需要访问的节点；
其中，步骤3.2中建立节点数组：将KD-Tree的节点进行编号，将指针的地址信息转化为节点的编号信息，使用一个节点数组保存所有的节点，节点内的所有指针信息均通过节点数组的下标索引来表示，所述指针信息包括左、右节点指针和每个叶子节点内的六个线索指针；
步骤3.2中建立面元编号数组：采用一个新的一维面元编号数组来存储节点包围盒内部的三角面元编号信息，在节点内保存有三角面元个数和节点内第一个三角形面元在面元编号数组中的索引，每个节点内所包含的三角形面元编号信息在面元编号数组中连续存储，通过三角面元个数和节点内第一个三角形面元在面元编号数组中的索引便可以访问任何一个节点内部三角面元的编号；
步骤3.2中进行射线与三角面元的可见性探查测试，根据叶子节点中保存的线索指针来搜索遍历，具体步骤如下：
5.1、根据AIM区域RWG基函数和PO区域三角形面元的位置坐标确定射线ray；
5.2、将射线ray和KD-Tree根节点进行相交测试，得到入射位置和出射位置；
5.3、当入射位置距离射线源点比出射位置距离源点近时，计算得到射线ray和节点包围盒的交点，如果当前交点所在的节点不是叶子节点，则根据交点的位置寻找叶子节点，如果交点与左或右子节点相交，则测试射线ray和左子节点或右子节点包围盒，更新入射位置、出射位置以及交点位置，并根据更新后的交点位置信息向下一层寻找子节点，直到找到叶子节点为止；
5.4、找到叶子节点后，需要将射线ray和该叶子节点内的所有三角面元进行相交测试，如果相交，计算相交位置，并更新最近位置λmin以及最近位置对应的三角面元编号Tmin，如果未相交，则不需要更新最近位置λmin；
5.5、相交测试完毕后，利用出射位置判断出射面，并根据叶子节点的线索指针定位下一个需要搜索的节点，如果计算出的下一个节点的索引合法，则将射线ray和下一个节点包围盒进行相交测试，重新计算入射位置和出射位置，返回步骤5.3；如果下一个节点的索引不是一个合法的索引，则表明遍历结束；将得到的最近位置对应的三角面元编号Tmin和PO区域当前进行可见面判断的三角面元的编号进行比较，如果相等，则表明该三角面元对于射线ray可见，位于亮区，否则表明未与该三角面元相交或者被其他的三角面元遮挡，位于暗区。
6.根据权利要求5所述的基于GPU的AIM-PO混合算法，其特征在于：对于AIM区域的每个RWG基函数，在GPU上构建数量与PO区域的三角面元数量一致的一维线程，所述线程可以并行进行射线与三角面元的可见性探查。
7.根据权利要求1所述的基于GPU的AIM-PO混合算法，其特征在于：
所述步骤4中，PO区域电流系数的求解需要计算矩阵A*和AIM区域电流系数IAIM的乘积，采用矩阵分块并行求解A*和IAIM的乘积，所述A*为一个K×N的矩阵，其中，N为AIM区域RWG基函数数目，K为PO区域三角面元数量，具体包括以下步骤：
4.1、将矩阵A*按行进行划分得：
<mrow>
<msup>
<mi>A</mi>
<mo>*</mo>
</msup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mi>m</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
上式中，A*被按行划分成了m+1个子矩阵Ai，该子矩阵的规模根据GPU全局存储区剩余空间来确定；
4.2、对A*的子矩阵Ai，设其维度为p×K，对Ai中各个元素的填充过程相互独立，在GPU中分配线程网格空间(p/t+1,K/t+1,1)，其中块内的线程为线程空间(t,t,1)，t为线程块的维度，每个线程计算一个矩阵元素，计算完毕后，便可计算Ai和IAIM的乘积，最终将得到的结果组合得PO区域表面电流系数IPO：
<mrow>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
<mo>=</mo>
<msup>
<mi>A</mi>
<mo>*</mo>
</msup>
<msup>
<mi>I</mi>
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>M</mi>
</mrow>
</msup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>A</mi>
<mn>0</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>M</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>M</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>A</mi>
<mn>2</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>M</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>A</mi>
<mi>m</mi>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>M</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
</mrow>
所述步骤5中，附加激励向量ΔV的求解需要计算矩阵Z*和PO区域电流系数IPO的乘积，矩阵Z*是一个K×N的矩阵，其规模和矩阵A*相同，具体包括以下步骤：
a、采取分块计算：将矩阵Z*进行分块：
<mrow>
<msup>
<mi>Z</mi>
<mo>*</mo>
</msup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>Z</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>Z</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>Z</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>Z</mi>
<mi>m</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
式中，矩阵Z*被按行划分成了m+1个子矩阵Zi；
b、对Zi并行填充后，可得Zi和IPO的乘积，附加激励向量ΔV被表示为：
<mrow>
<mi>&Delta;</mi>
<mi>V</mi>
<mo>=</mo>
<msup>
<mi>Z</mi>
<mo>*</mo>
</msup>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Z</mi>
<mn>0</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Z</mi>
<mn>1</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Z</mi>
<mn>2</mn>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Z</mi>
<mi>m</mi>
</msub>
<msup>
<mi>I</mi>
<mrow>
<mi>P</mi>
<mi>O</mi>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>.</mo>
</mrow>
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710638068.6A CN107562690A (zh) | 2017-07-31 | 2017-07-31 | 一种基于gpu的aim‑po混合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710638068.6A CN107562690A (zh) | 2017-07-31 | 2017-07-31 | 一种基于gpu的aim‑po混合方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107562690A true CN107562690A (zh) | 2018-01-09 |
Family
ID=60973694
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710638068.6A Withdrawn CN107562690A (zh) | 2017-07-31 | 2017-07-31 | 一种基于gpu的aim‑po混合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107562690A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110867218A (zh) * | 2018-08-27 | 2020-03-06 | 中国石油化工股份有限公司 | 一种分子电子能量计算方法及系统 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101794355A (zh) * | 2010-03-26 | 2010-08-04 | 中国人民解放军空军工程大学 | 电大物体电磁辐射和散射的计算机核外并行计算方法 |
| CN103888205A (zh) * | 2014-03-24 | 2014-06-25 | 上海华为技术有限公司 | 一种电磁波传播预测方法和装置 |
| CN106529082A (zh) * | 2016-12-02 | 2017-03-22 | 上海无线电设备研究所 | 一种快速计算电大尺寸目标电磁散射特征的方法 |
-
2017
- 2017-07-31 CN CN201710638068.6A patent/CN107562690A/zh not_active Withdrawn
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101794355A (zh) * | 2010-03-26 | 2010-08-04 | 中国人民解放军空军工程大学 | 电大物体电磁辐射和散射的计算机核外并行计算方法 |
| CN103888205A (zh) * | 2014-03-24 | 2014-06-25 | 上海华为技术有限公司 | 一种电磁波传播预测方法和装置 |
| CN106529082A (zh) * | 2016-12-02 | 2017-03-22 | 上海无线电设备研究所 | 一种快速计算电大尺寸目标电磁散射特征的方法 |
Non-Patent Citations (3)
| Title |
|---|
| XING WANG ET AL.: "Efficient Analysis of Antennas Mounted on Large-Scale Complex Platforms Using Hybrid AIM-PO Technique", 《IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION》 * |
| 张春杨: "基于KD-tree的目标电磁散射快速算法", 《中国优秀硕士学位论文全文数据库-基础科学辑》 * |
| 陈金鑫 等: "GPU加速的自适应积分法研究", 《微波学报》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110867218A (zh) * | 2018-08-27 | 2020-03-06 | 中国石油化工股份有限公司 | 一种分子电子能量计算方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Elsen et al. | Large calculation of the flow over a hypersonic vehicle using a GPU | |
| Stantchev et al. | Fast parallel particle-to-grid interpolation for plasma PIC simulations on the GPU | |
| Naumov et al. | AmgX: A library for GPU accelerated algebraic multigrid and preconditioned iterative methods | |
| Huthwaite | Accelerated finite element elastodynamic simulations using the GPU | |
| Castro et al. | GPU computing for shallow water flow simulation based on finite volume schemes | |
| Ergul et al. | A hierarchical partitioning strategy for an efficient parallelization of the multilevel fast multipole algorithm | |
| CN110275733A (zh) | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 | |
| Lastra et al. | Simulation of shallow-water systems using graphics processing units | |
| CN112733364B (zh) | 一种基于阻抗矩阵分块的箔条云散射快速计算方法 | |
| Dang et al. | The sliced coo format for sparse matrix-vector multiplication on cuda-enabled gpus | |
| Zubair et al. | An optimized multicolor point-implicit solver for unstructured grid applications on graphics processing units | |
| KR20120057274A (ko) | Ｉｄ－ｆｄｔｄ 방법을 이용한 전자기파 수치해석 방법 | |
| Stefanski | Implementation of FDTD-compatible Green's function on heterogeneous CPU-GPU parallel processing system | |
| Gao et al. | Fast RCS prediction using multiresolution shooting and bouncing ray method on the GPU | |
| CN107562690A (zh) | 一种基于gpu的aim‑po混合方法 | |
| Palmesi et al. | Highly parallel demagnetization field calculation using the fast multipole method on tetrahedral meshes with continuous sources | |
| Siebenborn et al. | A curved-element unstructured discontinuous Galerkin method on GPUs for the Euler equations | |
| Göddeke et al. | Finite and spectral element methods on unstructured grids for flow and wave propagation problems | |
| WO2007081303A1 (en) | Applications of interval arithmetic for reduction of number of computations in ray tracing problems | |
| Marek et al. | An efficient parallelization strategy for the adaptive integral method based on graph partitioning | |
| Hallock et al. | Improving reaction kernel performance in lattice microbes: particle-wise propensities and run-time generated code | |
| Kiss et al. | Fast analysis of metallic antennas by parallel moment method implemented on CUDA | |
| Barnes | Improved C1 shape functions for simplex meshes | |
| Zhang et al. | Parallel SVD-CBFM based on GPU acceleration | |
| Bessonov et al. | Parallelization of the preconditioned IDR solver for modern multicore computer systems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WW01 | Invention patent application withdrawn after publication | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20180109 |