平面二維非均勻泥沙OpenMP并行計算模型

于守兵

(黃河水利科學研究院,河南鄭州 450003)

平面二維水流泥沙模型在研究水流和泥沙運動中發揮著越來越重要的作用,當模型范圍較大、模擬精度要求較高時,計算網格經常達到幾萬甚至上百萬,另外,地形沖淤計算要求的時間尺度一般為幾個月甚至幾年,這些都導致了計算時間的劇增。為了提高計算效率,研究者提出多種并行算法。

常見的并行程序開發模式有兩種:一種是消息傳遞模式,如MPI,在標準的串行程序中加入消息傳遞函數,這種模式需要明確地劃分數據結構和重新編寫源程序,程序實現較為復雜;另一種是共享內存模式[1],如OpenMP,無需作較大修改即可嵌入源程序,簡單靈活。朱星明等[2]構建了基于網絡的大型高性能并行計算共享平臺。Mahinthakumar等[3-4]采用有限單元法結合MPI技術模擬地下水流運動。余欣等[5-6]實現了基于MPI消息傳遞模式的黃河下游二維水流泥沙數學模型并行計算。左一鳴等[7]自主開發了并行通訊平臺,解決了MPI不能實現進程遷移的問題。歐劍等[8]實現了一維河網非恒定流數學模型的并行計算。李 來等[9]采用OpenMP技術對二維水動力有限差分數學模型進行了并行優化試驗。

基于三角形-四邊形混合網格的有限體積水流泥沙模型,對復雜區域邊界具有良好的擬合性,能保持整個計算區域的物理守恒性,方便局部區域加密,求解過程簡單。控制方程的顯式求解特性和分組泥沙之間的獨立性尤其適合采用OpenMP技術進行并行算法設計。本文嘗試建立基于三角形-四邊形混合網格的平面二維非均勻泥沙有限體積并行計算模型。

1 模型控制方程

平面二維非均勻泥沙模型運動控制方程在很多數學模型文獻中都有提及,現采用的泥沙模型只考慮懸沙輸移,控制方程具體形式見文獻[10]。模型采用有限體積法進行離散求解,計算區域采用三角形和四邊形混合網格進行剖分。水流運動方程對流項采用Roe格式求解,詳細過程參見文獻[11]。泥沙運動方程對流項采用迎風格式求解。非均勻泥沙挾沙能力計算采用Hec-6模型的處理方法,床沙級配調整采用文獻[12]的處理方法。

2 并行計算模型設計

2.1 并行計算模型思路

基于混合網格的Roe格式有限體積法采用時間顯式求解,當前時間步水流變量的求解只與上一時間步有關。即每個單元的求解具有獨立性,各單元之間不存在相互影響。非均勻泥沙的通量計算以及分組含沙量Si和分組沖淤厚度Δ zbi的計算也存在這種類似的單元獨立性。因此,算法設計時可以將全體單元按編號順序自動分成若干并行區域,其數目等于設定的線程數目,具體的計算流程見圖1。

圖1 平面二維非均勻泥沙模型并行計算流程

2.2 OpenMP并行指令

OpenMP通過在串行Fortran源程序中加入以！＄OMP或！＄開頭的若干指令實現并行計算。下面主要介紹最常用的兩條指令。

2.2.1循環并行指令

！＄OMP parallel do和！＄OMP end parallel do 構造復合并行工作區,對do循環區域進行并行計算,也即將循環的次數按照線程數進行分解。例如:

！＄OMP parallel do private(rATmp)

do i=1,miNumC

rATmp=cG*rRoughness**2*rUBed(i)/(depNow

(i)**0.333)

uFlux(i)=uFlux(i)-uNow(i)*rATmp+0.0000895*depNow(i)*vNow(i)

vFlux(i)=vFlux(i)-vNow(i)*rATmp-0.0000895*depNow(i)*uNow(i)

uNow(i)=(rTFlowStep*uFlux(i)+depNow(i)*uNow(i))/depNew(i)

vNow(i)=(rTFlowStep*vFlux(i)+depNow(i)*vNow(i))/depNew(i)

end do

！＄OMP end parallel do

上述程序中用到中間變量rATmp,必須加上private屬性,也即每個線程中的同名中間變量的內存都是獨立的,否則會造成多個線程同時調用同一內存,導致結果出錯。

2.2.2數組并行指令

！＄OMP parallel workshare和！＄OMP end parallel workshare實現沒有顯式do循環的數組并行處理。例如:

！＄OMP parallel workshare

rUBed(1:miNumC)=sqrt(uNow(1:miNumC)**2+vNow(1:miNumC)**2)

shrU(1:miNumC)=c G*rUBed(1:miNumC)*rRoughness/depNow(1:miNumC)**0.1667

！＄OMP end parallelworkshare

在！＄OMP標志區域內的每條語句自動按線程數進行計算任務分解。

2.3 并行計算模型代碼實例

下面以平面二維非均勻泥沙數學模型中的水流運動對流項計算模塊和非均勻泥沙計算模塊為例進行說明。

2.3.1水流運動對流項計算模塊

a.采用Roe格式計算通過模型區域內每條邊(界面)的對流通量代碼:

！＄OMP parallel do private(iCelLR,rHUVZbL,rHUVZbR)

do i=1,miNumS

iCelLR(:)=atGrdS(i)%c(:)

rHUVZbL=(/depNow(iCelLR(1)),uNow(iCelLR

(1)),vNow(iCelLR(1)),zb(iCelLR(1))/)

rHUVZbR=(/depNow(iCelLR(2)),uNow(iCelLR

(2)),vNow(iCelLR(2)),zb(iCelLR(2))/)

call SlvFlwByRoe(sidFlux(1:3,i),rHUVZbL,rHUVZbR,atGrdS(i)%norm)

end do

！＄OMP end parallel do

b.根據前面得到的每條邊的通量計算每個控制體的對流通量代碼:

！＄OMP parallel do private(rATmp,j)

do i=1,miNumC

rATmp(1:3)=0.0

do j=1,4

rATmp(1:3)=rATmp(1:3)+arCof(j,i)*sidFlux(:,atGrdC(i)%s(j))*atGrdC(i)%l(j)

end do

hFlux(i)=-rATmp(1)

uFlux(i)=-rATmp(2)

vFlux(i)=-rATmp(3)

enddo

！＄OMP end parallel do

2.3.2泥沙運動計算模塊

a.計算通過每個控制體的分組泥沙通量代碼:

do i=1,miNumSGrp

call CmpCnvcUpwind(SLFlux(:,i),SLG(:,i))

end do

循環變量i表示泥沙分組編號。對于每個分組泥沙,采用CmpCnvcUpwind()函數計算通量時按控制體單元編號進行并行計算。

b.計算每個控制體的挾沙能力和分組泥沙對應的地形沖淤以及床沙級配調整代碼:

！＄OMP parallel do

do i=1,miNumC

call CmpSLTransYu(i,depNow(i),rUBed(i),

depNew(i))

zb(i)=zb(i)+atSdmt(i)%rWhlZbDlt

call UpdateGraduation(i)

end do

！＄OMP end parallel do

c.更新控制體單元水深代碼:

！＄OMP parallel workshare

depNew(1:miNumC)=wlNow(1:miNumC)+zb(1:miNumC)

！＄OMP end parallel workshare

3 算例計算

3.1 運行環境

并行程序在Dell Precision T1500工作站上運行,處理器為Intel(R)Core(TM)i7 CPU四核八線程,主頻為2.8GHz,內存為4G。操作系統為Microsoft WindowsWinXP Professional版本 2002 Service Pack3。采用的編譯器為Intel(R)Visual Fortran 11.1.3471,編譯環境為Microsoft Visual Studio 2008。

3.2 算例簡介

模型區域位于黃河上游包頭段(圖2),區域上游約5km處有昭君墳水文站,河道長約22km,主槽寬約600m,灘地寬2～5km。采用SMS軟件,以20m尺寸進行三角形-四邊形混合網格剖分,共得到21496個三角形單元,59 103個四邊形單元,共計80599個單元。

圖2 算例計算區域

模型上游邊界取鄰近的昭君墳水文站1981年8月22日至10月26日實測流量含沙量過程(圖3)。由于昭君墳水文站缺乏實測懸移質級配資料,為方便起見,將泥沙粒徑按照懸沙中值粒徑和床沙中值粒徑分為兩組。參考距離模型區域最近的上游巴彥高勒水文站和下游頭道拐水文站實測資料,懸沙中值粒徑取0.02mm,床沙中值粒徑取0.20mm。相應的,上游邊界懸沙來沙級配分別為100%和0%,模型區域河床級配分別為0%和100%。下游邊界根據水流流量關系曲線給定水位,其他按自由邊界處理。初始水位給定為模型起算時刻水位,初始流速給定為零,初始含沙量給定為模型起算時刻含沙量。

圖3 算例上游邊界流量含沙量過程

3.3 效果分析

并行計算的效果采用并行加速比[7]衡量。并行加速比SP定義為串行計算耗費時間TS與并行計算耗費時間TP比值。

并行計算中,系統線程數目不同,相應的加速比也不同。理論上講,通過調用OpenMP線程設置函數OMP_set-num_threads()可以設置任意個數的并行線程數。考慮到工作站固有線程數為8,為分析加速比隨線程數目的變化關系,本次數值試驗設置的最大線程數為8。

不同線程數目計算得到的沖淤量都是相同的,主槽為沖刷761萬m3,灘地淤積532萬m3。這與相同模型區域內的實體模型試驗結果主槽為沖刷711萬m3,灘地淤積604萬m3基本一致[13]。不同線程數目完成計算任務的運行時間和相應的加速比見表1。

表1 并行計算模型不同線程數完成計算任務的運行時間和加速比

從表1可以看出:①加速比隨著線程數目的增加而增加,在線程數目等于8時,也即計算機硬件固有的線程數目時,加速比達到最大值1.55。②線程數目增加,加速比的提高幅度是遞減的。在線程數目等于3時,加速比已經達到最大值的92%,而之后再增加5個線程數目,加速比只增加最大值的8%。

4 討 論

理想情況下并行計算模型的加速比SP應等于系統的線程數目。然而,由于平面二維非均勻泥沙計算模型的內在要求使得整個求解過程是多個并行計算模塊按順序組成的串行過程,另外還存在一些只能用串行進行的計算,此外,還存在一些其他并行運行開銷,這些方面共同導致實際的運行效率低于理想情況。

文獻[9]采用二維水動力學OpenMP并行模型模擬丁壩繞流案例,計算網格單元數為67334個,雙線程情況下的平均加速比為1.56。本文模型采用的是非結構網格,而且考慮了非均勻泥沙計算,程序各個模塊之間的聯系開銷花費相對較多。

總之,OpenMP技術通過在串行程序中加上若干并行指令即可在單機上實現并行計算,而且可以獲得一定的計算效率。隨著今后處理器的內部優化和集成多核技術的發展,OpenMP技術能夠進一步地發揮其功能。

5 結 語

基于三角形-四邊形混合網格的平面二維非均勻泥沙Roe格式有限體積模型,由于水流和泥沙運動方程都是顯式求解,每個單元變量計算都只與上一時刻的變量有關,這種網格單元的獨立性尤其適合并行計算。利用OpenMP循環并行指令和數組并行指令,可以很容易地將Fortran串行程序改造為并行程序。案例結果顯示,當線程數目等于計算機固有線程數目時,并行加速比達到最大值1.55。

OpenMP技術很適合單臺計算機共享內存結構上的并行計算,由于使用線程間共享內存的方式協調并行計算,在多核計算機上運行的效率很高,內存開銷小,編程語句簡潔直觀,很容易實現。隨著計算機硬件性能的提高,OpenMP技術將得到更廣泛的應用。

[1]ANANTH G,ANSHUL G,GEORGE K,et al.Introduction to the parallel computing[M].北京:機械工業出版社,2005.

[2]朱星明,涂彬,陳煜,等.水利科學計算并行計算平臺構建及算法實踐[J].水利水電技術,2006,37(8):121-125.

[3]MAHINTHAKUMARG,SAIED F.Implementationand performance analysis of a parallelmulticomponent groundwater transport code [C]Proceedings of theNinth SIAM Conference on Parallel Processing for Scientific Computing.San Antonio:Society for Industrial and Applied Mathematicians,1999.

[4]MAHINTHAKUMAR G,SAIED F.A hybrid MPI-OpenMP implementation of an implicit finite-element code on parallel architectures[J].International Journal of High Performance Computing Applications,2002,16(4):371-393.

[5]余欣,楊明,王敏.基于MPI的黃河下游二維水沙數學模型并行計算研究[J].人民黃河,2005,27(3):49-53.

[6]楊明,余欣,姜愷,等.水動力學數學模型并行計算技術研究及實現[J].泥沙研究,2007(3):1-3.

[7]左一鳴,崔廣柏.二維水動力模型的并行計算研究[J].水科學進展,2008,19(6):846-850.

[8]歐劍,張行南,左一鳴,等.一維河網非恒定流數學模型的并行計算[J].江蘇大學學報:自然科學版,2009,30(5):518-522.

[9]李 來,徐學軍,陳黎明,等.OpenMP在水動力數學模型并行計算中的應用[J].海洋工程,2010,28(3):112-117.

[10]王萬戰.黃河河口數學模擬系統關鍵技術研究[R].鄭州:黃河水利科學研究院,2009.

[11]于守兵.淹沒丁壩對水流結構的調整作用研究[D].南京:南京水利科學研究院,2010.

[12]韋直林,趙良奎,付小平.黃河泥沙數學模型研究[J].武漢水利電力大學學報,1999,30(5):21-25.

[13]田勇,于守兵,顧志剛.包神鐵路黃河大橋改擴建工程實體模型試驗[R].鄭州:黃河水利科學研究院,2011.