基于協同計算的FDTD區域分解并行算法研究

馮 圓,唐曉斌,龔曉燕

1.中國電子科學研究院,北京 100041

2.空軍預警學院,湖北 武漢 430019

3.火箭軍指揮學院,湖北 武漢 430013

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

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

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

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

1 協同計算平臺架構

1.1 協同系統設計

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

1.2 系統的架構

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

圖1 計算任務協同平臺整體架構Fig.1 Overall framework of computing tasks on the collaboration platform

1.3 硬件資源

本平臺的物理設備主要由HP刀片服務器集群、曙光刀片服務器集群以及部分機架式服務器以及若干本地/異地PC機組成，平臺的網絡拓撲結構如圖2所示。

圖2 網絡拓撲結構Fig.2 The topology structure on network

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

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

2 融合區域分解方法的算法原理

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

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

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

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

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

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

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

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

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

3 測試比對與分析

算例1：算法正確性檢驗

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

圖4 立方體與入射平面波的關系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并行算法在協同計算平臺上開展該天線的輻射特性計算。通過深入探討在不同虛擬拓撲結構下FDTD算法的并行性能。單元和陣列模型示意圖如圖6所示。

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

仿真參數為：三個方向的大小為0.736 m×1.1725 m×0.004 m。共有兩層介質板，每層的厚度均為0.002 m，底層為PEC反射板，底層上方為泡沫支撐結構，介電常數為1，不計損耗；泡沫上方為偶極子輻射單元，夾層為介質材料，介電常數為2.2，不計損耗；組陣：橫向16單元，縱向18單元。選擇計算頻率f=3 GHz，入射波為平面波，整個計算區域的網格大小取Δx=Δy=Δz=0.2 mm，總的迭代步數為10000步。選取的角度范圍0°≤φ≤360°，角度采樣間隔Δφ=0.5°，對于吸收邊界使用5層CPML，激勵使用Gauss基脈沖，計算的總網格數為1650×1960×630，總的網格數為203742000。計算結果如下：圖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并行算法計算的多單元陣列天線增益方向圖，并將該結果與商用軟件CST的計算結果進行對比。從phi=0°面和phi=90°面的天線增益方向圖可以看出，兩者吻合較好，表明本文的算法可以準確地解決復雜結構的電磁場問題。

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

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

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

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

考慮到計算量巨大和計算資源難以滿足計算的需求，難以進一步研究算法的并行性能，為此，本文將對三維的虛擬拓撲情況進行計算和分析。由于CPU核數選用32時，虛擬拓撲結構可以出現的組合形式較少，僅有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核、三維不同虛擬拓撲下的并行性能，從計算數據可以看出，虛擬拓撲的不同，計算結果也有較大的差異，計算效率也會有較大的區別。通過上面的并行計算測試可以得出：在多核計算系統中，實現MPI通信下的FDTD并行計算需要充分考慮虛擬拓撲的結構與子域網格的劃分方式，不同的分割方式會導致計算量的差異，影響計算的效率。

4 討論與總結

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

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

（2）在維數相同的虛擬拓撲結構中，拓撲結構的選取應該遵循一定的規律，即：拓撲分割和計算網格在形式上要一致，比如，計算的總網格數為1650×1960×630，拓撲結構x×y×z在分割時，最好保持y＞x＞z，這樣就有效地減少在子域交界面形成的網格數量，避免交界面上更多的數據通信，提高計算效率。

[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].北京航空航天大學學報,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]江樹剛,林中朝,張 玉,等.國產超級計算機實現10核FDTD并行計算[J].西安電子科技大學學報,2015,42(5):89-90