于守兵

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

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

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

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

1 模型控制方程

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

2 并行計算模型設(shè)計

2.1 并行計算模型思路

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

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

2.2 OpenMP并行指令

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

2.2.1循環(huán)并行指令

！＄OMP parallel do和！＄OMP end parallel do 構(gòu)造復(fù)合并行工作區(qū),對do循環(huán)區(qū)域進行并行計算,也即將循環(huán)的次數(shù)按照線程數(shù)進行分解。例如:

！＄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屬性,也即每個線程中的同名中間變量的內(nèi)存都是獨立的,否則會造成多個線程同時調(diào)用同一內(nèi)存,導(dǎo)致結(jié)果出錯。

2.2.2數(shù)組并行指令

！＄OMP parallel workshare和！＄OMP end parallel workshare實現(xiàn)沒有顯式do循環(huán)的數(shù)組并行處理。例如:

！＄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標志區(qū)域內(nèi)的每條語句自動按線程數(shù)進行計算任務(wù)分解。

2.3 并行計算模型代碼實例

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

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

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

！＄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.根據(jù)前面得到的每條邊的通量計算每個控制體的對流通量代碼:

！＄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

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

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

！＄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 運行環(huán)境

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

3.2 算例簡介

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

圖2 算例計算區(qū)域

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

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

3.3 效果分析

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

并行計算中,系統(tǒng)線程數(shù)目不同,相應(yīng)的加速比也不同。理論上講,通過調(diào)用OpenMP線程設(shè)置函數(shù)OMP_set-num_threads()可以設(shè)置任意個數(shù)的并行線程數(shù)?？紤]到工作站固有線程數(shù)為8,為分析加速比隨線程數(shù)目的變化關(guān)系,本次數(shù)值試驗設(shè)置的最大線程數(shù)為8。

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

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

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

4 討 論

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

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

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

5 結(jié) 語

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

OpenMP技術(shù)很適合單臺計算機共享內(nèi)存結(jié)構(gòu)上的并行計算,由于使用線程間共享內(nèi)存的方式協(xié)調(diào)并行計算,在多核計算機上運行的效率很高,內(nèi)存開銷小,編程語句簡潔直觀,很容易實現(xiàn)。隨著計算機硬件性能的提高,OpenMP技術(shù)將得到更廣泛的應(yīng)用。

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

[2]朱星明,涂彬,陳煜,等.水利科學計算并行計算平臺構(gòu)建及算法實踐[J].水利水電技術(shù),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的黃河下游二維水沙數(shù)學模型并行計算研究[J].人民黃河,2005,27(3):49-53.

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

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

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

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

[10]王萬戰(zhàn).黃河河口數(shù)學模擬系統(tǒng)關(guān)鍵技術(shù)研究[R].鄭州:黃河水利科學研究院,2009.

[11]于守兵.淹沒丁壩對水流結(jié)構(gòu)的調(diào)整作用研究[D].南京:南京水利科學研究院,2010.

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

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