二维风暴潮模式并行计算研究——修改稿.doc

二维风暴潮模式并行计算研究 作者简介：郑友华（1986-），男，重庆大足人，在读硕士生，主要从事海洋动力学与数值模拟研究。郑友华1 李训强1 朱首贤2 张文静1（1、解放军理工大学气象学院；2、河海大学港口海岸与近海工程学院；江苏南京，211100）摘 要：本文通过对串行程序的分析、探讨，提出了基于区域分解的二维风暴潮模式并行算法；通过生成派生数据类型的方式，解决了二维数组子块传送和不同类型数组对应并行化读写的问题；通过子程序全局通信的方式，解决了边界处理的复杂性，提高了模式的可扩展性。并在选取不同的网格规模、不同的进程数的条件下进行了测试，得到了较为理想的结果。关键词：风暴潮 区域分解 全局通信 派生数据类型风暴潮是我国沿海的主要气象海洋灾害之一。随着海岸、近岸和海洋开发越来越多，对提高风暴潮模拟精度以及时空分辨率的要求越来越迫切。目前，广泛用于风暴潮模拟与预报的二维模式大多数是基于串行计算技术的，而计算区域的广阔性、模拟精度和时空分辨率的要求必然使得计算时效性不能满足预报的需要，因此，为解决高分辨率模式所带来的大计算量的预报时效问题非常必要的。现实中所有的大型计算机都是基于并行计算技术的，即便是用于科研或气象海洋预报的PC机，也都是多核的，也可以用并行计算技术解决不能充分利用CPU的问题。因此，为提高二维风暴潮模式的计算速度，研究开发基于并行计算技术的二维风暴潮模式是最好的选择。本文通过对串行程序的分析，提出了基于区域分解的二维风暴潮并行算法，得到了较为理想的结果。加速比峰值接近11.0，并行程序与串行程序同步计算结果的相对误差小于0001。1 并行方法及环境介绍并行计算就是研究如何把一个需要非常巨大的计算能力才能解决的问题分成许多小的部分，然后把这些部分分配给许多计算机进行并行处理，最后将这些计算结果综合起来得到最终的结果。在当前并行机上，比较流行的并行编程环境可以分为三类：消息传递、共享存储和数据并行。其中，以MPI和PVM为代表的消息传递并行编程技术基于大粒度的进程级并行，具有最好的可移植性，几乎被当前流行的各类并行机所支持，具有很好的可扩展性14。完全可以满足二维风暴潮模式并行计算的需要。本文采用MPI消息传递并行编程技术，对二维风暴潮模式进行了并行化处理，并在某计算服务系统上进行模拟调试。该系统采用了基于 Infiniband DDR 互联的高性能并行文件系统，经测试聚合I/O 写带宽最大可达1.7GB/s，聚合读带宽可达140GB/s，完全消除了机群系统的I/O 瓶颈问题。2 串行代码分析将串行程序进行并行化改造前，首先需要对原串行程序进行分析，找出程序中的热点，也就是程序中计算最集中、花费时间最长的地方；然后通过对热点区域的数据相关性、算法的可并行性进行分析，将热点区域进行并行化，以达到最好的并行化效果8,9。2.1 差分方程李岩等11在史峰岩等5引入流速逆变张量建立的可隐式求解的自适应曲线网格二维流场模式的基础上加以改进，考虑更为合理的物理过程。模式增加了气压梯度力、球面曲率、水平湍流扩散和辐射应力物理过程，本文采用该模式。 （2.1.1） （2.1.2） （2.1.3）式中：t为时间；x,y为直角坐标系坐标；u,v分别为沿x,y方向的流速分量；h为海底到静止海面的距离（静水深）；为自静止海面向上起算的海面升高；g为重力加速度；T,P,M,F分别为辐射应力、卷破波切应力、涡动粘性、水底摩擦对单位质量水体产生的作用力，下标为在x,y方向的分量11。图1 变量空间配置示意图： ： :u :v如图1所示，变量空间配置为C网格。空间差分采用中央差格式，时间积分采用ADI方法。并引入广义曲线坐标变换：广义速度量为：，U,V也称为速度的逆变张量5。对（2.1.2）、（2.1.2）式进行离散，可得形如以下格式的方程组57： （2.1.4） （2.1.5）Ea、Eb、A(i,j)、C(i,j)的表达式和差分方程的离散很容易推导，这里不再详细给出，有兴趣的读者可参考文献6。显然，第j行各网格点形如（2.1.5）式的差分方程构成一个三对角矩阵的方程组，可以用追赶法6求解，得到，然后由（2.1.4）式得到。同理，第i列各网格点的差分方程构成一个三对角矩阵的方程组，可以求解得到，然后得到。2.2 时间诊断通过intel vtune性能分析工具对串行程序进行时间分析，分别以计算08年9月7日-10月7日（时间步长180s）金门海域（area1,网格规模为246238）和东山-大埕湾海域（area2,网格规模为418418）风暴增水的所用时间为样本，得知，时间积分计算广义速度变量U,V和水深H（H=h+）的时间分别约占90.73%、92.17%，风应力和压强梯度力的计算时间约占3.71%、3.87%，计算结果的输出约占1.87%、1.91%。3 并行方案设计从串行程序的分析特点来看，主要是对计算广义速度变量U,V和水深H的子程序以及计算风应力和压强梯度力的子程序以及结果输出进行并行处理。下面就以计算广义速度变量U,V和水深H的子程序的并行化处理为例，阐述本文的并行化方案。3.1 区域分解方法由于第j行各网格点的差分方程构成一个三对角矩阵的方程组，可以用追赶法求解，得到，然后得到。所以，在计算与时，第j行和第j+1行、第j-1行之间没有依赖关系；同理，在计算和时，第i列和第i+1列、第i-1列之间没有依赖关系。这为并行计算提供了便利。分解的目的是为了分配计算任务，让每个进程负责计算其中的一块子区域。为了提高并行计算的效率，充分利用计算资源，减少等待时间，在网格划分过程中应遵循并行计算的2个重要原则，即各处理器负载平衡和减少每个处理器间的通信8。为此，采用区域分解方法，将x方向上的各行均匀分配给各进程，计算出与，将所得结果带入y方向；同理，将y方向上的各列均匀分配给各进程，计算出和。下面以y方向的计算为例对区域分解方法在本文的实际应用以及通信处理进行阐述。Root 0Root 2Root 1Root R-1Root R-2收集Root R-1广播Root kK=0,1,2R-1图2、区域分解及其通信示意图如图2，假设计算格点共有N列，每一列有M个V的计算点（即有M行），程序共调用R个进程。由于在计算V时，网格的每列各组成一个封闭方程组，之间没有相关性。我们将N列按顺序均匀地分成R组，前R-1组每组计算T列数据，T为不小于N/R的最小整数，第R组有N-T(R-1)列，不难看出N-T(R-1)T。从进程0到进程R-1，每个进程按序计算一组。3.2 通信方法为了减少各处理器之间的通信量，采用同步通信方式，即除根进程外，其它各进程先将对应各组的计算结果发送给根进程，根进程集齐了一次迭代（即一个时间层）的计算结果后，再向其他所有进程广播，如图2。由于进程R-1所计算区域数据量不大于其他进程，因此，本文将根进程选为进程R-1。由于各进程计算出的和数据均为二维数组的一部分，类型均为实型（也可设定为双精度型），我们首先生成一个新的数据类型newtype，它是由T个数据块组成，每块包含M个实型数据，各块偏移M个元素，由于FORTRAN是列向量储存，传递一个newtype类型数就相当于传递一个MN数组其中的T列。这既避免了因编译器不同而使传送的二维数组数据产生紊乱，得到不可预知的结果，又避免了因根进程多次进行打解包操作而带来的长时间等待。其MPI+FORTRAN90的并行部分代码如下：call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)mpn=ceiling(n/ numprocs)call MPI_TYPE_VECTOR(mpn,m,m,MPI_REAL,NEWTYPE,ierr)call MPI_TYPE_COMMIT(NEWTYPE,ierr)do 100 mpjj=0,numprocs-1 if(myid.eq.mpjj)then do 130 i=mpm*mpjj+1,mpm*(mpjj+1)+1 if(i.gt.m)goto 100130 continue if(mpjj.ne.(numprocs-1)then call MPI_SEND(vv(1,mpm*mpjj+1),1,NEWTYPE,numprocs-1,mpjj,& &MPI_COMM_WORLD,ierror) Call MPI_SEND(z(1,mpm*mpjj+1),1,NEWTYPE,numprocs-1,numprocs+mpjj,&MPI_COMM_WORLD,ierror)endif endif100 continue if(myid.eq.(numprocs-1)then if(numprocs.gt.1)then do mpjj=0,numprocs-2 call MPI_RECV(vv(1,mpm*mpjj+1),1,NEWTYPE,MPI_ANY_SOURCE, mpjj,&MPI_COMM_WORLD,status,ierror)call MPI_RECV(z(1,mpm*mpjj+1),1,NEWTYPE,MPI_ANY_SOURCE,& &numprocs+mpjj,MPI_COMM_WORLD,status,ierror) enddo endif endif call MPI_BARRIER(MPI_COMM_WORLD,ierror) call MPI_BCAST(vv,m*n,MPI_REAL, numprocs-1 ,MPI_COMM_WORLD,ierror) call MPI_BCAST(z,m*n,MPI_REAL, numprocs-1 ,MPI_COMM_WORLD,ierror)3.3 边界处理为解决边界处理的复杂性，提高程序的可扩展性，本文不再引入虚拟区域从相邻进程接收边界值，而是采用将风场、气压等值以及U、和V、的计算放入不同的子程序中，每个子程序的结尾都采用上述通信方法将结果传送给每个进程。3.4 输入输出处理为了减少通信量，节约用于通信的时间，在输入部分每个进程分别读取。由于在并行子程序中通信已经完成，输出部分不存在通信，采用并行输出可以节省输出时间。本文采用两部非阻塞共享指针并行文件输出方式，为避免各格点数据紊乱，将各网格序号与增水等计算结果对应输出，同样，在输出前首先生成派生数据类型。4 结果检验以某计算服务系统为平台，对并行程序进行性能测试，测试环境如表1。测试区域分别为area1和area2，模拟时间为08年9月7日00：00-10月7日18：00（时间步长180s），风场资料采用非对称藤田模型10建立台风模型风场与背景风场的合成（由风场计算程序提供），水深资料采用海图水深与遥感资料的合成，开边界采用南海大区模拟结果。程序计算的测试结果如图34。计算机组成技术标准计算节点数8节点处理器2颗Intel Xeon 处理器（4 核），8GB内存操作系统Red Hat Enterprise 5.1网络Infiniband DDR聚合I/O 写带宽1.7GB/s聚合I/O 读带宽140GB/s节点之间的MPI带宽1.4GB/sMPI 延迟3.96us编译器Intel Fortran Compiler 64bitMPIIntel MPI for Intel Compiler表1 测试环境图3 加速比与进程数关系图图4 并行效率与进程数关系图从加速比与进程数关系图可以看出：随着进程数的增加，并行加速比逐渐增大，在进程数为18左右时时达到峰值，然后出现下降的趋势，说明进程间的数据交换和同步等待消耗时间逐渐增多，特别是进程数大于峰值进程数时，其通信时间甚至大于每个进程计算时间。并行效率反映加速比随进程数目的增加而增加的比例，往往随着节点数的增加而降低9。从并行效率与进程数关系图中可以看出：当进程数较少时，并行效率随进程数增加而下降较快，而后随进程数增加而缓慢下降，当进程数小于加速比峰值进程数18时，并行效率超过50%，是相当可观的。在进程数较少时，area2的并行效率略低于area1，随着进程数的增加，前者明显高于后者。说明当进程数较少时，通行时间主要受通信数据量影响，受到MPI带宽的限制；计算网格规模增大时，并行部分在整个程序计算中的比重增加。并行程序与串行程序同步计算结果的相对误差小于0001，由于差别微乎其微，这里不再给出。结论通过对串行程序的分析、探讨，提出了基于区域分解的二维风暴潮模式并行算法；并通过生成派生数据类型的方式，解决了二维数组子块传送和不同类型数组对应并行化读写的问题；通过子程序全局通信的方式，解决了边界处理的复杂性，提高了模式的可扩展性。通过测试结果可以看出，当网格较大时，加速比峰值接近11，优于基于矩阵并行分解的算法。该算法没有增加并行的复杂性，也没有降低程序可读性，可以在其它流场模式的并行研究和工程应用中推广。由于ADI方法的特性，U、V的计算必须交替进行，使得各进程在完成一组计算后必须进行同步操作，才能进行下一组的计算，这为动态实现各进程间负载平衡带来了困难，这也是作者今后要重点解决的问题。参考文献：1都志辉高性能计算并行编程技术MPI并行程序设计M北京：清华大学出版社，20012李永兵，并行编程环境与工具.长治学院学报，2009,26(2):41-43.3张林波，迟学斌，莫则尧，等并行计算导论M北京：清华大学出版社，20064陈国良并行计算结构算法编程M北京：高等教育出版社，20035史峰岩,孔亚珍.流速逆变张量隐式求解方法及其在航道港池流场计算中的应用.海洋学报，1998，20(4) :1724.6J.B.T.M.Willemse, G. S. Stelling, G. K. Verboom. Solving The Shallow Water Equations With an Orthogonal Coordinate Transformation. International Symposium on Computational Fluid Dynamics,Tokyo,sept.19857周朦,方国洪.二维长波方程的一个无条件稳定有限差分格式J. 海洋与湖沼, 1988, 19(2): 164172. 8张树生,陶建华.潮流场数学方程交错分区并行计算方法研究. 海洋技术,2005,24(1):72-769王建军,张明进.河道二维水沙数学模型并行计算技术研究. 水道港口,2009,30(3),222-22510朱首贤,沙文钰,丁平兴,陈希.近岸非对称型台风风场模型J. 华东师范大学学报，2002，3：66-71.11 李岩,沙文钰,杨支中,朱首贤. 一次登陆湛江台风风暴潮数值预报.海洋预报，2006，23（1）：27-32Parallel computer technology study on 2D numerical model of Storm TideZHENG Youhua1 LI Xunqiang1 ZHU Shouxian2 ZHANG Wenjing1（1 Instit