CN110543663A

CN110543663A - 一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法

Info

Publication number: CN110543663A
Application number: CN201910662152.0A
Authority: CN
Inventors: 张兴军; 邹年俊; 伍卫国; 李靖波; 周剑锋; 韩立; 董小社
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2019-07-22
Filing date: 2019-07-22
Publication date: 2019-12-06
Anticipated expiration: 2039-07-22
Also published as: CN110543663B

Abstract

一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，包括以下步骤：步骤1，构建目标网格所有的区域划分方式集合；步骤2，为每个区域划分方式Mx，My，Mz，threads四元组建立通信拓扑；步骤3，根据MPI通信面的大小对区域划分方式四元组进行过滤分类；步骤4，对S1的几何参数Mx，My，Mz执行面向cache‑aware的区域划分策略；步骤5，度量相同区域划分方式的不同进程、线程组合的性能差异；步骤6，计算出开销最小的区域划分方式四元组Mx，My，Mz，threads，将其作为输出的区域划分方式。本发明适用于粗粒度MPI+OpenMP混合并行方式。与其他的快速划分方式相比较，如通信量最小，往往能够找到更优的区域划分方式；与性能模拟器相比较，则具有高效的划分选择速度。

Description

一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法

技术领域

本发明涉及计算机程序的数据划分和布局优化领域，特别涉及一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法。

背景技术

计算流体动力学(Computational Fluid Dynamics，CFD)是一门交叉学科，它借助于数值计算理论和高性能计算技术对流体力学的问题进行数值模拟和分析。随着计算机领域相关技术的飞速发展，计算流体力学研究的问题越来越深入，导致求解所需要的计算规模、存储规模都远超单台计算机所的性能极限，因此当前的大规模数值模拟都往往运行在分布式集群上。在集群上运行CFD程序，面临的第一个问题就是如何将原本的网格数据进行合理的分解，以充分利用多台计算节点的全部性能，实现任务的高效并行。区域分解方法(Domain Decomposition Method，DDM)为计算流体力学模拟的并行求解提供了有效途径，对数据集的分解可以归结为对网格数据的区域分解。在相同的并行规模下，不同的区域分解方式虽然可以拥有相同规模的子求解域，但这些子求解域的数据布局不一样，可能引起计算过程的访存性能波动。同时，分布式CFD程序需要借助于MPI完成不同计算节点的任务之间的通信，数据布局引起通信面、消息大小和消息数目的改变，也会造成不同的通信开销。

为了提高CFD程序的求解效率，研究者们使用MPI、OpenMP、CUDA、OpenCL 等并行编程框架结合CPU、GPU、MIC等硬件设备进行了大量研究。但为了保证代码的可移植性，很多程序仅支持CPU的并行计算。对于这一类的代码，主要采用MPI和OpenMP 两种编程框架。对于这类并行程序，节点间并行采用消息通信，而节点内既可以采用消息通信，也可以采用共享内存。为了减少通信开销，可以采用粗粒度MPI+OpenMP混合并行计算。在这种并行方式中，每一个线程绑定一个CPU核，执行一个子网格块区域的求解，而MPI进程仅用于分配线程，管理资源和建立MPI通信。它理论上具有更好的可扩展性、易读性、可维护性。采用粗粒度MPI+OpenMP混合并行方式，则意味着每个网格块的通信面可以分类为进程间的MPI通信面和进程内的OpenMP共享内存通信面。不同的进程、线程组合带来不同的通信方式，从而有不同的通信开销。

另外，如何对网格进行划分也是一个需要仔细考虑的问题。本方法针对的是六面体结构网格，划分之后得到的子网格块是网格规模近似的六面体结构网格，以达到负载均衡的目的。因此至多存在六个通信面，可以分为三类，分别是X-Y面、X-Z面和Y-Z面。通常根据每个通信面的估计通信量来评判区域划分方式的优劣，但是这种方法较为粗糙，只能对通信量差别较大的区域划分方式起作用。也可以通过建立通信模型的方式，先获得程序运行的轨迹，再在模拟器上重现和放大来获得不同区域划分方式的性能差异，但这种方法时间开销比较大，对程序的插桩效果也受到平台的影响，得到的结果并非显著优于其他方法，即代价较大而收益不足。可以发现，通信面在通信时往往会经历数据打包、缓冲的过程，不同区域划分方式的数据布局存在不同的cahce失效和数据局部性。但是从何种角度来度量通信打包缓冲过程中的数据局部性并将其利用起来的研究比较少。这种建立性能模型的方法既考虑到了数据局部性，求解的时间开销也少，逐渐成为一种新的研究途径。此外，区域划分方式的研究往往都是针对纯MPI并行方式，极少有针对粗粒度 MPI+OpenMP混合并行方式的区域划分方式研究。

发明内容

本发明的目的在于提供一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，以解决上述问题。

为实现上述目的，本发明采用以下技术方案：

一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，包括以下步骤：

步骤1，构建目标网格所有的区域划分方式集合：子网格X方向数值M_x，子网格Y 方向数值M_y，子网格Z方向数值M_z，进程内的线程数目threads；定义集合为S；

步骤2，为每个区域划分方式M_x，M_y，M_z，threads四元组建立通信拓扑；

步骤3，根据MPI通信面的大小对区域划分方式四元组进行过滤分类，将进入下一步的区域划分方式集合记为S₁。

步骤4，对S₁的几何参数M_x，M_y，M_z执行面向cache-aware的区域划分策略；

步骤5，根据混合并行的通信过程和步骤2得到的通信拓扑建立最晚结束时刻模型，求出缩放因子，度量相同区域划分下不同进程、线程组合的通信性能差异；

步骤6，将步骤5得到的缩放因子与步骤4结合起来，计算出通信开销最小的区域划分方式四元组M_x，M_y，M_z，threads，将其作为输出的区域划分方式。

进一步的，步骤1具体包括：按照负载均衡的简单六面体区域分解方法来对原始网格进行区域分解，分解成若干个三维大小均为M_x×M_y×M_z的六面体网格；目标网格所有的区域划分方式M_x，M_y，M_z，threads四元组构成一个集合S。

进一步的，步骤2具体包括：按照区域分解时定义的规则，子网格块几何序号，进程号和线程号之间存在通过计算得到的映射关系，从而建立起某个区域划分方式的通信拓扑。

进一步的，步骤3具体包括：根据步骤2获得的通信拓扑，统计该区域划分方式产生的MPI通信面消息大小是否超过当前平台测得的通信消息阈值；如果所有的区域划分方式都存在MPI通信面消息大小超过了通信消息阈值的情况，则过滤失效，全部区域划分方式都进入步骤4；否则，只有MPI通信面消息大小未超过通信消息阈值的区域划分方式进入步骤4；将进入步骤4的区域划分方式集合记为S₁。

进一步的，步骤4具体包括：根据L2cahce对子网格块区域的容纳情况分成三种情况：

1)子网格块区域满足

X-Y面总的访存开销：

X-Z面总的访存开销：

Y-Z面总的访存开销：

α+(M_y×M_z-α)×ρ

总的访存开销：

2)子网格块区域满足：

X-Y面总的访存开销：

M_x×M_y

X-Z面总的访存开销：

Y-Z面总的访存开销：

αM_y+(M_y×M_z-αM_y)×ρ

总的访存开销：

3)子网格块区域满足： X-Y面总的访存开销：

M_x×M_y

X-Z面总的访存开销：

M_x×M_z

Y-Z面总的访存开销：

M_y×M_z

总的访存开销：

M_x×M_y+M_x×M_z+M_y×M_z；

根据这三种情况来对每种区域划分方式的M_x，M_y，M_z进行访存开销估算，获得S₁

中每种区域划分方式的性能优劣。

进一步的，步骤5具体包括：

模型反映MPI通信时间和OpenMP共享内存通信时间的重叠程度；首先求解结束OpenMP共享内存通信的时刻t_fixed；接着确定每个MPI通信结束的时刻t_{mpi_i}；然后比较每一个MPI非阻塞接收检查调用的开始时刻t_{ws_i}与该MPI通信接收真实结束时刻t_{mpi_i}的先后，以时刻更靠后的值作为接收结束时刻t_{we_i}；最后，当求出了最后一个MPI通信面的接收结束时刻t_we，再加上这个通信面的拷贝开销，便得到了最晚结束时刻t_final，从而求出缩放因子反映线程数为threads与线程数为1的性能比值。

与现有技术相比，本发明有以下技术效果：

本发明的目的在于通过对cache建立性能模型并度量混合并行时的通信重叠，从数据局部性优化的角度和计算与通信重叠的角度来研究结构网格在粗粒度MPI+OpenMP混合并行时的区域划分选择策略，从而达到求解时更好的点对点通信性能。为面向粗粒度MPI+OpenMP混合并行的结构网格CFD程序提供了高效的区域划分选择策略，同时从通信开销最小和数据局部性两个角度来度量点对点通信函数的时间开销，从而可以找到点对点通信开销较小的区域划分方式。首先，这个方法适用于粗粒度MPI+OpenMP混合并行方式。其次，与其他的快速划分方式相比较，如通信量最小，往往能够找到更优的区域划分方式；与性能模拟器相比较，则具有高效的划分选择速度。

附图说明

图1为面向混合并行的区域划分选择策略的执行流程；

图2为不同通信面的访问地址连续性；

图3为一个任务执行混合并行通信时间的时间顺序。

具体实施方式

以下结合附图对本发明进一步说明：

请参阅图1至图3，一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，包括以下步骤：

1)构建目标网格所有可能的区域划分方式集合(M_x，M_y，M_z，threads)。首先需要按照负载均衡的简单六面体区域分解方法来对原始网格进行区域分解，分解成数个子网格块，保证每个子网格块区域的网格数目相同，三维大小一致，均为M_x×M_y×M_z的六面体。由于同一个原始网格在三个方向上都具备多种划分选择，产生了多种可以表示为 (D_x，D_y，D_z)的划分组合。当限定了某一个并行规模D时，并行规模与划分组合的关系需要保证为

D＝D_x×D_y×D_z (1)

式(1)中：

D_i——i方向的划分数目；

而划分数目D_i与原始网格和划分后子网格块的数值大小关系为

式(2)中：

M_i——i方向的子网格块区域数值大小；

N_i——i方向的原始网格数值大小；

通过公式(1)枚举划分组合，并通过公式(2)对原始网格展开分解计算，可以得到某一个并行规模D下，所有可能的区域划分方式(M_x，M_y，M_z)的集合S`。由于采用了粗粒度MPI+OpenMP混合并行的方式，相同的并行规模D下面，可以存在不同的进程、线程组合，进程数P、每个进程内的线程数T和并行规模D三者之间的关系为

D＝P×T (3)

当选择不同的进程、线程组合时，节点间仍然保持MPI通信，节点内的进程间则采用 OpenMP共享内存通信，通信性能会发生差异。因此，目标网格所有可能的区域划分方式(M_x，M_y，M_z，threads)四元组可以构成一个集合S作为图1中的输入。S的参数中，前三个参数代表了划分后子网格的三维数值，第四个参数代表混合并行时每个进程内的粗粒度线程数目。

2)为每个区域划分方式(M_x，M_y，M_z，threads)四元组建立通信拓扑。由于采用的是静态的区域分解方法，一旦划分结束，每个子网格块之间的邻居关系便固定了。相邻的子网格块之间由于存在数据依赖，每轮迭代求解时都需要相互通信来完成边界数据的交换。因此子网格块几何序号，进程号和线程号之间存在关系为

ProNum＝blockNum/threads (4)

ThreNum＝blockNum％threads (5)

式(4)、(5)中：

ProNum——进程序号；

ThreNum——线程序号；

blockNum——子网格块几何序号；

threads——进程内线程数目；

因此，只要按照区域分解时定义的规则，知道每个邻居子网格块的几何序号，就能知道当前求解的子网格块与某个邻居之间的通信时采用跨进程通信(MPI)还是进程内通信(OpenMP)。从而建立起某个区域划分方式的通信拓扑。

3)根据MPI通信面的大小对区域划分方式四元组进行过滤分类。由于执行分布式并行计算任务时，MPI消息传输往往是瓶颈，在达到某个阈值之后，MPI的通信模式往往发生切换，从而引起更多的通信延迟。比如，在天河2号上测试得到，当MPI消息大小大于 64KB时，便会引起非阻塞通信开销的大幅上升。因此，按消息大小对集合S过滤分类有利于排除明显不高效的区域划分方式。由于子网格块区域的每个边界上都需要包含一个额外网格作为邻居子网格块区域的数据缓冲区(称之为ghost层)，子网格块区域的实际网格总数为(M_x+2)×(M_y+2)×(M_z+2)，因此三个方向的通信面大小分别为 (M_x+2)×(M_y+2)、(M_x+2)×(M_z+2)和(M_y+2)×(M_z+2)。根据步骤2) 获得的通信拓扑，统计该区域划分方式产生的MPI通信面是否超过当前平台测得的通信阈值。如果所有的区域划分方式都存在MPI通信面超过了通信阈值的情况，则过滤失效，全部区域划分方式都进入步骤4)；否则，只有MPI通信面未超过通信阈值的区域划分方式进入步骤4)。将进入步骤4)的区域划分方式集合记为S₁。

4)对S₁的几何参数(M_x，M_y，M_z)执行面向cache-aware的区域划分策略。如图2所示，不同的通信面在执行数据打包缓冲时，访存连续性不同。X-Y面的地址连续性好，通信面是完全连续的一片地址；Y-Z面的地址连续性差，通信面上每迭代一个元素，都要跨越 (M_x+2)个元素。受限于L2cache的容量和操作数组的数目，所能容纳的初始网格数目也不同。如果一开始元素便容纳在L2cache中，则访问开销假设为1,；如果元素不在L2 cache中，则访问开销假设为ρ，ρ一般取值大于9。由于cahce替换时是以cache line为单位，因此当一片连续的元素不在L2cache中时，这些元素的平均访问开销为

式(6)中：

cache_line_size——cache line的大小；

sizeof(FP)——每个元素所占字节数；

由于元素往往使用浮点数进行存取，每个浮点数是8个字节，而cache line往往是64字节。因此，可以根据L2cahce对子网格块区域的容纳情况分成三种情况进行讨论。

(1)子网格块区域满足

对于X-Y面来说，由于X方向的跨步为1，在一开始的αM_x次迭代时，其访存开销近似为1。之后由于cache line替换的原因，每8次迭代需要对更高层次的存储器件获取新的cache line，因此每8次的访存开销近似为7+ρ。对于X-Z面来说，由于X方向的跨步为1，因此一开始的M_x次迭代时，其访存开销近似为1。之后由于cache line替换的原因，以及X方向上跨步为1，每8次迭代需要对更高层次的存储器件获取新的cache line，因此每8次的访存开销近似为7+ρ。对于Y-Z面来说，由于每次跨步距离为M_x+2，除了一开始的α次的访存开销近似为1，之后每次的访存开销都近似为ρ。因此访存开销如表1所示。

表1访存模型在情况(1)的开销

(2)子网格块区域满足

对于X-Y面来说，由于都能够被L2cache容纳，因此每次的访存开销近似为1。对于X-Z面来说，由于X方向的跨步为1，因此一开始的αM_x次迭代时，其访存开销近似为 1。之后由于cache line替换的原因，每8次迭代需要对更高层次的存储器件获取新的 cacheline，因此每8次的访存开销近似为7+ρ。对于Y-Z面来说，由于每次跨步距离为 M_x+2，因此一开始的次迭代的访存开销近似为1，之后每次的访存开销都为ρ。因此访存开销如表2所示。

表2访存模型在情况(2)的开销

(3)子网格块区域满足

对于三个方向，由于数据都被L2cache容纳，因此其访存开销都近似为1。总的访存开销如表所示。此时，访存开销只与当前区域划分方式产生的三个方向的数值有关。在这种情况下，可能会出现多种划分方式求得的值相同。考虑到实际过程中往往不会出现完全被容纳的情况，因此优先考虑选择X方向数值更大的方式。

表3访存模型在情况(3)的开销

根据这三种情况来对每种区域划分方式的(M_x，M_y，M_z)进行估算，可以获得S₁中每种区域划分方式的性能优劣。

5)根据混合并行的通信过程和步骤2)得到的通信拓扑建立最晚结束时刻模型，度量相同区域划分方式的不同进程、线程组合的性能差异。如图3所示例，混合并行通信可以分为三个部分。第一部分在t₈时刻结束，两种通信面分别按照求解器的既定顺序依次从边界拷贝数据到MPI用户发送缓冲区并启动非阻塞发送和非阻塞接收函数，或者拷贝数据到OpenMP缓冲区并将互斥锁解锁。第二部分在t₈到t₁₂时刻之间，完成线程间共享缓冲区的OpenMP通信，并将数据拷贝到相应通信面的Ghost层。剩下是第三部分，结束MPI 的非阻塞通信，并将接收的数据拷贝到相应通信面的Ghost层。通过研究混合并行通信的执行过程，发现只需要比较图3中t₁₆时刻的数值，就能知道相同区域划分方式在不同进程、线程组合时的性能差异。这个模型实际上计算MPI通信时间和OpenMP通信时间的重叠程度。首先求解第二部分结束的时刻t_fixed如公式(7)：

式(7)中：

n——当前子网格区域的通信面的数目；

m——当前子网格区域的MPI通信面的数目；

k——通信面拷贝一个字节的平均开销；

s_i——第i个通信面总的字节数；

t_{lock_i}——对互斥锁进行一次加锁、解锁的开销。

接着确定某个MPI通信结束的时刻如公式(8)

式(8)中：

n′——该MPI通信面开始激活通道，发送、接收数据之前完成的所有通信面拷贝的数目；

m′——该MPI通信面开始激活通道，发送、接收数据之前总共完成的OpenMP通信面的数目；

T₀——每次MPI通信的固定开销；

s_i——该MPI通信面总的字节数；

t——MPI传输过程中每字节的平均时间开销

然后，比较每一个MPI非阻塞接收检查调用的开始时刻t_{ws_i}与该MPI通信接收真实结束时刻t_{mpi_i}的先后，以时刻更靠后的值作为接收结束时刻t_{we_i}的值如公式(9)

最后，求最晚结束时刻t_final。当求出了最后一个MPI通信面的接收结束时刻t_we，再加上这个通信面的拷贝开销，便得到了最晚结束时刻t_final如公式(10)

t_final＝t_we+k×s_i (10)

通过这个计算过程，可以比较相同的几何划分在不同进程、线程组合下的性能差异，通过缩放因子来表示，如公式(11)

式(11)中：

μ——缩放因子；

t_threads——采用某个线程组合时的最晚结束时刻；

t₁——采用单线程组合时的最晚结束时刻。

6)将步骤5)得到的缩放因子与步骤4)结合起来，便能够计算出开销最小的区域划分方式四元组(M_x，M_y，M_z，threads)，将其作为输出的区域划分方式。图1展示了上述几个步骤的执行流程。

Claims

1.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，包括以下步骤：

步骤1，构建目标网格所有的区域划分方式集合：子网格X方向数值M_x，子网格Y方向数值M_y，子网格Z方向数值M_z，进程内的线程数目threads；定义集合为S；

步骤3，根据MPI通信面的大小对区域划分方式四元组进行过滤分类，将进入下一步的区域划分方式集合记为S₁；

2.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，步骤1具体包括：按照负载均衡的简单六面体区域分解方法来对原始网格进行区域分解，分解成若干个三维大小均为M_x×M_y×M_z的六面体网格；目标网格所有的区域划分方式M_x，M_y，M_z，threads四元组构成一个集合S。

3.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，步骤2具体包括：按照区域分解时定义的规则，子网格块几何序号，进程号和线程号之间存在通过计算得到的映射关系，从而建立起某个区域划分方式的通信拓扑。

4.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，步骤3具体包括：根据步骤2获得的通信拓扑，统计该区域划分方式产生的MPI通信面消息大小是否超过当前平台测得的通信消息阈值；如果所有的区域划分方式都存在MPI通信面消息大小超过了通信消息阈值的情况，则过滤失效，全部区域划分方式都进入步骤4；否则，只有MPI通信面消息大小未超过通信消息阈值的区域划分方式进入步骤4；将进入步骤4的区域划分方式集合记为S₁。

5.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，步骤4具体包括：根据L2 cahce对子网格块区域的容纳情况分成三种情况：

1)子网格块区域满足

X-Y面总的访存开销：

X-Z面总的访存开销：

Y-Z面总的访存开销：

α+(M_y×M_z-α)×ρ

总的访存开销：

2)子网格块区域满足：

X-Y面总的访存开销：

M_x×M_y

X-Z面总的访存开销：

Y-Z面总的访存开销：

αM_y+(M_y×M_z-αM_y)×ρ

总的访存开销：

3)子网格块区域满足：

X-Y面总的访存开销：

M_x×M_y

X-Z面总的访存开销：

M_x×M_z

Y-Z面总的访存开销：

M_y×M_z

总的访存开销：

M_x×M_y+M_x×M_z+M_y×M_z；

根据这三种情况来对每种区域划分方式的M_x，M_y，M_z进行访存开销估算，获得S₁中每种区域划分方式的性能优劣。

6.一种面向粗粒度MPI+OpenMP混合并行的结构网格区域划分方法，其特征在于，步骤5具体包括：

模型反映MPI通信时间和OpenMP共享内存通信时间的重叠程度；首先求解结束OpenMP共享内存通信的时刻t_fixed；接着确定每个MPI通信结束的时刻t_{mpi_i}；然后比较每一个MPI非阻塞接收检查调用的开始时刻t_{ws_t}与该MPI通信接收真实结束时刻t_{mpi_i}的先后，以时刻更靠后的值作为接收结束时刻t_{we_i}；最后，当求出了最后一个MPI通信面的接收结束时刻t_we，再加上这个通信面的拷贝开销，便得到了最晚结束时刻t_final，从而求出缩放因子反映线程数为threads与线程数为1的性能比值。