馮 圓,唐曉斌,龔曉燕

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

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

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

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

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

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

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

1 協(xié)同計算平臺架構(gòu)

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

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

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

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

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

1.3 硬件資源

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

圖2 網(wǎng)絡(luò)拓撲結(jié)構(gòu)Fig.2 The topology structure on network

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

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

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

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

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

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

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

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

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

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

同理可以得到磁場的表達式。其磁場表達式為：

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

3 測試比對與分析

算例1：算法正確性檢驗

首先，檢驗上述算法的正確性。本文選取與文獻[6]一致的三維金屬導(dǎ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。采用平面波入射，計算頻率設(shè)置為f=3 GHz，取δ=Δx=Δy=Δz，將整個求解目標分解為5個子區(qū)域，再沿著x，y，z方向?qū)?dǎo)體均勻離散成10份。吸收邊界取在第10個網(wǎng)格外表明處。入射平面波方向見圖4所示，計算結(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：算法的并行性能測試

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

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

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

圖7 16×18天線陣圖中的電場分布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并行算法計算的多單元陣列天線增益方向圖，并將該結(jié)果與商用軟件CST的計算結(jié)果進行對比。從phi=0°面和phi=90°面的天線增益方向圖可以看出，兩者吻合較好，表明本文的算法可以準確地解決復(fù)雜結(jié)構(gòu)的電磁場問題。

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

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

表中虛擬拓撲是笛卡爾坐標的形式，即x×y×z，8×16×1表示的是在z方向上沒有拓撲結(jié)構(gòu)，虛擬拓撲的維數(shù)在現(xiàn)有數(shù)字上減1；16×16×2表示在三個方向上都有拓撲結(jié)構(gòu)，所以拓撲的維度是3。

從表2可以看出，隨著CPU核數(shù)的增加，虛擬的拓撲數(shù)目在增加，計算的時間大幅減少。拓撲的結(jié)構(gòu)對并行計算的影響情況將是下面討論的問題。本文分別對在相同CPU核數(shù)、不同維度虛擬拓撲和相同CPU核數(shù)、相同維度虛擬拓撲兩種情況下的并行性能進行測試，討論計算效率最優(yōu)的問題。

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

表4 CPU 64核、三維不同虛擬拓撲下的并行性能對比Table 4 Comparison of parallel performances between CPU 64 cores and 3D topology

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

4 討論與總結(jié)

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

（1）在進行MPI通信下的FDTD并行計算時，需要提前確定最優(yōu)的虛擬拓撲結(jié)構(gòu)。原則上是三維拓撲優(yōu)于二維拓撲，二維拓撲優(yōu)于一維拓撲，同時，在所有三維拓撲結(jié)構(gòu)中仍然需要挑選最優(yōu)拓撲方案；

（2）在維數(shù)相同的虛擬拓撲結(jié)構(gòu)中，拓撲結(jié)構(gòu)的選取應(yīng)該遵循一定的規(guī)律，即：拓撲分割和計算網(wǎng)格在形式上要一致，比如，計算的總網(wǎng)格數(shù)為1650×1960×630，拓撲結(jié)構(gòu)x×y×z在分割時，最好保持y＞x＞z，這樣就有效地減少在子域交界面形成的網(wǎng)格數(shù)量，避免交界面上更多的數(shù)據(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]馮 圓,代小霞,唐曉斌,等.基于分布式平臺的FDTD并行算法[J].北京航空航天大學(xué)學(xué)報,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]江樹剛,林中朝,張 玉,等.國產(chǎn)超級計算機實現(xiàn)10核FDTD并行計算[J].西安電子科技大學(xué)學(xué)報,2015,42(5):89-90