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

（1.安徽大學(xué)計算智能與信號處理教育部重點實驗室，安徽合肥 230039；2.安徽大學(xué)電氣工程與自動化學(xué)院，安徽合肥 230039）

0 引 言

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

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

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

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

1 SFDTD的基本理論

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

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

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

通過變分法，可寫成：

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

1.1.1 時間方向上的離散

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

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

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

當取m＝5且p＝4時，可得5級4階辛積分方法。

1.1.2 空間方向上的離散

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

網(wǎng)格節(jié)點與相應(yīng)的整數(shù)標號對應(yīng)，即

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

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

1.2 二維SFDTD算法

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

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

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

2 并行化算法

2.1 CUDA模型

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

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

圖1 CUDA程序的編程模型

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

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

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

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

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

3.1 二維SFDTD算法

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

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

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

圖3 x＝20時場值Ez分布圖

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

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

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

4 結(jié)束語

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

［1］ 黃志祥，吳先良.辛算法的穩(wěn)定性及數(shù)值色散性分析［J］.電子學(xué)報，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］ 呂英華.計算電磁學(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平臺的時域有限差分算法研究［J］.合肥工業(yè)大學(xué)學(xué)報：自然科學(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＿home＿new＿cn.html.

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