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

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

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

1.引 言

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

圖1 CPU與GPU的計(jì)算性能對比

在計(jì)算能力突出的同時，GPU與CPU相比還兼具高度并行、多線程多核心等特點(diǎn)。電磁場數(shù)值計(jì)算若能有效利用GPU的計(jì)算能力及結(jié)構(gòu)特點(diǎn)，無疑會極大提高運(yùn)算速度和效率。

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

2. GPU加速FDTD的原理

2.1 GPU硬件架構(gòu)特點(diǎn)

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

2.2 GPU程序執(zhí)行特點(diǎn)

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

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

2.3 GPU與FDTD運(yùn)算結(jié)合的優(yōu)勢

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

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

3. GPU加速三維FDTD的計(jì)算模型

3.1 網(wǎng)格對應(yīng)線程的劃分及場量更新機(jī)制

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

3.2 FDTD算法的選擇

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

(1)

(2)

3.3 電磁參數(shù)的輸入與優(yōu)化

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

3.4 運(yùn)算結(jié)果的存儲輸出

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

綜上提出的GPU運(yùn)算流程模型見圖2。

圖2 GPU運(yùn)算流程模型

4.算例驗(yàn)證與分析

為驗(yàn)證提出的GPU加速流程模型，GPU加速運(yùn)算被運(yùn)用于實(shí)際工程問題的FDTD求解，并與CPU運(yùn)算進(jìn)行結(jié)果及性能比較。

4.1 算例與計(jì)算平臺

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

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

4.2 計(jì)算結(jié)果的精確性

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

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

圖4 功率反射系數(shù)對比圖

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

4.3 加速性能分析

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

表1 GPU與CPU計(jì)算耗時比

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

5.總結(jié)與展望

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

在GPU加速三維FDTD算法的可行性與高速性得到驗(yàn)證的同時，還有許多方面值得后繼深入研究，例如FDTD網(wǎng)格與線程的對應(yīng)是否有更高效的方式，GPU加速FDTD運(yùn)算的誤差分析與控制機(jī)理等等，這也是未來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的光波導(dǎo)器件FDTD并行算法研究[D]. 山東大學(xué), 2007

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

[5] 劉 瑜. FDTD算法的網(wǎng)絡(luò)并行研究及其電磁應(yīng)用[D]. 電子科技大學(xué), 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ì)算[J]. 電波科學(xué)學(xué)報(bào), 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ì)算電磁問題上的發(fā)展[J]. 電波科學(xué)學(xué)報(bào), 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.