于守兵

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

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

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

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

1 模型控制方程

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

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

2.1 并行計(jì)算模型思路

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

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

2.2 OpenMP并行指令

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

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

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

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

2.2.2數(shù)組并行指令

！＄OMP parallel workshare和！＄OMP end parallel workshare實(shí)現(xiàn)沒(méi)有顯式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標(biāo)志區(qū)域內(nèi)的每條語(yǔ)句自動(dòng)按線程數(shù)進(jìn)行計(jì)算任務(wù)分解。

2.3 并行計(jì)算模型代碼實(shí)例

下面以平面二維非均勻泥沙數(shù)學(xué)模型中的水流運(yùn)動(dòng)對(duì)流項(xiàng)計(jì)算模塊和非均勻泥沙計(jì)算模塊為例進(jìn)行說(shuō)明。

2.3.1水流運(yùn)動(dòng)對(duì)流項(xiàng)計(jì)算模塊

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

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

！＄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泥沙運(yùn)動(dòng)計(jì)算模塊

a.計(jì)算通過(guò)每個(gè)控制體的分組泥沙通量代碼:

do i=1,miNumSGrp

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

end do

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

b.計(jì)算每個(gè)控制體的挾沙能力和分組泥沙對(duì)應(yīng)的地形沖淤以及床沙級(jí)配調(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 算例計(jì)算

3.1 運(yùn)行環(huán)境

并行程序在Dell Precision T1500工作站上運(yùn)行,處理器為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 算例簡(jiǎn)介

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

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

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

圖3 算例上游邊界流量含沙量過(guò)程

3.3 效果分析

并行計(jì)算的效果采用并行加速比[7]衡量。并行加速比SP定義為串行計(jì)算耗費(fèi)時(shí)間TS與并行計(jì)算耗費(fèi)時(shí)間TP比值。

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

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

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

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

4 討 論

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

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

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

5 結(jié) 語(yǔ)

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

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

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

[2]朱星明,涂彬,陳煜,等.水利科學(xué)計(jì)算并行計(jì)算平臺(tái)構(gòu)建及算法實(shí)踐[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ù)學(xué)模型并行計(jì)算研究[J].人民黃河,2005,27(3):49-53.

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

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

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

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

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

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

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

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