[FDTD]高性能计算机群系统并行FDTD算法研究

并行FDTD   MPI   高性能计算机群 gbNPD*7g9  
本文研究了基于高性能计算机群系统的以消息传递方式（MPI）执行的并行FDTD（时域有限差分）算法，具体描述了数据交换方法和三维方向切割时的新数据类型定义，完成3×3×3的并行FDTD算法研究，验证了算法的可行性和有效性，为扩展为任意处理结点计算奠定了基础。 p-n_ ">7  
.-[uQtyWW  
引言 !%b.k6%>w  
时域有限差分算法（Finite-Difference Time-Domain Method）作为一种电磁学领域的研究方法具有许多优点，但它的实现和优越性的发挥受到计算平台的极大限制，其算法稳定性要求空间网格的最大尺寸应小于波长的十分之一[1]，因此在高频或计算复杂目标电大尺寸物体的情况下，用FDTD方法计算需要定义庞大数量的网格数，而当前PC机计算速度和内存容量有限，限制了FDTD方法的发展。 NQefrof  
FDTD算法具有天然的并行性，其每一个网格上的电磁场值只与相邻网格上的电磁场值有关，而与远端网格的值无关[1]，因此适合于空间上划分计算区域，分配给每个计算结点单独的计算任务，只需要在结点间传递分割面上的某部分电场或磁场值即可。 
#W`
>vd}  
本文的计算方法是在由DELL PowerEdge 1850 计算服务器和高性能的Myrinet 2000数据交换网络组成的中国传媒大学高性能计算机群系统（峰值浮点计算能力达到7296亿次/秒）上，以Linux操作系统和MPI软件为基本平台，以FDTD为基本算法，通过使用不同的数据类型定义和最佳虚拟拓扑编写了3×3×3划分的三维并行算法，旨在对并行FDTD有个更加深入广泛的研究和应用。 K|$Dnma^n  
2 区域分解方法  e;8>/G  
为方便分析，先以一维FDTD运算为例，由图1所示。要计算n+1时间步z=k网格点的电场分量Exn+1(k)，需要上个时间步该点的电场分量Exn(k) 及上半个时间步的磁场分量Hyn+1/2(k-1/2), Hyn+1/2(k+1/2)，因此必须先完成Hyn+1/2(k＋1/2)从处理器1向处理器0的传递，而且每个时间步都必须传递一次。而要完成处理器1内磁场计算，就需要Exn+1(k)由处理器0向处理器1传递。由此可看出在FDTD递推过程中将有两次中断，一次为电场交换，一次为磁场交换，即每次电场递推完成后，程序暂停下来进行电场的数据交换，同样磁场递推完成后，程序也暂停下来进行磁场交换，这样程序多次暂停中断，影响了计算效率，下面介绍另一种方法。 Ep-{Ew{T_=  
如图2所示，在左右两个处理器中都计算k点的电场，处理器1计算k点电场时需要传递处理器0的磁场，处理器0计算k点电场时也需要传递处理器1的磁场，因此在相邻两个子区域间只需传递磁场而不需要电场传递，代价只是多计算一点的电场值，就可以使递推过程中每一时间步中断只有一次。数值实验[2]表明，第二种方法的计算效率比第一种高。 X;ef&n`U0  
w$Lpuun{  
图1 一维并行计算区域分解方法1 l/[0N@r~  
pC,MiV$c"  
r2?-QvQ  
图2 一维并行计算区域分解方法2 S^|Uzc  
:_v/a+\n  
如果并行计算区域是三维的，则可以按照需要对其进行三个方向的分割，所形成的区域为长方体，截断边界为平面，需要交换的数据是与分割面平行的磁场量。如图3所示，当对三维物体进行x方向分割时，需要传递的数据是与切割面平行的Hy、Hz分量，以此类推，当y方向切割，需要传递数据为Hx、Hz分量，当z方向分割，需要传递数据为Hx、Hy分量。 0JTDJZOz@#  
cF{5[?wS  
图3 三维模型分割方式 h_6QVab@  
3 虚拟进程拓扑及不连续数据发送 b%>vhj&F  
3．1 虚拟进程拓扑 AkdONKO8{  
在许多并行应用程序中，进程的线性排列不能充分地反映进程间在逻辑上的通信模型，因此进程被排列成二维或三维网格形式的拓扑模型，并且用一个图来描述逻辑进程排列，我们称这种逻辑进程排列为“虚拟拓扑”[3]，在引入虚拟进程拓扑后，不管具体并行机的处理器物理上是如何排列的，都可以实现所希望的并行计算模型。 hU=n>g>nx  
假设只有4个处理器可用，划分为2×2的计算区域，每个处理器计算一块区域，则引入虚拟拓扑后数据块和处理器进程阵列之间的对应关系，如图4所示。 H|>dF)%pj  
T)Z2=5V  
图4 引入虚拟进程后的FDTD数组 '/+l\.z"&  
假设2nx、2ny分别是单个处理器中x、y方向FDTD的网格数，考虑交换数据的存储，需要在串行程序数组的基础上向两边各扩展一个网格，按行列同时划分的数组需要和上下、左右的邻居同时通信，因此上下左右都需要预留出通信的部分，用来存放同各个方向邻居通信得到的数据，即在均匀划分后各块网格数大小为（nx+2）×（ny+2）。 ~?&;nTwHe  
在本文中用的是常用的笛卡尔（Cartesian）拓扑结构。在引入进程拓扑和虚拟进程后，当一个真实进程向一个虚拟进程发送数据或从一个虚拟进程接收数据时，该真实进程会立即正确返回，如同执行了一个空操作，因此可以把边界块和内部块统一看待，全部通信都用捆绑发送MPI_SENDRECV语句实现，不管计算区域在一个方向上划分为多少个进程，都可以用该语句实现数据的接收和发送，只要事先用MPI_CART_SHIFT语句获得进程间的关系即可方便实现上下、左右、前后的场值通信。 &$uQ$]&H  
3．2 不连续数据发送 hPtSY'
_@_  
如果按照通常的程序设计方式，会将每个分割面上的所有须传递的场量值利用一个循环语句逐个进行传递，这种方法很繁琐，最主要的问题在于其带来的通信开销很大，为了避免这种情况，在通信之前，可以将要进行通信的数据组合到一个数据块里面。为此，MPI具有处理不连续数据的功能，可以在发送方连续复制若干个旧数据类型形成一个“块”，再通过一定间隔复制该块形成新的数据类型，然后发送出去，这样原来多次的消息发送操作只需一次就可以完成，提高了并行计算的性能。 #UGtYD}"  
首先我们必须清楚数组在内存中的存储结构，在C语言中，数组在内存中是按z方向优先排列，因此z方向的数据是连续存放的，可以直接进行通信，但当上下左右各增加一个通信网格后，需要传递的数据个数只是nx×ny，而不是全部的（nx+2）×（ny+2）个网格，则连续的数据被割裂，要形成一个块来进行传递就需要定义新数据类型。 .um]1_= \  
如图5是以z方向分割时的切割面，在FDTD递推中需要传递的数据是Hx、Hy分量。在这里可以用MPI_TYPE_INDEXED来形成新的数据类型，它的功能是复制一个旧数据类型到多个块序列中，每个块可以包含不同的拷贝数目和具有不同的偏移。如图5所示，每个块只包含1个元素，共有nx×ny个块，元素之间的间隔是（nz+2），但是从（1，3，0）至（2，1，0）时却需要跨越3×（nz+2）的距离，因此必须根据不同的位置设置块之间不同的偏移量。 tK)E*!  
t{t*.{w  
图5 以z方向分割时的数据结构 {`fhcEC  
当以y方向分割时，在FDTD递推中需要传递的数据是Hx、Hz分量，如图6所示，阴影部分的几个数组在内存中是连续的，因此可以作为一个块，每块包含nz个元素，共有nx个块，每块第一个元素之间的间隔是等距的，即从（1，0，1）至（2，0，1）之间的距离，为（ny+2）×（nz+2），因此可以用更通用的生成器MPI_TYPE_VECTOR，它的功能是可以将旧数据类型等间隔的复制形成新的数据类型。 8mTjf Br  
sxM0c  
图6 以y方向分割时的数据结构 oeKI9p13\  
同理，当以x方向分割时，也可以用这个语句，只是参数与上不同，其每块包含nz个元素，共有ny个块，每块第一个元素之间的间隔为（nz+2）。按照不同方向切割用不同的方法来定义新数据类型，保证了程序的正确性和高效性。 ##VS%&{  
4 3×3×3三维划分并行设计 De`)`\U  
4．1 边界的并行化处理 UP=0>jjbn:  
要实现3×3×3的三维方向划分，只要在三个方向都划分为三个子区域，MPI虚拟拓扑沿着三个方向建立，并且根据不同的方向定义不同的新数据类型来进行分割面上的平行磁场值交换，就可以实现三维并行分割与计算。但是在集合FDTD计算时还必须考虑具体FDTD的吸收边界、连接边界以及输出边界的并行化处理[4]。 m#Rgelhk.  
在连接边界的并行化处理上，入射波的加入还是按照原来串行程序入射波的设置方法进行，只需要利用几个简单的线性关系表达式，将并行区域内的坐标转换为区域分割前串行的坐标值即可，不需要进行并行数据交换。 w^N xR,  
由于输出边界的实质是根据边界面上场源的数值求出远区场的数值，不涉及差分计算，故不需要区域间的数据传递，只需要根据区域的分割情况分别求出各区域内输出边界面上的场源数值，然后应用归约函数MPI_REDUCE来求和，将分散于各个计算区域内求得的每一点输出边界上场源值相加后集中归约到某一个区域内输出，即可得到与串行相同的结果。 >*rsRR  
本文所应用的FDTD采用的是各向异性介质PML吸收边界，层数设为10，如图7所示，若在x方向分割成三个区域，则左右两个区域的PML介质损耗值σx必须根据在该区域内具体的坐标值进行设置，而中间区域σx则为0，若在三个方向同时分割为3×3×3个区域，也根据此原则分别对σx、σy、σz进行设置即可。 *=i|E7Irg  
4．2 并行计算具体实现 qH'T~#S  
本文以微带线馈电的矩形贴片微带天线为例，天线结构如图8所示，贴片为12.45mm×16mm的长 .. zt{?Ntb  
[G:wPp.y  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

5 结论 o/273I  
本文对并行FDTD计算中的区域分割、虚拟拓扑、新数据类型定义以及FDTD的边界化处理进行了详细的描述，具体地讨论了并行FDTD的实现方法，对于PC网络或高性能计算机群中MPI编程环境下并行FDTD的研究工作具有现实的参考意义。 #7/;d=  
dH"wYMNL  
参考文献： t-_~jZ<  
[1]葛德彪，阎玉波. 电磁场时域有限差分方法[M] .西安：西安电子科技大学出版社， 2002 ``?]13XjK  
[2]余文华，苏涛，Raj Mittra，刘永峻. 并行时域有限差分[M] 北京：中国传媒大学出版社， 2005 _bB:1l?V  
[3]都志辉.  高性能计算并行编程技术——MPI并行程序设计[M]  北京：清华大学出版社， 2001 ximW!y7  
[4]冯峰、逯贵祯、关亚林.  二维FDTD算法的网络并行运算研究[J].  北京广播学院学报（自然科学版），2005，(2)：33－37