隱含變向時(shí)域有限差分方法的MPI實(shí)現(xiàn)

李太全，陳 威 (長(zhǎng)江大學(xué)物理科學(xué)與技術(shù)學(xué)院，湖北 荊州 434023)

隱含變向時(shí)域有限差分方法的MPI實(shí)現(xiàn)

李太全，陳 威 (長(zhǎng)江大學(xué)物理科學(xué)與技術(shù)學(xué)院，湖北 荊州 434023)

在MPI環(huán)境中，研究了適合于隱含變向時(shí)域有限差分算法(ADI-FDTD)的虛擬拓?fù)浜凸?jié)點(diǎn)間的數(shù)據(jù)通信。將并行對(duì)角占優(yōu)算法(PDD)應(yīng)用于ADI-FDTD，極大地減少了并行節(jié)點(diǎn)間的數(shù)據(jù)通信。進(jìn)一步分析數(shù)據(jù)通信的傳輸速率，提出了數(shù)據(jù)成批傳送方案，實(shí)現(xiàn)了ADI-FDTD的高效率并行計(jì)算。

MPI；隱含變向時(shí)域有限差分算法;三對(duì)角方程組;并行對(duì)角占優(yōu)算法

時(shí)域有限差分算法(FDTD)[1]在電磁場(chǎng)輻射和散射、微波電路以及電磁兼容等領(lǐng)域獲得了非常成功的應(yīng)用。然而，較長(zhǎng)的計(jì)算時(shí)間和較大的存儲(chǔ)空間是FDTD在PC系統(tǒng)上求解復(fù)雜電磁場(chǎng)問題的瓶頸。 解決這一問題的有效方法之一是采用并行FDTD方法[2-3]實(shí)現(xiàn)FDTD的分布運(yùn)算。隨著計(jì)算機(jī)網(wǎng)絡(luò)的迅速發(fā)展，基于消息傳遞接口(MPI)[4]的編程環(huán)境被廣泛采用，是目前最重要的并行編程工具，該系統(tǒng)具有移植性好、功能強(qiáng)大、效率高等諸多優(yōu)點(diǎn)，已經(jīng)成為國際上一種并行程序標(biāo)準(zhǔn)，被越來越多的硬件生產(chǎn)廠商和計(jì)算機(jī)用戶所接受。

顯式FDTD的時(shí)間步長(zhǎng)受Courant條件限制，在精細(xì)結(jié)構(gòu)的電磁分析中，由于空間步長(zhǎng)遠(yuǎn)小于電磁波波長(zhǎng)，其計(jì)算效率顯著降低。隱含變向時(shí)域有限差分算法(ADI-FDTD)[5-6]的提出，打破Courant條件限制，計(jì)算效率得到了明顯提高。然而，ADI-FDTD算法需要求解一組三對(duì)角方程，這組方程導(dǎo)致了并行ADI-FDTD方法復(fù)雜化。筆者將應(yīng)用三對(duì)角矩陣并行算法，實(shí)現(xiàn)ADI-FDTD的并行計(jì)算。

1 子區(qū)域劃分與虛擬拓?fù)?/h2>

將整個(gè)FDTD 計(jì)算區(qū)域劃分為若干個(gè)子區(qū)域，每個(gè)進(jìn)程(或線程)計(jì)算其中的一個(gè)子區(qū)域，各個(gè)進(jìn)程之間通過通信傳遞交界面上的電磁場(chǎng)量，以實(shí)現(xiàn)FDTD 中電磁場(chǎng)量的遞推。子區(qū)域的分解方法應(yīng)根據(jù)計(jì)算區(qū)域的尺度，在權(quán)衡運(yùn)算開銷和通信開銷的前提下選擇劃分方案。一般采用圖1所示的方法[2]，將計(jì)算區(qū)域沿著3個(gè)方向分解為子區(qū)域。MPI提供了迪卡兒拓?fù)浜蛨D拓?fù)?，根?jù)上述的子區(qū)域劃分，每個(gè)子區(qū)域?qū)?yīng)一個(gè)進(jìn)程，而每個(gè)進(jìn)程對(duì)應(yīng)拓?fù)渲械囊粋€(gè)節(jié)點(diǎn)，所以，應(yīng)選擇笛卡兒拓?fù)洹?/p>

2 ADI-FDTD的三對(duì)角方程組與PDD算法

隱含變向時(shí)域有限差分算法將一個(gè)時(shí)間步的電磁場(chǎng)量遞推分為2個(gè)亞時(shí)間步[4]，在每個(gè)亞時(shí)間步的遞推運(yùn)算中，6個(gè)電磁場(chǎng)分量有3個(gè)需要求解三對(duì)角方程，如果在第一個(gè)亞時(shí)間步選定Ex、Ey、Ez應(yīng)用方程組求解，則Hx、Hy、Hz可以直接計(jì)算。以Ex為例，對(duì)應(yīng)的三對(duì)角方程可表示為：

(1)

式中，α、β、γ、p、q是與空間步長(zhǎng)、時(shí)間步長(zhǎng)和介質(zhì)的電磁特性相關(guān)的系數(shù)。

第2個(gè)亞時(shí)間步計(jì)算式與式(1)類似。關(guān)于方程組(1)的求解，由于在z軸方向上Ex相互牽連，需要將圖1中沿z軸位于同一直線的子區(qū)域合并考慮，如圖2所示。

圖1 子區(qū)域分解 圖2 求解Ex相關(guān)聯(lián)的子區(qū)域

求解方程(1)可采用并行對(duì)角占優(yōu)算法(PDD)[7]。對(duì)確定的i、j，方程(1)可表示如下：

(2)

將矩陣A分解為：

V=[αmum,γm-1um-1,…,αm(Np-1)um(Np-1),γm(Np-1)um(Np-1)-1]

U=[um-1,um,…,um(Np-1)-1,um(Np-1)]

式中，A(p)為m×m階三對(duì)角矩陣，對(duì)應(yīng)圖2中的一個(gè)子區(qū)域；(p)表示子區(qū)域編號(hào)；um為Nz個(gè)元素的列向量，且除第n元素為1外，其他元素均為0。

式(2)的解可以寫成：

(3)

(4)

(5)

αk=-Δt/(μ(i+0.5,j,k-0.5)Δz2)

γk=-Δt(μ(i+0.5,j,k+0.5)Δz2)

(6)

1)確定當(dāng)前節(jié)點(diǎn)p對(duì)應(yīng)子區(qū)域中網(wǎng)格i、j對(duì)應(yīng)的A(p)和d(p)；

移動(dòng)網(wǎng)絡(luò)的網(wǎng)絡(luò)節(jié)點(diǎn)是人，而不是網(wǎng)頁，也就是說人人都是數(shù)據(jù)的產(chǎn)生著和攜帶者，無論是微信、照片、微博等任何網(wǎng)絡(luò)相聯(lián)系的，都已成為數(shù)據(jù)的產(chǎn)生者。所涉及到的資料總量是目前任何一個(gè)軟件都不可能接受、管理、處理、總結(jié)的。有國外機(jī)構(gòu)預(yù)測(cè)，全球的數(shù)據(jù)量將在2020年擴(kuò)大到目前的五十倍，這樣就意味著大數(shù)據(jù)時(shí)代的容量將會(huì)更多、更大、更難處理。

2)求解方程組：

(7)

4)求解方程組：

(8)

(9)

tC=2NxNytα+12NxNytβ

(10)

式中，tα為通信響應(yīng)時(shí)間；tβ為一字節(jié)數(shù)據(jù)的傳輸時(shí)間(浮點(diǎn)數(shù)為4字節(jié)單精度實(shí)數(shù))。對(duì)試驗(yàn)網(wǎng)絡(luò)測(cè)試，tα=26.3μs，tβ=0.0875μs。由于tα?tβ，將數(shù)據(jù)成批傳送會(huì)大大提高通信效率。為此，采用如下步驟求解方程(1)：

3)對(duì)i、j循環(huán)，求解方程(8)，并存儲(chǔ)y0，y1；

這樣， 原來的2×Ny×Ny次數(shù)據(jù)通信減少為2次，此時(shí)通信時(shí)間開銷為tC=2tα+12NxNytβ。

3 效率分析與數(shù)值試驗(yàn)

設(shè)子區(qū)域內(nèi)FDTD網(wǎng)格數(shù)為Nx×Ny×Nz，在一次循環(huán)迭代中，運(yùn)算開銷約為6(te+th)NxNyNz，其中te、th分別為一個(gè)電場(chǎng)、磁場(chǎng)計(jì)算所需時(shí)間，通信開銷約為18tα+28tβ(NxNy+NyNz+NzNy)，算法的效率可以表示為:

圖3 在Nx=Ny=Nz=N時(shí)，效率與網(wǎng)格數(shù)N的關(guān)系

可見，當(dāng)子區(qū)域的網(wǎng)格數(shù)Nx=Ny=Nz時(shí)，算法效率最高。

由個(gè)人計(jì)算機(jī)組成10個(gè)節(jié)點(diǎn)的集群系統(tǒng)，每個(gè)節(jié)點(diǎn)的CPU主頻2GHz，內(nèi)存1024MHz，100/1000M網(wǎng)卡，100M交換機(jī)。在該系統(tǒng)中測(cè)得，te=0.583μs，th=0.104μs，tα=26.3μs，tβ=0.0875μs。在該系統(tǒng)中，取Nx=Ny=Nz=N時(shí)，算法的效率如圖3所示，在N>25時(shí)，算法已具有較高的效率。進(jìn)一步檢驗(yàn)算法的效率，筆者計(jì)算了如圖4所示的環(huán)形天線。在相對(duì)介電常數(shù)為εr=4.5、邊長(zhǎng)為1000m、厚為h=50m的正方形基板上，有一個(gè)線寬b=12m、邊長(zhǎng)a=500m的正方形環(huán)帶狀線天線，激勵(lì)點(diǎn)位于底邊的中點(diǎn)E。激勵(lì)源輸出阻抗Z=100Ω，輸出電壓v(t)=exp(-4π(t-t0)2/τ2)，其中τ=100ps。將整個(gè)空間劃分為Nx×Ny×Nz=300×40×300網(wǎng)格，且Δx=Δy=4m，Δz=10m。取時(shí)間步長(zhǎng)Δt=2×Δx/c，將空間在XY面等分為2、4、9個(gè)子區(qū)域。對(duì)該系統(tǒng)進(jìn)行200時(shí)間步迭代，運(yùn)算時(shí)間如表1所示。

表1 不同區(qū)域劃分的運(yùn)算時(shí)間比較

圖5是在上述4種情況下，分別得到天線輸入口的S11，4種情況下的良好一致性說明了上述算法的正確性。

4 結(jié) 語

采用MPI編程環(huán)境，實(shí)現(xiàn)了ADI-FDTD的并行計(jì)算。將PDD算法應(yīng)用于并行ADI-FDTD，極大地減少了節(jié)點(diǎn)間的數(shù)據(jù)傳送量，減小了通信開銷。并采用交換信息的成批傳送，減少通信次數(shù)，雖然需要額外內(nèi)存開銷暫存臨時(shí)數(shù)據(jù)，但可以顯著提高了計(jì)算效率。

圖4 環(huán)形天線結(jié)構(gòu) 圖5 4種情況下天線的S11

[1]葛德彪, 閆玉波.電磁波時(shí)域有限差分方法[M].西安：西安電子科技大學(xué)出版社, 2002.

[2] 張玉, 李斌, 梁昌洪. PC集群系統(tǒng)中MPI并行FDTD 性能研究[J].電子學(xué)報(bào), 2005, 33(9):1694-1697.

[3] 楊利霞, 葛德彪, 鄭奎松，等.電各向異性介質(zhì)FDTD 并行算法的研究[J].電波科學(xué)學(xué)報(bào), 2006, 21(1):43-48.

[4] 都志輝.高性能計(jì)算并行編程技術(shù)[M].北京：清華大學(xué)出版社, 2001.

[5] Namiki T.3-D ADI-FDTD Method Unconditionally Stable Time-Domain Algorithm for Solving Full Vector Maxwells Equations[J].IEEE Transactions on Microwave Theory and Techniques, 2000, 48(10): 1743-1747.

[6] Namiki T.A New FDTD Algorithm Based on Alternating-Direction Implicit Method[J].IEEE Trans- actions on Microwave Theory and Techniques, 1999, 47(10): 2003-2007.

[7] Sun X H, Zhang H, Ni L.Efficient Tri- diagonal Solvers on Multicomputers[J].IEEE Trans-actions on Computers, 1992, 41(3): 286-296.

2012-10-15

國家自然科學(xué)基金項(xiàng)目(41140034)。

李太全(1961-)，男，博士，副教授，現(xiàn)主要從事電子技術(shù)方面的教學(xué)與研究工作。

TN911.2

A

1673-1409(2013)01-0006-04

[編輯] 洪云飛