基于協(xié)同計(jì)算的FDTD區(qū)域分解并行算法研究

馮 圓,唐曉斌,龔曉燕

1.中國(guó)電子科學(xué)研究院,北京 100041

2.空軍預(yù)警學(xué)院,湖北 武漢 430019

3.火箭軍指揮學(xué)院,湖北 武漢 430013

美國(guó)、西班牙、土耳其、比利時(shí)等國(guó)外學(xué)者對(duì)并行電磁算法的研究取得了卓著成果。例如多層快速多極子方法（MLFMA）、時(shí)域有限差分法（FDTD）、有限元法（FEM）等多種算法都有對(duì)應(yīng)的并行算法。目前，HFSS、CST、FEKO等常見電磁軟件都開發(fā)了并行計(jì)算模塊。各國(guó)已經(jīng)將并行電磁計(jì)算技術(shù)廣泛用于隱身/反隱身等電子裝備的分析、設(shè)計(jì)、研制、評(píng)估等各個(gè)方面，推動(dòng)了裝備性能的快速升級(jí)。但隨著科學(xué)技術(shù)的不斷發(fā)展和電子武器裝備應(yīng)用領(lǐng)域的不斷擴(kuò)大，越來(lái)越多的更復(fù)雜的電磁問(wèn)題被提出，例如THz頻段的隱身/反隱身設(shè)計(jì)?？梢?，電子信息系統(tǒng)仍然存在著大量挑戰(zhàn)性課題迫切需要解決，傳統(tǒng)的并行計(jì)算技術(shù)逐漸不能滿足日益發(fā)展的電磁計(jì)算需求。為此，研究性能更為強(qiáng)大的電磁協(xié)同計(jì)算是十分必要的。作為協(xié)同計(jì)算的最新發(fā)展，云計(jì)算應(yīng)用于高性能科學(xué)計(jì)算領(lǐng)域具有廣闊前景。

國(guó)內(nèi)的協(xié)同計(jì)算技術(shù)距離國(guó)外有一定的差距。與國(guó)外研究情況類似，目前國(guó)內(nèi)的大規(guī)模電磁計(jì)算主要采用消息傳遞接口模型（MPI），即傳統(tǒng)的并行計(jì)算。而電磁協(xié)同計(jì)算特別是電磁云計(jì)算的研究工作尚屬空白。

直到最近，云計(jì)算作為一種新型的協(xié)同計(jì)算模式，在我國(guó)的發(fā)展方興未艾。中國(guó)科學(xué)技術(shù)大學(xué)率先對(duì)云計(jì)算MapReduce模型在科學(xué)計(jì)算中的應(yīng)用進(jìn)行了探討，開發(fā)了面向高性能科學(xué)計(jì)算的MapReduce系統(tǒng)，該系統(tǒng)支持大規(guī)模計(jì)算的任務(wù)分配和自動(dòng)并行，實(shí)現(xiàn)了矩陣LU分解，計(jì)算效率接近于傳統(tǒng)MPI。遺憾的是，基于云計(jì)算的電磁計(jì)算理論與方法并沒(méi)有建立。國(guó)內(nèi)對(duì)并行電磁計(jì)算的研究主要集中于MoM、FDTD、FEM等幾種典型方法，存在計(jì)算效率不夠高、計(jì)算規(guī)模不夠大等突出問(wèn)題，無(wú)法充分滿足隱身/反隱身等國(guó)防電子裝備的迫切需求[1-3]。

本文基于協(xié)同平臺(tái)的FDTD區(qū)域分解并行算法，添加了區(qū)域分解方法，引入劃分子區(qū)域的虛擬邊界，使用Despres傳輸條件進(jìn)行劃分，對(duì)三維問(wèn)題進(jìn)行分解，通過(guò)求解整個(gè)區(qū)域的有限差分方程，實(shí)現(xiàn)對(duì)整個(gè)區(qū)域的求解。最后，針對(duì)本文提出的算法進(jìn)行算法驗(yàn)證和并行測(cè)試。

1 協(xié)同計(jì)算平臺(tái)架構(gòu)

1.1 協(xié)同系統(tǒng)設(shè)計(jì)

基于“主控中心—子集群”的兩層泛型結(jié)構(gòu)，為基于廣域網(wǎng)的電磁計(jì)算任務(wù)的執(zhí)行提供了一個(gè)高效的計(jì)算平臺(tái)，整個(gè)平臺(tái)通過(guò)瀏覽器（或標(biāo)準(zhǔn)API調(diào)用）為終端用戶提供服務(wù)。用戶在終端提交任務(wù)請(qǐng)求后，經(jīng)過(guò)廣域網(wǎng)傳遞到主控制節(jié)點(diǎn)，主控制節(jié)點(diǎn)通過(guò)業(yè)務(wù)管理模塊獲得任務(wù)信息后，調(diào)用主控中心調(diào)度器依據(jù)當(dāng)前平臺(tái)資源實(shí)時(shí)狀況做出調(diào)度，并將任務(wù)下發(fā)到各個(gè)子中心控制節(jié)點(diǎn)，利用子中心調(diào)度器得到調(diào)度結(jié)果。

1.2 系統(tǒng)的架構(gòu)

在兩層調(diào)度體系中，由主控中心首先對(duì)任務(wù)各屬性（計(jì)算復(fù)雜度、存儲(chǔ)容量要求、通信量、地理位置等）進(jìn)行分析，根據(jù)策略將任務(wù)進(jìn)一步分發(fā)至子集群的控制節(jié)點(diǎn)。主控制中心主要完成任務(wù)接收、任務(wù)解析、資源映射、任務(wù)調(diào)度或分發(fā)、結(jié)果反饋、資源監(jiān)控、用戶認(rèn)證、任務(wù)監(jiān)控與遷移等工作，并不進(jìn)行任務(wù)的具體計(jì)算，所有的計(jì)算工作都是在各個(gè)子計(jì)算中心完成，子計(jì)算中心是真正的電磁計(jì)算任務(wù)運(yùn)行承載者。當(dāng)子集群控制節(jié)點(diǎn)調(diào)度器收到任務(wù)請(qǐng)求后，會(huì)執(zhí)行一個(gè)算法的迭代，然后產(chǎn)生結(jié)果集，返回給用戶。

圖1 計(jì)算任務(wù)協(xié)同平臺(tái)整體架構(gòu)Fig.1 Overall framework of computing tasks on the collaboration platform

1.3 硬件資源

本平臺(tái)的物理設(shè)備主要由HP刀片服務(wù)器集群、曙光刀片服務(wù)器集群以及部分機(jī)架式服務(wù)器以及若干本地/異地PC機(jī)組成，平臺(tái)的網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)如圖2所示。

圖2 網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)Fig.2 The topology structure on network

各物理設(shè)備的型號(hào)、具體硬件配置、用途等如表1所示：

表1 各物理設(shè)備的具體配置信息Table 1 Specific configuration information for each physical device

2 融合區(qū)域分解方法的算法原理

文獻(xiàn)[4]已經(jīng)使用時(shí)域有限差分方法分析了一般三維電磁場(chǎng)問(wèn)題，考慮到針對(duì)大尺寸的三維問(wèn)題求解時(shí)，存儲(chǔ)量和計(jì)算量偏大等問(wèn)題，本文添加了區(qū)域分解方法，用以解決這一問(wèn)題。不失一般性考慮求解區(qū)域沿著y方向（別的方向一樣）進(jìn)行分解（見圖3），各子區(qū)域間無(wú)重疊。總的區(qū)域分為Ω1，Ω2，Ω3，…，Ωm，交界面為Γ1，Γ2，Γ3，…，Γm-1。

圖3 沿z方向區(qū)域分解示意圖Fig.3 Schematic domain decomposition at Z direction

引入劃分子區(qū)域的虛擬邊界，使用Despres傳輸條件進(jìn)行劃分，對(duì)三維問(wèn)題進(jìn)行分解。區(qū)域Ω1，Ω2，Ω3，…，Ωm滿足Maxwell方程，交界面Γ1，Γ2，Γ3，…，Γm-1滿足Despres傳輸條件，同時(shí)，截?cái)噙吔缟蠞M足CPML吸收邊界條件，Yee網(wǎng)格下在無(wú)源區(qū)域，并設(shè)定相對(duì)電導(dǎo)率為1，即μr=1，后面的迭代都將使用這個(gè)結(jié)果進(jìn)行計(jì)算。FDTD算法的處理過(guò)程同前所述，本節(jié)只討論實(shí)施區(qū)域分解的內(nèi)容。

相鄰兩個(gè)區(qū)域Ωm和Ωm+1間的虛擬邊界Γm,m+1，該條件為：

城市公共藝術(shù)是城市的思想，是一種當(dāng)代文化的形態(tài)，擁有良好城市公共藝術(shù)的城市，才是一座有感覺(jué)的城市。為此，我們應(yīng)該憑著對(duì)藝術(shù)的忠誠(chéng)，做好每一件與我們有緣的城市公共藝術(shù)作品，享受藝術(shù)創(chuàng)作帶來(lái)的樂(lè)趣，尋找屬于自己的藝術(shù)作品的符號(hào)，用城市公共藝術(shù)彰顯城市的文化和特質(zhì)，展示城市的形象與魅力。

考慮子區(qū)域Ωm和Ωm+1間的虛擬邊界Γm,m+1，顯然有，令并將其代入上式（1），得到Despres傳輸條件下的方程[5]，即：

兩個(gè)區(qū)域間的虛擬邊界上的Despres傳輸條件則為：

同理可以得到磁場(chǎng)的表達(dá)式。其磁場(chǎng)表達(dá)式為：

式（4）、（5）是子區(qū)域Ωm的上邊界和下邊界Despres傳輸條件的差分公式。通過(guò)這兩組公式有效地連接了上一個(gè)相鄰區(qū)域Ωm-1和下一個(gè)相鄰區(qū)域Ωm+1。采用類似的方法處理所有子區(qū)域Ω1，Ω2，Ω3，…，Ωm，就可以求解出整個(gè)區(qū)域的有限差分方程，實(shí)現(xiàn)對(duì)整個(gè)區(qū)域的求解。

3 測(cè)試比對(duì)與分析

算例1：算法正確性檢驗(yàn)

首先，檢驗(yàn)上述算法的正確性。本文選取與文獻(xiàn)[6]一致的三維金屬導(dǎo)體方塊，對(duì)其電磁散射特性進(jìn)行分析。目標(biāo)電尺寸為20 cm×20 cm×20 cm，介質(zhì)表面涂覆厚度d=2.1 mm，涂覆層的材料電磁參數(shù)為?r=4.8，σ=1.67×10-3s/m，σm=1.42×103Ω/m。采用平面波入射，計(jì)算頻率設(shè)置為f=3 GHz，取δ=Δx=Δy=Δz，將整個(gè)求解目標(biāo)分解為5個(gè)子區(qū)域，再沿著x，y，z方向?qū)?dǎo)體均勻離散成10份。吸收邊界取在第10個(gè)網(wǎng)格外表明處。入射平面波方向見圖4所示，計(jì)算結(jié)果見圖5所示。

圖4 立方體與入射平面波的關(guān)系Fig.4 Relation between the cube and the incident plane wave

圖5 立方體abcd折線上的電流分布Fig.5 Current distribution on a cubeABCD fold line

算例2：算法的并行性能測(cè)試

再以16× 18個(gè)傘形陣列天線為例，采用本文的FDTD并行算法在協(xié)同計(jì)算平臺(tái)上開展該天線的輻射特性計(jì)算。通過(guò)深入探討在不同虛擬拓?fù)浣Y(jié)構(gòu)下FDTD算法的并行性能。單元和陣列模型示意圖如圖6所示。

圖6 16×18陣列幾何模型Fig.6 16×18 array geometry model

仿真參數(shù)為：三個(gè)方向的大小為0.736 m×1.1725 m×0.004 m。共有兩層介質(zhì)板，每層的厚度均為0.002 m，底層為PEC反射板，底層上方為泡沫支撐結(jié)構(gòu)，介電常數(shù)為1，不計(jì)損耗；泡沫上方為偶極子輻射單元，夾層為介質(zhì)材料，介電常數(shù)為2.2，不計(jì)損耗；組陣：橫向16單元，縱向18單元。選擇計(jì)算頻率f=3 GHz，入射波為平面波，整個(gè)計(jì)算區(qū)域的網(wǎng)格大小取Δx=Δy=Δz=0.2 mm，總的迭代步數(shù)為10000步。選取的角度范圍0°≤φ≤360°，角度采樣間隔Δφ=0.5°，對(duì)于吸收邊界使用5層CPML，激勵(lì)使用Gauss基脈沖，計(jì)算的總網(wǎng)格數(shù)為1650×1960×630，總的網(wǎng)格數(shù)為203742000。計(jì)算結(jié)果如下：圖7顯示了288單元陣列天線的電場(chǎng)分布；圖8為288單元陣列天線增益方向圖3D顯示；圖9為288單元陣列天線在phi=0°和phi=90°兩個(gè)面上的輻射方向圖。

圖7 16×18天線陣圖中的電場(chǎng)分布Fig.7 The electric field distribution in 16×18 antenna array

圖8 288單元陣列天線增益3D方向圖Fig.8 Gain 3d pattern of 288 element array antenna

圖9 288單元陣列天線輻射方向圖Fig.9 Radiation pattern of 288 element array antenna

圖9是采用了FDTD并行算法計(jì)算的多單元陣列天線增益方向圖，并將該結(jié)果與商用軟件CST的計(jì)算結(jié)果進(jìn)行對(duì)比。從phi=0°面和phi=90°面的天線增益方向圖可以看出，兩者吻合較好，表明本文的算法可以準(zhǔn)確地解決復(fù)雜結(jié)構(gòu)的電磁場(chǎng)問(wèn)題。

針對(duì)算法的并行效率測(cè)試的參數(shù)設(shè)置及對(duì)應(yīng)的時(shí)間和并行效率值如表2所示。對(duì)比文獻(xiàn)[7]的結(jié)果，可以看出，對(duì)于CPU核數(shù)從32到1024的測(cè)試中，算法并行效率從100%降到了73.67%，雖然不及江樹剛等測(cè)試的結(jié)果，但是也在一定程度上說(shuō)明了本文算法的程序具有較好的并行性能。本文的結(jié)果之所以并行效率沒(méi)有達(dá)到參考文獻(xiàn)的值，可能是在虛擬拓?fù)涞倪x擇、計(jì)算平臺(tái)的硬件性能上的差異導(dǎo)致的。

表2 不同核數(shù)下的并行性能對(duì)比Table 2 Comparison of parallel performance under different cores

表中虛擬拓?fù)涫堑芽栕鴺?biāo)的形式，即x×y×z，8×16×1表示的是在z方向上沒(méi)有拓?fù)浣Y(jié)構(gòu)，虛擬拓?fù)涞木S數(shù)在現(xiàn)有數(shù)字上減1；16×16×2表示在三個(gè)方向上都有拓?fù)浣Y(jié)構(gòu)，所以拓?fù)涞木S度是3。

從表2可以看出，隨著CPU核數(shù)的增加，虛擬的拓?fù)鋽?shù)目在增加，計(jì)算的時(shí)間大幅減少。拓?fù)涞慕Y(jié)構(gòu)對(duì)并行計(jì)算的影響情況將是下面討論的問(wèn)題。本文分別對(duì)在相同CPU核數(shù)、不同維度虛擬拓?fù)浜拖嗤珻PU核數(shù)、相同維度虛擬拓?fù)鋬煞N情況下的并行性能進(jìn)行測(cè)試，討論計(jì)算效率最優(yōu)的問(wèn)題。

考慮到計(jì)算量巨大和計(jì)算資源難以滿足計(jì)算的需求，難以進(jìn)一步研究算法的并行性能，為此，本文將對(duì)三維的虛擬拓?fù)淝闆r進(jìn)行計(jì)算和分析。由于CPU核數(shù)選用32時(shí)，虛擬拓?fù)浣Y(jié)構(gòu)可以出現(xiàn)的組合形式較少，僅有4×4×2、2×4×4和4×2×4三種，所以，本文選用64核進(jìn)行測(cè)試。

表4 CPU 64核、三維不同虛擬拓?fù)湎碌牟⑿行阅軐?duì)比Table 4 Comparison of parallel performances between CPU 64 cores and 3D topology

表4討論了CPU 64核、三維不同虛擬拓?fù)湎碌牟⑿行阅?，從?jì)算數(shù)據(jù)可以看出，虛擬拓?fù)涞牟煌?，?jì)算結(jié)果也有較大的差異，計(jì)算效率也會(huì)有較大的區(qū)別。通過(guò)上面的并行計(jì)算測(cè)試可以得出：在多核計(jì)算系統(tǒng)中，實(shí)現(xiàn)MPI通信下的FDTD并行計(jì)算需要充分考慮虛擬拓?fù)涞慕Y(jié)構(gòu)與子域網(wǎng)格的劃分方式，不同的分割方式會(huì)導(dǎo)致計(jì)算量的差異，影響計(jì)算的效率。

4 討論與總結(jié)

本文基于協(xié)同平臺(tái)開展一種新的FDTD并行算法研究。針對(duì)電大尺寸的三維開域問(wèn)題求解時(shí)，存儲(chǔ)量和計(jì)算量偏大等問(wèn)題，添加了區(qū)域分解方法，引入分子區(qū)域的虛擬邊界，使用Despres傳輸條件進(jìn)行劃分，對(duì)三維問(wèn)題進(jìn)行分解。同時(shí)對(duì)CPU并行效率進(jìn)行了研究，根據(jù)相同CPU核數(shù)、不同維度虛擬拓?fù)浜拖嗤珻PU核數(shù)、相同維度虛擬拓?fù)鋬煞N情況開展并行性能測(cè)試，實(shí)現(xiàn)算法的并行效率的優(yōu)化。同時(shí)也得到：

（1）在進(jìn)行MPI通信下的FDTD并行計(jì)算時(shí)，需要提前確定最優(yōu)的虛擬拓?fù)浣Y(jié)構(gòu)。原則上是三維拓?fù)鋬?yōu)于二維拓?fù)洌S拓?fù)鋬?yōu)于一維拓?fù)?，同時(shí)，在所有三維拓?fù)浣Y(jié)構(gòu)中仍然需要挑選最優(yōu)拓?fù)浞桨福?/p>

（2）在維數(shù)相同的虛擬拓?fù)浣Y(jié)構(gòu)中，拓?fù)浣Y(jié)構(gòu)的選取應(yīng)該遵循一定的規(guī)律，即：拓?fù)浞指詈陀?jì)算網(wǎng)格在形式上要一致，比如，計(jì)算的總網(wǎng)格數(shù)為1650×1960×630，拓?fù)浣Y(jié)構(gòu)x×y×z在分割時(shí)，最好保持y＞x＞z，這樣就有效地減少在子域交界面形成的網(wǎng)格數(shù)量，避免交界面上更多的數(shù)據(jù)通信，提高計(jì)算效率。

[1]Kim KH,Park QH.Overlapping computation and communication of three-dimensional FDTD on a GPU cluster[J].Computer Physics Communications,2012,183(11):2364-2369

[2]Xu L,Xu Y,Jiang RL,et al.Implementation and optimization of three-dimensional UPML-FDTD algorithm on GPU cluster[J].Computer Engineering&Science,2013,35:160-162

[3]Yang CT,Huang CL,Lin CF.Hybrid CUDA,Open MP,and MPI parallel programming on multi core GPU clusters[J].Computer Physics Communications,2011,182(1):266-269

[4]馮 圓,代小霞,唐曉斌,等.基于分布式平臺(tái)的FDTD并行算法[J].北京航空航天大學(xué)學(xué)報(bào),2016,42(9):1874-1883

[5]Despres B.Domain decomposition method and the helmholts Problem[C].Strasbourg France:Proc.Int.Symp Mathemat.Numerical Aspects Wave Propagate Phenomena,1992:44-52

[6]Taflove A, Umashankar K. Radar cross section of general three-dimensional scatters[[J].IEEETransEMC,1983,25(4):433-440

[7]江樹剛,林中朝,張 玉,等.國(guó)產(chǎn)超級(jí)計(jì)算機(jī)實(shí)現(xiàn)10核FDTD并行計(jì)算[J].西安電子科技大學(xué)學(xué)報(bào),2015,42(5):89-90