基于GPU的高階辛FDTD算法的并行仿真研究

馬巍巍， 孫 冬， 吳先良， 孫兵兵

（1.安徽大學(xué)計(jì)算智能與信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230039；2.安徽大學(xué)電氣工程與自動(dòng)化學(xué)院，安徽合肥 230039）

0 引 言

時(shí)域有限差分算法（Finite－difference timedomain method，簡(jiǎn)稱FDTD）是計(jì)算電磁學(xué)的一種常用方法。但在電大尺寸或長(zhǎng)時(shí)間仿真中，由于數(shù)值色散性和穩(wěn)定性等的影響，數(shù)值誤差會(huì)隨時(shí)間累計(jì)，造成計(jì)算結(jié)果失真。高階辛?xí)r域有限差分算法（the high－order symplectic finite－difference time－domain method，簡(jiǎn)稱SFDTD）可以解決這一問題，其優(yōu)越性也已被證實(shí)［1－3］。

但SFDTD進(jìn)行電磁仿真相當(dāng)耗時(shí)，且占有更大的計(jì)算存儲(chǔ)空間［3］。本文提出采用SFDTD與并行算法相結(jié)合的方法來進(jìn)行電磁仿真，以縮短計(jì)算時(shí)間，取得了較好的結(jié)果。本文在時(shí)間方向上采用高階辛差分，在長(zhǎng)期仿真中保持麥克斯韋方程的辛結(jié)構(gòu)；空間方向上采用4階交錯(cuò)差分以減小數(shù)值色散。

并行化是數(shù)值計(jì)算的必然趨勢(shì)［4］。圖形處理器（graphics processing unit，簡(jiǎn)稱GPU）并行技術(shù)所需硬件成本較低，且具有較高的計(jì)算能力，在計(jì)算統(tǒng)一設(shè)備架構(gòu)（compute unified device architecture，簡(jiǎn)稱CUDA）出現(xiàn)后，解決了GPU計(jì)算模型中的很多限制，已成為并行計(jì)算領(lǐng)域里的熱點(diǎn)研究方向，在FDTD算法上的應(yīng)用也已有了很多的發(fā)展［5－6］。但到目前為止，對(duì)高階辛算法進(jìn)行并行的研究還比較少。

本文在分析SFDTD算法和CUDA架構(gòu)的基礎(chǔ)上，研究了基于GPU加速的二維電磁波的SFDTD并行算法。在最新的Fermi架構(gòu)下，實(shí)現(xiàn)了該算法，與單CPU相比，計(jì)算速度可提升數(shù)十倍，驗(yàn)證了所提算法的有效性。

1 SFDTD的基本理論

1.1 麥克斯韋方程的哈密頓表述

已知在無源、均勻、無耗媒質(zhì)中的麥克斯韋旋度方程為：

根據(jù)哈密頓系統(tǒng)理論［7］，其對(duì)應(yīng)的哈密頓函數(shù)的微分形式為：

通過變分法，可寫成：

ε和μ分別為媒質(zhì)的介電常數(shù)和磁導(dǎo)率。

1.1.1 時(shí)間方向上的離散

由（3）式從t＝0到Δt的精確解為：

為了近似時(shí)間演化矩陣exp（Δt（U＋V）），采用m級(jí)p階辛積分方法，即

其中，cl和dl為辛算子，具體系數(shù)見文獻(xiàn)［2］。

當(dāng)取m＝5且p＝4時(shí)，可得5級(jí)4階辛積分方法。

1.1.2 空間方向上的離散

在YEE網(wǎng)格下建立矩陣差分網(wǎng)格，任一場(chǎng)量函數(shù)F（x，y，z，t）在第n個(gè)時(shí)間步第l級(jí)步進(jìn)的離散表達(dá)式為：

網(wǎng)格節(jié)點(diǎn)與相應(yīng)的整數(shù)標(biāo)號(hào)對(duì)應(yīng)，即

對(duì)于普通的FDTD方法，每個(gè)時(shí)間步只需要1級(jí)時(shí)間步進(jìn)。本文用4階中心（交錯(cuò)）差分離散空間，1階導(dǎo)數(shù)為：

經(jīng)空間上的4階中心差分、時(shí)間上的4階離散后，建立了麥克斯韋方程的離散辛框架，即SFDTD（4，4）算法。

1.2 二維SFDTD算法

本文以二維TM波為例，采用5級(jí)辛積分方法和4階中心差分。辛算子及參數(shù)的具體取值見文獻(xiàn)［2］。為能使電場(chǎng)和磁場(chǎng)具有相同數(shù)量級(jí)，將方程式進(jìn)行歸一化處理，因此在無源的自由空間中歸一化TM波方程為：

由（9）式可以看出，本文所研究的SFDTD算法為4階中心差分，具有天然的并行性。其空間電磁場(chǎng)排布為：每1個(gè)電場(chǎng)分量被8個(gè)磁場(chǎng)分量環(huán)繞，每1個(gè)磁場(chǎng)分量被8個(gè)電場(chǎng)分量環(huán)繞。這種空間取樣方式同樣采用了YEE網(wǎng)格，只不過采用了雙環(huán)路的形式。

SFDTD電（磁）場(chǎng)的推進(jìn)僅需要周圍的磁（電）場(chǎng)的信息及上一步的本點(diǎn)場(chǎng)值，不需要考慮整個(gè)計(jì)算區(qū)域的場(chǎng)值分布，因此具有很好的空間并行優(yōu)勢(shì)。

2 并行化算法

2.1 CUDA模型

計(jì)算統(tǒng)一設(shè)備架構(gòu)（CUDA）是通用GPU計(jì)算的一種技術(shù)規(guī)范，由NVIDIA公司推出。它采用了類C語言進(jìn)行開發(fā)，寫出在顯示芯片上執(zhí)行的程序。在CUDA架構(gòu)下，一個(gè)程序分為在CPU上執(zhí)行的Host端和在顯示芯片上執(zhí)行的Device端。在Device中執(zhí)行的代碼稱為內(nèi)核（Kernal）。通常Host端程序會(huì)將數(shù)據(jù)準(zhǔn)備好后，復(fù)制到顯存中，再由GPU執(zhí)行Device端的程序，完成后，再由Host端程序?qū)⒔Y(jié)果從GPU中取回，編程模型如圖1所示。

CPU主要負(fù)責(zé)執(zhí)行串行部分的代碼，從圖1可以看出，在CPU中運(yùn)行的代碼只有一個(gè)線程，這部分通常是一些控制類的操作。接下來是并行執(zhí)行的部分，這部分都是在GPU中運(yùn)行的。顯示芯片執(zhí)行的最小單位是thread（線程），數(shù)個(gè)線程可組成一個(gè)block（塊）。數(shù)個(gè)塊可組成一個(gè)grid（網(wǎng)格）。關(guān)于CUDA編程模型的詳細(xì)介紹可參閱文獻(xiàn)［8］。

圖1 CUDA程序的編程模型

根據(jù)這一架構(gòu)的特點(diǎn)，在實(shí)現(xiàn)SFDTD并行算法時(shí)，由于SFDTD在空間上的并行性，計(jì)算區(qū)域被劃分成為多個(gè)塊，然后每個(gè)塊又劃分為多個(gè)線程，用每個(gè)線程來計(jì)算一個(gè)YEE元胞。具體實(shí)現(xiàn)的主要編程思路如圖2所示。

圖2 程序?qū)崿F(xiàn)流程圖

2.2 支持CUDA的圖形處理器及Fermi架構(gòu)

2007年發(fā)布的CUDA是世界上第1個(gè)針對(duì)圖形顯示器（GPU）的C語言開發(fā)環(huán)境。NVIDIA公司2007年后推出的顯存超過256 MB的GPU都可以運(yùn)行和開發(fā)基于CUDA的并行算法。而在短短的幾年間，支持CUDA的硬件有了多次重要的變化［9］?；贑UDA編寫的代碼具有很好的通用性，在任何一款支持CUDA的單機(jī)上都可以運(yùn)行。最新的Fermi架構(gòu)提升了通用計(jì)算的重要性，與之前架構(gòu)相比有很多的優(yōu)勢(shì)。例如，在G80／GT200系列中，F(xiàn)ermi構(gòu)架中每個(gè)線程塊支持高達(dá)1 536個(gè)線程；Fermi允許多個(gè)內(nèi)核同時(shí)占用一個(gè)SM的計(jì)算資源，從而有效提高了設(shè)置的利用率；在Fermi中，768 k B的二級(jí)緩存，可加速全局存儲(chǔ)器和紋理存儲(chǔ)器的讀取。

3 二維SFDTD算法的CUDA仿真實(shí)現(xiàn)

3.1 二維SFDTD算法

為驗(yàn)證本文算法的正確性和有效性，使用SFDTD算法計(jì)算帶有PML邊界的TM平面波點(diǎn)源輻射問題，分別使用CPU和GPU對(duì)其進(jìn)行仿真。本文的程序開發(fā)環(huán)境為VS 2008，軟件運(yùn)行環(huán)境為32位Windows 7，內(nèi)存1 GB。處理器型號(hào)：GPU為NVIDIA GeForce GT 440；顯存512 MB；96個(gè)流處理器。CPU英特爾Pentium（奔騰）4 3.00 GHz，內(nèi)存1 G。所選擇的激勵(lì)源為位于區(qū)域中心的高斯脈沖源。CUDA在執(zhí)行程序時(shí)，是以warp為單位。目前的CUDA裝置，一個(gè)warp中有32個(gè)線程，分成2組的halfwarp，一次至少執(zhí)行16個(gè)線程才能有效隱藏各種運(yùn)算的延遲。因此，取Block－size＝16，即每個(gè)塊中有16×16個(gè)線程，再建立每個(gè)網(wǎng)格中的塊，這樣，在程序運(yùn)行時(shí)每個(gè)線程塊中都會(huì)有256個(gè)線程在同時(shí)執(zhí)行內(nèi)核。由于時(shí)間上是高階辛積分，在具體編程時(shí)與FDTD方法不同，是雙重循環(huán)，故其系數(shù)需小心處理正確放入顯存后，才能得出正確的結(jié)果。

3.2 數(shù)值結(jié)果及討論

考慮設(shè)置為計(jì)算域中心的二維真空中隨時(shí)間為高斯源變化的點(diǎn)源輻射問題。計(jì)算結(jié)果與CPU的計(jì)算結(jié)果進(jìn)行比較，如圖3所示。圖3中沿x，y方向的離散點(diǎn)格數(shù)分別為GRIDx＝GRIDy＝64，取CFL＝0.5。邊界條件為分裂場(chǎng)PML，時(shí)間步為100步。

圖3 x＝20時(shí)場(chǎng)值Ez分布圖

由圖3可見，采用GPU和CPU的圖像是一致的，這證明了GPU編碼的正確性。對(duì)網(wǎng)格數(shù)取不同值得到的運(yùn)行結(jié)果如圖4所示。

圖4 CPU、GPU的計(jì)算時(shí)間、加速比與網(wǎng)格關(guān)系

從圖4可以看出，在場(chǎng)點(diǎn)較少時(shí)，GPU的優(yōu)勢(shì)并不明顯。這是由于此時(shí)可分配線程數(shù)太少，多流處理器并未充分利用，且較少的線程也無法隱藏?cái)?shù)據(jù)在內(nèi)存與顯存之間的傳遞時(shí)間。但隨著網(wǎng)格數(shù)的增加，GPU并行計(jì)算的優(yōu)勢(shì)就得到了明顯體現(xiàn)。如在處理640×640個(gè)網(wǎng)格時(shí)，速度提升了10倍。在仿真實(shí)際的物理問題時(shí)，YEE元胞的數(shù)量通常會(huì)在數(shù)萬，甚至百萬以上，正適合于將這些數(shù)據(jù)分配給GPU的數(shù)百個(gè)標(biāo)量處理器進(jìn)行并行計(jì)算。

4 結(jié)束語

本文基于圖形處理器，提出一種高階辛?xí)r域有限差分方法的并行算法，并將其應(yīng)用到二維散射場(chǎng)問題的分析中。結(jié)果表明，該算法與傳統(tǒng)的CPU中實(shí)現(xiàn)該算法相比，在精度相當(dāng)?shù)那闆r下可節(jié)省大量時(shí)間，有助于電磁場(chǎng)問題的應(yīng)用，具有極為重要的實(shí)用價(jià)值。

［1］ 黃志祥，吳先良.辛算法的穩(wěn)定性及數(shù)值色散性分析［J］.電子學(xué)報(bào)，2006，34（3）：535－538.

［2］ Hirono T，Seki S，Lui W，et al.A three－dimensional fourthorder finite－difference time－domain scheme using a symplectic integrator propagator［J］.IEEE Transactions on Microwave Theory and Techniques，2001，49（9）：1640－1648.

［3］ 吳 瓊，黃志祥，吳先良.基于高階辛算法求解Maxwell方程［J］.系統(tǒng)工程與電子技術(shù)，2006，28（2）：342－344.

［4］ 呂英華.計(jì)算電磁學(xué)的數(shù)值方法［M］.北京：清華大學(xué)出版社，2006：6－8.

［5］ Valcarce A，De La Roche G，Jie Z.A GPU approach to FDTD for radio coverage prediction［C］／／The 11th IEEE Singapore International Conference on Communication Systems，Guangzhou，2008：1585－1590.

［6］ 沈 琛，王 璐，胡玉娟，等.基于CUDA平臺(tái)的時(shí)域有限差分算法研究［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，35（5）：644－647.

［7］ 沙 威.高階辛算法的理論和應(yīng)用研究［D］.合肥：安徽大學(xué)，2008.

［8］ NVIDIA Corporation Technical Staff.NVIDIA CUDA programming guide 2.0［EB／OL］.［2011－12－05］.http：／／www.nvidia.cn／object／cuda＿h(yuǎn)ome＿new＿cn.html.

［9］ 濮元愷.改變翻天覆地史上最全Fermi架構(gòu)解讀［EB／OL］.［2010－03－26］.http：／／www.sina.com.cn／