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

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

1.引 言

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

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

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

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

2. GPU加速FDTD的原理

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

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

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

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

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

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

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

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

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

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

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

3.2 FDTD算法的選擇

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

(1)

(2)

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

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

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

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

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

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

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

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

4.1 算例與計(jì)算平臺(tái)

算例模型為一個(gè)由尺寸為1.27 mm×12.7 mm，x向與y向單元間距均為17.8 mm的無(wú)限薄金屬振子單元組成的頻率選擇表面。頻率選擇表面屬于周期性結(jié)構(gòu)，這里選用譜FDTD[9]配合周期性邊界的方法，通過(guò)對(duì)一個(gè)單元建模計(jì)算進(jìn)行分析?？臻g步長(zhǎng)設(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截面通過(guò)連接邊界引入入射波。時(shí)間步長(zhǎng)設(shè)置為5.295×10-13s，總共計(jì)算8192步。在計(jì)算過(guò)程中對(duì)z=85截面上場(chǎng)值進(jìn)行時(shí)間采樣，在全部時(shí)間步計(jì)算完成后，將采樣結(jié)果通過(guò)傅里葉變換轉(zhuǎn)換到頻域并通過(guò)Poyinting定理求得橫截面的功率函數(shù)，與入射波的對(duì)應(yīng)功率函數(shù)相比后得到頻率選擇表面的功率反射系數(shù)。

計(jì)算過(guò)程采用兩套GPU運(yùn)算平臺(tái)，GPU平臺(tái)1為Nvidia GT130M圖形處理器與512 MB顯存，GPU平臺(tái)2為T(mén)esla C2050圖形處理器與3 GB顯存，CPU對(duì)比平臺(tái)為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)算過(guò)程中抽取時(shí)間步t=1000時(shí)，y=30截面上Ey值進(jìn)行輸出比較，雙GPU平臺(tái)數(shù)值完全吻合，GPU平臺(tái)計(jì)算結(jié)果與CPU平臺(tái)符合很好如圖3所示。

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

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

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

4.3 加速性能分析

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

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

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

5.總結(jié)與展望

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

在GPU加速三維FDTD算法的可行性與高速性得到驗(yàn)證的同時(shí)，還有許多方面值得后繼深入研究，例如FDTD網(wǎng)格與線(xiàn)程的對(duì)應(yīng)是否有更高效的方式，GPU加速FDTD運(yùn)算的誤差分析與控制機(jī)理等等，這也是未來(lái)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)加速時(shí)域有限元的二維輻射計(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] 吳 霞, 周樂(lè)柱. 時(shí)域有限元法在計(jì)算電磁問(wèn)題上的發(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.