《二维风暴潮模式并行计算研究——修改稿》由会员分享，可在线阅读，更多相关《二维风暴潮模式并行计算研究——修改稿（8页珍藏版）》请在金锄头文库上搜索。
1、二维风暴潮模式并行计算研究 郑友华 1 李训强 1 朱首贤 2 张文静 1（1、解放军理工大学气象学院；2、河海大学港口海岸与近海工程学院；江苏南京，211100）摘 要：本文通过对串行程序的分析、探讨，提出了基于区域分解的二维风暴潮模式并行算法；通过生成派生数据类型的方式，解决了二维数组子块传送和不同类型数组对应并行化读写的问题；通过子程序全局通信的方式，解决了边界处理的复杂性，提高了模式的可扩展性。并在选取不同的网格规模、不同的进程数的条件下进行了测试，得到了较为理想的结果。关键词：风暴潮 区域分解 全局通信 派生数据类型风暴潮是我国沿海的主要气象海洋灾害之一。随着海岸、近岸和海洋开发越来
2、越多，对提高风暴潮模拟精度以及时空分辨率的要求越来越迫切。目前，广泛用于风暴潮模拟与预报的二维模式大多数是基于串行计算技术的，而计算区域的广阔性、模拟精度和时空分辨率的要求必然使得计算时效性不能满足预报的需要，因此，为解决高分辨率模式所带来的大计算量的预报时效问题非常必要的。现实中所有的大型计算机都是基于并行计算技术的，即便是用于科研或气象海洋预报的 PC 机，也都是多核的，也可以用并行计算技术解决不能充分利用 CPU 的问题。因此，为提高二维风暴潮模式的计算速度，研究开发基于并行计算技术的二维风暴潮模式是最好的选择。本文通过对串行程序的分析，提出了基于区域分解的二维风暴潮并行算法，得到了较为
3、理想的结果。加速比峰值接近 11.0，并行程序与串行程序同步计算结果的相对误差小于0001。1 并行方法及环境介绍并行计算就是研究如何把一个需要非常巨大的计算能力才能解决的问题分成许多小的部分，然后把这些部分分配给许多计算机进行并行处理，最后将这些计算结果综合起来得到最终的结果。在当前并行机上，比较流行的并行编程环境可以分为三类：消息传递、共享存储和数据并行。其中，以 MPI 和 PVM 为代表的消息传递并行编程技术基于大粒度的进程级并行，具有最好的可移植性，几乎被当前流行的各类并行机所支持，具有很好的可扩展性 14 。完全可以满足二维风暴潮模式并行计算的需要。 作者简介：郑友华（1986-）
4、 ，男，重庆大足人，在读硕士生，主要从事海洋动力学与数值模拟研究。本文采用 MPI 消息传递并行编程技术，对二维风暴潮模式进行了并行化处理，并在某计算服务系统上进行模拟调试。该系统采用了基于 Infiniband DDR 互联的高性能并行文件系统，经测试聚合 I/O 写带宽最大可达 1.7GB/s，聚合读带宽可达 140GB/s，完全消除了机群系统的I/O 瓶颈问题。2 串行代码分析将串行程序进行并行化改造前，首先需要对原串行程序进行分析，找出程序中的热点，也就是程序中计算最集中、花费时间最长的地方；然后通过对热点区域的数据相关性、算法的可并行性进行分析，将热点区域进行并行化，以达到最好的并行
5、化效果 8,9。2.1 差分方程李岩等 11在史峰岩等 5引入流速逆变张量建立的可隐式求解的自适应曲线网格二维流场模式的基础上加以改进，考虑更为合理的物理过程。模式增加了气压梯度力、球面曲率、水平湍流扩散和辐射应力物理过程，本文采用该模式。（2.1.1）()()0huhvtxyzz+=（2.1.2）xxuvgTPMF-+-（2.1.3）yytxyz+式中： t 为时间； x,y 为直角坐标系坐标； u,v 分别为沿 x,y 方向的流速分量； h 为海底到静止海面的距离（静水深） ； 为自静止海面向上起算的海面升高； g 为重力加速度； T,P,M,F分别为辐射应力、卷破波切应力、涡动粘性、水底
6、摩擦对单位质量水体产生的作用力，下标为在 x,y 方向的分量 11。图 1 变量空间配置示意图： ： :u :v如图 1 所示，变量空间配置为 C 网格。空间差分采用中央差格式，时间积分采用 ADI 方法。并引入广义曲线坐标变换： ),(),(yx广义速度量为： vxuyJvxudtU1tV， U,V 也称为速度的逆变张量 5。yJ对（2.1.2） 、 （2.1.2）式进行离散，可得形如以下格式的方程组 57 ：（2.1.4）EbaUnn2/12/1（2.1.5）DjiCjiAnjjij 2/1,/,/, ),(),( Ea、 Eb、 A(i,j)、 C(i,j)的表达式和差分方程的离散很容易
7、推导，这里不再详细给出，有兴趣的读者可参考文献6。显然，第 j 行各网格点形如（2.1.5）式的差分方程构成一个三对角矩阵的方程组，可以用追赶法 6求解，得到 ，然后由（2.1.4）式得到 。同理，第 i 列各网格点的差分方程2/1n2/1nU构成一个三对角矩阵的方程组，可以求解得到 ，然后得到 。V2.2 时间诊断通过 intel vtune 性能分析工具对串行程序进行时间分析，分别以计算 08 年 9 月 7 日-10月 7 日（时间步长 180s）金门海域（area1,网格规模为 246238）和东山-大埕湾海域（area2,网格规模为 418418）风暴增水的所用时间为样本，得知，时间
8、积分计算广义速度变量 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 行
9、和第 j+1 行、第 j-1 行之间没2/1n2/1nU2/1n/U有依赖关系；同理，在计算 和 时，第 i 列和第 i+1 列、第 i-1 列之间没有依赖关系。1nV这为并行计算提供了便利。分解的目的是为了分配计算任务，让每个进程负责计算其中的一块子区域。为了提高并行计算的效率，充分利用计算资源，减少等待时间，在网格划分过程中应遵循并行计算的 2 个重要原则，即各处理器负载平衡和减少每个处理器间的通信 8。为此，采用区域分解方法，将 x方向上的各行均匀分配给各进程，计算出 与 ，将所得结果带入 y 方向；同理，将 y2/1n/U方向上的各列均匀分配给各进程，计算出 和 。下面以 y 方向的计
10、算为例对区域分解方法V在本文的实际应用以及通信处理进行阐述。Root 0 Root 2Root 1 Root R-1Root R-2收 集Root R-1广播Root k K=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
11、-1，每个进程按序计算一组。3.2 通信方法为了减少各处理器之间的通信量，采用同步通信方式，即除根进程外，其它各进程先将对应各组的计算结果发送给根进程，根进程集齐了一次迭代（即一个时间层）的计算结果后，再向其他所有进程广播，如图 2。由于进程 R-1 所计算区域数据量不大于其他进程，因此，本文将根进程选为进程 R-1。由于各进程计算出的 和 数据均为二维数组的一部分，类型均为1nV实型（也可设定为双精度型） ，我们首先生成一个新的数据类型 newtype，它是由 T 个数据块组成，每块包含 M 个实型数据，各块偏移 M 个元素，由于 FORTRAN 是列向量储存，传递一个newtype 类型数
12、就相当于传递一个 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(
13、NEWTYPE,ierr)do 100 mpjj=0,numprocs-1if(myid.eq.mpjj)thendo 130 i=mpm*mpjj+1,mpm*(mpjj+1)+1if(i.gt.m)goto 100130 continueif(mpjj.ne.(numprocs-1)thencall 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_
14、COMM_WORLD,ierror)endifendif100 continueif(myid.eq.(numprocs-1)thenif(numprocs.gt.1)thendo mpjj=0,numprocs-2call 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)e
15、nddoendifendifcall 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
16、.4 输入输出处理为了减少通信量，节约用于通信的时间，在输入部分每个进程分别读取。由于在并行子程序中通信已经完成，输出部分不存在通信，采用并行输出可以节省输出时间。本文采用两部非阻塞共享指针并行文件输出方式，为避免各格点数据紊乱，将各网格序号与增水等计算结果对应输出，同样，在输出前首先生成派生数据类型。4 结果检验以某计算服务系统为平台，对并行程序进行性能测试，测试环境如表 1。测试区域分别为area1 和 area2，模拟时间为 08 年 9 月 7 日 00：00-10 月 7 日 18：00（时间步长 180s） ，风场资料采用非对称藤田模型 10建立台风模型风场与背景风场的合成（由风场计算程序提供） ，水深资料采用海图水深与遥感资料的合成，开边界采用南海大区模拟结果。程序计算的测试结果如图 34。计算机组成 技术
