MPI下三維FDTD并行運(yùn)算的分析與實(shí)現(xiàn)

宋永，楊闊

（阿壩師范高等專(zhuān)科學(xué)校 四川 汶川 623002）

1966年K.S.Yee首先提出的一種以Maxwell方程為基礎(chǔ)的解決電磁場(chǎng)問(wèn)題的數(shù)值計(jì)算方法[1]____時(shí)域有限差分法（Fnite Difference Time DomainFDTD）。FDTD算法將Maxwell方程中的兩個(gè)旋度方程直接轉(zhuǎn)化為差分形式，將電磁場(chǎng)進(jìn)行空間和時(shí)間上的離散化，得到電磁場(chǎng)演化的迭代方程組，實(shí)現(xiàn)在一定體積內(nèi)和一段時(shí)間上對(duì)連續(xù)電磁場(chǎng)的數(shù)據(jù)取樣壓縮。

在計(jì)算技術(shù)的迅速發(fā)展下，時(shí)域有限差分法被廣泛接受和應(yīng)用。但是隨著計(jì)算量的增大和問(wèn)題的復(fù)雜化，F(xiàn)DTD的實(shí)現(xiàn)和優(yōu)越性的發(fā)揮也受到硬件的極大限制。原因是其算法穩(wěn)定性要求空間網(wǎng)格的最大尺寸應(yīng)小于波長(zhǎng)的十分之一，計(jì)算時(shí)間又正比于網(wǎng)絡(luò)總數(shù)處理電大目標(biāo)電磁問(wèn)題時(shí)所需內(nèi)存很大。這使得在當(dāng)前的PC機(jī)計(jì)算能力下，F(xiàn)DTD算法僅在低頻范圍獲得了較多應(yīng)用，在分析電大目標(biāo)電磁問(wèn)題中受到很大限制。而采用并行FDTD算法是解決這一問(wèn)題的有效途徑之一，從串行計(jì)算到并行計(jì)算可以使計(jì)算能力大幅度提高。

1 FDTD并行計(jì)算環(huán)境簡(jiǎn)介

FDTD算法本身具有適合于開(kāi)展并行運(yùn)算的優(yōu)勢(shì)[2]：由于FDTD方法中每一網(wǎng)格點(diǎn)的電場(chǎng)（磁場(chǎng)）分量的迭代公式只與它本身上一時(shí)刻的場(chǎng)值和周?chē)W(wǎng)格磁場(chǎng)（電場(chǎng)）上半個(gè)時(shí)間步的值有關(guān)，而與計(jì)算區(qū)域內(nèi)其他遠(yuǎn)端網(wǎng)格點(diǎn)的場(chǎng)量沒(méi)有直接關(guān)系（如圖1）。因而利用FDTD的這種局限性可以將FDTD區(qū)域在空間上分割成若干個(gè)子區(qū)域分別進(jìn)行計(jì)算，這一點(diǎn)與并行運(yùn)算的思路是統(tǒng)一的。

圖1 FDTD算法中Yee氏網(wǎng)格Fig.1 Yee of FDTD arithmetic

1.1 并行計(jì)算機(jī)——曙光4000A

曙光4000系列是面向網(wǎng)格計(jì)算的高性能計(jì)算機(jī)系統(tǒng)，它采用網(wǎng)格技術(shù)、體系結(jié)構(gòu)以構(gòu)件性（Component）、標(biāo)準(zhǔn)性（Standard）、協(xié)作性 （Coorcdinate）為基準(zhǔn)采用了服務(wù)化（Service）、安全化 （Secnrit）、專(zhuān) 業(yè) 化 （Specialization）、智 能 化（Intelligence）的技術(shù)。曙光4000A是曙光系列中計(jì)算速度最快的一種系統(tǒng)[3]。

1.2 MPI消息傳遞機(jī)制

消息傳遞技術(shù)MPI（Message Passing Interface消息傳遞接口）

MPI是一種在專(zhuān)業(yè)計(jì)算機(jī)領(lǐng)域使用得比較廣泛的技術(shù)，是全球工業(yè)、政府和科研部門(mén)聯(lián)合推出的適合進(jìn)程間進(jìn)行標(biāo)準(zhǔn)消息傳遞的并行程序設(shè)計(jì)平臺(tái)；是目前應(yīng)用最廣泛的并行程序設(shè)計(jì)平臺(tái)。消息傳遞就是指并行執(zhí)行的部分之間通過(guò)消息傳遞來(lái)交換信息、協(xié)調(diào)步伐、控制執(zhí)行。消息傳遞一般是面向分布式存儲(chǔ)結(jié)構(gòu)但是它也可以適用于共享存儲(chǔ)的并行機(jī)。

1.3 MPI的基本函數(shù)

MPI有6個(gè)最基本的函數(shù)這6個(gè)基本函數(shù)構(gòu)成了一個(gè)完整的MPI程序所需例程的最小集。

1）MPI初始化：MPI-Init （&arge&argv） 該初始化函數(shù)在MPI應(yīng)用程序的開(kāi)始部分調(diào)用，它使應(yīng)用程序中各個(gè)并行執(zhí)行的任務(wù)獲得初始信息。

2）返回通信域包含的進(jìn)程個(gè)數(shù)：MPI-Comm-size（commsize）

3）返回當(dāng)前進(jìn)程標(biāo)識(shí)如果有n個(gè)進(jìn)程進(jìn)程的標(biāo)識(shí)號(hào)從0到n-1。

MPI-Comm-rank（comm.，rank）

4）消息發(fā)送：將首地址為buf的count個(gè)數(shù)據(jù)類(lèi)型為datatype的數(shù)據(jù)由當(dāng)前進(jìn)程發(fā)送給denst進(jìn)程。

intMPI-Send（void*bufintcountMPI-Datatypedatatypeintdest int tag MOI-Comm comm）

5）消息接收：從sounce進(jìn)程接收 count個(gè)數(shù)據(jù)類(lèi)型為datatpye的數(shù)據(jù)，并保存到首地址為buf的內(nèi)存空間中。

Int MPI-Recv（void*bufint countMPI-Datatype datatype int sourceint tag MPI-Comm commMP-Status*status）

6） MPI結(jié)束：MPI-Finalize（）.

1.4 MPI程序流程設(shè)計(jì)

MPI環(huán)境的初始化和結(jié)束流程：在調(diào)用MPI例程之前各個(gè)進(jìn)程都應(yīng)該執(zhí)MPI_Init進(jìn)行初始化，接著用MPI_Comm_size函數(shù)來(lái)獲取缺省組 （Group）的大小。調(diào)用MPI_Comm_Rank函數(shù)獲取調(diào)用進(jìn)程在缺省組中的邏輯編號(hào)（從0開(kāi)始），然后進(jìn)程可以根據(jù)需要向其他節(jié)點(diǎn)發(fā)送消息或接收其他節(jié)點(diǎn)的消息。經(jīng)常調(diào)用的函數(shù)是MPI_Send和MPI_Revc，最后當(dāng)不需要調(diào)用任何MPI進(jìn)程后調(diào)用MPI_Finalize函數(shù)消除MPI環(huán)境進(jìn)程，此時(shí)可以結(jié)束。

2 MPI下三維FDTD的并行分析和運(yùn)算實(shí)現(xiàn)

2.1 FDTD的三維基礎(chǔ)

在FDTD離散中，電場(chǎng)和磁場(chǎng)各節(jié)點(diǎn)的空間分布采用的是著名的Yee氏網(wǎng)格。這種電磁場(chǎng)分量的空間取樣方式不僅符合法拉第電磁感應(yīng)定律和安培環(huán)路定律的自然結(jié)構(gòu)，而且這種電磁場(chǎng)各分量的空間相對(duì)位置也適合于麥克斯韋方程的差分計(jì)算，能夠恰當(dāng)?shù)孛枋鲭姶艌?chǎng)的傳播特性。

在三維的PML介質(zhì)中每個(gè)場(chǎng)分量分解為兩個(gè)子分量。在笛卡爾坐標(biāo)系中6個(gè)分量則變?yōu)?2個(gè)子分量，記為Exy、Exz、Eyx、Eyz、Ezx、Ezy、Hxy、Hxz、Hyx、Hyz、Hzx、Hzy。 進(jìn)行三維區(qū)域分割后，相鄰空間是在一個(gè)面元上進(jìn)行數(shù)據(jù)變換，這樣數(shù)據(jù)的交換量便增大，同時(shí)時(shí)間并行也變得復(fù)雜。如圖2所示的沿二軸的分割圖，可以看出兩節(jié)點(diǎn)需要交換的數(shù)據(jù)是與分割面平行的場(chǎng)量 Ey、Ez、Ey、Ez。

圖2 三維空間區(qū)域分割后的場(chǎng)分布Fig.2 Regional segmentation of field distribution in three dimensional space after

2.2 用捆綁發(fā)送接收實(shí)現(xiàn)三維FDTD的應(yīng)用

捆綁發(fā)送和接收操作把發(fā)送一個(gè)消息到一個(gè)目的地和從另一個(gè)進(jìn)程接收一個(gè)消息合并到一個(gè)調(diào)用中，源和目的可以是相同的。捆綁發(fā)送、接收操作雖然在語(yǔ)句上等同于一個(gè)發(fā)送操作和一個(gè)接收操作的結(jié)合，但是它們可以有效地避免單獨(dú)書(shū)寫(xiě)發(fā)送或接收操作時(shí)由于次序的錯(cuò)誤而造成的死鎖。這是因?yàn)樵摬僮饔赏ㄐ畔到y(tǒng)來(lái)實(shí)現(xiàn)的，系統(tǒng)會(huì)優(yōu)化通信次序從而有效地避免不合理的通信次序，最大限度地避免死鎖的產(chǎn)生[4]。對(duì)于三維的FDTD并行程序，采用C++類(lèi)函數(shù)來(lái)編程。在編程的過(guò)程中，主要涉及到頭文件的設(shè)置，并行初始化激勵(lì)源的加載，數(shù)據(jù)的傳遞以及對(duì)邊界處的處理。對(duì)于吸收邊界采用的是PML吸收邊界，邊界處電導(dǎo)率的大小與所處位置的關(guān)系式是 σ=σmax（ρ/d）n。σmax為邊界處的最大電導(dǎo)率。ρ、d 分別是吸收邊界的位置以及吸收層的厚度。n是ρ和σmax之間所遵循的關(guān)系。n=1是線性關(guān)系，本程序中采用的n=2的拋物線關(guān)系的。對(duì)于三維FDTD電場(chǎng)量E，磁場(chǎng)量H各有6個(gè)分量分別 Exy、Exz、Eyx、Eyz、Ezx、Ezy、Hxy、Hxz、Hyx、Hyz、Hzx、Hzy。

2.3 算法檢驗(yàn)

在同一微分高斯脈沖激勵(lì)源的前提下，分別采用串行、并行計(jì)算對(duì)同一問(wèn)題。共計(jì)算500個(gè)時(shí)間步，可以得到每一點(diǎn)在不同時(shí)間步采用不同進(jìn)程的場(chǎng)值幅值。對(duì)于同一個(gè)位置不同的時(shí)間步通過(guò)串行、并行得到的數(shù)據(jù)進(jìn)行比較，采用了兩個(gè)點(diǎn)的比較，點(diǎn)（25，30，30）和點(diǎn)（25，50，50）所得到的結(jié)果如圖 3 所示。 左圖是點(diǎn)（25，30，30），右圖是點(diǎn)（25，50，50），由圖示結(jié)果可知該并行程序中傳遞函數(shù)是正確的，整個(gè)程序是合理的。交換邊界處的數(shù)據(jù)能進(jìn)行互相接收和發(fā)送。

2.4 并行性能分析

并行計(jì)算的目標(biāo)是使用p個(gè)處理器執(zhí)行程序并獲得比在單個(gè)處理器上高p倍的性能。順序執(zhí)行時(shí)間和并行執(zhí)行時(shí)間之比稱(chēng)作加速比[5-6]。

加速比（Sp）=串行執(zhí)行時(shí)間（Topt）/并行執(zhí)行時(shí)間（Tp）

并行計(jì)算的效率（也稱(chēng)作處理器使用率）是加速比除以處理器個(gè)數(shù)：

效率（Ep）=加速比（Sp）/p

圖3 同一點(diǎn)不同時(shí)間步的磁場(chǎng)幅度值Fig.3 At the same time steps of different magnetic field amplitude values

為了達(dá)到加速比Sp，并行執(zhí)行時(shí)間必須是串行執(zhí)行時(shí)間的l/p。因?yàn)閮H有p個(gè)處理器，這意味著每個(gè)處理器必須分配到同樣數(shù)量的工作。所有處理器在并行執(zhí)行全過(guò)程中都必須在執(zhí)行程序，而且算法在并行化過(guò)程中不能引入額外的操作。

在并行計(jì)算中，通信開(kāi)銷(xiāo)已經(jīng)成為影響計(jì)算機(jī)效率的一個(gè)原因，這是因?yàn)槊總€(gè)處理機(jī)要獨(dú)立管理自己的進(jìn)程。為了達(dá)到較高的加速比，并行程序面臨的一個(gè)很重要的問(wèn)題就是負(fù)載平衡。它是指合理的在處理器之間重新分配負(fù)載平衡，達(dá)到系統(tǒng)的綜合性能最優(yōu)。一般情況下，根據(jù)各節(jié)點(diǎn)的性能不同應(yīng)該承擔(dān)不同規(guī)模的工作。計(jì)算速度快及內(nèi)存大者應(yīng)承擔(dān)較大的任務(wù)量，反之則承擔(dān)較小的任務(wù)量，這樣可以避免快節(jié)點(diǎn)等待慢節(jié)點(diǎn)同步而造成的計(jì)算時(shí)間浪費(fèi)。負(fù)載平衡過(guò)程必須在平衡計(jì)算能力和通訊開(kāi)銷(xiāo)的基礎(chǔ)上平衡各節(jié)點(diǎn)的負(fù)載。以PML作為吸收邊界的FDTD并行算法則是一種典型的同步并行算法。每一個(gè)時(shí)間步即是控制全局的并行時(shí)鐘，不同的進(jìn)程在每一個(gè)時(shí)間步中都要嚴(yán)格的按照FDTD迭代方式進(jìn)行電場(chǎng)、磁場(chǎng)的數(shù)據(jù)交換。因此每個(gè)進(jìn)程所承擔(dān)的計(jì)算量負(fù)載的平衡是提高FDTD并行算法效率的關(guān)鍵之一。

2.5 并行效率的分析

對(duì)于三維的自由空間，點(diǎn)源輻射可以用數(shù)值比較解決問(wèn)題。文中采用的是計(jì)算500個(gè)時(shí)間步的方式，分別用串行和并行時(shí)間來(lái)比較其效率。消息的發(fā)送和接收采用MPI捆綁式的堵塞式且對(duì)不同的網(wǎng)格大小進(jìn)行計(jì)算，時(shí)間單位是s。其中np=2采用的是一維分割，np=4采用的是二維分割，np=8采用的是三維分割，計(jì)算結(jié)果如表1所示。

對(duì)表中的數(shù)據(jù)進(jìn)行分析：對(duì)于50×50×50串并行所用時(shí)間進(jìn)行比較，隨著進(jìn)程的增多所耗時(shí)間也是成倍的增加。加速比sp＜1，此時(shí)的網(wǎng)格大小只是適用于串行不適合并行計(jì)算；對(duì)于網(wǎng)格 100×100×100、200×200×100、200×200×200 采用串行和并行運(yùn)算。對(duì)于np=2的墻上時(shí)間進(jìn)行比較，可以得到此時(shí)均可以得到sp＞1，則該網(wǎng)格是適合并行算法的；對(duì)于100×100×100，采用不同進(jìn)程所需要的時(shí)間以及得到的加速比進(jìn)行比較可以得到：采用一維分割時(shí)加速比Sp＞1，但是對(duì)于二、三維分割方式，加速比是小于1的。由表1可知，對(duì)于200×200×100采用二維分割能得到最大的加速比。對(duì)于200×200×200的網(wǎng)格，三維分割才可以達(dá)到最大的加速比，其值高達(dá)2.214。這是因?yàn)閴ι蠒r(shí)間由兩部分組成：空間計(jì)算時(shí)間和通信時(shí)間。加速比的高低取決與這兩者的比例，因此可以得出：

表1 不同網(wǎng)格采用不同進(jìn)程運(yùn)行的結(jié)果Tab.1 Different grid using different process operation results

1）當(dāng)計(jì)算的網(wǎng)格數(shù)比較小的時(shí)候，比如50×50×50時(shí)，由于通信也需要一定的時(shí)間導(dǎo)致整個(gè)空間的計(jì)算時(shí)間比通信時(shí)間要小，采用并行算法時(shí)要比串行時(shí)間長(zhǎng)，得到的加速比小于1。采用的并行進(jìn)程越多通信量成倍的增加，則所得到的加速比也是成倍的降低效率的；

2）只有計(jì)算網(wǎng)格達(dá)到一定的數(shù)目，比如 100×100×100時(shí)，采用并行計(jì)算才能優(yōu)于串行計(jì)算達(dá)到加速比Sp＞1。此時(shí)在各個(gè)進(jìn)程之間的通信時(shí)間要小于計(jì)算時(shí)間，計(jì)算時(shí)間是主要的。所以并不是所有的并行計(jì)算都要比串行計(jì)算快的。

3）當(dāng)采用并行計(jì)算的時(shí)，空間的分割方式與網(wǎng)格劃分的多少是有關(guān)的。

①當(dāng)整個(gè)計(jì)算空間不是太大時(shí)，100×100×100通信時(shí)間在整個(gè)墻上時(shí)間占了一個(gè)較大的分量，那么此時(shí)采用最簡(jiǎn)單的一維分割采用較小數(shù)目的進(jìn)程則能減小通信時(shí)間。比如兩個(gè)進(jìn)程時(shí)，此時(shí)得到的加速比是要高十多個(gè)進(jìn)程的。如果此時(shí)采用多個(gè)進(jìn)程將導(dǎo)致加速比Sp＜1的情況出現(xiàn)。

②當(dāng)整個(gè)計(jì)算網(wǎng)格數(shù)大到一定的程度時(shí)，比如200×200×200，由于計(jì)算時(shí)間占了很大的一部分，相對(duì)通信時(shí)間的分量減小，無(wú)論采用幾維的分割方式，Sp＞1的趨勢(shì)是不會(huì)改變的。此時(shí)采用三維的分割方式多個(gè)進(jìn)程并行可以得到最大的加速比和效率。

③采用三維的分割，每維上的進(jìn)程數(shù)是要根據(jù)實(shí)際情況來(lái)選擇的，并不是越多越好，每增加一個(gè)進(jìn)程就是增加一倍的通信時(shí)間。對(duì)于實(shí)際計(jì)算中，計(jì)算步是成千上萬(wàn)的網(wǎng)格數(shù)，也是多于106個(gè)。此時(shí)采用合適的并行計(jì)算，根據(jù)模型采取合適的分割方式以及各個(gè)方向上的進(jìn)程數(shù)顯得尤為重要，恰當(dāng)劃分空間網(wǎng)格可以達(dá)到一個(gè)最大的加速比。

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

文中主要用數(shù)值說(shuō)明了基于MPI消息傳遞機(jī)制的FDTD并行算法分析網(wǎng)格數(shù)、進(jìn)程數(shù)、分割方式，三者對(duì)效率的影響，定量的給出三者的關(guān)系。網(wǎng)格數(shù)比較小的時(shí)候適合串行。只有網(wǎng)格數(shù)達(dá)到一定的數(shù)目，并行才能優(yōu)于串行計(jì)算，而此時(shí)采用不同維數(shù)并行計(jì)算對(duì)效率影響比較大，效率的提高依賴(lài)于網(wǎng)格劃分方式和進(jìn)程數(shù)。文中給出了網(wǎng)格數(shù)以及適合的劃分方式和進(jìn)程數(shù)。這將為以后進(jìn)行并行計(jì)算如何選擇一個(gè)合適的分割方式和進(jìn)程數(shù)提供一個(gè)參考。

[1]Yee K S.Numerical solution of initial boundary value problems involving Maxwell equa-tion in isltropic media[J].IEEE Trans.Antennas Propagat.，1966，14（3）:302-307.

[2]Shlager K L.，Schneider John B.Selective survey of the finitedifference time-domain literature[J].IEEE Antennas&Pro Pagation Magazine，1995，37（4）：53-57.

[3]都志輝.高性能計(jì)算并行編程技術(shù)-MPI并行程序設(shè)計(jì)[M].北京：清華大學(xué)出版社，2001.

[4]李太全.時(shí)域有限差分算法的MPI實(shí)現(xiàn)與效率分析 [J].長(zhǎng)江大學(xué)學(xué)報(bào)：自科版，2007，4（1）:32-34.

LI Tai-quan.Implementation of FDTD in MPI and analysis on efficiency[J].Journal of Yangtze University：Natural Science Edition，2007，4（1）:32-34.

[5]Michael J.Quinn.Parallel Programming in C with MPI and OpenMP[M].北京：清華大學(xué)出版社，2003.

[6]Quinn M J.MPI與OpenMP并行程序設(shè)計(jì) （C語(yǔ)言版）[M].北京：清華大學(xué)出版社，2004.