氧碘化学激光器数值模拟中的多块并行通信算法

郭　红　李　艳　安恒斌

1(北京应用物理与计算数学研究所　北京　100094)2(计算物理重点实验室(北京应用物理与计算数学研究所)　北京　100094)(guo_hong@iapcm.ac.cn)

氧碘化学激光器数值模拟中的多块并行通信算法

郭红1,2李艳1安恒斌1,2

1(北京应用物理与计算数学研究所北京100094)2(计算物理重点实验室(北京应用物理与计算数学研究所)北京100094)(guo_hong@iapcm.ac.cn)

摘要为实现氧碘化学激光器大规模数值模拟，基于JASMIN(J parallel adaptive structured mesh applications infrastructure)框架设计实现了氧碘化学激光器数值模拟的多块并行通信算法.该算法针对模拟中构造多块结构网格间通信关系困难、氧碘化学激光器数值模拟填充物理边界复杂以及网格块内块间通信调度策略不统一等通信问题，实现了块间关系自动识别算法用以自动识别多块结构网格间通信关系，构造了特殊物理边界数据结构帮助实现物理边界的填充，并采用了统一的通信调度策略以降低通信时间.数值结果表明：该算法解决了氧碘化学激光器大规模数值模拟中的通信问题，并很好地模拟了氧碘化学激光器装置，使用简单，可以扩展到上千个处理器核.

关键词氧碘化学激光器；大规模计算；多块结构网格；并行通信算法；JASMIN

激光器被广泛地应用于工业、医疗、科研等领域.氧碘化学激光是目前波长最短的化学激光，适合大气传输，原理上和实验上均已证明它发展到高功率和高光束质量是可行的.氧碘化学激光器(chemical oxygen iodine laser, COIL)成为最接近实际应用的激光器之一，其在工业应用方面越来越受到关注，应用范围涉及汽车制造业、造船业、石油及矿业开采业、水下切割和焊接作业等.

随着计算机技术和模拟方法的发展，数值模拟在氧碘化学激光器的研究过程中起了重要作用[1]，与改进硬件和发展新的诊断技术一起被看成改进COIL性能的重要研究手段[2-5].要准确模拟COIL，需要三维Navier-Stokes方程组和详细扩散模型[6].许多国家开展了COIL的三维数值模拟工作，包括俄罗斯[7-8]、日本[9]、美国[10-11]等，美国代表最高水平，其模拟程序COBALT60还在校验中.

氧碘化学激光器模拟是一种三维粘性可压缩化学反应流问题[12-13]，通常具有复杂的几何外形.为适应氧碘化学激光器的复杂几何外形、反映真实流动特征，我国利用多块结构网格(块与块之间交界面上的网格结点一一对应)来离散计算区域，用VICON程序(由俄罗斯激光器件公司开发的计算流体力学程序)对COIL喷管冷流场进行了三维模拟，得出的结果较准确地反映了模拟装置的几何变化[13].然而，要准确模拟氧碘混合喷管内流场，模拟区域需要取得足够大，这就需要大规模并行计算.在氧碘化学激光器的大规模多块并行数值模拟中，将面临3个困难：

1) 构造多块结构网格间通信关系困难.在多块结构网格应用程序中，由于各块网格可采用任意的拓扑结构，这就有可能导致网格块与块间的索引空间不统一，使得相邻面的消息传递异常困难.目前国际上先进的支撑多块结构网格拼接的并行软件，包括SAMRAI[14]，CHOMBO[15]等，根据网格块间旋转平移关系来描述网格块间拼接关系，计算块间拼接索引范围.这种方法扩展性强，适合大规模并行计算.然而，其实现需要用户自行输入每对网格块间旋转平移关系，输入数据通常超过上百行，使用非常困难.尤其对于多块结构网格中出现的非结构拼接，即奇异拼接，实现其通信关系描述更加困难.

2) 填充物理边界复杂.在氧碘化学激光器数值模拟中，由于物理边界的填充与计算方向、流体方向相关，其通信方式与网格内部不同，数据通信需要进行单独处理.因此，物理边界的处理也增加了实现多块结构网格拼接并行程序的复杂性，并且容易影响程序通信性能.

3) 通信调度策略不统一.通信调度策略不统一包括网格块内、块间通信调度策略不统一，以及结构拼接与奇异拼接通信调度策略不统一.对于前者，并行软件中通常采用的通信调度策略局限于单块结构网格，没有对多块情形形成统一模板，即每块网格先各自实现块内通信，然后再进行块间通信.如此一来，通信次数正比于相邻网格块对的个数，影响通信性能.对于后者，由于结构拼接与奇异拼接的通信关系描述不统一，导致2种通信分别实现，同样会增加处理器核等待的时间.

针对上述问题，本文将面向氧碘化学激光器大规模数值模拟，基于JASMIN(J parallel adaptive structured mesh applications infrastructure)框架[16-17]设计多块并行通信算法，基于该算法自动计算多块结构网格块间通信关系，并支撑氧碘化学激光器大规模数值模拟并行程序扩展到上千核.

1控制方程及数值算法

1.1控制方程

控制方程包括质量方程、动量方程、能量方程、组分方程和湍流方程.基于笛卡儿坐标系N种组分控制方程的积分守恒形式为

(1)

其中t为时间.

U=(ρ,ρux,ρuy,ρuz,ρe,ρf1,ρf2,…,ρfN-1)T

(2)

为列守恒变量，ρ为密度，ux,uy,uz为速度分量，fi(i=1,2,…,N-1)为组分i的质量分数，e为单位质量的总能，V为任意气体体积，dS为垂直于包围体积V的边界的外法线向量，F为对流项通量，G为粘性项通量，

H=(0,0,0,0,Q,S1,S2,…,SN-1)

(3)

为源项.Q为能量源项，Si(i=1,2,…,N-1)为组分质量源项.

对流项通量F的具体表达式如下：

F=Fxi+Fyj+Fzk,

(4)

其中：

(ρe+p)ux,ρuxf1,ρuxf2,…,ρuxfN-1)；

(ρe+p)uy,ρuyf1,ρuyf2,…,ρuyfN-1)；

(ρe+p)uz,ρuzf1,ρuzf2,…,ρuzfN-1).

粘性项通量G的表达式为

G=Gxi+Gyj+Gzk,

(5)

其中：

Gx=(0,-σxx,-τxy,-τxz,-uxσxx-

uyτxy-uzτxz+qx,J1x,J2x,…,J(N-1)x)；

Gy=(0,-τxy,-σyy,-τyz,-uxτxy-

uyσyy-uzτyz+qy,J1y,J2y,…,J(N-1)y)；

Gz=(0,-τxz,-τyz,-σzz,-uxτxz-

uyτyz-uzσzz+qz,J1z,J2z,…,J(N-1)z).

其中，σxx,σyy,σzz,τyz,τxy,τxz为应力分量，q为热通量，J为组分扩散通量.

1.2数值算法

氧碘化学激光器数值模拟采用的计算方法是欧拉方法与物理过程分裂[13].在时间步[t,t+Δt]内，整个问题分为3个相继的问题：

(6)

其中,U*(t)=U(t)；

(7)

其中,U**(t)=U*(t+Δt)；

(8)

其中,U***(t)=U**(t+Δt).

U(t+Δt)=U***(t+Δt).

(9)

对于式(6)描述的非粘性问题，采用二阶Roe-ENO显格式[18]或隐式格式逼近.对于式(7)描述的粘性问题，采用二阶加权隐式差分格式逼近.式(8)为化学反应问题，采用显式六阶龙格——库塔方法求解.

物理边界条件采用无滑固壁，即高密度条件(低Knudsen数)下的固体不渗透表面.由于在计算方法中采用物理过程分裂后整个问题分解为3个独立的问题，相应的边界条件也分别给定.

粘性部分给定的物理边界条件是：速度为0，温度为用户给定的壁面温度.

ux=uy=uz=0，

T=TW，

fi=(fi)n,i=1,2,…,N-1，

(10)

其中，下标n指内网格，TW是用户给定的壁面温度，若TW≤0.0，是绝热壁，则T=Tn.

非粘性部分给定的物理边界条件是：速度的法向分量为0，切向分量等于最近的内网格的值.

化学反应过程是在每个网格内各组分之间的反应过程，不涉及物理边界条件.

2多块并行通信算法

本文基于JASMIN框架设计实现了氧碘化学激光器数值模拟的多块并行算法.为了帮助提高Cache命中率，进而提升访存性能，每块网格都被划分为若干独立的小计算区域，这些小计算区域被称为网格片.网格片是数值计算的基本单位，数据通信也是基于网格片实现.

在氧碘化学激光器数值模拟所采用的计算方法中，除化学反应过程是单点计算，非粘性和粘性气流问题均需借助网格片周围网格单元完成数值计算.这需要填充网格片的影像区.根据几何位置的不同，影像区分为3种：网格块内影像区、网格块间影像区以及物理边界影像区.

网格片的块内影像区由同一网格块中相邻网格片进行填充，此时该网格片与相邻网格片的索引空间相同，影像区填充算法同单块结构网格算法，容易实现.而网格块间影像区和物理边界影像区的填充则相对复杂，下面将分别讨论.

2.1网格块间影像区填充方法

本文基于JASMIN框架中用于多块结构网格拼接的块间关系描述方法[19]对网格块间拼接关系进行描述.该方法通过块间索引空间的旋转关系、平移关系来计算网格块间拼接索引范围.块间索引空间的旋转关系是一对邻居网格中索引空间相互关系.图1列出了二维4种右手系索引空间的16种邻居组合，三维则包括24种右手系索引空间的576种组合.平移关系则是由2块网格中空间坐标相同的2个网格结点索引计算得到.

Fig. 1　Translations of two blocks in 2-dimensional space.图1　二维空间旋转关系

通过块间关系描述方法可以计算出2块网格的相交索引范围，进而将通信数据的索引空间由源网格块转换至目的网格块，再填充到相应的索引位置上，实现网格块间拼接.这种方法避免了传统算法中记录块间拼接索引范围的过程，并且可以直接用于描述任意2个网格片的拼接关系，扩展性强，适合大规模并行计算.然而，基于块间关系描述方法直接进行计算，依靠手工输入邻居信息非常复杂，耗时耗力且不易查找错误，离实际应用还有一定差距.采用JASMIN框架中的块间关系自动识别工具可以解决上述问题.块间关系自动识别工具根据结构拼接计算算法实现，利用网格几何坐标计算网格结构拼接的块对之间的旋转关系、平移关系，而对于图2所示的奇异拼接却不适用.

Fig.2　Singularity.图2　奇异拼接

由于在奇异拼接处的数值计算的特殊性，多数并行软件采用了对奇异拼接专门描述的方法，顺序记录奇异拼接所涉及的所有网格块块号，并统一描述其通信关系.奇异拼接的通信关系描述同样非常复杂，依靠手工输入非常困难.针对氧碘化学激光器数值模拟中出现的奇异拼接，在结构拼接计算算法的基础上，设计和实现了块间奇异拼接计算算法.块间奇异拼接计算算法与块间结构拼接计算算法一起组成完备的块间关系自动识别算法.下面以三维算法为例阐述块间关系自动识别算法(算法1).

算法1. 块间关系自动识别算法.

输入: 存储网格结点坐标的文件cfile；

输出: 网格块间旋转、平移关系及奇异拼接关系.

Step1. 读取网格结点坐标并存储到数组

coords中；

if (blocki∩blockj≠∅)

将blocki,blockj的ID插入到

neighbor_list中；

其中，i,j∈{1,2,…,block_num},i≠j；

end if

forneighbor_list{

沿3个坐标轴方向{I,J,K}分别计算blocki与blockj的拼接范围；

if (沿某个坐标轴方向拼接点数量≥2) 记为Xc,其中X∈{I,J,K}；

if (Xc_num==2) {

blocki与blockj面拼接；

任意选取3个不在同一直线的交点计算旋转关系；

} else if (Xc_num==1) {

blocki与blockj棱拼接；

任意选取2个交点计算拼接轴方向的旋转关系；

根据拼接棱所处的坐标面计算另外2个方向的旋转关系；

}else if (Xc_num==0) {

blocki与blockj点拼接；

寻找blocki与blockj共同的邻居blockk,且blockk与blocki,blockj至少1个面拼接；

通过blockk与blocki,blockj的旋转关系计算blocki与blockj的旋转关系；

} end if

选取一个交点计算blocki与blockj的平移关系；

} end for

forneighbor_list{

获取包含同一拼接棱的所有邻居对blocki与blockj的ID，将ID插入到block_list；

清除block_list中相同block的ID；

ifblock_list中网格块数量≠4

奇异拼接；

根据拼接顺序对block_list中ID排序，同时计算拼接范围.

} end for

基于块间关系自动识别算法，用户只需输入存储网格结点坐标的文件名，即可自动生成块间关系，并直接应用于多块结构网格拼接的计算中.存储网格结点坐标的文件，可由大多数网格生成软件得到，例如GRIDGEN[20].

2.2物理边界影像区填充方法

Fig. 3　The description of physical boundaries.图3物理边界描述

Fig. 4　The management of physical boundary conditions.图4物理边界处理方法

物理边界影像区填充方法可以满足氧碘化学激光器多块结构网格数值模拟中物理边界填充的多样化需求.

2.3统一的通信调度策略

由于数据通信基于网格片实现，网格片间通信是并行应用程序的主要时间开销.在氧碘化学激光器多块并行数值模拟中，网格片间通信包括块内网格片间通信、块间网格片间通信以及物理边界网格片间通信.其中，块间网格片间通信又包括结构拼接和奇异拼接通信.

网格片间通信调度策略是实现多块结构网格拼接所面临的主要问题.通信调度不统一，会增加处理器核同步等待时间，从而导致通信性能差.JASMIN框架中的基于网格片的统一模板的通信调度策略[19]，通过将块间旋转关系、平移关系引入通信调度策略，创建块内网格片间及结构拼接块间网格片间统一模板的通信调度；通过构造临时数据结构，使得块间奇异拼接及物理边界填充可以间接利用上述通信调度，从而构造出4种网格片间通信模式的统一通信调度策略.基于该通信调度策略，可以统一实现网格块内块间通信，有效地减少通信中核同步等待时间.

3数值实验

本文模拟了向碘蒸气注入氦气作为稀释气体的COIL混合喷管的流动问题.其中，网格块总数为66，如图5所示，网格单元总数450万.在该算例中，共包括221组结构拼接网格块邻居以及12组奇异拼接邻居.通过块间关系自动识别算法在单个处理器核上计算网格块间通信关系仅需要不到10 min.

Fig. 5　The blocks’ distribution.图5　网格块分布图

图6所示为增益分布图,由图6可以看出，增益最大值出现在中心位置，形成了合适的增益区，延长了气流方向的增益区，为光腔内激光提取做好了准备.

Fig. 6　The distribution of gain of secondary flow.图6　副气流增益分布图

表1所示为固定问题规模450万网格单元，处理器核数在1～2 048范围内变化，运行氧碘化学激光器数值模拟程序10个时间步的数值计算与通信时间总和及加速比.从表1可知，相对于单处理器核，在2 048个处理器核上模拟程序可以得到420倍多的加速.

Table 1　Parallel Speedup

4结论

针对传统的多块并行算法在氧碘化学激光器数值模拟中构造多块结构网格间通信关系困难、网格块内块间通信调度策略不统一以及填充物理边界的特殊性等问题，基于JASMIN框架设计实现了氧碘化学激光器数值模拟的多块并行通信算法.该算法有效地解决了传统多块并行算法中的使用困难及通信性能瓶颈问题.数值实验表明，该算法解决了多块结构网格并行程序实现困难的问题，很好地模拟了氧碘化学激光器装置，在2 048处理器核上的加速比超过420，可以扩展到上千个处理器核.

参考文献

[1]Koop G, Hartlove J, Clendening C, et al. Airborne laser flight-weighted laser module (FLM) and coil moduling support[C]Proc of the 31st AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 2000: Article 2421

[2]Truesdell K A, Helms C A, Hager G D. Coil development in the USA[C]Proc of the 25th AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 1994: Article 2421

[3]Duff E A, Truesdell K A. Chemical oxygen iodine laser (COIL) technology and development[C]Proc of SPIE 5414. Bellingham, WA: SPIE, 2004: 52-68

[4]Deng Liezheng, Shi Wenbo, Du Shuyan, et al. Removal of water vapor content in singlet oxygen flow using cold-jet method[J]. High Power Laser and Particle Beams, 2012, 24(7): 1539-1543 (in Chinese)(邓列征, 石文波, 杜姝彦, 等. 冷射流方法脱除单重态氧气流中的水汽[J]. 强激光与粒子束, 2012, 24(7): 1539-1543)

[5]Zhang Lin, Xian Yuqiang. Design and simulation of supporting stabilizing system for separated COIL resonator cavity[J]. High Power Laser and Particle Beams, 2013, 25(S0): 29-32 (in Chinese)(张林, 鲜玉强. 分体型COIL谐振腔支撑致稳系统设计与仿真分析[J]. 强激光与粒子束, 2013, 25(S0): 29-32)

[6]Madden T J, Solomon W C. A detailed comparison of a computational fluid dynanic simulation and a laboratory experiment for a COIL laser[C]Proc of the 28th AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 1997: Article 2387

[7]Boreysho A S. Integration of COIL in high-power laser systems[C]Proc of SPIE 5777. Bellingham, WA: SPIE, 2005: 16-31

[8]Boreysho A S, Savin A V, Malkov V M. Problems and solutions in COIL gas dynamics[C]Proc of SPIE 5777. Bellingham, WA: SPIE, 2005: 142-148

[9]Sugimoto D, Tei K, Ito T, et al. CFD based optimization of iodine mixing for supersonic COIL[C]Proc of SPIE 5120. Bellingham, WA: SPIE, 2003: 446-451

[10]Madden T J, Miller J H. An analysis of mechanisms of flow unsteadiness in chemical oxygen-iodine laser flowfields[C]Proc of the 35th AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 2004: Article 2728

[11]Miller J H, Shang J S. Validation of parallel Navier-Stokes codes for chemical oxygeniodine laser flowfields[C]Proc of the 33rd AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 2002: Article 1090

[12]Buggeln R C, Shamroth S, Lampson A I, et al. Three-dimensional (3-D) Navier-Stokes analysis of the mixing and power extraction in a supersonic chemical oxygen iodine laser (COIL) with transverseI2injection[C]Proc of the 25th AIAA Plasmadynamics and Lasers Conf. Reston, VA: AIAA, 1994: Article 2435

[13]Li Yan, Li Shouxian, Yu Zhen. Feasibility of numerical simulation of flow through COIL 3-D nozzle with VICON code[J]. High Power Laser and Particle Beams, 2006, 18(5): 736-740 (in Chinese)(李艳, 李守先, 余真. VICON程序模拟氧碘化学激光器3维喷管流场的可行性探讨[J]. 强激光与粒子束, 2006, 18(5): 736-740)

[14]Hornung R D, Wissink A M, Kohn S R. Managing complex data and geometry in parallel structured AMR applications[J]. Engineering with Computers, 2006, 22: 181-195

[15]Colella P, Graves D T, Ligocki T J, et al. Chombo software package for AMR applications design document, LBNL-6616E[R]. Berkeley, CA: Lawrence Berkeley National Labotatory, 2014

[16]Mo Zeyao, Zhang Aiqing, Cao Xiaolin, et al. JASMIN: A parallel software infrasture for scientific computing[J]. Frontiers of Computer Science in China, 2010, 4(4): 480-488

[17]Mo Zeyao, Zhang Aiqing, Cao Xiaolin, et al. User’s guide for JASMIN (2.0 version), T09-JMJL-01[R]. Beijing: Institute of Applied Physics and Computational Mathematics, 2011 (in Chinese)(莫则尧, 张爱清, 曹小林, 等. 并行自适应结构网格应用支撑软件框架(JASMIN 2.0版)用户指南, T09-JMJL-01[R]. 北京: 北京应用物理与计算数学研究所, 2011)

[18]Yang J Y, Hsu C A. High-resolution, nonoscillatory schemes for unsteady compressible flows[J]. AIAA Journal, 1992, 30(6): 1570-1575

[19]Guo Hong, Mo Zeyao, Zhang Aiqing. A parallel module for multiblock structured mesh in JASMIN and its applications. Computer Engineering & Science, 2012, 34(8): 69-74 (in Chinese)(郭红, 莫则尧, 张爱清. JASMIN框架中多块结构网格拼接并行计算及其应用[J]. 计算机工程与科学, 2012, 34(8): 69-74)

[20]Steinbrenner J P, Chawner J R, Fouts C L. The gridgen 3D multiple block grid generation system, WRDC-TR-90-3022[R]. Dayton, OH: Wright Patterson Air Force Base, 1990

Guo Hong, born in 1978. PhD. Member of China Computer Federation. Associate professor in Institute of Applied Physics and Computational Mathematics, Beijing, China. Her research interests include parallel numerical algorithms, large scale communi-cation algorithms, and parallel software infrastures.

Li Yan, born in 1975. Master. Associate professor. Her main research interests include numerical simulation of chemical oxygen-iodine laser and domain decom-position methods.

An Hengbin, born in 1974. PhD. Professor in Institute of Applied Physics and Com-putational Mathematics, Beijing, China. His main research interests include parallel algorithms and applications, domain specific parallel software infrastures.

A Parallel Communication Algorithm in Supersonic COIL’s Calculations Using Multiblock Mesh

Guo Hong1,2, Li Yan1, and An Hengbin1,2

1(InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094)2(LaboratoryofComputationalPhysics(InstituteofAppliedPhysicsandComputationalMathematics),Beijing100094)

AbstractIn this paper, a parallel communication algorithm of supersonic chemical oxygen iodine laser (COIL)’s calculation using multiblock structured mesh has been designed and implemented. This communication algorithm that is designed for large scale supersonic chemical oxygen iodine laser’s simulation is based on JASMIN(J parallel adaptive structured mesh applications infrastructure) infrastructure. There are several communication problems in the large scale supersonic chemical oxygen iodine laser’s simulation such as the complexity of blocks’ connecting relationship description, the complexity of management of boundary conditions and the nonuniform of communication schedule. The communication algorithm in this paper includes a blocks’ relationship recognition algorithm to compute blocks’ connecting relationship automatically, a special data structure to help managing boundary conditions and a unified communication schedule which can reduce communication time. According to our calculation results, the communication problems can be resolved by the communication algorithm in this paper and supersonic chemical oxygen iodine laser based on this communication algorithm can be realized easily and simulated regularly. The simulation with 4.5 million mesh cells in this paper can run efficiently on thousands of processor cores.

Key wordschemical oxygen iodine laser; large scale calculation; multiblock structured mesh; parallel communication algorithm; J parallel adaptive structured mesh applications infrastructure (JASMIN)

收稿日期：2014-12-30；修回日期：2015-06-09

基金项目：国家自然科学基金项目(61370066,11171039)；国家“九七三”重点基础研究发展计划基金项目(2011CB309702)；国家“八六三”高技术研究发展计划基金项目(2012AA01A309)

中图法分类号TP301

This work was supported by the National Natural Science Foundation of China (61370066,11171039), the National Basic Research Program of China (973 Program) (2011CB309702), and the National High Technology Research and Development Program of China (863 Program) (2012AA01A309).