时域有限差分法的图形处理单元的加速

摘 要：时域有限差分法，即FDTD（Finite Difference Time Domain），是计算电磁学的一种重要方法。作为一种天然的并行算法，它的计算过程可以划分为多个同时进行相似计算的子计算。这个方法主要是把麦克斯韦方程在时间上和空间上进行差分化，并且通过时间领域上的更新来模仿电磁场的变化来计算问题，因而有利于解决很多电磁场问题。而图形处理单元即GPU（Graphic Processing Unit）相对于CPU的高性能计算速度以及NVIDA公司生产的GPU特有的高并行结构，为时域有限差分的加速提供了可能。

关键字：时域有限差分法；图形处理单元；麦克斯韦方程；并行算法

1 FDTD的基本原理

FDTD算法是1966年K.S.Yee发表在AP上的一篇论文建立起来的，后被称为Yee网格空间离散方式。核心思想是把带时间变量的麦克斯韦旋度方程转化为差分形式，故中文称之为“时有限差分法”。麦克斯韦方程如下：

其中H是磁场强度，E是电场强度，D是电位移，B是磁通量密度。

上述两个矢量方程描述了麦克斯韦方程中磁场与电场复杂交错的关系。由于在三维空间中每个矢量方程又可以分解为三个标量方程，因此该方程组可以化为6个标量方程如下所示：

这样的话便将问题的几何空间离散为空间网格，电场和磁场的分量便被置于空间离散的网格点上。而这是FDTD计算的前提。我们再用差分近似替代麦克斯韦这6个标量方程中的时间和空间导数，构造一系列方程，均以前一时间步电磁场瞬时值来“预测”后一时间步电磁场的瞬时值，由此构造时间不断向前推近的算法，来模拟时域中的电磁场变化过程。

1966年，Yee首次给出了麦克斯韦旋度方程的一组差分形式，这组方程在空间和时间上一离散的形式给出，使用的是中心差分法。对空间（X轴方向）的中心差分法离散公式如下：

对Y轴，Z轴方向的中心差分离散以此类推。

对时间的中心差分离散公示如下：

图1为Yee单元网格的结构：

我们根据6个麦克斯韦标量方程，结合上述中心差分法，便可以得到三维问题下的FDTD更新方程如下：

其余的五个方程也如法可以写出，因此任何时刻可一次算出一个点，并行算法可计算出多个点。通过这些运算可以交替算出电场磁场在各个时间步的值。

上述方法运行时忽略了两个很重要的问题，一是数值色散，就是该方法会产生一定误差，如何维持该方法可行的精度问题，二是计算机中模拟的无限过程与计算机本身储存内存有限相矛盾的问题。

解决数值色散的关键是空间网格大小的选择，时间步长△t，空间步长△x， △y， △z必须满足一定的关系，否则就使得数值表现不稳定，表现为：随着计算步数的增加，计算場量的数值会无限的增大，这种增大不是由于误差积累造成的，而是由于电磁波的传播关系被破坏造成的。所以△t， △x， △y， △z必须满足一定的关系以保证稳定性，其数值稳定条件为如下：

我们可以软件设置吸收边界来解决第二个问题。时域有限差分网格将在某处被截断。这要求在网格截断处不能引起波的明显反射。不完善的问题空间截断会引起数值反射，在计算空间经一定的模拟时间后会恶化计算结果。所以添加吸收边界及如何添加合适的吸收边界也已成为FDTD算法研究中的热点问题。

综上所述其程序流程图算法如图2。

2 基于matlab的FDTD算法仿真

基于上述算法流程图，采用matlab语言来进行FDTD算法的仿真。其部分代码如下：

%一维的情况

%电场只有Ex，磁场只有Hy

%每个网格用了2个时间步

%采用了简单的ABC；

%划分了200个网格；

clc；clear；

low1=0；

low2=0；

high1=0；

high2=0；

%figure（1）； %line（'erasemode'，'xor'）；

KE=200；

for i=1：KE； % 电场初始化

e（i）=0；

end；

for i=1：KE； % 磁场初始化

h（i）=0；

end；

npml=10；

for k=1：KE；

g2（k）=1；

g3（k）=1；

f2（k）=1；

f3（k）=1；

end

for k=1：npml；

xn=0.3*（（npml-k+1）/npml）^3；

g2（k）=1/（1+xn）；

g2（KE-k）=1/（1+xn）；

g3（k）=（1-xn）/（1+xn）；

g3（KE-k）=（1-xn）/（1+xn）；

xn=0.3*（（npml-k+1-0.5）/npml）^3；

f2（k）=1/（1+xn）；

f2（KE-k）=1/（1+xn）；

f3（k）=（1-xn）/（1+xn）；

f3（KE-k）=（1-xn）/（1+xn）；

end

for t=0：2300；

for k=2：KE；

e（k）=g3（k）*e（k）+ g2（k）*0.5*（h（k-1）-h（k））； % 电场迭代

end；

% e（100）=e（100）+exp（（-0.5）*（（40-t）/12）^2）； % 软源

% e（101）=exp（（-0.5）*（（40-t）/12）^2）； % 硬源

%e（101）=cos（0.2*t）； % 正弦源

% e（1）=low2； % 吸收边界

% low2=low1；

% low1=e（2）；

% e（201）=high2；

% high2=high1；

% high1=e（200）；

% e（1）=0； %金属边界

% e（201）=0；

for k=1：KE-1；

h（k）=f3（k）*h（k）+f2（k）*0.5*（e（k）-e（k+1））； % 磁場迭代

end；

% plot（h）；

plot（e）；

axis（[1 201 -1.0 1.0]）； % 吸收边界显示需要

drawnow；

pause（0.001）；

end；

部分运行结果如图3：（此结果由上述代码运行产生）

图3中间为正弦激励源，随着时间不断把电磁场向前推进，由于计算机存储空间限制，在两边设置吸收边界，抵达吸收边界的波形会被吸收掉。 （下转第238页）

（上接第226页）

3 CPU与GPU加速的选择

CPU的知名度远大于GPU，然而此次FDTD算法的加速却更加青睐GPU。众所周知，CPU（中央处理单元）是计算机的核心设备，用来完成各种复杂的计算任务，应用范围广泛。而GPU（图形处理单元）的工作目的只有一个，就是为计算机提供图像处理所需指令和数据。由于其执行任务单一，设计周期相对较短，技术突破十分迅速。不仅如此低成本高性能的GPU是一种高度并行的数据流处理器，显性结构的数据并行是GPU不在需要大量复杂的逻辑控制。数据流经过高速内存接口时不再需要大量缓存。这样就可以空出大量存储空间用于计算单元。对于解决FDTD这样需要大规模计算的算法十分有利。虽然GPU的设计主要是进行图像处理，并不是用于计算。但图像处理过程与矢量运算极为相似，为GPU进行电磁场矢量计算奠定了基础。下表给出两者之间的性能比较。

图4为GPU和CPU最大理论处理速度：（单位：/GFlops）

然而在众多的GPU中，进行FDTD算法研究的学者一直选择了与FDTD天然并行结构不谋而合的NVIDIA显卡。要想在这种类型上的GPU进行FDTD计算的编程，必须首先熟悉它的硬件操作语言CUDA。因为要在GPU上进行数学运算，需要先对GPU中的着色单元进行编程。CUDA语言是一门低级语言，对于掌握了C语言和Matlab这些高级语言的人来说低级语言要复杂的多，这也阻碍了GPU的进一步应用。总之，GPU在FDTD计算方面有着巨大的应用前景。

4 结论

本文简明扼要的介绍了FDTD算法的由来，基本原理，数值色散和吸收边界问题。并且给出了基于Matlab语言的仿真过程。进一步分析了GPU和CPU加速FDTD算法的优劣情况，指明了基于GPU加速的可行性。FDTD算法作为一种重要的电磁场计算方法，需要占据大量储存空间并且有这天然的并行结构。而GPU这个本来运用于图像处理的显卡，由于其并行数据流结构可以提供大量计算空间，使得计算速度相对于CPU有了极大提升。相信GPU还会在未来的FDTD计算道路上发挥自己更大的作用。

参考文献：

[1]葛德彪，闫玉波.电磁波时域有限差分方法[M].西安：西安电子科技大学出版社，2005.

[2]沈琛.基于GPU加速的FDTD算法对电磁辐射与散射问题的研究[D].安徽大学硕士学位论文，2010（04）.

[3]（美）AtefElsherbeni，VeyselDemir . Matlab模拟的电磁学时域有限差分法[M].北京：国防工业出版社，2013.

[4]杜文革.基于多GPU的FDTD并行算法及其在电磁仿真中的应用[D].山东大学博士学位论文，2011（03）.

[5]张波.基于图形处理器的时域有限差分算法研究[J].电波科学学报，2011（12）.

[6]马巍巍，孙冬.基于 GPU的高阶辛 FDTD 算法的并行仿真研究[J].合肥工业大学学报（自然科学版），2012（07）.

项目名称：安徽大学 大学生科研训练计划项目 “基于GPU的FDTD并行算法”。项目编号：KYXL2014073

作者简介：潘东旭，男，中共党员，安徽大学2012级电气工程系本科在读。