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

張立紅 余文華 楊小玲

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

引 言

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

1.理論分析

1.1 FDTD方法

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

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

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

1.2 SSE指令集

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

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

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

2.數(shù)值實驗及結(jié)果

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

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

第一級數(shù)據(jù)并行基于MPI庫。把要進行仿真的區(qū)域進行區(qū)域分解，各子區(qū)域的電磁場值分別獨立計算，負責計算各子域的進程之間通過MPI庫的消息傳遞函數(shù)進行通信。

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

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

2.2 SSE加速實現(xiàn)

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

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

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

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

4）把原來的最內(nèi)層循環(huán)展開，循環(huán)次數(shù)變?yōu)樵瓉淼乃姆种唬ㄟ@就是SSE指令集對FDTD仿真進行加速的原理）。

5）磁場值的遞推計算。

計算Hz的偽代碼如下:

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

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

2.3 數(shù)值實驗結(jié)果

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

表1 計算時間及加速比

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

3.結(jié) 論

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

［1］余文華，楊小玲，劉永俊，等.并行FDTD和IBM BlueGene／L巨型計算機結(jié)合求解電大尺寸的電磁問題［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算法實現(xiàn)研究［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，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］葛德彪，閆玉波.電磁場時域有限差分方法［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］都志輝.高性能計算并行編程技術(shù)-MPI并行程序設(shè)計［M］.北京:清華大學出版社，2001.

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