基于圖形處理器的時域有限差分算法硬件加速

張 波 薛正輝 任 武 李偉明 盛新慶

(北京理工大學信息與電子學院，北京 100081)

1.引 言

長久以來，作為微型計算機的核心硬件之一，CPU一直承擔著絕大多數運算任務，也是傳統數值計算軟件面對的對象。而近年來，隨著用戶對于實時高解析率的圖像處理需求的急劇增加，原本僅負責圖形處理的圖形處理器(GPU)呈現出驚人的發展趨勢。發展至今，GPU的浮點處理能力及帶寬已全面超越同期CPU，圖1為CPU和使用開放圖形程序接口(OpenGL)與計算統一設備架構(CUDA)GPU的掃描性能對比，可見GPU的浮點運算能力相較同期CPU具有突出的優勢。

圖1 CPU與GPU的計算性能對比

在計算能力突出的同時，GPU與CPU相比還兼具高度并行、多線程多核心等特點。電磁場數值計算若能有效利用GPU的計算能力及結構特點，無疑會極大提高運算速度和效率。

隨著英偉達(Nvidia)公司在2002年提出可編程流處理器的概念以及在2006年推出CUDA運算架構，GPU運算能力的通用化在近年來取得了長足的發展。在計算電磁學領域，2004年Krakiwsky等人實現了GPU加速二維電磁場FDTD運算[1]，其后美國斯坦福大學和萊斯大學、加拿大卡爾加里大學的研究者廣泛開展了利用GPU進行電磁場數值計算的研究[2-3]。2007年山東大學韓林等人開展了利用GPU結合網絡并行運算技術的二維FDTD算法研究，以光波導器件為分析對象進行了探索[4]。2008年，電子科技大學劉瑜等人對GPU加速二維ADI-FDTD算法進行了研究[5]。同期的西南交通大學劉昆等人開展了GPU加速時域有限元的二維輻射計算研究[6-7]。然而迄今為止，GPU加速并未廣泛運用于FDTD運算中，原因在于工程計算多面向三維空間，而三維FDTD計算的加速相較二維情況，在算法優化等方面提出了諸多新要求，對加速效率的要求也更加嚴格。如何針對GPU硬件架構選擇適宜的FDTD算法進行優化，并且將三維空間中的激勵源引入和吸收邊界設置等流程高效地在GPU上得以實現，從而使得GPU加速FDTD運算真正進入工程應用階段，是本文研究的問題。

2. GPU加速FDTD的原理

2.1 GPU硬件架構特點

GPU架構的突出特點是采用高度并行化的流處理器群作為運算單元，以Nvidia GT130M GPU為例，該GPU包含32個流處理器，以8個一組為單位組成4個流處理器群，處理器核心頻率為600 MHz.與CPU運算體系中內存對應，GPU運算主要存儲器為顯存，按功能分為可被全部流處理器讀寫的全局顯存、可被同一線程塊內所有線程讀寫的線程共享顯存以及單線程獨享的線程顯存。

2.2 GPU程序執行特點

與硬件架構相適應，GPU程序執行方式相比CPU程序并行化程度更高，更深入硬件層次。具體執行方式是以線程(Thread)為單位并行執行開發者編寫的“Kernel”函數。詳細流程是開發者編寫執行計算任務的Kernel函數，規定函數執行的總線程數以及劃分出線程網格(Grid)，線程網格會被細分為線程塊(Block)后輸送給GPU.GPU接到線程塊以及Kernel函數后，自動將每個線程塊以32個線程為一個單位(即一個Warp)進行劃分，分配給多個流處理器群并行執行。

盡管GPU同時運行的線程數受其硬件能力制約，但Kernel函數的并發總線程數卻可以遠超出這一限制，這得益于GPU的Warp執行機制，即多個Warp可以輪流運行于多個流處理器群，甚至當某個流處理器群中的Warp處于等待狀態時，它可以臨時運行其他閑置的Warp.多線程多層面的并行機制極大提高了GPU的硬件利用率，同時GPU還提供了同一線程塊中的線程間同步功能，可以有效地進行程序流程控制。

2.3 GPU與FDTD運算結合的優勢

FDTD算法的場量更新機制非常適宜并行化。以FDTD算法中電場更新為例，求解某網格第n+1步電場E時，需讀取的場值僅包括該網格第n步電場E及第n+1/2步時相鄰網格的磁場H，這意味著場量更新過程對于各網格的空間更新次序不存在任何要求。這使得編寫Kernel程序以及分割線程網格時，可以大幅偏重于線程執行效率的提高，而非保證線程的執行次序。

FDTD算法的算術指令符合流處理器的運算能力。GPU相較CPU而言，其運算能力的優勢來源于多流處理器設計以及并行化的線程架構，就單個流處理器而言，它的單精度運算能力及精度并未達到同級別CPU的高度。相較其他電磁場數值算法，FDTD算法不采用矩陣運算，基本不涉及求冪求積分等計算，對GPU運算效率的提高和誤差控制無疑十分有利。

3. GPU加速三維FDTD的計算模型

3.1 網格對應線程的劃分及場量更新機制

以使用Nvidia Gt130M顯卡計算56×56×120尺寸空間為例，由于顯卡只有32個流處理器，4個流處理器群，而空間網格共有376320個，因此為每個網格分配一個線程既無必要也影響效率。在實際計算時，采取了一種“XY平面各點并行，Z向循環推進更新”的模式：將線程塊大小設置為32×8，即Gt130M顯卡所能支持的最大線程塊大小，將XY平面的XY向網格數除以線程塊對應維度大小并向上取整，得到線程網格大小為2×7，至此XY平面網格已一一映射為線程。在計算時，首先將所有線程指針指向Z=1時XY平面各點對應存儲空間，在所有線程運算完成后，將每個線程指針指向Z向下一個網格進行新一輪計算，循環往復直至全空間更新完成。這種更新模式在提高多Warp執行機制效率的同時節省了運算資源。由于實際參與迭代運算的總線程數(32×8×2×7=3584)大于XY平面的網格數(56×56=3136)，每個線程在進行場量運算時，需要判斷其對應網格是否處于求解范圍內，若不處于求解范圍內則不參與運算，否則顯存讀取地址會發生沖突進而導致計算結果出錯。

3.2 FDTD算法的選擇

對于并行FDTD算法，所有網格的場量更新采取同一類型運算公式有利于提高程序運行效率和進行誤差控制分析。因此算法模型選擇完全匹配層(UPML)作為吸收邊界，并在包括UPML層的全空間內統一采取Taflove所提出的兩步迭代法進行運算[8]，以Hx為例的場量更新公式見下式。

(1)

(2)

3.3 電磁參數的輸入與優化

對于入射波的引入采取了總場/散射場法，當線程判斷到其對應網格處于連接邊界上時，會進行入射波引入或消去處理，這與CPU計算類似，因此不再贅述。

3.4 運算結果的存儲輸出

由于顯存與內存之間的數據傳輸會占用可觀的運算資源與時間，而FDTD計算結果往往只需要經空間與時間采樣后的場量數據。因此GPU運算時盡量將采樣后的數據存放于顯存中，等待全部計算過程完畢后再進行輸出，若采樣時間點過多或采樣面過大，會依據可用顯存空間進行分次存儲和輸出，總體原則遵循在顯存容量允許的情況下，盡可能減少GPU對內存空間的訪問。

綜上提出的GPU運算流程模型見圖2。

圖2 GPU運算流程模型

4.算例驗證與分析

為驗證提出的GPU加速流程模型，GPU加速運算被運用于實際工程問題的FDTD求解，并與CPU運算進行結果及性能比較。

4.1 算例與計算平臺

算例模型為一個由尺寸為1.27 mm×12.7 mm，x向與y向單元間距均為17.8 mm的無限薄金屬振子單元組成的頻率選擇表面。頻率選擇表面屬于周期性結構，這里選用譜FDTD[9]配合周期性邊界的方法，通過對一個單元建模計算進行分析。空間步長設置為0.318 mm，網格大小為56×56×120，x向與y向網格四周采取周期性邊界條件進行處理，z向網格兩端各設置10層UPML層作為吸收邊界，單元模型放置于z=45截面中心并在z=75截面通過連接邊界引入入射波。時間步長設置為5.295×10-13s，總共計算8192步。在計算過程中對z=85截面上場值進行時間采樣，在全部時間步計算完成后，將采樣結果通過傅里葉變換轉換到頻域并通過Poyinting定理求得橫截面的功率函數，與入射波的對應功率函數相比后得到頻率選擇表面的功率反射系數。

計算過程采用兩套GPU運算平臺，GPU平臺1為Nvidia GT130M圖形處理器與512 MB顯存，GPU平臺2為Tesla C2050圖形處理器與3 GB顯存，CPU對比平臺為Intel Core2 T6500處理器與3 GB內存，CPU與GPU運算的模型參數以及運算參數完全一致。

4.2 計算結果的精確性

為驗證GPU加速的數值精確性，運算過程中抽取時間步t=1000時，y=30截面上Ey值進行輸出比較，雙GPU平臺數值完全吻合，GPU平臺計算結果與CPU平臺符合很好如圖3所示。

在全部時間步計算完成后，分別將GPU與CPU計算的結果進行后期處理，求出對應的功率反射系數，并與商業軟件CST 2010計算所得結果進行比較。雙GPU平臺數值完全吻合，GPU平臺與CPU平臺運算結果對比見圖4。

圖4 功率反射系數對比圖

由圖4可見，GPU與CPU的計算結果吻合程度極佳，與CST軟件所得結果也符合較好。通過后續計算比較，GPU計算所得功率反射系數與CPU計算結果之間差異僅為0.44%左右，足以滿足絕大多數情況下工程需要。

4.3 加速性能分析

在計算過程中，將CPU平臺與雙GPU平臺計算至同一時間步時所耗時間記錄于表1。

表1 GPU與CPU計算耗時比

由表1可見，在整個運算過程中GPU平臺的運算性能保持穩定，GPU平臺1相對CPU平臺加速比穩定在23倍以上，GPU平臺2由于科學計算專用GPU的使用，加速比達到了174倍以上，有力證明了GPU加速FDTD運算的高效性。

5.總結與展望

通過對GPU加速FDTD算法的原理探討與算例分析可以看到，計算流程中網格與線程的映射、算法的選擇與優化、模型與入射波的輸入、運算結果的存儲等方面都在本文所提出的加速流程模型中得以高效實現。遵循這一流程的GPU加速FDTD運算，在滿足運算精度需要的前提下，相較傳統CPU運算，其運算速度大幅提高，不僅可以勝任大規模工程計算，適應性也極為優秀。在計算平臺成本方面，算例中GPU平臺1的芯片只是Nvidia顯卡的中低端型號，目前主流Nvidia顯卡都已經具備CUDA運算功能，這無疑大大降低了GPU加速FDTD的門檻，當需要進行大規模高性能科學計算時，可以使用GPU平臺2中的專業級GPU芯片，相應加速比也更加優秀。

在GPU加速三維FDTD算法的可行性與高速性得到驗證的同時，還有許多方面值得后繼深入研究，例如FDTD網格與線程的對應是否有更高效的方式，GPU加速FDTD運算的誤差分析與控制機理等等，這也是未來GPU加速FDTD研究的方向。

[1] KRAKIWSKY S E, TUMER L E, OKONIEWSKI M M. Acceleration of finite-difference time-domain (FDTD) using graphics processor units (GPU)[C]//IEEE MTTS International Microwave Symposium Digest, 2004, 2: 1033-1036.

[2] STEFANSKI T P, DRYSDALE T D. Acceleration of the 3D ADI-FDTD method using graphics processor units[J]. IEEE MTT-S International Microwave Symposium Digest (MTT), 2009, 1: 241-244.

[3] PRICE D K, HUMPHREY J R, KELMELIS E J. GPU-Based Accelerated 2D and 3D FDTD Solvers[J]. Physics and Simulation of Optoelectronic Devices XV, 2007, 6468(1): 22-25.

[4] 韓 林. 基于GPU的光波導器件FDTD并行算法研究[D]. 山東大學, 2007

LIN Han. GPU Based Optical Waveguide FDTD Parallel Research[D]. Shandong University, 2007. (in Chinese)

[5] 劉 瑜. FDTD算法的網絡并行研究及其電磁應用[D]. 電子科技大學, 2008

LIU Yu. Parallel FDTD Algorithm Based on Network and Applications in Electromagnetic Problems[D]. University of Electronic Science and Technology of China, 2008. (in Chinese)

[6] 劉 昆, 王曉斌, 廖 成. 圖形處理器(GPU)加速時域有限元的二維輻射計算[J]. 電波科學學報, 2008, 23(1): 111-114.

LIU Kun, WANG Xiaobin, LIAO Cheng. Acceleration of time-domain finite element 2-D radiation using graphics processor units(GPU)[J]. Chinese Journal of Radio Science, 2008, 23(1): 111-114. (in Chinese)

[7] 吳 霞, 周樂柱. 時域有限元法在計算電磁問題上的發展[J]. 電波科學學報, 2008, 23(6): 1208-1216.

WU Xia, ZHOU Lezhu. Application and development of time-domain finite element method on EM analysis[J]. Chinese Journal of Radio Science, 2008, 23(1): 1208-1216. (in Chinese)

[8] TAFLOVE A. Computational Electrodynamics: The Finite-Difference Time-Domain Method Third Edition[M]. Norwood MA: Artech House, 2005.

[9]AMINIAN A and RAHMAT-SAMII Y. Spectral FDTD: a novel technique for the analysis of oblique incident plane wave on periodic structures[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(6): 1818-1525.