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

宋永，楊闊

（阿壩師范高等專科學校 四川 汶川 623002）

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

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

1 FDTD并行計算環(huán)境簡介

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

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

1.1 并行計算機——曙光4000A

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

1.2 MPI消息傳遞機制

消息傳遞技術MPI（Message Passing Interface消息傳遞接口）

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

1.3 MPI的基本函數(shù)

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

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

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

3）返回當前進程標識如果有n個進程進程的標識號從0到n-1。

MPI-Comm-rank（comm.，rank）

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

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

5）消息接收：從sounce進程接收 count個數(shù)據(jù)類型為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程序流程設計

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

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

2.1 FDTD的三維基礎

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

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

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

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

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

2.3 算法檢驗

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

2.4 并行性能分析

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

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

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

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

圖3 同一點不同時間步的磁場幅度值Fig.3 At the same time steps of different magnetic field amplitude values

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

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

2.5 并行效率的分析

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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

[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]都志輝.高性能計算并行編程技術-MPI并行程序設計[M].北京：清華大學出版社，2001.

[4]李太全.時域有限差分算法的MPI實現(xiàn)與效率分析 [J].長江大學學報：自科版，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].北京：清華大學出版社，2003.

[6]Quinn M J.MPI與OpenMP并行程序設計 （C語言版）[M].北京：清華大學出版社，2004.