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

(西安電子科技大學(xué) 陜西省超大規(guī)模電磁計(jì)算重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

隨著無(wú)線電技術(shù)和微電子技術(shù)的飛速發(fā)展，大型設(shè)備和電子電路的電磁兼容問(wèn)題變得越來(lái)越普遍且重要。以電磁場(chǎng)仿真算法為基礎(chǔ)的電磁場(chǎng)仿真軟件作為天線設(shè)計(jì)、電磁兼容等工程問(wèn)題的基礎(chǔ)工具，具有非常重要的實(shí)際意義。在計(jì)算電磁學(xué)領(lǐng)域的眾多電磁算法中，時(shí)域有限差分(Finite Difference Time Domain，F(xiàn)DTD)方法、有限元(Finite Element Method，F(xiàn)EM)方法和矩量法(Method of Moments，MoM)是3種精度較高的經(jīng)典算法。而時(shí)域間斷伽略金(Discontinuous Galerkin Time Domain，DGTD)方法兼具時(shí)域有限差分方法的時(shí)域?qū)掝l帶特性和時(shí)域有限元(Finite Element Time Domain，F(xiàn)ETD)方法的非結(jié)構(gòu)網(wǎng)格特性，具有很高的理論精度，尤其在計(jì)算寬帶場(chǎng)路耦合問(wèn)題和復(fù)雜媒質(zhì)問(wèn)題時(shí)更具優(yōu)勢(shì)[1-7]。

不同于時(shí)域有限元方法求解強(qiáng)耦合單元的特性，時(shí)域間斷伽略金方法使用間斷伽略金法降低了網(wǎng)格間的耦合性，求解過(guò)程中計(jì)算各單元的獨(dú)立矩陣，避免了計(jì)算大型稀疏矩陣，計(jì)算復(fù)雜度遠(yuǎn)小于時(shí)域有限元方法[8]。時(shí)域間斷伽略金方法中常用的顯式迭代方案有l(wèi)eapfrog和ADI等[9-10]，在時(shí)域間斷伽略金方法的顯式迭代并行計(jì)算過(guò)程中，各個(gè)單元在迭代更新時(shí)所需的相鄰單元信息為前一時(shí)刻信息，因此需要與該單元的相鄰單元進(jìn)行信息交換。在并行迭代計(jì)算過(guò)程中，這一信息交換過(guò)程是時(shí)域間斷伽略金方法的主要通信熱點(diǎn)。相鄰單元的數(shù)量會(huì)明顯影響并行計(jì)算過(guò)程中的通信量，即通信負(fù)載。將單元?jiǎng)澐纸o各并行進(jìn)程的過(guò)程決定了各個(gè)進(jìn)程的計(jì)算量，即各進(jìn)程的計(jì)算負(fù)載。通信負(fù)載與計(jì)算負(fù)載共同決定了算法的并行效率。受制于區(qū)域劃分方法的局限性[11-12]，與采用結(jié)構(gòu)性網(wǎng)格進(jìn)行離散的時(shí)域有限差分[13-14]方法相比，在時(shí)域間斷伽略金方法并行計(jì)算的過(guò)程中需要在計(jì)算負(fù)載和通信負(fù)載間做取舍[15-16]。同時(shí)在并行迭代過(guò)程中，進(jìn)程間的數(shù)據(jù)收發(fā)過(guò)程非常頻繁，因此在大規(guī)模并行狀態(tài)下難以取得很高的并行效率[15]。筆者將給出一種適用于大規(guī)模并行的時(shí)域間斷伽略金通信策略。

1 并行DGTD原理與通信策略

無(wú)源麥克斯韋方程組的微分形式中，電場(chǎng)E和磁場(chǎng)H表示為

(1)

(2)

(3)

其中，Z為單元i的阻抗，Y為單元i的導(dǎo)納；Z+為單元i相鄰單元的阻抗，Y+為單元i相鄰單元的導(dǎo)納；E+、H+為單元i相鄰單元的電場(chǎng)和磁場(chǎng)。由于數(shù)值通量的存在，求解式(2)的過(guò)程中需要相鄰單元的電磁場(chǎng)值。對(duì)式(2)進(jìn)行伽略金檢驗(yàn)后得到[9]：

(4)

(5)

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

(6)

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

(7)

其中，αcpn為相鄰進(jìn)程數(shù)，對(duì)于立方體均勻分割的情況，αcpn可視為6。γ為每個(gè)相鄰網(wǎng)格所需的通信量，β為各個(gè)網(wǎng)格每次迭代時(shí)的計(jì)算量。式(7)表明，對(duì)每個(gè)進(jìn)程來(lái)說(shuō)，其通信量與進(jìn)程數(shù)正相關(guān)，計(jì)算量與進(jìn)程數(shù)負(fù)相關(guān)，當(dāng)計(jì)算規(guī)模(進(jìn)程數(shù))擴(kuò)大時(shí)，程序的計(jì)算效率會(huì)因?yàn)橥ㄐ帕康奶嵘档汀Ｋ惴看蔚牡臅r(shí)間滿足

(8)

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

(9)

其中，Tc為基準(zhǔn)核數(shù)的計(jì)算時(shí)間，c為基準(zhǔn)核數(shù)。視Sm、γp、Sc、β為常數(shù)，可以得到：

(10)

(11)

(12)

式(10)表明，在均勻劃分的理想情況下，點(diǎn)對(duì)點(diǎn)通信的DGTD方法迭代過(guò)程的并行效率與并行規(guī)模和通信周期呈負(fù)相關(guān)。其中，Cm是一次迭代中方形結(jié)構(gòu)一個(gè)面上所有網(wǎng)格的通信時(shí)間，Cc是單進(jìn)程情況下整個(gè)模型的總迭代計(jì)算時(shí)間。顯然，在大規(guī)模并行時(shí)想要取得較好的并行效率，就需要減小每次通信的通信周期。

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

實(shí)際應(yīng)用過(guò)程中，一個(gè)進(jìn)程往往與多個(gè)進(jìn)程相關(guān)聯(lián)，順序通信方式下總通信周期數(shù)近似為進(jìn)程總數(shù)，每個(gè)通信周期內(nèi)的通信量遠(yuǎn)小于通信緩沖區(qū)大小，造成了通信資源的浪費(fèi)。

當(dāng)一對(duì)相鄰進(jìn)程進(jìn)行通信時(shí)，與這兩個(gè)進(jìn)程沒(méi)有相鄰關(guān)系的另一對(duì)相鄰進(jìn)程稱為這兩個(gè)進(jìn)程的無(wú)關(guān)通信組。在同一個(gè)通信周期內(nèi)，幾個(gè)無(wú)關(guān)通信組能夠同時(shí)進(jìn)行通信操作，從而避免通信帶寬的浪費(fèi)。如圖3(a)所示，按進(jìn)程號(hào)從小到大的順序，依次分析每個(gè)通信周期內(nèi)的無(wú)關(guān)通信組，能夠得到通信周期最小的通信拓?fù)?，稱為最小周期通信拓?fù)?。以下給出獲取最小周期通信拓?fù)涞膫未a。

最小周期通信拓?fù)鋫未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。

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

2 數(shù)值算例

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

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

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

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

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

3 結(jié)束語(yǔ)

在并行規(guī)模較大時(shí)，點(diǎn)對(duì)點(diǎn)通信的DGTD方法的并行效率主要取決于每次迭代時(shí)通信過(guò)程的通信周期數(shù)。筆者采用最小周期通信策略能夠有效降低通信時(shí)的通信周期，使算法在進(jìn)程間進(jìn)行點(diǎn)對(duì)點(diǎn)通信時(shí)取得更好的并行效率，從而提升DGTD方法在大規(guī)模并行計(jì)算時(shí)的計(jì)算效率，進(jìn)而使其更高效地應(yīng)用于工程問(wèn)題的計(jì)算中。