CN110633149A - 均衡非结构网格单元计算量的并行负载均衡方法 - Google Patents

Abstract

本发明公开了一种均衡非结构网格单元计算量的并行负载均衡方法，目的是提高非结构网格应用的并行计算负载均衡性能和并行计算速度。技术方案是先生成二维非结构网格，然后使用非结构网格结构化处理方法对非结构网格相关文件进行局部化，生成非结构网格局部方向数组和关键元素数组，接着基于非结构网格局部方向数组，采用递归选择模板的方式进行网格单元之间计算量层面的负载均衡，最后采用网格剖分方法进行处理器之间网格量层面的负载均衡，采用子网格区域分配的方式将计算量分配给处理器。本发明在处理器之间网格量负载均衡的基础上实现了网格单元之间计算量的均衡，使处理器的计算量与其计算能力相匹配，提高了负载均衡性能和并行计算速度。

Description

均衡非结构网格单元计算量的并行负载均衡方法

技术领域

本发明属于高性能计算领域负载均衡方法，具体涉及一种均衡非结构网格单元计算量的并行负载均衡方法。

背景技术

工程计算可解决无法进行实验研究或无法进行复杂理论分析的问题，现已广泛应用于航空航天、汽车船舶制造、桥梁建筑设计、气象预报分析、环境工程等众多工程应用领域。一般来说，工程计算的主要过程是对网格上的物理量进行迭代计算。研究对象的复杂度增加时，网格量与计算量越来越大，为了减少研发的时间成本，常使用高性能计算机进行大规模并行计算。

负载均衡是影响并行计算效率的关键因素之一。负载均衡是指将计算量按计算能力分配到各个处理器上并行执行，使得所有处理器的计算能力得以充分发挥。计算量均衡时，所有处理器同时结束计算，整体计算时间最短，并行效率和计算效率最优。计算量若不均衡，部分处理器将空闲等待，造成计算资源浪费，同时整体计算时间也将延长，导致工程应用的研发成本增加。因此，采用负载均衡方法为各处理器分配均衡的计算量是大规模并行计算中不可或缺的环节。

负载均衡的方法与采用的网格类型有关。工程应用中，网格主要分为结构网格和非结构网格两类。结构网格沿着笛卡尔坐标系的坐标轴方向规则排列，网格单元的形状为四边形或六面体，易获得高阶精度和高计算效率，但其网格生成等处理困难且耗时，需要大量人工干预。非结构网格没有排列和形状的限制，可采用商业软件自动实现网格生成、网格变形、网格剖分等处理，更适用于外形复杂的工程应用。

目前，非结构网格工程的负载均衡方法(称为背景技术1)是通过网格剖分实现处理器之间网格量层面的负载均衡，即根据处理器信息将非结构网格进行区域分解与分配，使得处理器获得与之计算能力相匹配的网格量。背景技术1的流程图如图1所示，包括以下步骤：

第一步，生成二维非结构网格。

使用专业软件将非结构网格工程应用需要计算的二维空间进行网格化处理，生成二维非结构网格并输出二维非结构网格相关文件。常用的专业软件包括ANSYS公司的Fluent软件[http://www.ansys.com，版本18.0以上]等商业软件、美国国家航空航天局的FUN3D软件[Biedron R等，FUN3D Manual:13.2，2017，版本13.0以上]等In-house软件、以及OpenFOAM[https://www.openfoam.com.版本3.0以上]、SU2[https://su2code.github.io，版本5.0以上]等开源软件。非结构网格相关文件包括网格单元文件、网格点文件、网格边文件、网格单元的邻接单元文件、网格点的邻接单元文件、网格边的邻接单元文件等。

第二步，根据处理器信息和非结构网格相关文件，采用网格剖分方法进行处理器之间网格量层面的负载均衡。

2.1根据处理器总数及各处理器计算能力，初始化处理器计算能力数组Ps，

其中N_p为处理器总数，P_i为第i个处理器的计算能力；

2.2根据第一步生成的非结构网格相关文件初始化网格相关数据结构。具体步骤如下：

2.2.1令网格单元文件为网格单元数组cs，

其中N_c为网格单元总数，C_kc为第kc个网格单元。网格单元采用网格点编号数组表示，C_kc＝{kp₁,...,kp_w,...,kp_W}表示第kc个网格单元由第kp₁,...,kp_w,...,kp_W个网格点构成，是一个W边形网格单元，W为正整数。如C₁＝{1,4,6}表示第1个网格单元由第1、4、6个网格点构成，是一个三角形网格单元；

2.2.2令网格点文件为网格点数组pts，

其中N_pt为网格点总数，Pt_kp为第kp个网格点。网格点采用笛卡尔坐标数组表示，Pt_kp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数，如Pt₁＝{1.0,2.0}表示二维网格中第1个网格点的坐标为(1.0,2.0)；

2.2.3令网格边文件为网格边数组eds，

其中N_e为网格边总数，E_ke为第ke条网格边。网格边采用顶点(即网格边两端的网格点)编号数组表示，E_ke＝{kp₁,kp₂}表示第ke条网格边的两个顶点分别为第kp₁,kp₂个网格点，如E₁＝{1,4}表示第1条网格边的两个顶点分别为第1、4个网格点；

2.3根据第一步生成的非结构网格相关文件初始化网格邻接关系数据结构。

具体步骤如下：

2.3.1令网格单元的邻接单元文件为数组cCells，

其中CC_kc表示网格单元C_kc的邻接网格单元数组。如CC₁＝{2,3,5,7,11}表示与网格单元C₁相邻接的网格单元为C₂,C₃,C₅,C₇,C₁₁；

2.3.2令网格点的邻接单元文件为数组pCells，

其中PC_kp表示网格点Pt_kp的邻接网格单元数组。如PC₁＝{1,2,4,5,7,8}表示与网格点Pt₁相邻接的单元为C₁,C₂,C₄,C₅,C₇,C₈；

2.3.3令网格边的邻接单元文件为数组eCells，

其中EC_ke表示网格边E_ke的邻接网格单元数组。如EC₁＝{1,2}表示与边E₁相邻接的网格单元为C₁,C₂；

2.4根据处理器计算能力数组Ps及N_p,N_c，计算各处理器的最优网格量。处理器的最优网格量与其计算能力成正比，以第i个处理器为例，其最优网格量M_i可表示为：

其中1≤k≤N_p且k为整数；

2.5使用商业软件(如Fluent、OpenFOAM等)剖分非结构网格。在商业软件中选定剖分算法，按照第2.4步计算得到的最优网格量M_i进行剖分，将非结构网格单元剖分成N_p个子网格区域，第i个子网格区域的网格量等于最优网格量M_i，其中1≤i≤N_p。剖分后生成子网格区域的非结构网格相关文件，包含相邻子网格区域的边界信息；

2.6采用子网格区域分配的方式，将计算量分配给处理器。将第i个子网格区域的网格量分配给第i个处理器，其中1≤i≤N_p，并生成子网格区域对应的非结构网格相关文件。网格量分配后，第i个处理器将获得M_i个网格单元，即M_i个网格单元上的计算量分配给了第i个处理器。

上述方法实现了网格量在处理器间的分配与均衡，但其并不一定能获得较好的负载均衡性能。因为单个非结构网格单元的计算量不一定相同，背景技术1为处理器分配了均衡的网格单元个数，并不等价于为处理器分配了均衡的计算量。非结构网格工程应用中，网格单元C_kc的计算量正比于C_kc模板单元的个数，C_kc更新及计算时所需使用的网格单元称为C_kc的模板单元。一般情况下，为保证非结构网格应用的计算精度和鲁棒性，常采用共点单元作为模板单元，即采用C_kc的邻接网格单元数组CC_kc中的网格单元作为C_kc的模板单元。由于非结构网格排列无序且不规则，各网格单元的共点单元个数常不相等，导致各网格单元的计算量不均衡。尤其对于复杂的工程应用，网格单元之间计算量的差异更大。此时背景技术1采用网格量分配的方式，不能实现计算量在处理器间的均衡分配，影响了并行负载均衡性能。

因此在处理器之间网格量均衡的基础上，进一步研究网格单元计算量均衡的负载分配方法，对提高非结构网格工程的并行计算负载均衡性能具有重要意义。

发明内容

本发明要解决的技术问题是：针对二维工程应用中使用最广泛的三角形、四边形及三角形-四边形混合非结构网格，在现有网格量均衡分配的基础上，提出一种均衡非结构网格单元计算量的并行负载均衡方法。通过为非结构网格单元重新选择模板单元，使得每一个非结构网格单元的计算量均相等，进而使得处理器在获得均衡网格量的同时，也获得均衡的计算量，真正提高非结构网格应用的并行计算负载均衡性能和并行计算速度。

本发明技术方案如下：

第一步，生成二维非结构网格并初始化非结构网格相关数据结构。

1.1使用专业软件将非结构网格工程应用需要计算的二维空间进行网格化处理，生成非结构网格并输出非结构网格相关文件。此步与技术背景1第一步相同。

1.2根据非结构网格相关文件初始化相关数据结构，方法是：

1.2.1令网格单元文件为网格单元数组cs，其中N_c为网格单元总数，C_kc为第kc个网格单元；网格单元采用网格点编号数组表示，C_kc＝{kp₁,...,kp_w,...,kp_W}表示第kc个网格单元由第kp₁,...,kp_w,...,kp_W个网格点构成，是一个W边形网格单元，W为正整数；

1.2.2令网格点文件为网格点数组pts，

其中N_pt为网格点总数，Pt_kp为第kp个网格点；网格点采用笛卡尔坐标数组表示，Pt_kp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数；

1.2.3令网格单元的邻接单元文件为数组cCells，其中CC_kc表示网格单元C_kc的邻接网格单元编号组成的数组；

1.2.4令网格点的邻接单元文件为数组pCells，

其中PC_kp表示网格点Pt_kp的邻接网格单元编号组成的数组；

第二步，使用非结构网格结构化处理方法对非结构网格相关文件进行局部化，生成非结构网格局部方向数组和关键元素数组，方法是：

2.1根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和非结构网格局部方向数组，方法是：

2.1.1定义关键元素数组

其中El_kc为C_kc的关键元素编号组成的数组。本发明中C_kc的关键元素是用于确定局部方向的C_kc的网格边或网格点。具体而言，三角形网格单元的关键元素为3条网格边和1个网格点，如当C₁为三角形网格且El₁＝{2,4,7,9}时，表示C₁的4个关键元素分别为C₁的网格边E₂、E₄、E₇以及C₁的网格点Pt₉；四边形网格单元的关键元素为4条网格边，如当C₁为四边形网格且El₁＝{2,4,7,9}时，表示C₁的4个关键元素分别为构成C₁的网格边E₂、E₄、E₇、E₉。初始化

均为{0,0,0,0}，从而与有意义的单元编号区别以避免歧义；

2.1.2根据

定义非结构网格局部方向数组locDir，

其中LD_kc表示网格单元C_kc的局部方向向量组成的数组。LD_kc包含4个数据项，分别表示C_kc的4个局部方向向量，由El_kc确定：

·LD_kc[1]为El_kc[3]边中点指向El_kc[1]边中点构成的向量；

·若C_kc为四边形，LD_kc[2]为El_kc[4]边中点指向El_kc[2]边中点构成的向量；若C_kc为三角形，LD_kc[2]为El_kc[4]网格点指向El_kc[2]边中点构成的向量；

·LD_kc[3]为El_kc[1]边中点指向El_kc[3]边中点构成的向量；

·若C_kc为四边形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]边中点构成的向量；若C_kc为三角形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]网格点构成的向量。

由上述定义易知LD_kc[1]＝-LD_kc[3],LD_kc[2]＝-LD_kc[4],1≤kc≤N_c。由于局部方向需要在关键元素确定后才能确定，故

的每一项都初始化为零向量；

2.2采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法是：

2.2.1定义并初始化阵面推进相关数据结构，方法是：

2.2.1.1定义并初始化阵面数组fronts。fronts为动态数组，数据项个数记为N_f，每一项表示已进行访问的网格点的编号。fronts初始化为位于壁面边界(壁面边界指物体表面，如火箭表面等)的网格点的编号数组，N_f初始化为位于壁面边界的网格点的总数；

2.2.1.2定义并初始化下一阵面数组nextFronts。nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为N_nf。nextFronts初始化为空集，N_nf初始化为0；

2.2.1.3初始化阵面层数组frontsLevel，

其中Fl_kp表示网格点Pt_kp的阵面层值。若kp在数组fronts中，初始化Fl_kp＝0，否则令Fl_kp＝-1；

2.2.1.4定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；

2.2.1.5初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝N_f、令最终尾索引finalTail＝N_f；

2.2.2以head为fronts索引，获取网格点的编号front[head]，为方便表述，记front[head]＝j,1≤j≤N_pt；

2.2.3以j为网格点的邻接单元数组pCells的索引，

获取所有与网格点Pt_j相邻的网格单元编号组成的数组PC_j，记PC_j的数据项总数为PCnum；

2.2.4初始化网格单元索引i＝1；

2.2.5以i为PC_j的索引，获取第i个与Pt_j相邻的网格单元编号PC_j[i]，为方便表述，令o＝PC_j[i],1≤o≤N_c。若El_o[1]＝0，则C_o的关键元素和局部方向还没有确定，转2.2.6；若El_o[1]≠0，转2.2.9；

2.2.6根据网格单元数组cs判断C_o的形状，若C_o包含三个数据项，说明C_o为三角形，转2.2.7；否则C_o包含四个数据项，转2.2.8；

2.2.7确定三角形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，具体流程如下：

2.2.7.1根据网格单元数组cs获得C_o三个顶点的编号C_o[1]、C_o[2]、C_o[3]；

2.2.7.2以C_o[1]、C_o[2]、C_o[3]为阵面层数组frontsLevel索引，获得对应阵面层值

2.2.7.3对比

按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]重排序，重排序后的C_o[1]、C_o[2]、C_o[3]记为kp₁、kp₂、kp₃，即

2.2.7.4若

转2.2.7.5；否则

转2.2.7.6；

2.2.7.5更新El_o[4]＝kp₃，El_o[2]为kp₁、kp₂组成的网格边的编号，El_o[1]为kp₁、kp₃组成的网格边的编号，El_o[3]为kp₂、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.6更新El_o[4]＝kp₁，El_o[2]为kp₂、kp₃组成的网格边的编号，El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₁、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.7根据El_o，按照第2.1.2步中局部方向的定义，计算得到LD_o，即：

·LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

·LD_o[2]为El_o[4]网格点指向El_o[2]边中点构成的向量；

·LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

·LD_o[4]为El_o[2]边中点指向El_o[4]网格点构成的向量。

2.2.7.8更新nextFronts，具体步骤为：

2.2.7.8.1若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.2若

且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.3若

且kp₃不在nextFronts中，将kp₃加入nextFronts，N_nf＝N_nf+1，转2.2.9；

2.2.8确定四边形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，具体流程如下：

2.2.8.1根据网格单元数组cs获得C_o四个顶点的编号C_o[1]、C_o[2]、C_o[3]、C_o[4]；

2.2.8.2以C_o[1]、C_o[2]、C_o[3]、C_o[4]为frontsLevel索引，获得对应阵面层值

2.2.8.3对比

按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]、C_o[4]重排序，将重排序后的C_o[1]、C_o[2]、C_o[3]、C_o[4]记为kp₁、kp₂、kp₃、kp₄，即

2.2.8.4更新El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₃、kp₄组成的网格边的编号。若kp₁、kp₃构成C_o的网格边，El_o[2]更新为kp₁、kp₃组成的网格边的编号，El_o[4]更新为kp₂、kp₄组成的网格边的编号；否则若kp₁、kp₄构成C_o的网格边，El_o[2]更新为kp₁、kp₄组成的网格边的编号，El_o[4]更新为kp₂、kp₃组成的网格边的编号；

2.2.8.5根据El_o，按照第2.1.2步中局部方向的定义，计算得到LD_o，即

·LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

·LD_o[2]为El_o[4]边中点指向El_o[2]边中点构成的向量；

·LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

·LD_o[4]为El_o[2]边中点指向El_o[4]边中点构成的向量。

2.2.8.6更新nextFronts，具体步骤为：

2.2.8.6.1若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.8.6.2若

且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1，转2.2.9；

2.2.9令i＝i+1，若i＜PCnum+1，转2.2.5；否则转2.2.10。

2.2.10令head＝head+1，若head＜tail+1，转2.2.2；否则转2.2.11；

2.2.11根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：

2.2.11.1将nextFronts的数据项加入到fronts中，更新N_f＝N_f+N_nf；

2.2.11.2更新frontsLevel中数据项Fl_kp＝n+1，其中kp为nextFronts的数据项；

2.2.11.3更新n＝n+1；

2.2.11.4更新nextFronts为空集，令N_nf＝0；

2.2.11.5更新tail＝finalTail＝N_f；

2.2.12若head＜finalTail+1，转2.2.2；否则所有网格单元的关键元素和局部方向都已确定，得到关键元素数组keyEl和局部方向数组locDir，转第三步；

第三步，基于非结构网格局部方向数组，采用递归选择模板的方式进行网格单元之间计算量层面的负载均衡。

3.1根据第一步生成的非结构网格相关文件初始化网格边数组。令网格边文件为网格边数组eds，

其中N_e为网格边总数，E_ke为第ke条网格边。网格边采用顶点(即网格边两端的网格点)编号数组表示，E_ke＝{kp₁,kp₂}表示第ke条网格边的两个顶点分别为第kp₁,kp₂个网格点，如E₁＝{1,4}表示第1条网格边的两个顶点分别为第1、4个网格点；

3.2根据第一步生成的非结构网格相关文件初始化网格边的邻接单元数组。令网格边的邻接单元文件为数组eCells，

其中EC_ke表示网格边E_ke的邻接网格单元数组。如EC₁＝{1,2}表示与边E₁相邻接的网格单元为C₁,C₂；

3.3定义并初始化模板单元个数N。N为常数，表示需要为每个网格单元选择的模板单元个数，也为4个方向的总模板单元个数，N与网格单元的计算量成正比。N的取值与具体应用相关，N为正整数，一般为4的倍数；

3.4定义并初始化模板单元数组stencil，

其中S_kc是C_kc的模板单元编号组成的数组。由于本发明将为每个网格单元选择N个模板单元，故S_kc有N个数据项，S_kc的每个数据项为C_kc的模板单元编号。S_kc的每个数据项都初始化为0，以避免与有意义的单元编号重复而引起歧义。

3.5初始化网格单元递归的单元索引kc＝1；

3.6以kc为局部方向数组locDir和关键元素数组keyEl的索引，分别获取C_kc的4个局部方向向量：LD_kc[1]、LD_kc[2]、LD_kc[3]、LD_kc[4]，以及4个局部方向对应的4个关键元素：El_kc[1]、El_kc[2]、El_kc[3]、El_kc[4]；

3.7初始化局部方向索引j＝1；

3.8若C_kc为三角形网格且j＝4，LD_kc[j]对应的关键元素El_kc[4]为网格点，转3.9进行按点模板选择；若C_kc不为三角形网格或j≠4，LD_kc[j]对应的关键元素El_kc[j]为网格边，转3.10进入按边模板选择；

3.9沿局部方向LD_kc[j]进行按点模板选择。具体步骤为：

3.9.1初始化递归次数n＝N/4。沿4个方向的总模板单元个数为N，则每个局部方向需要选择的模板个数为n＝N/4；

3.9.2定义并初始化递归相关参数。递归过程和模板单元选择过程需要记录和更新一些参数，称为递归相关参数。对按点选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，点参数poiA＝El_kc[j]；

3.9.3根据递归相关参数dirA,cellA,poiA以及网格点的邻接单元数组pCells进行模板单元选择。具体步骤如下：

3.9.3.1以poiA为pCells的索引，获取与网格点Pt_poiA邻接的单元的编号数组pCells[poiA]；

3.9.3.2从数组pCells[poiA]的网格单元中，选出C_kc的模板单元C_tmp(其中tmp是pCells[poiA]的数据项)，使得向量

与dirA的夹角最小，

指以C_cellA中心为起点，C_tmp中心为终点的向量，如图3所示；

3.9.3.3根据C_tmp更新模板单元数组stencil，

方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.9.4令n＝n-1，若n＞0，表示还需递归n次，转3.9.5；若n≤0，沿LD_kc[j]方向的递归结束，转3.11；

3.9.5根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数，并确定下一次递归的模板选择方式。三个递归相关参数更新具体如下：

3.9.5.1令cellA＝tmp，即将C_cellA更新为C_tmp以保证递归沿局部方向LD_kc[j]进行；

3.9.5.2若poiA＝El_tmp[4]，转3.9.5.3；若poiA≠El_tmp[4]，转3.9.5.4；

3.9.5.3优先按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式。令dirA＝LD_tmp[2]。更新后的dirA指向的关键元素是C_tmp的网格边，故下一次递归为按边选择，需要初始化新的递归相关参数：令边参数edgA＝El_tmp[2]，转3.10.3按边选择；

3.9.5.4按网格拓扑信息更新dirA，并确定下一次递归的模板选择方式。以tmp为网格单元数组cs的索引，获取C_tmp的三个顶点，三个顶点分别为Pt_poiA、Pt_kp1和Pt_kp2，其中1≤kp1≤N_pt,1≤kp2≤N_pt。从Pt_kp1、Pt_kp2以及它们构成的网格边E_ke(1≤ke≤N_e)的中点midPt这三个点中选出点selPt(其中selPt＝Pt_kp1或Pt_kp2或midPt)，使得向量

与dirA的夹角最小，是以Pt_poiA为起点，以selPt为终点的向量，如图4和图5所示。基于这样的网格拓扑信息，令

若selPt＝Pt_kp1或Pt_kp2，则令poiA＝kp1或kp2(如图4中selPt＝Pt_kp2，令poiA＝kp2)，转3.9.3按点选择；否则selPt＝midPt，如图5所示，

指向的关键元素为网格边E_ke的中点，则初始化新的递归相关参数：令边参数edgA＝ke，转3.10.3按边选择；

3.10沿局部方向LD_kc[j]进行按边模板选择。具体步骤为：

3.10.1初始化递归次数n＝N/4；

3.10.2定义并初始化递归相关参数。对按边选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，边参数edgA＝El_kc[j]；

3.10.3根据递归相关参数dirA,cellA,edgA以及网格边的邻接单元数组eCells进行模板单元选择。具体步骤如下：

3.10.3.1以edgA为eCells的索引，获取网格边E_edgA的邻接单元数组eCells[edgA]；

3.10.3.2根据eCells[edgA]选择按边选择的模板单元。数组eCells[edgA]包含两个网格单元，C_cellA和C_tmp，其中1≤tmp≤N_c。对于三角形和四边形网格而言，网格边的邻接单元只有两个，即C_cellA和C_tmp。选择C_tmp为本次按边选择的C_kc的模板；

3.10.3.3根据C_tmp更新模板单元数组stencil，

方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.10.4令n＝n-1，若n＞0，则还需递归n次，转3.10.5；若n≤0，沿LD_kc[j]方向的递归结束，转3.11；

3.10.5根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数并确定其模板选择方式。三个递归相关参数更新方法如下：

3.10.5.1令cellA＝tmp；

3.10.5.2按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式。根据第2.3步可知，edgA一定是C_tmp的关键元素，即edgA＝El_tmp[1]或El_tmp[2]或El_tmp[3]。根据第2.3步局部方向与关键元素的对应关系对递归参数进行更新，具体而言：

若edgA＝El_tmp[1]则更新dirA＝LD_tmp[3]，更新边参数edgA＝El_tmp[3]，转3.10.3按边选择；

若edgA＝El_tmp[3]则更新dirA＝LD_tmp[1]，更新边参数edgA＝El_tmp[1]，转3.10.3按边选择；

若edgA＝El_tmp[2]则更新dirA＝LD_tmp[4]，需要初始化新的递归相关参数：点参数poiA＝El_tmp[4]，转3.9.3按点选择；

3.11令j＝j+1，若j＞4，则C_kc的4个局部方向都已进行模板选择，转3.12；若j≤4，转3.8执行下一个局部方向的模板选择；

3.12令kc＝kc+1，若kc＞N_c，则所有网格单元都已选择N个模板单元，即模板单元数组stencil中每个数据项都已生成完毕，转3.13；若kc≤N_c，转3.6执行下一个网格单元的模板选择；

3.13将模板单元数组stencil输出到模板单元文件，通过为每个非结构网格单元分配相同个数模板单元的方式，为每个非结构网格单元分配相同的计算量。模板单元数组

其中数组S_kc中包含N个数据项，即每个网格单元都具有N个模板单元。模板单元个数与计算量成正比，因此每个非结构网格单元分配到了相同的计算量。

第四步，根据处理器信息、模板单元文件、网格相关数据结构以及网格邻接关系数据结构，采用网格剖分方法进行处理器之间网格量层面的负载均衡，生成子网格区域的非结构网格相关文件和子网格区域的模板单元文件，采用子网格区域分配的方式将计算量分配给处理器，方法是：

4.1根据处理器总数及各处理器计算能力，初始化处理器计算能力数组Ps，

其中N_p为处理器总数，P_i为第i个处理器的计算能力，与背景技术1的第2.1步方法相同；

4.2根据处理器计算能力数组Ps及N_p,N_c，计算各处理器的最优网格量。以第i个处理器为例，其最优网格量M_i可表示为：其中1≤k≤N_p且k为整数，与背景技术1的第2.4步方法相同；

4.3使用商业软件(如版本18.0以上的Fluent、版本3.0以上的OpenFOAM等)剖分非结构网格。采用商业软件按照第4.2步计算得到的最优网格量M_i将第一步生成的非结构网格剖分成N_p个子网格区域，第i个子网格区域的网格量等于最优网格量M_i，其中1≤i≤N_p。剖分后生成子网格区域的非结构网格相关文件，与背景技术1的第2.5步方法相同。剖分后生成子网格区域的模板单元文件，子网格区域的模板单元文件是第3.13步生成的模板单元文件的子集。具体而言，子网格区域的模板单元文件只包含子网格区域网格单元对应的模板单元信息，而第3.13步生成的模板单元文件包含所有网格单元的模板单元信息。因此子网格区域的模板单元文件中每个网格单元的计算量仍然相等，记为N_com；

4.4采用子网格区域分配的方式将计算量分配给处理器：将第i个子网格区域的网格量及模板单元文件分配给第i个处理器，其中1≤i≤N_p，第i个处理器获得与之计算能力匹配的计算量。第4.2步和第4.3步保证了第i个子网格区域的网格量等于最优网格量M_i，与第i个处理器的计算能力成正比。同时第i个子网格区域的模板单元文件保证了第i个子网格区域中每个网格单元的计算量均为常数N_com。因此将第i个子网格区域的网格量及模板单元文件分配给第i个处理器后，第i个处理器可获得M_i个网格单元，即获得的计算量为N_com·M_i，与第i个处理器的计算能力成正比。

采用本发明可以达到以下技术效果：

与现有技术相比，本发明在处理器之间网格量负载均衡的基础上，进一步实现了网格单元之间计算量的均衡，使得非结构网格并行应用中，处理器的计算量与其计算能力相匹配，从而提高负载均衡性能和并行计算速度。具体原因为：

1)本发明中每个处理器具有与之计算能力匹配的网格单元量。在第4.3步中结合了网格剖分算法，保证了第i个子网格区域的网格量等于第i个处理器的最优网格量M_i，故每个处理器可以分配与之计算能力相匹配的网格单元量；

2)本发明中每个网格单元具有相同的计算量。本发明在第三步为每个网格单元选择了相同个数的模板单元，模板单元个数与计算量成正比，因此相比于背景技术1，本发明可实现每个网格单元都具有相同的计算量。

综合上述两个层面可知，本发明通过两个层面的负载均衡，使得每个处理器分配到与之计算能力匹配的计算量，提高了整体的负载均衡性能和并行计算速度。

附图说明

图1是背景技术公布的背景技术1的总体流程图。

图2是本发明的总体流程图。

图3是本发明第3.9.3.2步按点选择递归中模板单元选择过程示意图。

图4是本发明第3.9.5.4步按点选择递归中下一次递归方向向量参数更新示意图。

图5是本发明第3.9.5.4步按点选择递归中下一次递归方向向量参数更新示意图。

图6是本发明第三步的流程图。

具体实施方式

图2是本发明的总体流程图。如图2所示，本发明包括以下步骤：

第一步，生成二维非结构网格并初始化非结构网格相关数据结构。

1.1使用专业软件将非结构网格工程应用需要计算的二维空间进行网格化处理，生成非结构网格并输出非结构网格相关文件。

1.2根据非结构网格相关文件初始化相关数据结构，方法是：

1.2.1令网格单元文件为网格单元数组cs，

其中N_c为网格单元总数，C_kc为第kc个网格单元；网格单元采用网格点编号数组表示，C_kc＝{kp₁,...,kp_w,...,kp_W}表示第kc个网格单元由第kp₁,...,kp_w,...,kp_W个网格点构成，是一个W边形网格单元，W为正整数；

1.2.2令网格点文件为网格点数组pts，

其中N_pt为网格点总数，Pt_kp为第kp个网格点；网格点采用笛卡尔坐标数组表示，Pt_kp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数；

1.2.3令网格单元的邻接单元文件为数组cCells，

其中CC_kc表示网格单元C_kc的邻接网格单元编号组成的数组；

1.2.4令网格点的邻接单元文件为数组pCells，

其中PC_kp表示网格点Pt_kp的邻接网格单元编号组成的数组；

第二步，使用非结构网格结构化处理方法对非结构网格相关文件进行局部化，生成非结构网格局部方向数组和关键元素数组，方法是：

2.1根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和非结构网格局部方向数组，方法是：

2.1.1定义关键元素数组其中El_kc为C_kc的关键元素编号组成的数组。初始化

均为{0,0,0,0}；

2.1.2根据

定义非结构网格局部方向数组locDir，

其中LD_kc表示网格单元C_kc的局部方向向量组成的数组。LD_kc包含4个数据项，分别表示C_kc的4个局部方向向量，由El_kc确定：

·LD_kc[1]为El_kc[3]边中点指向El_kc[1]边中点构成的向量；

·若C_kc为四边形，LD_kc[2]为El_kc[4]边中点指向El_kc[2]边中点构成的向量；若C_kc为三角形，LD_kc[2]为El_kc[4]网格点指向El_kc[2]边中点构成的向量；

·LD_kc[3]为El_kc[1]边中点指向El_kc[3]边中点构成的向量；

·若C_kc为四边形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]边中点构成的向量；若C_kc为三角形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]网格点构成的向量。

的每一项都初始化为零向量；

2.2采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法是：

2.2.1定义并初始化阵面推进相关数据结构，方法是：

2.2.1.1定义并初始化阵面数组fronts。fronts为动态数组，数据项个数记为N_f，每一项表示已进行访问的网格点的编号。fronts初始化为位于壁面边界的网格点的编号数组，N_f初始化为位于壁面边界的网格点的总数；

2.2.1.2定义并初始化下一阵面数组nextFronts。nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为N_nf。nextFronts初始化为空集，N_nf初始化为0；

2.2.1.3初始化阵面层数组frontsLevel，其中Fl_kp表示网格点Pt_kp的阵面层值。若kp在数组fronts中，初始化Fl_kp＝0，否则令Fl_kp＝-1；

2.2.1.4定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；

2.2.1.5初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝N_f、令最终尾索引finalTail＝N_f；

2.2.2以head为fronts索引，获取网格点的编号front[head]，为方便表述，记front[head]＝j,1≤j≤N_pt；

2.2.3以j为网格点的邻接单元数组pCells的索引，

获取所有与网格点Pt_j相邻的网格单元编号组成的数组PC_j，记PC_j的数据项总数为PCnum；

2.2.4初始化网格单元索引i＝1；

2.2.5以i为PC_j的索引，获取第i个与Pt_j相邻的网格单元编号PC_j[i]，为方便表述，令o＝PC_j[i],1≤o≤N_c。若El_o[1]＝0，则C_o的关键元素和局部方向还没有确定，转2.2.6；若El_o[1]≠0，转2.2.9；

2.2.6根据网格单元数组cs判断C_o的形状，

若C_o包含三个数据项，说明C_o为三角形，转2.2.7；否则C_o包含四个数据项，转2.2.8；

2.2.7确定三角形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，具体流程如下：

2.2.7.1根据网格单元数组cs获得C_o三个顶点的编号C_o[1]、C_o[2]、C_o[3]；

2.2.7.2以C_o[1]、C_o[2]、C_o[3]为阵面层数组frontsLevel索引，获得对应阵面层值

2.2.7.3对比

按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]重排序，重排序后的C_o[1]、C_o[2]、C_o[3]记为kp₁、kp₂、kp₃，即

2.2.7.4若

转2.2.7.5；否则

转2.2.7.6；

2.2.7.5更新El_o[4]＝kp₃，El_o[2]为kp₁、kp₂组成的网格边的编号，El_o[1]为kp₁、kp₃组成的网格边的编号，El_o[3]为kp₂、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.6更新El_o[4]＝kp₁，El_o[2]为kp₂、kp₃组成的网格边的编号，El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₁、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.7根据El_o，按照第2.1.2步中局部方向的定义，计算得到LD_o，即：

·LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

·LD_o[2]为El_o[4]网格点指向El_o[2]边中点构成的向量；

·LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

·LD_o[4]为El_o[2]边中点指向El_o[4]网格点构成的向量。

2.2.7.8更新nextFronts，具体步骤为：

2.2.7.8.1若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.2若

且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.3若且kp₃不在nextFronts中，将kp₃加入nextFronts，N_nf＝N_nf+1，转2.2.9；

2.2.8确定四边形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，具体流程如下：

2.2.8.1根据网格单元数组cs获得C_o四个顶点的编号C_o[1]、C_o[2]、C_o[3]、C_o[4]；

2.2.8.2以C_o[1]、C_o[2]、C_o[3]、C_o[4]为frontsLevel索引，获得对应阵面层值

2.2.8.3对比

按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]、C_o[4]重排序，将重排序后的C_o[1]、C_o[2]、C_o[3]、C_o[4]记为kp₁、kp₂、kp₃、kp₄，即

2.2.8.4更新El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₃、kp₄组成的网格边的编号。若kp₁、kp₃构成C_o的网格边，El_o[2]更新为kp₁、kp₃组成的网格边的编号，El_o[4]更新为kp₂、kp₄组成的网格边的编号；否则若kp₁、kp₄构成C_o的网格边，El_o[2]更新为kp₁、kp₄组成的网格边的编号，El_o[4]更新为kp₂、kp₃组成的网格边的编号；

2.2.8.5根据El_o，按照第2.1.2步中局部方向的定义，计算得到LD_o，即

·LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

·LD_o[2]为El_o[4]边中点指向El_o[2]边中点构成的向量；

·LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

·LD_o[4]为El_o[2]边中点指向El_o[4]边中点构成的向量。

2.2.8.6更新nextFronts，具体步骤为：

2.2.8.6.1若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.8.6.2若且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1，转2.2.9；

2.2.9令i＝i+1，若i＜PCnum+1，转2.2.5；否则转2.2.10。

2.2.10令head＝head+1，若head＜tail+1，转2.2.2；否则转2.2.11；

2.2.11根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：

2.2.11.1将nextFronts的数据项加入到fronts中，更新N_f＝N_f+N_nf；

2.2.11.2更新frontsLevel中数据项Fl_kp＝n+1，其中kp为nextFronts的数据项；

2.2.11.3更新n＝n+1；

2.2.11.4更新nextFronts为空集，令N_nf＝0；

2.2.11.5更新tail＝finalTail＝N_f；

2.2.12若head＜finalTail+1，转2.2.2；否则所有网格单元的关键元素和局部方向都已确定，得到关键元素数组keyEl和局部方向数组locDir，转第三步；

第三步，基于非结构网格局部方向数组，采用递归选择模板的方式进行网格单元之间计算量层面的负载均衡。如图6所示，包括以下步骤：

3.1根据第一步生成的非结构网格相关文件初始化网格边数组。令网格边文件为网格边数组eds，

其中N_e为网格边总数，E_ke为第ke条网格边。网格边采用顶点(即网格边两端的网格点)编号数组表示，E_ke＝{kp₁,kp₂}表示第ke条网格边的两个顶点分别为第kp₁,kp₂个网格点；

3.2根据第一步生成的非结构网格相关文件初始化网格边的邻接单元数组。令网格边的邻接单元文件为数组eCells，

其中EC_ke表示网格边E_ke的邻接网格单元数组。如EC₁＝{1,2}表示与边E₁相邻接的网格单元为C₁,C₂；

3.3定义并初始化模板单元个数N。N为常数，表示需要为每个网格单元选择的模板单元个数，也为4个方向的总模板单元个数，N与网格单元的计算量成正比。N的取值与具体应用相关，N为正整数，一般为4的倍数；

3.4定义并初始化模板单元数组stencil，

其中S_kc是C_kc的模板单元编号组成的数组。由于本发明将为每个网格单元选择N个模板单元，故S_kc有N个数据项，S_kc的每个数据项为C_kc的模板单元编号。S_kc的每个数据项都初始化为0。

3.5初始化网格单元递归的单元索引kc＝1；

3.6以kc为局部方向数组locDir和关键元素数组keyEl的索引，分别获取C_kc的4个局部方向向量：LD_kc[1]、LD_kc[2]、LD_kc[3]、LD_kc[4]，以及4个局部方向对应的4个关键元素：El_kc[1]、El_kc[2]、El_kc[3]、El_kc[4]；

3.7初始化局部方向索引j＝1；

3.8若C_kc为三角形网格且j＝4，LD_kc[j]对应的关键元素El_kc[4]为网格点，转3.9进行按点模板选择；若C_kc不为三角形网格或j≠4，LD_kc[j]对应的关键元素El_kc[j]为网格边，转3.10进入按边模板选择；

3.9沿局部方向LD_kc[j]进行按点模板选择。具体步骤为：

3.9.1初始化递归次数n＝N/4。沿4个方向的总模板单元个数为N，则每个局部方向需要选择的模板个数为n＝N/4；

3.9.2定义并初始化递归相关参数。递归过程和模板单元选择过程需要记录和更新一些参数，称为递归相关参数。对按点选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，点参数poiA＝El_kc[j]；

3.9.3根据递归相关参数dirA,cellA,poiA以及网格点的邻接单元数组pCells进行模板单元选择。具体步骤如下：

3.9.3.1以poiA为pCells的索引，获取与网格点Pt_poiA邻接的单元的编号数组pCells[poiA]；

3.9.3.2从数组pCells[poiA]的网格单元中，选出C_kc的模板单元C_tmp(其中tmp是pCells[poiA]的数据项)，使得向量与dirA的夹角最小，

指以C_cellA中心为起点，C_tmp中心为终点的向量，如图3所示；图3是本发明第3.9.3.2步按点选择递归中模板单元选择过程示意图。图3中C_cellA、C_tmp为网格单元，箭头所示dirA表示方向向量，虚线箭头为单元中心连线构成的向量，Pt_poiA表示网格点。图3描述从网格单元C_cellA出发，沿着方向向量dirA进行按点(Pt_poiA)选择模板单元的过程，其中Pt_poiA为dirA与C_cellA的交点。除C_cellA外，与点Pt_poiA相邻的网格单元共有5个(C₁,C₂,C₃,C₄,C₅)。以C_cellA中心为起点，分别以C₁,C₂,C₃,C₄,C₅中心为终点，可以构成5个向量，即

由图可看出，与方向向量dirA夹角最小，故将选出C₃作为此次按点选择的模板单元C_tmp。

3.9.3.3根据C_tmp更新模板单元数组stencil，

方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.9.4令n＝n-1，若n＞0，表示还需递归n次，转3.9.5；若n≤0，沿LD_kc[j]方向的递归结束，转3.11；

3.9.5根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数，并确定下一次递归的模板选择方式。三个递归相关参数更新具体如下：

3.9.5.1令cellA＝tmp，即将C_cellA更新为C_tmp以保证递归沿局部方向LD_kc[j]进行；

3.9.5.2若poiA＝El_tmp[4]，转3.9.5.3；若poiA≠El_tmp[4]，转3.9.5.4；

3.9.5.3优先按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式。令dirA＝LD_tmp[2]。更新后的dirA指向的关键元素是C_tmp的网格边，故下一次递归为按边选择，需要初始化新的递归相关参数：令边参数edgA＝El_tmp[2]，转3.10.3按边选择；

3.9.5.4按网格拓扑信息更新dirA，并确定下一次递归的模板选择方式。以tmp为网格单元数组cs的索引，获取C_tmp的三个顶点，三个顶点分别为Pt_poiA、Pt_kp1和Pt_kp2，其中1≤kp1≤N_pt,1≤kp2≤N_pt。从Pt_kp1、Pt_kp2以及它们构成的网格边E_ke(1≤ke≤N_e)的中点midPt这三个点中选出点selPt(其中selPt＝Pt_kp1或Pt_kp2或midPt)，使得向量

与dirA的夹角最小，是以Pt_poiA为起点，以selPt为终点的向量，如图4和图5所示。

图4是本发明第3.9.5.4步按点选择递归中下一次递归方向向量参数更新示意图。图4中点Pt_kp1、Pt_kp2、Pt_poiA为网格单元C_tmp的网格点，E_ke为Pt_kp1与Pt_kp2构成的边，midPt为E_ke的中点，箭头所示dirA表示方向向量，虚线箭头为Pt_poiA与Pt_kp1、Pt_kp2、midPt构成的向量。图4描述沿dirA按点(Pt_poiA)选择出模板单元C_tmp后，确定下一次递归为按点(Pt_kp2)选择的过程。以Pt_poiA为起点，分别以点Pt_kp1、点Pt_kp2、边中点midPt为终点可以构成3个向量，即

由图可看出，向量

与方向向量dirA夹角最小，选择Pt_kp2作为selPt，下一次递归为按点(Pt_kp2)选择。

图5是本发明第3.9.5.4步按点选择递归中下一次递归方向向量参数更新示意图。图5中点Pt_kp1、Pt_kp2、Pt_poiA为网格单元C_tmp的网格点，E_ke为Pt_kp1与Pt_kp2构成的边，midPt为E_ke的中点，箭头所示dirA表示方向向量，虚线箭头为Pt_poiA与Pt_kp1、Pt_kp2、midPt构成的向量。图5描述沿dirA按点(Pt_poiA)选择出模板单元C_tmp后，确定下一次递归为按边(E_ke)选择。以Pt_poiA为起点，分别以点Pt_kp1、点Pt_kp2、边中点midPt为终点可以构成3个向量，即

由图可看出，向量

与方向向量dirA夹角最小，选择midPt作为selPt，下一次递归为按点(E_ke)选择。

基于这样的网格拓扑信息，令

若selPt＝Pt_kp1或Pt_kp2，则令poiA＝kp1或kp2(如图4中selPt＝Pt_kp2，令poiA＝kp2)，转3.9.3按点选择；否则selPt＝midPt，如图5所示，

指向的关键元素为网格边E_ke的中点，则初始化新的递归相关参数：令边参数edgA＝ke，转3.10.3按边选择；

3.10沿局部方向LD_kc[j]进行按边模板选择。具体步骤为：

3.10.1初始化递归次数n＝N/4；

3.10.2定义并初始化递归相关参数。对按边选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，边参数edgA＝El_kc[j]；

3.10.3根据递归相关参数dirA,cellA,edgA以及网格边的邻接单元数组eCells进行模板单元选择。具体步骤如下：

3.10.3.1以edgA为eCells的索引，获取网格边E_edgA的邻接单元数组eCells[edgA]；

3.10.3.2根据eCells[edgA]选择按边选择的模板单元。数组eCells[edgA]包含两个网格单元，C_cellA和C_tmp，其中1≤tmp≤N_c。对于三角形和四边形网格而言，网格边的邻接单元只有两个，即C_cellA和C_tmp。选择C_tmp为本次按边选择的C_kc的模板；

3.10.3.3根据C_tmp更新模板单元数组stencil，

方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.10.4令n＝n-1，若n＞0，则还需递归n次，转3.10.5；若n≤0，沿LD_kc[j]方向的递归结束，转3.11；

3.10.5根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数并确定其模板选择方式。三个递归相关参数更新方法如下：

3.10.5.1令cellA＝tmp；

3.10.5.2按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式。根据第2.3步可知，edgA一定是C_tmp的关键元素，即edgA＝El_tmp[1]或El_tmp[2]或El_tmp[3]。根据第2.3步局部方向与关键元素的对应关系对递归参数进行更新，具体而言：

若edgA＝El_tmp[1]则更新dirA＝LD_tmp[3]，更新边参数edgA＝El_tmp[3]，转3.10.3按边选择；

若edgA＝El_tmp[3]则更新dirA＝LD_tmp[1]，更新边参数edgA＝El_tmp[1]，转3.10.3按边选择；

若edgA＝El_tmp[2]则更新dirA＝LD_tmp[4]，需要初始化新的递归相关参数：点参数poiA＝El_tmp[4]，转3.9.3按点选择；

3.11令j＝j+1，若j＞4，则C_kc的4个局部方向都已进行模板选择，转3.12；若j≤4，转3.8执行下一个局部方向的模板选择；

3.12令kc＝kc+1，若kc＞N_c，则所有网格单元都已选择N个模板单元，即模板单元数组stencil中每个数据项都已生成完毕，转3.13；若kc≤N_c，转3.6执行下一个网格单元的模板选择；

3.13将模板单元数组stencil输出到模板单元文件，通过为每个非结构网格单元分配相同个数模板单元的方式，为每个非结构网格单元分配相同的计算量。模板单元数组

其中数组S_kc中包含N个数据项，即每个网格单元都具有N个模板单元。模板单元个数与计算量成正比，因此每个非结构网格单元分配到了相同的计算量。

第四步，根据处理器信息、模板单元文件、网格相关数据结构以及网格邻接关系数据结构，采用网格剖分方法进行处理器之间网格量层面的负载均衡，生成子网格区域的非结构网格相关文件和子网格区域的模板单元文件，采用子网格区域分配的方式将计算量分配给处理器，方法是：

4.1根据处理器总数及各处理器计算能力，初始化处理器计算能力数组Ps，其中N_p为处理器总数，P_i为第i个处理器的计算能力，与背景技术1的第2.1步方法相同；

4.2根据处理器计算能力数组Ps及N_p,N_c，计算各处理器的最优网格量。以第i个处理器为例，其最优网格量M_i可表示为：

其中1≤k≤N_p且k为整数，与背景技术1的第2.4步方法相同；

4.3使用商业软件(如版本18.0以上的Fluent、版本3.0以上的OpenFOAM等)剖分非结构网格。采用商业软件按照第4.2步计算得到的最优网格量M_i将第一步生成的非结构网格剖分成N_p个子网格区域，第i个子网格区域的网格量等于最优网格量M_i，其中1≤i≤N_p。剖分后生成子网格区域的非结构网格相关文件，与背景技术1的第2.5步方法相同。剖分后生成子网格区域的模板单元文件，子网格区域的模板单元文件是第3.13步生成的模板单元文件的子集。具体而言，子网格区域的模板单元文件只包含子网格区域网格单元对应的模板单元信息，而第3.13步生成的模板单元文件包含所有网格单元的模板单元信息。因此子网格区域的模板单元文件中每个网格单元的计算量仍然相等，记为N_com；

4.4采用子网格区域分配的方式将计算量分配给处理器：将第i个子网格区域的网格量及模板单元文件分配给第i个处理器，其中1≤i≤N_p，第i个处理器获得与之计算能力匹配的计算量。第4.2步和第4.3步保证了第i个子网格区域的网格量等于最优网格量M_i，与第i个处理器的计算能力成正比。同时第i个子网格区域的模板单元文件保证了第i个子网格区域中每个网格单元的计算量均为常数N_com。因此将第i个子网格区域的网格量及模板单元文件分配给第i个处理器后，第i个处理器可获得M_i个网格单元，即获得的计算量为N_com·M_i，与第i个处理器的计算能力成正比。

Claims (9)

1.一种均衡非结构网格单元计算量的并行负载均衡方法，其特征在于包括以下步骤：

第一步，生成二维非结构网格，方法是：

1.1 将非结构网格工程应用需要计算的二维空间进行网格化处理，生成非结构网格并输出二维非结构网格相关文件；二维非结构网格相关文件包括网格单元文件、网格点文件、网格边文件、网格单元的邻接单元文件、网格点的邻接单元文件、网格边的邻接单元文件；

1.2 根据非结构网格相关文件初始化相关数据结构，方法是：

1.2.1 令网格单元文件为网格单元数组cs，其中N_c为网格单元总数，C_kc为第kc个网格单元；网格单元采用网格点编号数组表示，C_kc＝{kp₁,...,kp_w,...,kp_W}表示第kc个网格单元由第kp₁,...,kp_w,...,kp_W个网格点构成，是一个W边形网格单元，W为正整数；

1.2.2 令网格点文件为网格点数组pts，其中N_pt为网格点总数，Pt_kp为第kp个网格点；网格点采用笛卡尔坐标数组表示，Pt_kp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数；

1.2.3 令网格单元的邻接单元文件为数组cCells，其中CC_kc表示网格单元C_kc的邻接网格单元编号组成的数组；

1.2.4 令网格点的邻接单元文件为数组pCells，

其中PC_kp表示网格点Pt_kp的邻接网格单元编号组成的数组；

第二步，使用非结构网格结构化处理方法对非结构网格相关文件进行局部化，生成非结构网格局部方向数组和关键元素数组，方法是：

2.1 根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和非结构网格局部方向数组，方法是：

2.1.1 定义关键元素数组其中El_kc为C_kc的关键元素编号组成的数组；三角形网格单元的关键元素为3条网格边和1个网格点，四边形网格单元的关键元素为4条网格边，初始化

均为{0,0,0,0}；

2.1.2 根据keyEl定义非结构网格局部方向数组locDir，

其中LD_kc表示网格单元C_kc的局部方向向量组成的数组，LD_kc包含4个数据项，分别表示C_kc的4个局部方向向量，由El_kc确定：

LD_kc[1]为El_kc[3]边中点指向El_kc[1]边中点构成的向量；

若C_kc为四边形，LD_kc[2]为El_kc[4]边中点指向El_kc[2]边中点构成的向量；若C_kc为三角形，LD_kc[2]为El_kc[4]网格点指向El_kc[2]边中点构成的向量；

LD_kc[3]为El_kc[1]边中点指向El_kc[3]边中点构成的向量；

若C_kc为四边形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]边中点构成的向量；若C_kc为三角形，LD_kc[4]为El_kc[2]边中点指向El_kc[4]网格点构成的向量；

将

的每一项都初始化为零向量；

2.2 采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法是：

2.2.1 定义并初始化阵面推进相关数据结构，方法是：

2.2.1.1 定义并初始化阵面数组fronts，fronts为动态数组，数据项个数记为N_f，每一项表示已进行访问的网格点的编号；fronts初始化为位于壁面边界即物体表面的网格点的编号数组，N_f初始化为位于壁面边界的网格点的总数；

2.2.1.2 定义并初始化下一阵面数组nextFronts，nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为N_nf；nextFronts初始化为空集，N_nf初始化为0；

2.2.1.3 初始化阵面层数组frontsLevel，

其中Fl_kp表示网格点Pt_kp的阵面层值；若kp在数组fronts中，初始化Fl_kp＝0，否则令Fl_kp＝-1；

2.2.1.4 定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；

2.2.1.5 初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝N_f、令最终尾索引finalTail＝N_f；

2.2.2 以head为fronts索引，获取网格点的编号front[head]，记front[head]＝j,1≤j≤N_pt；

2.2.3 以j为网格点的邻接单元数组pCells的索引，

获取所有与网格点Pt_j相邻的网格单元编号组成的数组PC_j，记PC_j的数据项总数为PCnum；

2.2.4 初始化网格单元索引i＝1；

2.2.5 以i为PC_j的索引，获取第i个与Pt_j相邻的网格单元编号PC_j[i]，令o＝PC_j[i],1≤o≤N_c；若El_o[1]＝0，则C_o的关键元素和局部方向还没有确定，转2.2.6；若El_o[1]≠0，转2.2.9；

2.2.6 根据网格单元数组cs判断C_o的形状，

若C_o包含三个数据项，说明C_o为三角形，转2.2.7；否则C_o包含四个数据项，转2.2.8；

2.2.7 确定三角形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，转2.2.9；

2.2.8 确定四边形网格单元C_o的关键元素，得到El_o，根据El_o计算得到局部向量LD_o，并更新nextFronts，转2.2.9；

2.2.9 令i＝i+1，若i<PCnum+1，转2.2.5；否则转2.2.10；

2.2.10 令head＝head+1，若head<tail+1，转2.2.2；否则转2.2.11；

2.2.11 根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：

2.2.11.1 将nextFronts的数据项加入到fronts中，更新N_f＝N_f+N_nf；

2.2.11.2 更新frontsLevel中数据项Fl_kp＝n+1，其中kp为nextFronts的数据项；

2.2.11.3 更新n＝n+1；

2.2.11.4 更新nextFronts为空集，令N_nf＝0；

2.2.11.5 更新tail＝finalTail＝N_f；

2.2.12 若head<finalTail+1，转2.2.2；否则所有网格单元的关键元素和局部方向都已确定，得到关键元素数组keyEl和局部方向数组locDir，转第三步；

第三步，基于非结构网格局部方向数组，采用递归选择模板的方式进行网格单元之间计算量层面的负载均衡，方法是：

3.1 根据第一步生成的非结构网格相关文件初始化网格边数组，令网格边文件为网格边数组eds，其中N_e为网格边总数，E_ke为第ke条网格边；网格边采用顶点即网格边两端的网格点编号数组表示，E_ke＝{kp₁,kp₂}表示第ke条网格边的两个顶点分别为第kp₁,kp₂个网格点；

3.2 根据第一步生成的非结构网格相关文件初始化网格边的邻接单元数组，令网格边的邻接单元文件为数组eCells，其中EC_ke表示网格边E_ke的邻接网格单元数组；

3.3 定义并初始化模板单元个数N，N为正整数，表示需要为每个网格单元选择的模板单元个数，也为4个方向的总模板单元个数，N与网格单元的计算量成正比；

3.4 定义并初始化模板单元数组stencil，

其中S_kc是C_kc的模板单元编号组成的数组，S_kc有N个数据项；将S_kc的每个数据项都初始化为0；

3.5 初始化网格单元递归的单元索引kc＝1；

3.6 以kc为局部方向数组locDir和关键元素数组keyEl的索引，分别获取C_kc的4个局部方向向量：LD_kc[1]、LD_kc[2]、LD_kc[3]、LD_kc[4]，以及4个局部方向对应的4个关键元素：El_kc[1]、El_kc[2]、El_kc[3]、El_kc[4]；

3.7 初始化局部方向索引j＝1；

3.8 若C_kc为三角形网格且j＝4，LD_kc[j]对应的关键元素El_kc[4]为网格点，转3.9进行按点模板选择；若C_kc不为三角形网格或j≠4，LD_kc[j]对应的关键元素El_kc[j]为网格边，转3.10进入按边模板选择；

3.9 沿局部方向LD_kc[j]进行按点模板选择，具体步骤为：

3.9.1 初始化递归次数n＝N/4；

3.9.2 定义并初始化递归相关参数，对按点选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，点参数poiA＝El_kc[j]；

3.9.3 根据递归相关参数dirA,cellA,poiA以及网格点的邻接单元数组pCells进行模板单元选择，具体步骤如下：

3.9.3.1 以poiA为pCells的索引，获取与网格点Pt_poiA邻接的单元的编号数组pCells[poiA]；

3.9.3.2 从数组pCells[poiA]的网格单元中，选出C_kc的模板单元C_tmp，其中tmp是pCells[poiA]的数据项，使得向量与dirA的夹角最小，

指以C_cellA中心为起点，C_tmp中心为终点的向量；

3.9.3.3 根据C_tmp更新模板单元数组stencil，

方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.9.4 令n＝n-1，若n>0，转3.9.5；若n≤0，转3.11；

3.9.5 根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数，并确定下一次递归的模板选择方式，dirA,cellA,poiA更新方法是：

3.9.5.1 令cellA＝tmp，即将C_cellA更新为C_tmp以保证递归沿局部方向LD_kc[j]进行；

3.9.5.2 若poiA＝El_tmp[4]，转3.9.5.3；若poiA≠El_tmp[4]，转3.9.5.4；

3.9.5.3 优先按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式，方法是令dirA＝LD_tmp[2]；令边参数edgA＝El_tmp[2]，转3.10.3；

3.9.5.4 按网格拓扑信息更新dirA，并确定下一次递归的模板选择方式：以tmp为网格单元数组cs的索引，获取C_tmp的三个顶点，三个顶点分别为Pt_poiA、Pt_kp1和Pt_kp2，其中1≤kp1≤N_pt,1≤kp2≤N_pt；从Pt_kp1、Pt_kp2以及它们构成的网格边E_ke的中点midPt这三个点中选出点selPt，使得向量

与dirA的夹角最小，

是以Pt_poiA为起点，以selPt为终点的向量；基于这样的网格拓扑信息，令

若selPt＝Pt_kp1或Pt_kp2，则令poiA＝kp1或kp2，转3.9.3按点选择；否则selPt＝midPt，

指向的关键元素为网格边E_ke的中点，则初始化新的递归相关参数：令边参数edgA＝ke，转3.10.3按边选择；

3.10 沿局部方向LD_kc[j]进行按边模板选择，方法为：

3.10.1 初始化递归次数n＝N/4；

3.10.2 定义并初始化递归相关参数，对按边选择的递归参数进行初始化：令方向向量参数dirA＝LD_kc[j]，单元参数cellA＝kc，边参数edgA＝El_kc[j]；

3.10.3 根据递归相关参数dirA,cellA,edgA以及网格边的邻接单元数组eCells进行模板单元选择，具体方法为：

3.10.3.1 以edgA为eCells的索引，获取网格边E_edgA的邻接单元数组eCells[edgA]；

3.10.3.2 根据eCells[edgA]选择按边选择的模板单元：选择C_tmp为本次按边选择的C_kc的模板；

3.10.3.3 根据C_tmp更新模板单元数组stencil，方法是将数组S_kc中第一个为0的数据项更新为C_tmp；

3.10.4 令n＝n-1，若n>0，转3.10.5；若n≤0，转3.11；

3.10.5 根据递归相关参数dirA,cellA,poiA和选出的模板单元C_tmp更新下一次递归的相关参数并确定其模板选择方式，dirA,cellA,poiA更新方法如下：

3.10.5.1 令cellA＝tmp；

3.10.5.2 按C_tmp的局部方向更新dirA，并确定下一次递归的模板选择方式：

若edgA＝El_tmp[1]则更新dirA＝LD_tmp[3]，更新边参数edgA＝El_tmp[3]，转3.10.3按边选择；

若edgA＝El_tmp[3]则更新dirA＝LD_tmp[1]，更新边参数edgA＝El_tmp[1]，转3.10.3按边选择；

若edgA＝El_tmp[2]则更新dirA＝LD_tmp[4]，点参数poiA＝El_tmp[4]，转3.9.3按点选择；

3.11 令j＝j+1，若j>4，则C_kc的4个局部方向都已进行模板选择，转3.12；若j≤4，转3.8；

3.12 令kc＝kc+1，若kc>N_c，则所有网格单元都已选择N个模板单元，即模板单元数组stencil中每个数据项都已生成完毕，转3.13；若kc≤N_c，转3.6；

3.13 将模板单元数组stencil输出到模板单元文件，通过为每个非结构网格单元分配相同个数模板单元的方式，为每个非结构网格单元分配相同的计算量；模板单元数组

其中数组S_kc中包含N个数据项，即每个网格单元都具有N个模板单元；

第四步，根据处理器信息、模板单元文件、网格相关数据结构以及网格邻接关系数据结构，采用网格剖分方法进行处理器之间网格量层面的负载均衡，生成子网格区域的非结构网格相关文件和子网格区域的模板单元文件，采用子网格区域分配的方式将计算量分配给处理器，方法是：

4.1 根据处理器总数及各处理器计算能力，初始化处理器计算能力数组Ps，

其中N_p为处理器总数，P_i为第i个处理器的计算能力；

4.2 根据处理器计算能力数组Ps及N_p,N_c，计算各处理器的最优网格量：第i个处理器的最优网格量

其中1≤k≤N_p且k为整数；

4.3 剖分非结构网格：按照M_i将第一步生成的非结构网格剖分成N_p个子网格区域，第i个子网格区域的网格量等于最优网格量M_i，其中1≤i≤N_p；剖分后生成子网格区域的非结构网格相关文件，剖分后还生成子网格区域的模板单元文件，子网格区域的模板单元文件是第3.13步生成的模板单元文件的子集；子网格区域的模板单元文件中每个网格单元的计算量仍然相等，记为N_com；

4.4 采用子网格区域分配的方式将计算量分配给处理器：将第i个子网格区域的网格量及模板单元文件分配给第i个处理器，其中1≤i≤N_p，第i个处理器获得M_i个网格单元，即获得的计算量为N_com·M_i。

2.如权利要求1所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于1.1步所述将工程应用需要计算的二维空间进行网格化处理的方法是使用专业软件，专业软件指版本18.0以上的Fluent软件、版本13.0以上的FUN3D软件、版本3.0以上的OpenFOAM、版本5.0以上的SU2。

3.如权利要求1所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于2.2.7步所述确定三角形网格单元C_o的关键元素El_o和局部方向LD_o，并更新nextFronts的流程为：

2.2.7.1 根据网格单元数组cs获得C_o三个顶点的编号C_o[1]、C_o[2]、C_o[3]；

2.2.7.2 以C_o[1]、C_o[2]、C_o[3]为阵面层数组frontsLevel索引，获得对应阵面层值

2.2.7.3 对比按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]重排序，重排序后的C_o[1]、C_o[2]、C_o[3]记为kp₁、kp₂、kp₃，即

2.2.7.4 若

转2.2.7.5；否则

转2.2.7.6；

2.2.7.5 更新El_o[4]＝kp₃，El_o[2]为kp₁、kp₂组成的网格边的编号，El_o[1]为kp₁、kp₃组成的网格边的编号，El_o[3]为kp₂、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.6 更新El_o[4]＝kp₁，El_o[2]为kp₂、kp₃组成的网格边的编号，El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₁、kp₃组成的网格边的编号，转2.2.7.7；

2.2.7.7 根据El_o，按照第2.1.2步中局部方向的定义计算得到LD_o；

2.2.7.8 更新nextFronts，具体步骤为：

2.2.7.8.1 若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.2 若

且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1；

2.2.7.8.3 若

且kp₃不在nextFronts中，将kp₃加入nextFronts，N_nf＝N_nf+1，转2.2.9。

4.如权利要求3所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于2.2.7.7步所述根据El_o计算得到LD_o的方法是：

LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

LD_o[2]为El_o[4]网格点指向El_o[2]边中点构成的向量；

LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

LD_o[4]为El_o[2]边中点指向El_o[4]网格点构成的向量。

5.如权利要求1所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于2.2.8步所述确定四边形网格单元C_o的关键元素El_o和局部方向LD_o，并更新nextFronts的流程如下：

2.2.8.1 根据网格单元数组cs获得C_o四个顶点的编号C_o[1]、C_o[2]、C_o[3]、C_o[4]；

2.2.8.2 以C_o[1]、C_o[2]、C_o[3]、C_o[4]为frontsLevel索引，获得对应阵面层值

2.2.8.3 对比

按照阵面层值从小到大的顺序，将C_o[1]、C_o[2]、C_o[3]、C_o[4]重排序，将重排序后的C_o[1]、C_o[2]、C_o[3]、C_o[4]记为kp₁、kp₂、kp₃、kp₄，即

2.2.8.4 更新El_o[1]为kp₁、kp₂组成的网格边的编号，El_o[3]为kp₃、kp₄组成的网格边的编号；若kp₁、kp₃构成C_o的网格边，El_o[2]更新为kp₁、kp₃组成的网格边的编号，El_o[4]更新为kp₂、kp₄组成的网格边的编号；否则若kp₁、kp₄构成C_o的网格边，El_o[2]更新为kp₁、kp₄组成的网格边的编号，El_o[4]更新为kp₂、kp₃组成的网格边的编号；

2.2.8.5 根据El_o，按照第2.1.2步中局部方向的定义，计算得到LD_o；

2.2.8.6 更新nextFronts，具体步骤为：

2.2.8.6.1 若

且kp₁不在nextFronts中，将kp₁加入nextFronts，N_nf＝N_nf+1；

2.2.8.6.2 若

且kp₂不在nextFronts中，将kp₂加入nextFronts，N_nf＝N_nf+1。

6.如权利要求5所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于2.2.8.5步所述根据El_o计算得到LD_o的方法是：

LD_o[1]为El_o[3]边中点指向El_o[1]边中点构成的向量；

LD_o[2]为El_o[4]边中点指向El_o[2]边中点构成的向量；

LD_o[3]为El_o[1]边中点指向El_o[3]边中点构成的向量；

LD_o[4]为El_o[2]边中点指向El_o[4]边中点构成的向量。

7.如权利要求1所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于3.3步所述N的取值为4的倍数。

8.如权利要求1所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于4.3步剖分非结构网格的方法是采用含剖分算法的软件。

9.如权利要求8所述的均衡非结构网格单元计算量的并行负载均衡方法，其特征在于所述含剖分算法的软件指版本18.0以上的Fluent软件、版本3.0以上的OpenFOAM。

Images