加速并行时域有限差分仿真的新方法

张立红 余文华 杨小玲

（1.中国传媒大学信息工程学院，北京 100024；2.中国人民武装警察部队学院，河北 廊坊 065000；3.Penn State University，USA PA 16802）

引 言

时域有限差分（FDTD）法最早由K.S.Yee在1966年提出，经过几十年的发展，FDTD已经形成了一套比较完善的方法体系，相对于其他的计算电磁学方法，FDTD因其简单灵活而受到广大电磁计算研究者的欢迎，但它的实现却面临着一些问题，如庞大的计算量是普通PC机所不能满足的，因此，科学家、电磁工作者等提出了各种各样的并行算法来解决这些问题，如基于消息传递接口（MPI）的并行技术、基于OpenMP的共享存储编程技术以及基于映射文件的技术等［1-4］。最近，也有一些文献提出用图形处理单元（GPU）对FDTD算法进行加速［5］。文章提出了一种利用单指令多数据流式扩展（SSE）指令集来加速并行FDTD仿真的新方法，用C语言开发了基于MPI库、OpenMP和SSE指令集的三维并行FDTD代码，最后以具体的电磁仿真实例验证了新方法的可行性和加速效率，并将其与普通并行FDTD仿真方法进行了对比。

1.理论分析

1.1 FDTD方法

在FDTD方法中，电磁波传播以及电磁波与物质的相互作用是通过电场和磁场在空间和时间上的差分递推实现的，空间某处的电场值可以由该处上一时间步的电场值和其周围上半个时间步的四个磁场值计算得到，而空间某处的磁场值可以由该处上一时间步的磁场值和其周围上半个时间步的四个电场值计算得到。在FDTD方法中，电磁场值的位置和递推关系可以用式（1）和图1表示。公式式（1）表示的是磁场沿z轴方向分量的递推公式，其他两个方向的分量以及电场分量的递推公式与（1）完全相似［6］。

图1 电磁场值关系图

从递推公式（1）和图1可以看出，FDTD方法具有与生俱来的并行性:FDTD方法中，每一个网格点的磁场（电场）分量的迭代公式只与它自己上一时间步的值和它周围网格点电场（磁场）上半个时间步的值有关，而与计算区域内其他网格点的场值没有直接关系，非常适合并行计算［7］；而且，递推公式中，几乎所有的计算都是对一组数据进行相同的加、减或乘法操作，非常适合单指令多数据（SIMD）模式的并行处理。

1.2 SSE指令集

SSE指令集是Intel在其芯片中实现了基于寄存器的SIMD架构之后提供的指令集，它使用8个独立的128位寄存器，允许SIMD计算同时作用于4个紧缩的单精度浮点数据单元。SSE指令集包括70条指令，其中包含提高3D图形运算效率的50条SIMD浮点运算指令、12条多媒体扩展（MMX）整数运算增强指令和8条优化内存中连续数据块传输指令［8］。AMD处理器也加入了对SSE指令集的支持。现在市场上能够买到的处理器大都支持SSE指令集。

图2 单指令多数据操作

因为SSE指令集是单指令多数据操作，因此，可以通过循环展开来减少运算时间，从而提高运算速度。SSE指令集要求它的操作数是一种新的紧缩类型，对于float类型的数据，SSE指令集的操作数是把4个float标量数据压缩成一个类型的向量数据。图2是一个典型的单指令多数据的操作。图2中的a和b都是类型的向量数据，都是由4个float类型的标量数据压缩而成的，SSE指令对a和b进行加法运算操作，得到一个类型的向量数据，存放在Result中。

2.数值实验及结果

2.1 三级数据并行结构

普通的基于MPI或OpenMP的并行FDTD算法都是一级或两级并行结构［9］，把基于SSE指令集的新加速方法引入到FDTD仿真后，程序形成了三级数据并行结构。

第一级数据并行基于MPI库。把要进行仿真的区域进行区域分解，各子区域的电磁场值分别独立计算，负责计算各子域的进程之间通过MPI库的消息传递函数进行通信。

第二级数据并行采用OpenMP共享存储编程实现。首先利用OpenMP生成多个线程，然后将每个子区域的计算再分配给各个线程并行执行，算法实现框架如下:

第三级数据并行利用SSE指令集实现。对于单精度浮点运算，普通的运算操作一次得到一个计算结果，而使用SSE指令集，一个运算可以得到四个计算结果，从而实现细粒度数据并行，加快了计算速度。

2.2 SSE加速实现

在利用SSE加速并行FDTD算法时，仅对电磁场的递推部分进行了加速，其中包括整个计算区域的电磁场的递推过程和卷积完全匹配层（CPML）吸收边界［10］的处理过程。先讨论整个计算区域的电磁场的递推情况。以计算磁场沿z轴方向的分量Hz为例，可按照以下步骤实现:

1）定义SSE所需的类型的变量，并为其赋值（作为SSE指令的操作数）。

2）把电磁场递推公式中所需的系数加载到SSE寄存器中。

3）把float类型的指针变量转换成SSE所需的类型的指针变量。

4）把原来的最内层循环展开，循环次数变为原来的四分之一（这就是SSE指令集对FDTD仿真进行加速的原理）。

5）磁场值的递推计算。

计算Hz的伪代码如下:

CPML吸收边界的处理方法与电磁场值递推过程类似，例如，在计算沿y轴方向的CPML区域的磁场时，可以参照前面的步骤，实现伪代码如下:

为了优化程序，提高缓存命中率，还可以把计算电磁场值的基本递推过程和CPML吸收边界的处理过程合并起来，通过判断j的取值是否在CPML吸收边界区域内来决定是否进行电磁场值边界的更新，实现框架如下:

2.3 数值实验结果

为了验证新方法的加速效率，文章进行了实验测试，分别计算了40×40×40、80×80×80和160×160×160个均匀网格的真空中电磁波的传播，其中，激励源为高斯脉冲源，放置在立方体计算区域的正中心，电磁场初始值均设为0，CPML吸收边界为6层。实验平台是PC机，CPU为Intel的T2300（双核），1.66GHz，时间步为400，实验结果如表1所示。从测试结果可以看出，使用了SSE指令集加速的代码比普通的并行代码所需的计算时间大大减少。

表1 计算时间及加速比

普通并行代码一个指令进行一次运算操作，得到一个结果值，而SSE代码一个指令进行一次运算，得到四个结果值，因此，理论上使用SSE指令集加速时，最理想情况是加速比等于4，实验在160×160×160均匀网格时间步为400的情况下，得到的加速比为2.62，加速效果比较好。

3.结 论

提出了利用SSE指令集来加速并行FDTD算法的方法，开发了三级数据并行结构的三维FDTD仿真代码，在Intel T2300的PC机上实现了对基于MPI和OpenMP的三维并行FDTD仿真的加速，得到的加速比为2.62。使用SSE指令集加速，无需额外购买任何硬件，只需改变部分并行代码即可实现，因此，使用SSE指令集来加快运算速度，从而减少运行时间是一种高效、经济的新途径。

［1］余文华，杨小玲，刘永俊，等.并行FDTD和IBM BlueGene／L巨型计算机结合求解电大尺寸的电磁问题［J］.电波科学学报，2006，21（4）:562-566.YU Wenhua，YANG Xiaoling，LIU Yongjun，et al.Solving electrically large EM problems using parallel FDTD and IBM BlueGene／L supercomputer［J］.Chinese Journal of Radio Science，2006，21（4）:562-566.（in Chinese）

［2］雷继兆，梁昌洪，丁 伟，等.机载天线辐射特性的并行FDTD分析［J］.电波科学学报，2008，23（6）:1139-1143.LEI Jizhao，LIANG Changhong，DING Wei，et al.A-nalysis of radiation characters of airborne antennas with parallel FDTD［J］.Chinese Journal of Radio Science，2008，23（6）:1139-1143.（in Chinese）

［3］YU W，LIU Y，SU T，et al.A robust parallel conformal FDTD processing package using the MPI library［J］.IEEE Ant.and Prop.Mag.，2005，47（3）:39-59.

［4］刘 瑜，梁 正，杨梓强.基于映射文件的电磁并行FDTD算法实现研究［J］.电波科学学报，2008，23（4）:634-639.LIU Yu，LIANG Zheng，YANG Ziqiang.Implementation of parallel FDTD algorithm based on mapped file［J］.Chinese Journal of Radio Science，2008，23（4）:643-639.（in Chinese）

［5］Xu K，Fan Z H，Ding D Z，et al.GPU accelerated unconditionally stable crank-Nicolson FDTD method for the analysis of three-dimensional microwave circuits［J／OL］.Progress in Electromagnetics Research（PIER），2010，102:381-395［2011-03-25］.http:／／www.jpier.org／PIER／pier.php paper=10020606.

［6］葛德彪，闫玉波.电磁场时域有限差分方法［M］.西安电子科技大学出版社，2002.

［7］余文华，苏 涛，Raj Mittra，等.并行时域有限差分［M］.北京:中国传媒大学出版社，2005.

［8］Intel corporation.Intel Architecture Optimization Reference Manual［M／OL］.USA:Intel corporation，1999［2011-03-25］.http:／／www.intel.com／design／pentiumii／manuals／245127.htm.

［9］都志辉.高性能计算并行编程技术-MPI并行程序设计［M］.北京:清华大学出版社，2001.

［10］TAFLOVE A，HAGNESS S C.Computational Electrodynamics the Finite-Difference Time Domain Method［M］.London:Artech House，2005.