大规模CFD多区结构网格任务负载平衡算法

唐 波，王勇献

(国防科学技术大学计算机学院，湖南 长沙 410073)

1 引言

计算流体力学CFD(Computational Fluid Dynamics)是一门利用计算机对流体的复杂流动规律进行数值模拟的学科。在CFD数值模拟中, 首先要把计算区域空间离散化为网格, 生成模拟计算过程中所需要的网格文件。网格主要分为结构网格和非结构网格, 本文主要关注结构网格。随着科学技术的发展, CFD数值模拟精度要求越来越高, 网格规模越来越大, 要满足这种需求, 唯一的可行手段就是发展大规模并行计算。

当采用结构网格时, 计算区域被离散为多个大小不等的网格区块, 在模拟计算过程中, 这些区块被分配到不同进程上进行流体流动规律的数值模拟。如何使得数值模拟过程中各个进程上任务负载平衡, 是提高并行CFD程序性能的关键技术之一。对于当前并行CFD数值模拟应用, 任务负载平衡问题已经变得越来越复杂, 随着并行计算规模的加大, 进程间通信开销的大小对并行性能的影响在加大; 网格块数和进程数的不匹配, 网格块的计算量与进程的计算能力不匹配, 使得传统的分割或组合策略不能很好满足负载平衡的需求; 随着近年来异构并行机的出现, 不同类型处理器进程上的计算能力不一样, 不同类型处理器进程间的通信能力不一样, 使得传统的负载平衡算法已不能满足负载平衡要求。因此，有必要研究一种新的负载平衡算法来满足当前的并行CFD任务负载平衡需求。

目前, 国际上已有很多并行CFD任务负载平衡软件, 如MB-split; 也有多种负载平衡算法, 如Streng M等[1]和Ytterstrom A[2]提出的三种经典的负载平衡算法:递归对分法、贪婪负载平衡算法、增量贪婪负载平衡算法, Saule E提出的一种新的启发式算法等[3]。国内也有相应的研究, 如张娟等人[4]提出的基于Block的递归二分对无向图剖分算法, 郑秋亚等人[5]提出的一种类似贪婪负载平衡算法的方法, 李桂波等人[6]提出的基于遗传算法的任务负载平衡算法等。上述软件和算法都只适应特定的任务分配情形, 不能很好地满足当前并行CFD任务负载平衡的需求。

2 问题描述

2.1 计算量和通信开销的度量

计算量与通信开销是衡量大规模并行计算负载平衡的两个重要因素, 它们都与具体的应用问题紧密相关。在典型的多区结构网格CFD应用中, 根据所采用的求解器类型的不同, 计算量跟离散网格单元呈不同的关系。例如, 在经典的显式Runge-Kutta求解器及隐式LUSGS求解器中, 每个网格单元上执行的计算量大体相等, 各网格块的计算量与网格块中网格单元数成正比。本文重点关注这两类求解器, 因此计算量可以简化为网格单元数。

Figure1 Parallel computer architecture图1 图1 并行计算机体系结构

并行计算的通信开销则依赖于高性能计算平台的硬件体系结构及并行计算组织方式。当前的主流并行计算机的体系结构如图1所示, 一台并行计算机内有多个结节, 各个结节的配置一样, 每个结节内可以包含多种类型的处理器, 每个处理器中可以创建多个进程。其中, 不同处理器类型进程的计算能力不一样, 不同进程之间的通信开销随着进程所在处理器类型不一样和是否在同一结节内而不同。如P1与P1、P1与P2、P2与P2之间的通信开销随着它们是否在同一结节内而不同。

因此,如何对进程间的通信开销进行准确建模与度量将变得异常重要。在以往CFD任务负载平衡的研究中, 通常用通信量代替通信开销, 当通信开销与通信量二者不成正比时, 这种代替不能准确地反映进程间通信所占时间。经过测试和总结, 我们用式(1)的分段线性函数来刻画进程间的通信开销与通信量的关系:

C(x)=ax+b,x∈[x1,x2]

(1)

其中，x代表通信量(单位:KB),C(x)为x通信量下的通信开销(单位:ms),a、b是与机器硬件相关的参数。

2.2 负载平衡问题的数学描述

在并行CFD数值模拟计算中, 假定整个CFD流场区域由Nz个网格块组成, 而高性能计算资源包含(图1中)结节Nnode个, 每个结节内创建a1个P1类型进程, 创建a2个P2类型进程, 整个计算的总进程数为Np=Nnode×(a1+a2)。CFD任务负载平衡需要将Nz个网格块分配到Np个进程上, 使得其达到任务负载平衡。下面建立其数学描述。

CFD任务负载平衡问题的数学描述:定义进程i上期望负载计算量为Wi, ave;每给定一种任务分配方式M,我们用Wi(M)表示进程i上的实际计算量,Ci(M)表示进程i与其余进程间的通信总开销。定义进程i上的计算负载平衡因子gi(M)如式(2)所示,它刻画了进程i上实际计算负载与期望计算负载间的相对偏差量。定义任务分配方案M下的最大通信开销f(M)如式(3)所示, 它表示通信负载最大的进程上的通信开销量：

(2)

(3)

本文的CFD任务负载平衡模型旨在寻找一个分配方案使得以下约束优化问题成立：

(1)目标函数：minf(M)；

(2)约束条件：gi(M)=0,i∈{0, 1, 2, …,Np-1}。

其中，M为不同的分配方案。

3 算法设计

在真实的并行CFD任务负载平衡应用中, 我们通常会遇到进程数远大于网格块数或者网格块计算量大于进程所期望负载计算量的情况, 这时我们需要对网格块进行分割，以满足负载平衡的需求; 也会遇到进程数远小于网格块数的情况, 这时我们需要对网格块进行组合映射，以达到各个进程上的负载平衡。为了适应各种应用情形, 本文提出的负载平衡算法分为三步：首先，通过一定原则对原始网格块进行分割, 使其块数大于进程数, 并使得体积最大网格块的计算量不超过进程所期望负载计算量; 然后，对第一步分割后的网格块进行组合分配, 这是一个NP完全问题, 本文采用遗传算法来解决这一组合映射问题; 当经过前两步处理仍存在计算负载不平衡的进程时, 则启动第三步处理, 对所有这些进程上的计算量进行增减微调, 直至所有进程都达到负载平衡需求。下面对这三步进行详细介绍。

3.1 网格块循环分割算法及网格块分割策略

对于给定的进程数目及原始网格分布状态,可计算出每个进程上计算的平均网格量,以此为基准,当原始网格中存在体积超大的网格块时,为了完成负载平衡,必须对这样的网格块进行分割。这时可实施以下的网格块循环分割算法, 使每个网格块的体积大小不超过单个进程的平均网格量,从而可使用第二步完成网格块向进程上的组合映射。网格块循环分割算法的核心是维护一个网格块列表和一个进程列表,通过反复迭代,逐次完成网格块的分割。

3.1.1 网格块循环分割算法

输入原始网格块, 找出未标记体积最大网格块, 记该网格块为Zunsign, 其计算量为Nmax; 计算进程的集合, 找出未标记计算能力最大的进程, 如果所有进程都已标记则选出计算能力最小的进程, 记该进程为Pexpect, 其期望负载计算量记为Wexpect。算法如下：

1. while(存在Zunsign)

2. if(Nmax

3. 标记Zunsign和Pexpect;

4. else{

5. 分割Zunsign, 标记计算量为Wexpect的子网格块和Pexpect;

6. 其余子网格块添加到原始网格未标记网格块的序列中}

7. end if

8. end while

输出适合向进程上映射的网格块集合, 其表现在以下两个方面:网格块数大于进程数; 网格块的计算量不超过进程所期望负载的计算量。

3.1.2 网格块分割策略

CFD数值模拟中, 网格块网格单元层数太少会对模拟计算的精度产生很大影响, 甚至导致计算结果不收敛。在并行CFD计算过程中, 对接边界面积与通信量成正比。因此，在网格块分割过程中要避免出现过“薄”的网格块, 同时使得分割后子网格块的表面积较小。为此我们提出了网格块的多维分割方法来满足这种需求。

本文多维分割方法分三种情况对网格块进行不同维数的分割:

(1)如果该网格块的形状接近正方体, 则对其进行三维分割, 分割过程如图2所示。图2中标记“r”的网格块为分割所需的网格块, 各子网格块的大小不一样, 为了方便画图, 此处子网格块大小取为一样。图3和图4作相同解释。

Figure 2 Three-dimensional split example图2 三维分割示例

(2)如果该网格块的形状为“扁平”状, 则对其进行二维分割, 分割过程如图3所示。

Figure 3 Two-dimensional split example图3 二维分割示例

(3)如果该网格块的形状为“长条”状, 则对其进行一维分割, 分割过程如图4所示。

Figure 4 One-dimensional split example图4 一维分割示例

上述多维分割算法, 本文称为非智能多维分割算法。该算法在实际分割过程中, 容易产生过小的网格块, 如所需子网格块的大小和需分割父网格块的大小相差不大, 这时对父网格块进行三维分割, 将会产生体积很小的子网格块。因此, 本文在非智能多维分割算法的基础上, 提出了智能多维分割算法, 其过程如下:

输入:待分割的最大网格块, 其计算量为Nmax; 当前进程期望负载的计算量, 记为Wexpect。

输出:分割所需的子网格块和剩余的子网格块。

1. if(Wexpect≤1/8*Nmax)

2. if(是否适合三维分割)

3. 三维分割;

4. else if(是否适合二维分割)

5. 二维分割;

6. else 一维分割;

7. end if

8. else if(Wexpect>1/8*Nmax&&Wexpect≤1/4*Nmax)

9. if(是否适合二维分割)

10. 二维分割;

室间隔缺损是儿科常见先天性心脏病，约20%，其中缺损口径小、分流量少等症状较轻者一般临床症状并不典型，症状严重者会出现呼吸窘迫和左心衰竭等症状，由于左心室与右心室存在压差，左向右分流时可导致肺血多引发充血性心力衰竭、肺动脉高压等症状，给患儿的生命安全带来了极大的威胁[1-2]。研究表明，在室间隔缺损并发症出现之前予以针对性的措施进行干预，是可以痊愈的，其中治疗方式选择至关重要[3-4]。本研究探讨介入治疗、外科手术在室间隔缺损患儿的应用价值，现报道如下。

11. else 一维分割;

12. end if

13. else if(Wexpect>1/4*Nmax)

14. 一维分割;

15. end if

两种多维分割算法, 在相同的条件下分割网格块的效果将在实验部分给出。

3.2 网格块向进程上的分配算法

网格块向进程上的分配问题是一个NP完全问题, 随着CFD数值模拟规模的扩大, 问题解的搜索空间急剧增加, 普通的搜索算法已不能满足CFD任务分配的需求, 本文采用遗传算法来解决该问题。

3.2.1 遗传算法的基本流程

Figure 5 Flow chart of genetic algorithm图5 遗传算法流程图

由于CFD任务分配问题的特殊性, 在应用遗传算法解决该问题的过程中, 为了减小遗传算法的搜索空间、提高遗传算法的搜索能力、加快遗传算法的收敛速度, 我们对遗传算法的以下操作做了相应改进, 以适应CFD任务分配问题。

3.2.2 染色体编码和种群初始化

在遗传算法中, 一个染色体对应问题的一个解。常规遗传算法的染色体编码一般采用二进制编码, 对于本文的任务分配问题是不适用的, 本文染色体x的编码形式如式(4)所示:

(4)

本文用随机算法按照式(4)的编码方式生成遗传算法的初始种群。随机生成的编码,可能产生部分进程分配的网格块为空的个体,称为无效个体。对于这种情况,首先找出网格块数为0的进程,然后从网格块数最多的进程中随机分出一个网格块给此空进程,如果存在多个空进程,循环执行上述操作,直至空进程数为0。

3.2.3 选择操作、交叉操作、变异操作

本文中选择操作采用经典的轮盘赌方法和最佳个体保存方法;变异操作采用传统的单点变异操作,采用3.2.2节中提到的方法处理变异操作过程中产生的无效个体。

为了避免因传统交叉操作产生的无效个体,本文引用并改进了由Goldberg D E和Lingle R[7]于1985年提出的PMX操作,使其能够满足该问题的需求。例如,对于两个父个体,随机产生两个位串交叉点“|”,位串交叉点之间的区域为交叉区域,如表1所示。首先交换父个体的两个交叉区域,得到表1中的中间体;为了避免出现空进程,进行如下操作,中间体的交叉区域映射关系有:

2→1, 1→2, 4→3

根据映射关系,对M1中间体交叉区域外的部分进行映射交换。如第一对映射关系2→1。找出M1中间体交叉区域外的第一个2所在位置,使其变换为1,然后依次根据映射关系进行类似变换,最终得到M1的子个体如表1所示。同样对M2中间体进行映射交换,可得最终的子个体。

Table 1 Crossover operation表1 交叉操作

3.2.4 适应度函数设计

我们采用带惩罚函数[8]的适应度函数设计, 如式(5)所示:

(5)

其中，M为分配方案(个体)；F(M)为个体M下的适应度值;f(M)刻画了个体M下进程间的最大通信开销;g(M)=max(|gi(M)|), 刻画了个体M下进程上负载平衡因子绝对值的最大值;φ为惩罚函数,r1为权重系数,r2为惩罚系数。

对于不考虑通信开销的情形, 可令r1为0。

惩罚函数有很多确定方法, 为了扩大约束条件g(M)对适应度的影响, 我们设定:

φ(g(M))=g2(M)

(6)

对于惩罚系数r2, 取个体M下的最大通信开销f(M)的值, 此时适应度函数为:

(7)

其中r1为可调节的参数。当r1较大时, 为弱约束条件;当r1较小时, 为强约束条件。

3.3 进程上计算量二次优化算法

经过3.2节处理后, 原始网格块在空间上分为了Np个区域, 每个区域对应一个进程。如图6所示，实心黑圈代表网格块, 连线代表网格块之间的相邻关系, 网格块数目为11; 虚线划分的区域与进程对应, 进程数为3。由于3.2节中的最终网格块分配方案是由带惩罚函数的适应度函数选出的最优个体, 不一定能很好地满足2.2节中提出的约束条件, 也就是各进程上的计算量不一定能达到负载平衡要求。对于没有达到负载平衡的进程我们需要对其计算量进行二次优化, 使其达到计算量的负载平衡。

Figure 6 Grids partition example图6 网格区域划分示例

该算法主要分为三个步骤:(1) 当前分配方案, 若需要二次优化, 选出需要二次优化的进程P; (2) 从进程集合中找出合适的进程Psuit, 使其与P共同完成二次优化; (3) 从进程P中分离出所需的网格块Zsuit分配给Psuit。下面分别描述上面三个步骤的细节:

第1步:在当前的分配方案中, 若进程最大的负载平衡因子大于所期望的值(此处取为0.05) , 则该进程P需要二次优化。

第2步:

(1)列出在空间上与P相邻的进程链表Plist,Plist中负载平衡因子最小的进程为Plmin;

(2)if(Plmin的负载平衡因子< -0.05)

Plmin→Psuit;

(3)else 找出进程集合中负载平衡因子最小的进程Pmin,Pmin→Psuit;

(4)end if

第3步:

(1)定义对象为网格块的链表Zlist;

(2)if(P与Psuit相邻)

(3) 列出P中与Psuit有相邻关系的网格块, 添加给Zlist;

(4)else列出P中的网格块, 添加给Zlist;

(5)end if

(6)计算Psuit所需的计算量Nrequire;

(7)if(Zlist中是否存在网格块Zsuit的计算量N, 使得Nrequire*(1-0.05)

(8)Zsuit分配给Psuit, 从P中删除Zsuit;

(9)else 使用分割策略分割Zlist中计算量最大的网格块, 分割出所需的网格块Zsuit;

(10)Zsuit分配给Psuit, 分割剩余的子网格块分配给P, 从P中删除被分割的网格块;

(11)end if

4 数值实验结果

4.1 通信模型中参数的确定

本文的数值实验平台为国家超级计算长沙中心的高性能计算系统, 每个结节由两个6核的Intel Xeon X5670的CPU组成, 主频为2.93 GHz, 双精度浮点峰值性能是70.32 GFlops。

在上述平台上, 经过简单的实验测试与数据拟合, 结果表明，通信开销与通信量之间均呈分段线性函数关系, 但结节内与结节间两种情形有区别, 结果如表2所示。

Table 2 Test results of communication overhead model表2 通信开销模型测试结果

4.2 CFD算例负载平衡性能测试

4.2.1 同构平台的测试

在4.1节通信开销模型测试的基础上, 不同的CFD算例在不同的策略下进行负载平衡性能测试。根据网格块数与进程数目的相对多少，将算例分为两种情形, 下文分别称为“情形A”和“情形B”。共测试了四种负载平衡策略, 分别用a、b、c、d表示(如图7和8所示), 其中a为文献[1]提到的贪婪负载平衡算法,b、c、d为本文算法。a采用一维分割策略,b采用智能多维分割策略,c采用非智能多维分割策略,d采用一维分割策略。我们将以文献中的策略a作为性能比较的基准。

情形A:网格块数较少情形的测试结果(原始网格包含4个网格块, 12亿网格单元)。

Figure 7 Test results of load balancing in less zone number图7 网格块数较少情形的负载平衡测试结果

测试结果显示, 在此种情况下, 本文算法b、c、d与传统算法a在性能比较上具有很大的优势。随着进程数的增加, 算法a的最大负载平衡因子绝对值急剧增加, 显现出了各进程上计算任务的不平衡性, 这远离了负载平衡的约束条件, 因此本文算法优于传统算法。此处a方法64进程后的情况不予测试。对于b与c和d之间的比较, 就进程间最大通信开销、对接边界面数目、对接边界面面积三方面而言,b具有明显优势。实验结果表明了b算法能够使得分割负载后网格块之间的对接边界面数较少, 对接边界面面积最少, 从而在各进程计算任务负载平衡的基础上使得进程间的最大通信开销最小。综上所述, 在网格块数较少的情况下,b为综合性能最好的方法。

情形B:网格块数较多情形的测试结果(原始网格包含265个网格块, 15 000 000网格单元)。

Figure 8 Test results of load balancing in more zone number图8 网格块数较多情形的负载平衡测试结果

测试结果显示, 在此种情况下, 本文算法b、c、d与传统算法a在性能比较上具有很大的优势。在进程数较少的时候, 都能满足最大负载平衡因子绝对值小于所期望的值(本文为0.05), 但在进程间通信开销方面本文算法具有明显优势。随着进程数的增加, 算法a的最大负载平衡因子绝对值急剧增加, 显现出了各进程上计算任务的不平衡性, 这远离了负载平衡的约束条件, 综合分析本文算法优于传统算法。当进程数小于网格块数或相差不大时, 进程间最大通信开销、对接边界面数目、对接边界面面积三方面而言,b较c具有明显优势，b与d相当。随着进程数的继续增加, 当进程数远远大于网格块数时, 情形变换为A, 如图7所示, 此时在进程间最大通信开销、对接边界面数目、对接边界面面积三方面的比较中,b具有明显优势。实验结果表明，在网格块数较多的情况下,b能使空间上具有相邻关系的网格块分配到同一进程上, 从而减少各进程间的通信开销。在此种情况下，算法b能使分割后网格中的对接边界面数较少, 对接边界面面积最小, 从而在各进程上计算任务负载平衡的基础上使得进程间的最大通信开销最小。综上所述, 在网格块数较多的情况下, 算法b是综合性能最好的方法。

4.2.2 异构平台的测试

如图1所示并行计算机, 假设一个结节包含第一类处理器两个, 第二类处理器一个, 每个处理器中创建两个进程。假设第二类处理器的计算能力是第一类的两倍, 本文异构平台下只作理论测试, 所以进程间的通信开销简化为4.1节测试的通信开销模型。测试结果如图9所示,a′、b′、c′为本文算法，a′采用智能多维分割策略,b′采用非智能多维分割策略,c′采用一维分割策略(原始网格包含4个网格块, 12亿网格单元)。结果显示, 本文算法a′、b′、c′同样适合异构平台, 对于a′与b′和c′之间的比较, 就最大负载平衡因子绝对值来说, 三种方法都能满足要求；就进程间最大通信开销、对接边界面数目、对接边界面面积三方面而言,a′具有明显优势。综上所述, 本文算法同样适合异构平台, 采用智能多维分割方法的本文算法a′为综合性能最好的方法。

Figure 9 Test results of load balancing on heterogeneous platform图9 异构平台上的负载平衡测试结果

4.3 CFD实例负载平衡测试

在4.1节中提到的实验平台下进行CFD实例负载平衡测试, 其中网格包含八个等体积的网格块, 800万网格单元, 采用消息传递并行编程模型分别在20、40、60个CPU进程上进行并行模拟。表3列出了4.2.1节中四种方法的实际负载平衡性能(时间为CFD并行计算迭代20步所用的时间)。

对所列四种负载平衡指标对CFD并行测试结果进行分析, 结果表明:随着进程数的增加, 方法b在并行CFD任务负载平衡方面具有最优的稳定效果, 并且具有很好的可扩展性(例如，进程上的最大计算时间随着进程数的增加线性减少), 而且进程间最大通信开销也是四种方法中最小的。

综合上述所有测试, 智能多维分割算法能够使负载平衡后的最大负载平衡因子满足要求, 使得进程间的最大通信开销最小, 对接边界面数目较小, 对接边界面面积最小。

Table 3 Test result of load balancing of CFD application example表3 CFD实例负载平衡测试结果

5 结束语

本文提出的大规模CFD多区结构网格任务负载平衡算法, 集网格块的分割、网格块的组合、进程上计算量的优化于一体, 解决了传统负载平衡方法适应性不强的问题。该算法既适应同构计算平台, 也适应新型异构平台; 既适应网格块较少的情形, 也适应网格块较多的情形。该算法对通信开销方面做了很好的度量。实验结果表明, 在适当的分割策略下本文算法能够取得很好的并行CFD任务负载平衡性能, 减少并行CFD程序的运行时间。随着并行规模的扩大, 本文所用遗传算法的搜索空间会逐渐增加, 负载平衡处理时间会渐渐加大, 后续工作将研究如何对大规模负载平衡问题中的遗传算法进行加速。

[1] Streng M. Load balancing for computational fluid dynamics calculationsp[M]∥High Performance Computing in Fluid Dynamics Netherlands,Netherlands:Springer, 1996.

[2] Ytterström A. A tool for partitioning structured multiblock meshes for parallel computational mechanics[J]. International Journal of High Performance Computing Applications, 1997, 11(4):336-343.

[3] Saule E, Bas E Ö, Catalyürek I V. Load-balancing spatially located computations using rectangular partitions[J]. Journal of Parallel and Distributed Computing,2012,70(10):1201-1214.

[4] Zhang Juan,Lu Lin-sheng.Automatic partition algorithm based on multi-region and multi-code problem[J]. Computer Engineering, 2010, 36(9):73-76.(in Chinese)

[5] Zhang Qiu-ya, Liu San-yang, Zuo Da-hai. CFD parallel computing and load balancing research using multi-block structured grids[J]. Chinese Journal of Engineering Mathematics, 2010, 27(2):219-224.(in Chinese)

[6] Li Gui-bo, Yang Guo-wei. Study on parallel computation and load balance strategy based on multiblock structured grid[J]. Journal of Astronautics, 2011, 32(6):1224-1230.(in Chinese)

[7] Goldberg D E, Lingle R.Alleles, Loci, and the travelling salesman problem[C]∥Proc of International Conference on Genetic Algorithms and Their Applications, 1985:154-159.

[8] Chen Guo-liang, Wang Xu-fa, Zhuang Zhen-quan, et al. Genetic algorithm and its application[M]. Beijing:People’s Posts and Telecommunications Publishing House, 2001.(in Chinese)

附中文参考文献：

[4] 张娟, 陆林生.基于多区域多代码问题的自动分块算法[J].计算机工程, 2010, 36(9):73-76.

[5] 郑秋亚, 刘三阳, 左大海. 多块结构化网格CFD并行计算和负载平衡研究[J].工程数学学报, 2010, 27(2):219-224.

[6] 李桂波, 杨国伟. 基于多块结构网格的并行计算及负载平衡研究[J]. 宇航学报, 2011, 32(6):1224-1230.

[8] 陈国良, 王煦法, 庄镇泉, 等. 遗传算法及其应用[M]. 北京:人民邮电出版社, 2001.