平面二维非均匀泥沙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.