非結構網格瞬態電磁場計算中的高效通信方法

李岷軒，江樹剛，吳慶愷，林中朝

(西安電子科技大學 陜西省超大規模電磁計算重點實驗室，陜西 西安 710071)

隨著無線電技術和微電子技術的飛速發展，大型設備和電子電路的電磁兼容問題變得越來越普遍且重要。以電磁場仿真算法為基礎的電磁場仿真軟件作為天線設計、電磁兼容等工程問題的基礎工具，具有非常重要的實際意義。在計算電磁學領域的眾多電磁算法中，時域有限差分(Finite Difference Time Domain，FDTD)方法、有限元(Finite Element Method，FEM)方法和矩量法(Method of Moments，MoM)是3種精度較高的經典算法。而時域間斷伽略金(Discontinuous Galerkin Time Domain，DGTD)方法兼具時域有限差分方法的時域寬頻帶特性和時域有限元(Finite Element Time Domain，FETD)方法的非結構網格特性，具有很高的理論精度，尤其在計算寬帶場路耦合問題和復雜媒質問題時更具優勢[1-7]。

不同于時域有限元方法求解強耦合單元的特性，時域間斷伽略金方法使用間斷伽略金法降低了網格間的耦合性，求解過程中計算各單元的獨立矩陣，避免了計算大型稀疏矩陣，計算復雜度遠小于時域有限元方法[8]。時域間斷伽略金方法中常用的顯式迭代方案有leapfrog和ADI等[9-10]，在時域間斷伽略金方法的顯式迭代并行計算過程中，各個單元在迭代更新時所需的相鄰單元信息為前一時刻信息，因此需要與該單元的相鄰單元進行信息交換。在并行迭代計算過程中，這一信息交換過程是時域間斷伽略金方法的主要通信熱點。相鄰單元的數量會明顯影響并行計算過程中的通信量，即通信負載。將單元劃分給各并行進程的過程決定了各個進程的計算量，即各進程的計算負載。通信負載與計算負載共同決定了算法的并行效率。受制于區域劃分方法的局限性[11-12]，與采用結構性網格進行離散的時域有限差分[13-14]方法相比，在時域間斷伽略金方法并行計算的過程中需要在計算負載和通信負載間做取舍[15-16]。同時在并行迭代過程中，進程間的數據收發過程非常頻繁，因此在大規模并行狀態下難以取得很高的并行效率[15]。筆者將給出一種適用于大規模并行的時域間斷伽略金通信策略。

1 并行DGTD原理與通信策略

無源麥克斯韋方程組的微分形式中，電場E和磁場H表示為

(1)

(2)

(3)

其中，Z為單元i的阻抗，Y為單元i的導納；Z+為單元i相鄰單元的阻抗，Y+為單元i相鄰單元的導納；E+、H+為單元i相鄰單元的電場和磁場。由于數值通量的存在，求解式(2)的過程中需要相鄰單元的電磁場值。對式(2)進行伽略金檢驗后得到[9]：

(4)

(5)

其中，Γlp為與單元i屬于相同進程的相鄰單元面，Γcp為與單元i屬于不同進程的相鄰單元面。以式(5)為基礎，采用不同時間差分近似方法可以得到關于ei和hi的并行迭代方程。同時，在不同進程的交界面單元處，各進程的邊界單元需要與相鄰進程單元交換信息。

(6)

除邊緣進程外(被分配區域處于模型邊緣的進程)，各進程的通信量Pmessage與計算量Pamount之比為

(7)

其中，αcpn為相鄰進程數，對于立方體均勻分割的情況，αcpn可視為6。γ為每個相鄰網格所需的通信量，β為各個網格每次迭代時的計算量。式(7)表明，對每個進程來說，其通信量與進程數正相關，計算量與進程數負相關，當計算規模(進程數)擴大時，程序的計算效率會因為通信量的提升而降低。算法每次迭代所消耗的時間滿足

(8)

其中，T為每次迭代的總時間，Tm為每次迭代的通信時間，Tc為每次迭代的計算時間，Sm為計算平臺每個通信周期的通信速度，γp為各進程在每個通信周期中的平均通信量，Sc為計算平臺的運算速度，N為整個通信過程的非自鎖點對點通信次數(通信周期)。顯然，要減少算法的迭代時間，需要提高算法的計算速度和通信速度。使用圖形處理器(Graphic Processor Unit，GPU)和多核處理器(Many Integerated Core，MIC)卡等異構加速技術和中央處理器(Central Processing Unit，CPU)多線程技術能夠直接而有效地提升算法的計算速度[15-17]；筆者著重討論提升算法通信速度的方式。根據文獻[10]中并行效率的定義，n個進程計算時的并行效率Rn可以表示為

(9)

其中，Tc為基準核數的計算時間，c為基準核數。視Sm、γp、Sc、β為常數，可以得到：

(10)

(11)

(12)

式(10)表明，在均勻劃分的理想情況下，點對點通信的DGTD方法迭代過程的并行效率與并行規模和通信周期呈負相關。其中，Cm是一次迭代中方形結構一個面上所有網格的通信時間，Cc是單進程情況下整個模型的總迭代計算時間。顯然，在大規模并行時想要取得較好的并行效率，就需要減小每次通信的通信周期。

基于消息傳遞接口(Message Passing Interface，MPI)中的MPI_ISEND和MPI_IRECV兩個接口，可以方便地完成相鄰進程的信息交換[12]。為了降低算法在通信時的通信周期，需要在通信過程開始之前規定好各進程點對點通信的順序。最簡單的通信順序為進程號由小到大依次通信的順序通信方式。以圖1(a)所示的進程關聯情況為例，其進程間的相鄰關系可以用圖1(b)所示的關聯矩陣表示。圖1(b)中列號為本地進程號，行號為相鄰進程號，進程間有相鄰關系的元素為1，無相鄰關系的元素為0。采用順序通信方式后，各個進程的通信順序即通信拓撲如圖1(c)所示，各個矩陣元素值代表對應進程之間所屬的通信周期數；該策略通信過程中的實際通信情況如圖2所示，圖中由虛線箭頭連接的方框表示一組收發過程，橫軸為通信周期數。

實際應用過程中，一個進程往往與多個進程相關聯，順序通信方式下總通信周期數近似為進程總數，每個通信周期內的通信量遠小于通信緩沖區大小，造成了通信資源的浪費。

當一對相鄰進程進行通信時，與這兩個進程沒有相鄰關系的另一對相鄰進程稱為這兩個進程的無關通信組。在同一個通信周期內，幾個無關通信組能夠同時進行通信操作，從而避免通信帶寬的浪費。如圖3(a)所示，按進程號從小到大的順序，依次分析每個通信周期內的無關通信組，能夠得到通信周期最小的通信拓撲，稱為最小周期通信拓撲。以下給出獲取最小周期通信拓撲的偽代碼。

最小周期通信拓撲偽代碼：

DO WHILE NotAllCommunicationFinished

NowPeriod++；

DO send_rank=0，mpi_size

IF comm_check(send_rank).equal..false. THEN

DO recv_rank=ConnectProcessor0，ConnectProcessorN；

IF comm_check(recv_rank).equal..false. THEN

comm_check(send_rank)=.true.；CommMatrix(send_rank，recv_rank)=NowPeriod；

comm_check(recv_rank)=.true.；CommMatrix(recv_rank，send_rank)=NowPeriod；exit；

ENDIF

ENDDO

ENDIF

ENDDO

ENDDO。

最小周期通信的實際通信過程如圖3(b)所示。相對于點對點通信中簡單的順序通信方式，最小周期通信方式的總通信周期更少；根據式(10)，其理論并行效率也更高。

2 數值算例

以圖4(a)所示的矩形波導縫隙天線輻射特性計算為例，驗證方法的有效性。矩形波導尺寸為710 mm×13.2 mm×6.5 mm，使用同軸線在波導中心處饋電，采用頻帶為14～20 GHz的調制高斯脈沖激勵源激勵，計算時使用由Intel Xeon Gold 6140 CPU @ 2.30 GHz處理器組成的計算集群；該集群共16節點，每節點含兩塊CPU處理器及72 GB內存；每個處理器核心數為18，最大線程數為36。計算時，給每個處理器核心分配兩個進程。其中，縫隙天線16 GHz的xOz面和yOz面增益如圖4(b)和圖4(c)所示。表1給出了使用順序通信和最小周期通信的計算情況。可以看出，在不同計算核心數并行規模的情況下，順序通信的通信周期遠大于最小周期通信，導致并行效率大大降低，這一點在大規模計算時更為突出。

表1 矩形波導縫隙天線計算情況(迭代663 217步)

文獻[15]采用非阻塞的順序通信策略，在計算1.6×106網格量的模型時，由6進程擴展至30進程(5倍)時效率降至90%以下。文中算例網格量為1.1×106，與文獻[15]中算例有一定相似性。文中算例在核心數為936的計算過程中，每個進程平均計算的網格量為1 176，每對進程的單次通信的網格量約為121個單元的信息。以120進程為基準，采用最小周期通信策略擴展8倍的并行效率達到了約89%，同時其通信周期約為順序通信的3%，這表明，采用最小周期通信策略能有效降低程序的通信周期，從而提升算法的并行效率。

以圖5(a)所示的汽車模型為例，驗證并行規模較大時，最小周期通信策略的并行效率。汽車尺寸為4.8 m×1.8 m×1.45 m，計算頻率為1 GHz，計算平臺為神威·太湖之光超級計算機，激勵源為水平極化平面波，入射方向沿z軸負向，極化方向沿x軸正向，計算θ=90°，φ由0°到360°的xOy面雙站散射截面(Radar Cross Section，RCS)。需要注意的是，神威·太湖之光超級計算機采用由申威26010異構眾核處理器組成的計算節點，該處理器由1個主核和64個從核組成，稱為一個核組，每個節點包含4個核組[18]。計算過程中，每個核組作為一個進程共同參與計算。DGTD方法在申威26010異構眾核處理器上運行時需要對主核和從核分別進行編程[13]，而筆者主要討論核組間的通信過程，對此不多敘述。

圖5(b)展示了DGTD與MoM的RCS計算結果。該算例迭代10 000步的計算時間和并行效率如圖6所示。以800進程為基準，當并行規模擴展到8 000進程(并行規模擴展10倍)時，并行效率約為70.38%。根據式(8)，圖6中的通信時間受通信周期數和進程數增加的共同影響而逐漸上升，這是影響其并行效率的主要因素。

3 結束語

在并行規模較大時，點對點通信的DGTD方法的并行效率主要取決于每次迭代時通信過程的通信周期數。筆者采用最小周期通信策略能夠有效降低通信時的通信周期，使算法在進程間進行點對點通信時取得更好的并行效率，從而提升DGTD方法在大規模并行計算時的計算效率，進而使其更高效地應用于工程問題的計算中。