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

張立紅 余文華 楊小玲

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

引 言

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

1.理論分析

1.1 FDTD方法

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

圖1 電磁場值關系圖

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

1.2 SSE指令集

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

圖2 單指令多數據操作

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

2.數值實驗及結果

2.1 三級數據并行結構

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

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

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

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

2.2 SSE加速實現

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

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

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

3）把float類型的指針變量轉換成SSE所需的類型的指針變量。

4）把原來的最內層循環展開，循環次數變為原來的四分之一（這就是SSE指令集對FDTD仿真進行加速的原理）。

5）磁場值的遞推計算。

計算Hz的偽代碼如下:

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

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

2.3 數值實驗結果

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

表1 計算時間及加速比

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

3.結 論

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

［1］余文華，楊小玲，劉永俊，等.并行FDTD和IBM BlueGene／L巨型計算機結合求解電大尺寸的電磁問題［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算法實現研究［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，Fan 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］都志輝.高性能計算并行編程技術-MPI并行程序設計［M］.北京:清華大學出版社，2001.

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