時域有限差分法的圖形處理單元的加速

摘 要：時域有限差分法，即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級電氣工程系本科在讀。