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

摘 要：時(shí)域有限差分法，即FDTD（Finite Difference Time Domain），是計(jì)算電磁學(xué)的一種重要方法。作為一種天然的并行算法，它的計(jì)算過程可以劃分為多個(gè)同時(shí)進(jìn)行相似計(jì)算的子計(jì)算。這個(gè)方法主要是把麥克斯韋方程在時(shí)間上和空間上進(jìn)行差分化，并且通過時(shí)間領(lǐng)域上的更新來模仿電磁場(chǎng)的變化來計(jì)算問題，因而有利于解決很多電磁場(chǎng)問題。而圖形處理單元即GPU（Graphic Processing Unit）相對(duì)于CPU的高性能計(jì)算速度以及NVIDA公司生產(chǎn)的GPU特有的高并行結(jié)構(gòu)，為時(shí)域有限差分的加速提供了可能。

關(guān)鍵字：時(shí)域有限差分法；圖形處理單元；麥克斯韋方程；并行算法

1 FDTD的基本原理

FDTD算法是1966年K.S.Yee發(fā)表在AP上的一篇論文建立起來的，后被稱為Yee網(wǎng)格空間離散方式。核心思想是把帶時(shí)間變量的麥克斯韋旋度方程轉(zhuǎn)化為差分形式，故中文稱之為“時(shí)有限差分法”。麥克斯韋方程如下：

其中H是磁場(chǎng)強(qiáng)度，E是電場(chǎng)強(qiáng)度，D是電位移，B是磁通量密度。

上述兩個(gè)矢量方程描述了麥克斯韋方程中磁場(chǎng)與電場(chǎng)復(fù)雜交錯(cuò)的關(guān)系。由于在三維空間中每個(gè)矢量方程又可以分解為三個(gè)標(biāo)量方程，因此該方程組可以化為6個(gè)標(biāo)量方程如下所示：

這樣的話便將問題的幾何空間離散為空間網(wǎng)格，電場(chǎng)和磁場(chǎng)的分量便被置于空間離散的網(wǎng)格點(diǎn)上。而這是FDTD計(jì)算的前提。我們?cè)儆貌罘纸铺娲溈怂鬼f這6個(gè)標(biāo)量方程中的時(shí)間和空間導(dǎo)數(shù)，構(gòu)造一系列方程，均以前一時(shí)間步電磁場(chǎng)瞬時(shí)值來“預(yù)測(cè)”后一時(shí)間步電磁場(chǎng)的瞬時(shí)值，由此構(gòu)造時(shí)間不斷向前推近的算法，來模擬時(shí)域中的電磁場(chǎng)變化過程。

1966年，Yee首次給出了麥克斯韋旋度方程的一組差分形式，這組方程在空間和時(shí)間上一離散的形式給出，使用的是中心差分法。對(duì)空間（X軸方向）的中心差分法離散公式如下：

對(duì)Y軸，Z軸方向的中心差分離散以此類推。

對(duì)時(shí)間的中心差分離散公示如下：

圖1為Yee單元網(wǎng)格的結(jié)構(gòu)：

我們根據(jù)6個(gè)麥克斯韋標(biāo)量方程，結(jié)合上述中心差分法，便可以得到三維問題下的FDTD更新方程如下：

其余的五個(gè)方程也如法可以寫出，因此任何時(shí)刻可一次算出一個(gè)點(diǎn)，并行算法可計(jì)算出多個(gè)點(diǎn)。通過這些運(yùn)算可以交替算出電場(chǎng)磁場(chǎng)在各個(gè)時(shí)間步的值。

上述方法運(yùn)行時(shí)忽略了兩個(gè)很重要的問題，一是數(shù)值色散，就是該方法會(huì)產(chǎn)生一定誤差，如何維持該方法可行的精度問題，二是計(jì)算機(jī)中模擬的無(wú)限過程與計(jì)算機(jī)本身儲(chǔ)存內(nèi)存有限相矛盾的問題。

解決數(shù)值色散的關(guān)鍵是空間網(wǎng)格大小的選擇，時(shí)間步長(zhǎng)△t，空間步長(zhǎng)△x， △y， △z必須滿足一定的關(guān)系，否則就使得數(shù)值表現(xiàn)不穩(wěn)定，表現(xiàn)為：隨著計(jì)算步數(shù)的增加，計(jì)算場(chǎng)量的數(shù)值會(huì)無(wú)限的增大，這種增大不是由于誤差積累造成的，而是由于電磁波的傳播關(guān)系被破壞造成的。所以△t， △x， △y， △z必須滿足一定的關(guān)系以保證穩(wěn)定性，其數(shù)值穩(wěn)定條件為如下：

我們可以軟件設(shè)置吸收邊界來解決第二個(gè)問題。時(shí)域有限差分網(wǎng)格將在某處被截?cái)唷＿@要求在網(wǎng)格截?cái)嗵幉荒芤鸩ǖ拿黠@反射。不完善的問題空間截?cái)鄷?huì)引起數(shù)值反射，在計(jì)算空間經(jīng)一定的模擬時(shí)間后會(huì)惡化計(jì)算結(jié)果。所以添加吸收邊界及如何添加合適的吸收邊界也已成為FDTD算法研究中的熱點(diǎn)問題。

綜上所述其程序流程圖算法如圖2。

2 基于matlab的FDTD算法仿真

基于上述算法流程圖，采用matlab語(yǔ)言來進(jìn)行FDTD算法的仿真。其部分代碼如下：

%一維的情況

%電場(chǎng)只有Ex，磁場(chǎng)只有Hy

%每個(gè)網(wǎng)格用了2個(gè)時(shí)間步

%采用了簡(jiǎn)單的ABC；

%劃分了200個(gè)網(wǎng)格；

clc；clear；

low1=0；

low2=0；

high1=0；

high2=0；

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

KE=200；

for i=1：KE； % 電場(chǎng)初始化

e（i）=0；

end；

for i=1：KE； % 磁場(chǎng)初始化

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））； % 電場(chǎng)迭代

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））； % 磁場(chǎng)迭代

end；

% plot（h）；

plot（e）；

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

drawnow；

pause（0.001）；

end；

部分運(yùn)行結(jié)果如圖3：（此結(jié)果由上述代碼運(yùn)行產(chǎn)生）

圖3中間為正弦激勵(lì)源，隨著時(shí)間不斷把電磁場(chǎng)向前推進(jìn)，由于計(jì)算機(jī)存儲(chǔ)空間限制，在兩邊設(shè)置吸收邊界，抵達(dá)吸收邊界的波形會(huì)被吸收掉。 （下轉(zhuǎn)第238頁(yè)）

（上接第226頁(yè)）

3 CPU與GPU加速的選擇

CPU的知名度遠(yuǎn)大于GPU，然而此次FDTD算法的加速卻更加青睞GPU。眾所周知，CPU（中央處理單元）是計(jì)算機(jī)的核心設(shè)備，用來完成各種復(fù)雜的計(jì)算任務(wù)，應(yīng)用范圍廣泛。而GPU（圖形處理單元）的工作目的只有一個(gè)，就是為計(jì)算機(jī)提供圖像處理所需指令和數(shù)據(jù)。由于其執(zhí)行任務(wù)單一，設(shè)計(jì)周期相對(duì)較短，技術(shù)突破十分迅速。不僅如此低成本高性能的GPU是一種高度并行的數(shù)據(jù)流處理器，顯性結(jié)構(gòu)的數(shù)據(jù)并行是GPU不在需要大量復(fù)雜的邏輯控制。數(shù)據(jù)流經(jīng)過高速內(nèi)存接口時(shí)不再需要大量緩存。這樣就可以空出大量存儲(chǔ)空間用于計(jì)算單元。對(duì)于解決FDTD這樣需要大規(guī)模計(jì)算的算法十分有利。雖然GPU的設(shè)計(jì)主要是進(jìn)行圖像處理，并不是用于計(jì)算。但圖像處理過程與矢量運(yùn)算極為相似，為GPU進(jìn)行電磁場(chǎng)矢量計(jì)算奠定了基礎(chǔ)。下表給出兩者之間的性能比較。

圖4為GPU和CPU最大理論處理速度：（單位：/GFlops）

然而在眾多的GPU中，進(jìn)行FDTD算法研究的學(xué)者一直選擇了與FDTD天然并行結(jié)構(gòu)不謀而合的NVIDIA顯卡。要想在這種類型上的GPU進(jìn)行FDTD計(jì)算的編程，必須首先熟悉它的硬件操作語(yǔ)言CUDA。因?yàn)橐贕PU上進(jìn)行數(shù)學(xué)運(yùn)算，需要先對(duì)GPU中的著色單元進(jìn)行編程。CUDA語(yǔ)言是一門低級(jí)語(yǔ)言，對(duì)于掌握了C語(yǔ)言和Matlab這些高級(jí)語(yǔ)言的人來說低級(jí)語(yǔ)言要復(fù)雜的多，這也阻礙了GPU的進(jìn)一步應(yīng)用。總之，GPU在FDTD計(jì)算方面有著巨大的應(yīng)用前景。

4 結(jié)論

本文簡(jiǎn)明扼要的介紹了FDTD算法的由來，基本原理，數(shù)值色散和吸收邊界問題。并且給出了基于Matlab語(yǔ)言的仿真過程。進(jìn)一步分析了GPU和CPU加速FDTD算法的優(yōu)劣情況，指明了基于GPU加速的可行性。FDTD算法作為一種重要的電磁場(chǎng)計(jì)算方法，需要占據(jù)大量?jī)?chǔ)存空間并且有這天然的并行結(jié)構(gòu)。而GPU這個(gè)本來運(yùn)用于圖像處理的顯卡，由于其并行數(shù)據(jù)流結(jié)構(gòu)可以提供大量計(jì)算空間，使得計(jì)算速度相對(duì)于CPU有了極大提升。相信GPU還會(huì)在未來的FDTD計(jì)算道路上發(fā)揮自己更大的作用。

參考文獻(xiàn)：

[1]葛德彪，閆玉波.電磁波時(shí)域有限差分方法[M].西安：西安電子科技大學(xué)出版社，2005.

[2]沈琛.基于GPU加速的FDTD算法對(duì)電磁輻射與散射問題的研究[D].安徽大學(xué)碩士學(xué)位論文，2010（04）.

[3]（美）AtefElsherbeni，VeyselDemir . Matlab模擬的電磁學(xué)時(shí)域有限差分法[M].北京：國(guó)防工業(yè)出版社，2013.

[4]杜文革.基于多GPU的FDTD并行算法及其在電磁仿真中的應(yīng)用[D].山東大學(xué)博士學(xué)位論文，2011（03）.

[5]張波.基于圖形處理器的時(shí)域有限差分算法研究[J].電波科學(xué)學(xué)報(bào)，2011（12）.

[6]馬巍巍，孫冬.基于 GPU的高階辛 FDTD 算法的并行仿真研究[J].合肥工業(yè)大學(xué)學(xué)報(bào)（自然科學(xué)版），2012（07）.

項(xiàng)目名稱：安徽大學(xué) 大學(xué)生科研訓(xùn)練計(jì)劃項(xiàng)目 “基于GPU的FDTD并行算法”。項(xiàng)目編號(hào)：KYXL2014073

作者簡(jiǎn)介：潘東旭，男，中共黨員，安徽大學(xué)2012級(jí)電氣工程系本科在讀。