徐奇澎，郭裕順

(杭州電子科技大學(xué)電子信息學(xué)院，杭州 310018)

FDTD 法即時(shí)域有限差分方法(Finite Difference in Time Domain)，是把Maxwell 方程式在時(shí)間和空間域進(jìn)行差分，通過(guò)電場(chǎng)和磁場(chǎng)的交替計(jì)算，來(lái)模擬電磁波傳播的一種數(shù)值計(jì)算方法［1］。自上世紀(jì)六十年代出現(xiàn)這種算法以來(lái)，人們已對(duì)它作了大量研究，是目前電磁波數(shù)值仿真的主要方法，被廣泛地應(yīng)用于科學(xué)與工程中各種電磁場(chǎng)與電磁波問(wèn)題。

FDTD 算法應(yīng)用中存在的一個(gè)主要問(wèn)題是:當(dāng)計(jì)算區(qū)域較大，即空間差分產(chǎn)生較多網(wǎng)格時(shí)，需要耗費(fèi)大量的計(jì)算資源與時(shí)間［2］。為此降低計(jì)算代價(jià)，減少計(jì)算時(shí)間，一直是重要的研究課題。以往采取的措施不外乎兩個(gè)方面:一是改進(jìn)算法，設(shè)法提高計(jì)算效率;二是在通常的計(jì)算平臺(tái)上采用并行計(jì)算。

FPGA 的發(fā)展為解決這類(lèi)大規(guī)模計(jì)算問(wèn)題提供了一種新途徑。FDTD 計(jì)算的特點(diǎn)是在固定的一個(gè)空間網(wǎng)格中進(jìn)行重復(fù)的算術(shù)運(yùn)算，數(shù)據(jù)流規(guī)則，控制流較為簡(jiǎn)單，因此可利用FPGA 可編程的特點(diǎn)，設(shè)計(jì)專(zhuān)門(mén)的算術(shù)運(yùn)算電路，實(shí)現(xiàn)FDTD 的計(jì)算。一塊FPGA 芯片可同時(shí)實(shí)現(xiàn)多個(gè)乘法器、加法器，輕易實(shí)現(xiàn)并行運(yùn)算。這樣針對(duì)FDTD 算法專(zhuān)門(mén)設(shè)計(jì)的運(yùn)算電路，可大大提高計(jì)算速度。最早用FPGA 實(shí)現(xiàn)FDTD 算法的是Schneider，他用定點(diǎn)整數(shù)格式實(shí)現(xiàn)了FDTD 算法的一維運(yùn)算，但網(wǎng)格單元僅10個(gè)［3］。Gandhi、Durbano 分別用浮點(diǎn)數(shù)實(shí)現(xiàn)了FDTD 算法的二維、三維運(yùn)算，但受限于當(dāng)時(shí)FPGA 的性能，這些設(shè)計(jì)工作頻率較低，性能不高，未能體現(xiàn)硬件計(jì)算的優(yōu)勢(shì)［4］。后來(lái)Wu Shuguang 在Gandhi 的研究基礎(chǔ)上把數(shù)據(jù)存儲(chǔ)位置由主機(jī)系統(tǒng)內(nèi)存改為FPGA 片上SRAM，大大減少了數(shù)據(jù)傳輸耗費(fèi)的時(shí)間［5］。最近Hideki Kawaguchi 等采用分配式并行存儲(chǔ)使數(shù)據(jù)讀寫(xiě)速度進(jìn)一步加快，當(dāng)FPGA 工作頻率在50 MHz時(shí)，計(jì)算速度比軟件工作在2.8 GHz 的個(gè)人電腦上快了50倍［6－7］。因此可以預(yù)計(jì)隨著FPGA 的進(jìn)一步發(fā)展，用FPGA 實(shí)現(xiàn)FDTD 算法的優(yōu)勢(shì)將更加明顯。

1 二維FDTD 算法

本文考慮二維FDTD 算法的實(shí)現(xiàn)。二維情形下所有電磁量與Z 坐標(biāo)無(wú)關(guān)［8］，即?/?z=0。這時(shí)電磁場(chǎng)的直角分量可以分為獨(dú)立的兩組，即以Hx、Hy、Ez為一組的TM 波，及以Ex、Ey、Hz為一組的TE波［9］。以TM 波為例，二維TM 波的Yee 元胞如圖1所示，電場(chǎng)值根據(jù)周?chē)?個(gè)磁場(chǎng)值完成更新，磁場(chǎng)值根據(jù)周?chē)?個(gè)電場(chǎng)值完成更新，更新關(guān)系滿(mǎn)足右手螺旋法則。

圖1 二維TM 波的Yee 元胞

為了便于硬件實(shí)現(xiàn)，對(duì)上述關(guān)系中的1/2 按照一定規(guī)則做出處理［10］。處理后TM 波的計(jì)算公式如下:

下面介紹用FPGA 來(lái)完成上述計(jì)算的方法。

2 算法的硬件實(shí)現(xiàn)

用硬件實(shí)現(xiàn)式(1)～式(3)的計(jì)算首先要確定數(shù)字的表示方法。本文采用符合IEEE-754 32 位標(biāo)準(zhǔn)的單精度浮點(diǎn)數(shù)格式，即將一個(gè)數(shù)表示為V=(－1)S×1.M×2E-BIAS，其中最高位S 是符號(hào)位，隨后的8 位E是指數(shù)位，后面的23 位M 是尾數(shù)位，BIAS=28－1=127。

圖2 電路結(jié)構(gòu)

圖3 計(jì)算電路

實(shí)現(xiàn)一次FDTD 計(jì)算的工作流程如圖4所示，分為以下幾步:

(1)初始化:確定要計(jì)算的網(wǎng)格大小，所要計(jì)算的時(shí)間步，計(jì)算式系數(shù)。

(2)加入激勵(lì)后，整個(gè)計(jì)算電路開(kāi)始工作。本文中因?yàn)榧?lì)信號(hào)是添加在RAM Ez中，所以先計(jì)算Ez值。根據(jù)計(jì)數(shù)器計(jì)算出相應(yīng)的地址值，從RAM 中取得相鄰兩個(gè)單元的磁場(chǎng)值，Ez的計(jì)算電路完成一個(gè)點(diǎn)的Ez值計(jì)算，E值計(jì)數(shù)器會(huì)加1，然后檢查其是否數(shù)到最大值。如果計(jì)數(shù)器沒(méi)有達(dá)到最大值則再次計(jì)算Ez值，進(jìn)行RAM Ez中下一個(gè)點(diǎn)的更新。如果是計(jì)數(shù)器達(dá)到最大值則表明RAM Ez中所有的點(diǎn)都完成一個(gè)時(shí)間步的更新，轉(zhuǎn)換更新算式開(kāi)始進(jìn)行磁場(chǎng)值的更新。

(3)Hx和Hy的更新同時(shí)進(jìn)行，當(dāng)一個(gè)點(diǎn)的Hx和Hy值更新完成之后，H值計(jì)數(shù)器會(huì)加1，然后檢查其是否數(shù)到最大值。如果計(jì)數(shù)器沒(méi)有達(dá)到最大值則再次計(jì)算Hx和Hy值，進(jìn)行RAM Hx和RAM Hy中下一個(gè)點(diǎn)的更新。如果計(jì)數(shù)器達(dá)到最大值則表明RAM Hx和RAM Hy中所有的點(diǎn)都完成一個(gè)時(shí)間步的更新，同時(shí)也意味著網(wǎng)格中所有的點(diǎn)都完成了一個(gè)時(shí)間步的電場(chǎng)值和磁場(chǎng)值的更新。

(4)完成一個(gè)時(shí)間步更新后，時(shí)間步計(jì)數(shù)器T加1，然后檢查T(mén) 是否小于nstep，如果是則表示還未完成設(shè)定的時(shí)間步，回到算法第3 步，開(kāi)始網(wǎng)格內(nèi)所有節(jié)點(diǎn)下一個(gè)時(shí)間步的更新。如果不是則表示計(jì)算完成設(shè)定的時(shí)間步，整個(gè)算法完成。

圖4 算法流程

在設(shè)計(jì)過(guò)程中我們可以根據(jù)并行原理［11－12］和FPGA 架構(gòu)特點(diǎn)，設(shè)計(jì)多個(gè)計(jì)算電路，讓這些計(jì)算電路同時(shí)運(yùn)算，可以大大減少運(yùn)算所需要的時(shí)間，提高算法的效率。

3 實(shí)驗(yàn)結(jié)果

根據(jù)上述方案，用Verilog HDL 編寫(xiě)了代碼，以ALTERA 公司的EP3C55F484C6為目標(biāo)器件進(jìn)行了編譯，用ModelSim 進(jìn)行了仿真。設(shè)置計(jì)算網(wǎng)格大小為60 ×60，在網(wǎng)格中間加一個(gè)高斯脈沖g(t)=exp(－0.5(t－t0)2/δ2)充當(dāng)激勵(lì)源，觀察3 600個(gè)網(wǎng)格點(diǎn)經(jīng)過(guò)60個(gè)步長(zhǎng)時(shí)間的運(yùn)算結(jié)果和運(yùn)算所需要的時(shí)間。圖5 就是仿真所得部分結(jié)果，端口data_out 輸出的是EZ RAM 中的值，圖中可以看到輸出的值是符合IEEE 754 標(biāo)準(zhǔn)的32 位單精度浮點(diǎn)數(shù)，由于采用了流水線(xiàn)設(shè)計(jì)，每個(gè)時(shí)鐘輸出一個(gè)運(yùn)算結(jié)果。圖6 是3 600個(gè)點(diǎn)在經(jīng)過(guò)60 步后仿真結(jié)果。

為了比較FPGA 計(jì)算的結(jié)果，我們對(duì)同樣計(jì)算用Visual C++6.0 編寫(xiě)了軟件程序，結(jié)果如圖7，兩者誤差不超過(guò)0.11%。FPGA 設(shè)計(jì)的最高工作頻率可達(dá)118.39 MHz，完成整個(gè)運(yùn)算過(guò)程需670 000個(gè)時(shí)鐘周期。假設(shè)FPGA 工作頻率為10 MHz，則完成全部運(yùn)算需要0.067 s，而軟件程序在主頻2.2 GHz的個(gè)人電腦上運(yùn)行時(shí)間是1.48 s，因此FPGA 計(jì)算要快幾乎22倍。當(dāng)采用2倍并行計(jì)算時(shí)，完成全部運(yùn)算需368 310個(gè)時(shí)鐘周期，加速比可達(dá)40倍。因此FPGA 計(jì)算對(duì)提高FDTD 算法速度的效果是十分明顯的。

圖5 ModelSim 仿真結(jié)果

圖6 FPGA 運(yùn)算結(jié)果

圖7 C 語(yǔ)言運(yùn)算結(jié)果

4 結(jié)論

本文對(duì)FDTD 算法的FPGA 實(shí)現(xiàn)進(jìn)行了初步研究。以?xún)删S情形為例，編寫(xiě)了Verilog HDL 實(shí)現(xiàn)代碼，在Altera FPGA 上進(jìn)行了編譯仿真，獲得了初步試驗(yàn)結(jié)果。在設(shè)計(jì)中采用了流水線(xiàn)技術(shù)、并行運(yùn)算等手段，同時(shí)采用了FPGA 片內(nèi)RAM 作為存儲(chǔ)單元，取得了較好的加速效果。隨著制造工藝的發(fā)展，F(xiàn)PGA 片內(nèi)資源會(huì)越來(lái)越豐富，基于FPGA 的計(jì)算可能成為提高FDTD 算法性能的一條有效途徑。

［1］Yee K S.Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equation in Isotropic Media［J］.IEEE Trans Antennas Propagate，1966，(14):302－307.

［2］丁偉.時(shí)域有限差分法關(guān)鍵技術(shù)及其應(yīng)用研究［D］.西安電子科技大學(xué)，2007.

［3］Schneider R N，Turner L E，Okoniewski M M.Application of FPGA Technology to Accelerate the Finite-Difference Time-Domain(FDTD)Method［C］//The 2002 ACM International Symposium on Field-Programmable Gate Arrays(FPGA'02)，2002:24－26.

［4］Durbano J P，Ortiz F E，Humphrey J R.Hardware Implementation of a Three-Dimensional Finite-Difference Time-Domain Algorithm［C］//IEEE Antennas and Wireless Propagation Letters，2003:2，54－57.

［5］Wu S.An FPGA Implementation of FDTD Codes for Reconfigurable High Performance Computing［D］.M.S.Thesis，Dept.of Electrical and Computer Engineering and Computer Science，University of Cincinnati，2005.

［6］Kawaguchi H，F(xiàn)ujita Y，F(xiàn)ujishima Y.Improved Architecture of FDTD/FIT Dedicated Computer for Higher Performance Computation［J］.IEEE Transactions on Magnetics，2008，44(6):1226－1229.

［7］Fujita Y，Kawaguchi H.Full Custom PCB Implementation of FDTD/FIT Dedicated Computer［J］.IEEE Transactions on Magnetics，2009，45(3):1100－1103.

［8］張煊，張明德.二維光子晶體波導(dǎo)及耦合器的特性模擬［J］.電子器件，2005，28(6):63－67.

［9］葛德彪，閏玉波.電磁波時(shí)域有限差分方法［M］.西安電子科技大學(xué)出版社，2002.

［10］Dennis M Sullivan.Electromagnetic Simulation Using the FDTD Method［M］.IEEE Microwave Theory and Techniques Society，2000.

［11］雷繼兆.PC和服務(wù)器集群下的并行FDTD 算法及其應(yīng)用研究［D］.西安電子科技大學(xué)，2009.

［12］劉瑜.FDTD 算法的網(wǎng)絡(luò)并行研究及其電磁應(yīng)用［D］.電子科技大學(xué)，2008.