一种数值模拟计算中的负载均衡方法

摘要

本发明公开了一种数值模拟计算中的负载均衡方法，涉及稀薄气体动力学数值模拟领域，本方法利用网格剖分METIS串行库实现MPI并行DSMC方法的基于模拟粒子数分布为权重的计算网格重新剖分，并根据剖分前后网格单元的全局编号不变这一关系完成流场信息的整体映射，最终实现了DSMC方法数值模拟计算中的负载均衡。

说明书

技术领域

本发明涉及稀薄气体动力学数值模拟领域，具体地，涉及一种数值模拟计算中的负载均衡方法。

背景技术

直接模拟蒙特卡洛(Direct Simulation Monte Carlo, DSMC)方法是求解稀薄气体动力学的主要工具。这一方法直接从流动物理出发，利用大量的模拟分子来仿真真实的气体流动。与一般的计算流体力学数值格式相比，DSMC方法计算量较大，尤其是在流动的稀薄程度不是很高的情况。将这一方法应用于复杂外形时通常需要采用并行计算。常见的并行方式有MPI、OpenMP和GPU等，其中MPI的使用较为普遍。DSMC方法MPI并行程序一般采用SPMD（Single-Program Multuple-Data）的方式来进行。利用网格剖分工具METIS，将计算网格根据网格数目划分为若干个大致相等的分区，然而将这些分区加载到不同的计算进程上，网格之间的信息传递通过MPI消息来完成。这样可大大提升DSMC方法的计算效率，拓展其应用领域和范围。

然而，DSMC方法的计算负载并不由计算网格直接决定，而是主要受模拟粒子数目控制。随着计算的进行，计算区域内的模拟粒子会不断地运动、碰撞，流场结构也会动态变化直至达到定常状态。模拟粒子在高密度区域集中，而在低密度区域比较稀疏，从而造成计算负载的不平衡，影响并行计算效率。根据网格数目在计算开始时刻进行的静态分区，不能满足高效并行计算的需求。

发明内容

本发明目的是为了满足DSMC方法的高效并行计算的需求，有效的实现DSMC方法的计算负载均衡，提高并行计算的效率。

为实现上述目的，本发明提供了一种数值模拟计算中的负载均衡方法，所述方法包括：

在根进程内读入网格文件，对所述网格文件进行第一次剖分，获得若干个第一分区网格文件，将若干个所述第一分区网格文件分配至各处理器进程；

各处理器进程分别读入对应的第一分区网格文件，开始气体动力数值模拟并行计算；

判断流场内的气体流动是否达到定常，当气体流动达到定常时，统计各处理器进程内网格单元的全局编号以及网格单元内的模拟粒子数目；

在根进程内，将网格单元的全局编号与网格单元内粒子数目进行匹配，获得匹配结果；

在根进程内，基于上述匹配结果以网格单元内模拟粒子数目为权重对所述网格文件进行第二次剖分，获得若干个第二分区网格文件；

统计各处理器进程内的模拟粒子信息，将所述模拟粒子信息存入第一数组；

根据所述第二次剖分的分区结果，按照各处理器进程的编号信息对所述第一数组进行排序，获得排序结果；

基于所述排序结果，将流场内的模拟粒子信息映射到对应的所述第二分区网格文件中；

各处理器进程分别读入对应的第二分区网格文件，进行气体动力数值模拟并行计算。

其中，本方法的原理为：本申请发明人研究发现DSMC方法的计算负载并不由计算网格直接决定，而是主要受模拟粒子数目控制。随着计算的进行，计算区域内的模拟粒子会不断地运动、碰撞，流场结构也会动态变化直至达到定常状态。模拟粒子在高密度区域集中，而在低密度区域比较稀疏，从而造成计算负载的不平衡，影响并行计算效率。根据网格数目在计算开始时刻进行的静态分区，不能满足高效并行计算的需求。

而网格剖分工具METIS在执行时可以指定节点权重。如果将网格单元内的模拟粒子数目作为权重进行分区，那么得到的各分区将保持模拟粒子总数大致相等，这意味着在之后的续算中计算负载基本平衡，并行效率将显著提升。因此，本方法利用网格剖分METIS串行库实现MPI并行DSMC方法的基于模拟粒子数分布为权重的计算网格重新剖分，并根据剖分前后网格单元的全局编号不变这一关系完成流场信息的整体映射，最终实现了DSMC方法数值模拟计算中的负载均衡。

优选的，所述方法还包括将所述全局编号以及所述模拟粒子数目收集至主进程的相应数组中。因现在调用的分区程序是在主进程中串行执行的，所以将这些作为权重信息的数据收集至主进程中。

优选的，本方法通过MPI全局函数MPI_Gatherv()将所述全局编号以及所述模拟粒子数目收集至主进程的相应数组中。其中，MPI_Gatherv()函数用于长度不一样的数据收集，MPI_Gather()函数用于长度一样的数据收集。

优选的，所述方法具体包括：

统计各处理器进程中的网格单元数nCellInProc，利用MPI中的MPI_Gather()函数将各处理器进程的网格单元数信息收集至根进程中的nCellInProcArray数组中；

在根进程内计算相对位移数组displs；

基于所述网格单元数nCellInProc和所述数组displs，调用所述MPI全局函数MPI_Gatherv()将所述全局编号以及所述模拟粒子数目收集至主进程的相应数组中。

优选的，本方法中网格单元的全局编号与网格单元内粒子数目的匹配方式为：基于网格单元的全局编号对网格单元内粒子数信息进行排序匹配。

其中，网格单元内的粒子数将作为权重信息传递给METIS分区程序，排序匹配的作用是使得模拟粒子数与网格单元编号正确匹配。

优选的，本方法采用网格剖分工具METIS对所述网格文件进行剖分。

优选的，所述模拟粒子信息包括：模拟粒子的所在单元编号、模拟粒子所属的组份编号、模拟粒子的空间三维坐标、模拟粒子在三维方向的速度大小、模拟粒子转动能和模拟粒子振动能。

优选的，本方法按照各处理器进程编号由小到大的顺序对所述第一数组进行排序。其中，这样设计的目的是MPI发送时是根据进程编号从小到大的顺序进行的。

优选的，本方法还包括：

自定义获得模拟粒子信息的第一数据结构；

基于所述第一数据结构，获得MPI数据结构；

基于所述MPI数据结构，基于MPI函数将流场内的模拟粒子信息映射到对应的所述第二分区网格文件中；其中，所述第一数据结构为：

procID，zoneID，cellID，specieID，XYZ(3)，UVW(3)，erot，evib；

其中，procID是模拟粒子发往的进程编号，zoneID是模拟粒子所在网格分区编号，cellID是模拟粒子所在网格单元编号，specieID是模拟粒子组分编号，XYZ(3)是模拟粒子的三维空间坐标，UVW(3)是模拟粒子的速度三分量，erot是模拟粒子的转动能，evib是模拟粒子的振动能。

其中，本发明根据模拟粒子的数据结构定义了相应的MPI自定义数据结构。模拟粒子信息中既包含整型数据，也包含浮点型数据，且数据项较多。MPI的内置数据类型不能直接发送接收粒子信息。如果单独发送，则较为繁琐。本发明根据MPI的相关语法，定义了模拟粒子的专属数据结构，即第一数据结构，使得可以直接发送整个粒子的所有信息。

优选的，所述基于所述第一数据结构，获得MPI数据结构，具体包括：

定义blocklengths(0) = 4，blocklengths(1) = 8，oldtypes(0) = MPI_INTEGER，oldtypes(1) = MPI_DOUBLE_PRECISION；

调用MPI_GET_ADDRESS函数计算出偏移量offsets(0)和offsets(1)

调用MPI_TYPE_CREATE_STRUCT函数和MPI_TYPE_COMMIT函数完成moletype的定义和注册，其中，blocklengths(0)为第一个数据块的长度，blocklengths(1)第二个数据块的长度，oldtypes(0)为第一个数据块的数据类型，oldtypes(1)为第二个数据块的数据类型，moletype为新定义的MPI数据结构。

网格文件第二次剖分后，所有模拟粒子需要从原有的网格体系映射到新的网格体系下，信息发送属于N对N操作，本发明基于自定义的粒子数据结构，利用MPI全局函数MPI_ALLTOALLV()实现了原有流场的成功映射，使得计算可以在定常结果的基础上继续进行。

本发明提供的一个或多个技术方案，至少具有如下技术效果或优点：

与不考虑计算负载平衡调整的MPI并行方案相比，本发明以网格单元内的模拟粒子数作为权重进行重分区操作，并完成流场的映射，使得计算负载在各计算进程上平衡分布，可大大提高并行的计算效率。

附图说明

此处所说明的附图用来提供对本发明实施例的进一步理解，构成本发明的一部分，并不构成对本发明实施例的限定；

图1为数值模拟计算中的负载均衡方法的流程示意图；

图2为重分区前计算负载的分布示意图；

图3为重分区后计算负载的分布示意图；

图4为重分区操作前的分区示意图；

图5为重分区操作后的分区示意图。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在相互不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述范围内的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。

实施例

请参考图1，图1为数值模拟计算中的负载均衡方法的流程示意图，本发明实施例提供了一种数值模拟计算中的负载均衡方法，所述方法包括：

在根进程内读入网格文件，对所述网格文件进行第一次剖分，获得若干个第一分区网格文件，将若干个所述第一分区网格文件分配至各处理器进程；

各处理器进程分别读入对应的第一分区网格文件，开始气体动力数值模拟并行计算；

判断流场内的气体流动是否达到定常，当气体流动达到定常时，统计各处理器进程内网格单元的全局编号以及网格单元内的模拟粒子数目；

在根进程内，将网格单元的全局编号与网格单元内粒子数目进行匹配，获得匹配结果；

在根进程内，基于上述匹配结果以网格单元内模拟粒子数目为权重对所述网格文件进行第二次剖分，获得若干个第二分区网格文件；

统计各处理器进程内的模拟粒子信息，将所述模拟粒子信息存入第一数组；

根据所述第二次剖分的分区结果，按照各处理器进程的编号信息对所述第一数组进行排序，获得排序结果；

基于所述排序结果，将流场内的模拟粒子信息映射到对应的所述第二分区网格文件中；

各处理器进程分别读入对应的第二分区网格文件，进行气体动力数值模拟并行计算。

本申请利用网格剖分METIS串行库实现MPI并行DSMC方法的基于模拟粒子数分布为权重的计算网格重新剖分，并根据剖分前后网格单元的全局编号不变这一关系完成流场信息的整体映射。具体步骤如下：

在根进程内读入网格文件，利用METIS完成无权重条件下的网格剖分（此时各分区内网格单元数基本相同），并将各分区网格分配至各处理器进程；其中，分区网格与处理器进程可以是一一对应关系，也可以是一个处理器进程对应多个分区网格；

各处理器进程分别读入上述步骤中生成的分区网格，开始DSMC并行计算（即开始气体动力数值模拟并行计算）；

判断流场内的气体流动是否达到定常，当流场内的气体流动达到定常时（一般认为，计算步数=5~10为气流横穿计算域全场所需步数，气流横穿计算域全场所需步数=计算域流向长度/（气体流向速度*时间步长）），开始进行计算负载平衡操作；

统计各处理器进程内网格单元的全局编号cellNoInProcArray以及网格单元内的模拟粒子数目nMoleInCellProcArray，并通过MPI全局函数MPI_Gatherv()将它们收集至主进程的相应数组中。在根进程内，按照网格单元的全局编号对单元内模拟粒子数信息进行排序；（统计各进程中的网格单元数nCellInProc，利用MPI中的MPI_Gather()函数将各进程的网格单元数信息收集至根进程中的nCellInProcArray数组中，在根进程内计算相对位移数组displs；网格单元内的模拟粒子数将作为权重信息传递给METIS分区程序，排序的作用是使得模拟粒子数与网格单元编号正确匹配，其中，MPI_Gatherv()程序需要用到nCellInProcArray和displs等数组信息）

在根进程内调用METIS分区程序，进行考虑以网格单元内模拟粒子数目为权重的重分区，即对原始网格文件进行重新分区；

统计各处理器进程内的模拟粒子信息，将其存入moleArray数组，模拟粒子信息包括模拟粒子的所在单元编号、模拟粒子所属组份编号、模拟粒子的空间三维坐标、模拟粒子在三维方向上的速度大小、模拟粒子的转动能、模拟粒子振动能等。根据新的分区关系，按照进程号从小到大的顺序，对moleArray数组进行排序。其中，模拟粒子所属的组分编号指的是氮气分子、氧气分子等，在程序中它们有一个相应的数字编号。

通过MPI全局函数MPI_ALLTOALLV()将流场内的模拟粒子信息映射到新的分区体系中。在信息传递过程中，为了方便，可以根据模拟粒子数据类型定义新的MPI数据类型。

各处理器进程分别读入新的分区网格，重新初始化部分数组（记录达到定常状态的流场数据的数组），在原有流场基础上进行续算。

本发明根据模拟粒子的数据结构定义了相应的MPI自定义数据结构。模拟粒子信息中既包含整型数据，也包含浮点型数据，且数据项较多。MPI的内置数据类型不能直接发送接收粒子信息。如果单独发送，则较为繁琐。本实施例根据MPI的相关语法，定义了模拟粒子的专属数据结构，使得可以直接发送整个粒子的所有信息。

模拟粒子的数据结构用fortran语言type自定义数据格式定义，具体如下：

procID，zoneID，cellID，specieID，XYZ(3)，UVW(3)，erot，evib；

其中，procID是模拟粒子将发往的进程编号（整型标量），zoneID是模拟粒子所在分区编号（整型标量），cellID是模拟粒子所在单元编号（整型标量），specieID是模拟粒子组分编号（整型标量），XYZ(3)是模拟粒子的三维空间坐标（浮点型数组），UVW(3)是模拟粒子的速度三分量（浮点型数组），erot是模拟粒子的转动能（浮点型标量），evib是模拟粒子的振动能（浮点型标量）。

定义MPI自定义数据类型moletype如下：

procID，zoneID，cellID，specieID，XYZ(3)，UVW(3)，erot，evib；

其中，procID是模拟粒子将发往的进程编号（整型标量），zoneID是模拟粒子所在分区编号（整型标量），cellID是模拟粒子所在单元编号（整型标量），specieID是模拟粒子组分编号（整型标量），XYZ(3)是模拟粒子的三维空间坐标（浮点型数组），UVW(3)是模拟粒子的速度三分量（浮点型数组），erot是模拟粒子的转动能（浮点型标量），evib是模拟粒子的振动能（浮点型标量）。

根据上述数据结构，可知数据类型有两块，分别是整型和浮点型，其中整型数据有4个，浮点型数据有8个。首先定义blocklengths(0) = 4，blocklengths(1) = 8，oldtypes(0) = MPI_INTEGER，oldtypes(1) = MPI_DOUBLE_PRECISION。接着调用MPI_GET_ADDRESS函数计算出偏移量offsets(0)和offsets(1)。最后调用MPI_TYPE_CREATE_STRUCT函数和MPI_TYPE_COMMIT函数即可完成moletype的定义和注册。

重新分区后，所有模拟粒子需要从原有的网格体系映射到新的网格体系下，信息发送属于N对N操作，本发明基于自定义的粒子数据结构，利用MPI全局函数MPI_ALLTOALLV()实现了原有流场的成功映射，使得计算可以在定常结果的基础上继续进行。

图1给出了负载平衡操作的流程图。接下来给出一个圆柱绕流算例的测试结果。来流Ma数=10，来流Knudsen数=0.25，来流为氩气。

图2给出了重分区前计算负载的分布，图3给出了重分区后计算负载的分布，其中，DLB操作为重分区操作，横坐标代表分区号，纵坐标代表该分区的粒子数，从图2-3中可以看出，经过重分区操作后，计算负载的不平衡性得到极大改善。每千步的运行时间从9m58s降到2m7s，耗时减少约79%。

图4给出了重分区操作前的分区示意图，图5给出了重分区操作后的分区示意图，其中，横坐标代表空间x坐标，纵坐标代表空间y坐标，从图4-5中可以看出，经过重分区操作后，高密度区域的分区数目明显增多。

与不考虑计算负载平衡调整的MPI并行方案相比，本发明以网格单元内的模拟粒子数作为权重进行重分区操作，并完成流场的映射，使得计算负载在各计算进程上平衡分布，可大大提高并行的计算效率。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。