加速并行時(shí)域有限差分仿真的新方法

張立紅 余文華 楊小玲

（1.中國(guó)傳媒大學(xué)信息工程學(xué)院，北京 100024；2.中國(guó)人民武裝警察部隊(duì)學(xué)院，河北 廊坊 065000；3.Penn State University，USA PA 16802）

引 言

時(shí)域有限差分（FDTD）法最早由K.S.Yee在1966年提出，經(jīng)過(guò)幾十年的發(fā)展，F(xiàn)DTD已經(jīng)形成了一套比較完善的方法體系，相對(duì)于其他的計(jì)算電磁學(xué)方法，F(xiàn)DTD因其簡(jiǎn)單靈活而受到廣大電磁計(jì)算研究者的歡迎，但它的實(shí)現(xiàn)卻面臨著一些問(wèn)題，如龐大的計(jì)算量是普通PC機(jī)所不能滿足的，因此，科學(xué)家、電磁工作者等提出了各種各樣的并行算法來(lái)解決這些問(wèn)題，如基于消息傳遞接口（MPI）的并行技術(shù)、基于OpenMP的共享存儲(chǔ)編程技術(shù)以及基于映射文件的技術(shù)等［1-4］。最近，也有一些文獻(xiàn)提出用圖形處理單元（GPU）對(duì)FDTD算法進(jìn)行加速［5］。文章提出了一種利用單指令多數(shù)據(jù)流式擴(kuò)展（SSE）指令集來(lái)加速并行FDTD仿真的新方法，用C語(yǔ)言開(kāi)發(fā)了基于MPI庫(kù)、OpenMP和SSE指令集的三維并行FDTD代碼，最后以具體的電磁仿真實(shí)例驗(yàn)證了新方法的可行性和加速效率，并將其與普通并行FDTD仿真方法進(jìn)行了對(duì)比。

1.理論分析

1.1 FDTD方法

在FDTD方法中，電磁波傳播以及電磁波與物質(zhì)的相互作用是通過(guò)電場(chǎng)和磁場(chǎng)在空間和時(shí)間上的差分遞推實(shí)現(xiàn)的，空間某處的電場(chǎng)值可以由該處上一時(shí)間步的電場(chǎng)值和其周圍上半個(gè)時(shí)間步的四個(gè)磁場(chǎng)值計(jì)算得到，而空間某處的磁場(chǎng)值可以由該處上一時(shí)間步的磁場(chǎng)值和其周圍上半個(gè)時(shí)間步的四個(gè)電場(chǎng)值計(jì)算得到。在FDTD方法中，電磁場(chǎng)值的位置和遞推關(guān)系可以用式（1）和圖1表示。公式式（1）表示的是磁場(chǎng)沿z軸方向分量的遞推公式，其他兩個(gè)方向的分量以及電場(chǎng)分量的遞推公式與（1）完全相似［6］。

圖1 電磁場(chǎng)值關(guān)系圖

從遞推公式（1）和圖1可以看出，F(xiàn)DTD方法具有與生俱來(lái)的并行性:FDTD方法中，每一個(gè)網(wǎng)格點(diǎn)的磁場(chǎng)（電場(chǎng)）分量的迭代公式只與它自己上一時(shí)間步的值和它周圍網(wǎng)格點(diǎn)電場(chǎng)（磁場(chǎng)）上半個(gè)時(shí)間步的值有關(guān)，而與計(jì)算區(qū)域內(nèi)其他網(wǎng)格點(diǎn)的場(chǎng)值沒(méi)有直接關(guān)系，非常適合并行計(jì)算［7］；而且，遞推公式中，幾乎所有的計(jì)算都是對(duì)一組數(shù)據(jù)進(jìn)行相同的加、減或乘法操作，非常適合單指令多數(shù)據(jù)（SIMD）模式的并行處理。

1.2 SSE指令集

SSE指令集是Intel在其芯片中實(shí)現(xiàn)了基于寄存器的SIMD架構(gòu)之后提供的指令集，它使用8個(gè)獨(dú)立的128位寄存器，允許SIMD計(jì)算同時(shí)作用于4個(gè)緊縮的單精度浮點(diǎn)數(shù)據(jù)單元。SSE指令集包括70條指令，其中包含提高3D圖形運(yùn)算效率的50條SIMD浮點(diǎn)運(yùn)算指令、12條多媒體擴(kuò)展（MMX）整數(shù)運(yùn)算增強(qiáng)指令和8條優(yōu)化內(nèi)存中連續(xù)數(shù)據(jù)塊傳輸指令［8］。AMD處理器也加入了對(duì)SSE指令集的支持?，F(xiàn)在市場(chǎng)上能夠買到的處理器大都支持SSE指令集。

圖2 單指令多數(shù)據(jù)操作

因?yàn)镾SE指令集是單指令多數(shù)據(jù)操作，因此，可以通過(guò)循環(huán)展開(kāi)來(lái)減少運(yùn)算時(shí)間，從而提高運(yùn)算速度。SSE指令集要求它的操作數(shù)是一種新的緊縮類型，對(duì)于float類型的數(shù)據(jù)，SSE指令集的操作數(shù)是把4個(gè)float標(biāo)量數(shù)據(jù)壓縮成一個(gè)類型的向量數(shù)據(jù)。圖2是一個(gè)典型的單指令多數(shù)據(jù)的操作。圖2中的a和b都是類型的向量數(shù)據(jù)，都是由4個(gè)float類型的標(biāo)量數(shù)據(jù)壓縮而成的，SSE指令對(duì)a和b進(jìn)行加法運(yùn)算操作，得到一個(gè)類型的向量數(shù)據(jù)，存放在Result中。

2.數(shù)值實(shí)驗(yàn)及結(jié)果

2.1 三級(jí)數(shù)據(jù)并行結(jié)構(gòu)

普通的基于MPI或OpenMP的并行FDTD算法都是一級(jí)或兩級(jí)并行結(jié)構(gòu)［9］，把基于SSE指令集的新加速方法引入到FDTD仿真后，程序形成了三級(jí)數(shù)據(jù)并行結(jié)構(gòu)。

第一級(jí)數(shù)據(jù)并行基于MPI庫(kù)。把要進(jìn)行仿真的區(qū)域進(jìn)行區(qū)域分解，各子區(qū)域的電磁場(chǎng)值分別獨(dú)立計(jì)算，負(fù)責(zé)計(jì)算各子域的進(jìn)程之間通過(guò)MPI庫(kù)的消息傳遞函數(shù)進(jìn)行通信。

第二級(jí)數(shù)據(jù)并行采用OpenMP共享存儲(chǔ)編程實(shí)現(xiàn)。首先利用OpenMP生成多個(gè)線程，然后將每個(gè)子區(qū)域的計(jì)算再分配給各個(gè)線程并行執(zhí)行，算法實(shí)現(xiàn)框架如下:

第三級(jí)數(shù)據(jù)并行利用SSE指令集實(shí)現(xiàn)。對(duì)于單精度浮點(diǎn)運(yùn)算，普通的運(yùn)算操作一次得到一個(gè)計(jì)算結(jié)果，而使用SSE指令集，一個(gè)運(yùn)算可以得到四個(gè)計(jì)算結(jié)果，從而實(shí)現(xiàn)細(xì)粒度數(shù)據(jù)并行，加快了計(jì)算速度。

2.2 SSE加速實(shí)現(xiàn)

在利用SSE加速并行FDTD算法時(shí)，僅對(duì)電磁場(chǎng)的遞推部分進(jìn)行了加速，其中包括整個(gè)計(jì)算區(qū)域的電磁場(chǎng)的遞推過(guò)程和卷積完全匹配層（CPML）吸收邊界［10］的處理過(guò)程。先討論整個(gè)計(jì)算區(qū)域的電磁場(chǎng)的遞推情況。以計(jì)算磁場(chǎng)沿z軸方向的分量Hz為例，可按照以下步驟實(shí)現(xiàn):

1）定義SSE所需的類型的變量，并為其賦值（作為SSE指令的操作數(shù)）。

2）把電磁場(chǎng)遞推公式中所需的系數(shù)加載到SSE寄存器中。

3）把float類型的指針變量轉(zhuǎn)換成SSE所需的類型的指針變量。

4）把原來(lái)的最內(nèi)層循環(huán)展開(kāi)，循環(huán)次數(shù)變?yōu)樵瓉?lái)的四分之一（這就是SSE指令集對(duì)FDTD仿真進(jìn)行加速的原理）。

5）磁場(chǎng)值的遞推計(jì)算。

計(jì)算Hz的偽代碼如下:

CPML吸收邊界的處理方法與電磁場(chǎng)值遞推過(guò)程類似，例如，在計(jì)算沿y軸方向的CPML區(qū)域的磁場(chǎng)時(shí)，可以參照前面的步驟，實(shí)現(xiàn)偽代碼如下:

為了優(yōu)化程序，提高緩存命中率，還可以把計(jì)算電磁場(chǎng)值的基本遞推過(guò)程和CPML吸收邊界的處理過(guò)程合并起來(lái)，通過(guò)判斷j的取值是否在CPML吸收邊界區(qū)域內(nèi)來(lái)決定是否進(jìn)行電磁場(chǎng)值邊界的更新，實(shí)現(xiàn)框架如下:

2.3 數(shù)值實(shí)驗(yàn)結(jié)果

為了驗(yàn)證新方法的加速效率，文章進(jìn)行了實(shí)驗(yàn)測(cè)試，分別計(jì)算了40×40×40、80×80×80和160×160×160個(gè)均勻網(wǎng)格的真空中電磁波的傳播，其中，激勵(lì)源為高斯脈沖源，放置在立方體計(jì)算區(qū)域的正中心，電磁場(chǎng)初始值均設(shè)為0，CPML吸收邊界為6層。實(shí)驗(yàn)平臺(tái)是PC機(jī)，CPU為Intel的T2300（雙核），1.66GHz，時(shí)間步為400，實(shí)驗(yàn)結(jié)果如表1所示。從測(cè)試結(jié)果可以看出，使用了SSE指令集加速的代碼比普通的并行代碼所需的計(jì)算時(shí)間大大減少。

表1 計(jì)算時(shí)間及加速比

普通并行代碼一個(gè)指令進(jìn)行一次運(yùn)算操作，得到一個(gè)結(jié)果值，而SSE代碼一個(gè)指令進(jìn)行一次運(yùn)算，得到四個(gè)結(jié)果值，因此，理論上使用SSE指令集加速時(shí)，最理想情況是加速比等于4，實(shí)驗(yàn)在160×160×160均勻網(wǎng)格時(shí)間步為400的情況下，得到的加速比為2.62，加速效果比較好。

3.結(jié) 論

提出了利用SSE指令集來(lái)加速并行FDTD算法的方法，開(kāi)發(fā)了三級(jí)數(shù)據(jù)并行結(jié)構(gòu)的三維FDTD仿真代碼，在Intel T2300的PC機(jī)上實(shí)現(xiàn)了對(duì)基于MPI和OpenMP的三維并行FDTD仿真的加速，得到的加速比為2.62。使用SSE指令集加速，無(wú)需額外購(gòu)買任何硬件，只需改變部分并行代碼即可實(shí)現(xiàn)，因此，使用SSE指令集來(lái)加快運(yùn)算速度，從而減少運(yùn)行時(shí)間是一種高效、經(jīng)濟(jì)的新途徑。

［1］余文華，楊小玲，劉永俊，等.并行FDTD和IBM BlueGene／L巨型計(jì)算機(jī)結(jié)合求解電大尺寸的電磁問(wèn)題［J］.電波科學(xué)學(xué)報(bào)，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］雷繼兆，梁昌洪，丁 偉，等.機(jī)載天線輻射特性的并行FDTD分析［J］.電波科學(xué)學(xué)報(bào)，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］劉 瑜，梁 正，楊梓強(qiáng).基于映射文件的電磁并行FDTD算法實(shí)現(xiàn)研究［J］.電波科學(xué)學(xué)報(bào)，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，F(xiàn)an 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］葛德彪，閆玉波.電磁場(chǎng)時(shí)域有限差分方法［M］.西安電子科技大學(xué)出版社，2002.

［7］余文華，蘇 濤，Raj Mittra，等.并行時(shí)域有限差分［M］.北京:中國(guó)傳媒大學(xué)出版社，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］都志輝.高性能計(jì)算并行編程技術(shù)-MPI并行程序設(shè)計(jì)［M］.北京:清華大學(xué)出版社，2001.

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