CN112765871B

CN112765871B - 一种基于曲线坐标的并行粒子追踪方法和装置

Info

Publication number: CN112765871B
Application number: CN202110370636.5A
Authority: CN
Inventors: 关清帝; 梁剑寒; 张�林; 陈文武; 王玉奇
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2021-04-07
Filing date: 2021-04-07
Publication date: 2021-06-18
Anticipated expiration: 2041-04-07
Also published as: CN112765871A

Abstract

本申请涉及一种基于曲线坐标的并行粒子追踪方法、装置、计算机设备和存储介质。所述方法包括：建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为局部网格结构、全局网格非结构的多块网格，所述局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同。在局部网格内基于曲线坐标直接实现定位与追踪，在多个局部网格间基于块与块之间的拓扑信息实现粒子的坐标转换与通信，最终实现全局的粒子追踪。采用本方法能够实现块结构网格下相对简洁的数据结构和粒子追踪实现过程，同时也很好的处理了分区间复杂拓扑问题，数值实现简洁，计算效率高，十分适合拉格朗日方法的大规模并行计算。

Description

一种基于曲线坐标的并行粒子追踪方法和装置

技术领域

本申请涉及计算流体力学领域，特别是涉及一种基于曲线坐标的并行粒子追踪方法、装置、计算机设备和存储介质。

背景技术

在计算流体力学中，拉格朗日粒子方法是研究流体力学问题的一种重要手段，广泛应用在两相流、光滑粒子方法以及拉格朗日随机方法等问题的研究中。在欧拉网格中，通过拉格朗日粒子处理固（液）相是计算流体力学中常见的问题。拉格朗日粒子可以精确处理湍流扩散等空间输运问题，消除了网格方法引起的数值扩散问题。虽然粒子方法提供了许多潜在的优势，但这种方法也产生了许多需要解决的问题。其中比较核心的问题就是如何实现网格中粒子的定位和追踪。目前的网格-粒子方法中，根据采用的网格和粒子追踪策略，可以简单分为两类，一类是基于结构网格的粒子追踪方法，另一类是基于非结构网格的粒子追踪方法。

基于结构网格的搜索-定位过程复杂，需要重复多次的向量计算来判断粒子是否在当前网格中，而且在并行计算中，由于分区间的网格信息需要进行通信，因此结构网格块与块之间的搜索过程更加复杂且难以应用到复杂的网格拓扑中。非结构网格下的粒子追踪存在的问题是，需要设计基于面的数据结构，其在追踪过程中需要粒子到各个边（面）的垂直距离，数值实现相对复杂，而且每个粒子都要进行大量的定位计算，计算效率较低，非结构网格另一个问题就是计算精度也相对较低。传统的粒子追踪方法，不论是基于结构还是非结构网格的追踪方法，均涉及粒子物理坐标与网格物理坐标间复杂的向量计算，尤其在三维计算中，粒子追踪过程更加复杂耗时。

因此，现有的粒子追踪方法存在导致拉格朗日方法的并行效率较低的问题。

发明内容

基于此，有必要针对上述技术问题，提供一种提高拉格朗日方法并行计算效率的基于曲线坐标的并行粒子追踪方法、装置、计算机设备和存储介质。

一种基于曲线坐标的并行粒子追踪方法，所述方法包括：

建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为局部网格结构、全局网格非结构的多块网格，所述局部网格为结构化的曲线网格，每个所述局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同；

获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据所述第一曲线坐标得到待追踪粒子在所述第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据所述第一曲线坐标，以及所述第一局部网格和所述第二局部网格对应的曲线坐标系的转换关系，得到所述待追踪粒子在所述第二局部网格中的第二曲线坐标；

根据所述第二曲线坐标，得到所述待追踪粒子所在的网格索引。

在其中一个实施例中，还包括：在所述获取待追踪粒子在第一局部网格中的第一曲线坐标之后，获取追踪时间；

计算所述待追踪粒子到所述第一局部网格和所述第二局部网格的边界的距离，以及到达所述边界的穿越时间；

当所述穿越时间大于等于所述跟踪时间，待追踪粒子不会穿越所述第一局部网格的边界；

当所述穿越时间小于所述跟踪时间，待追踪粒子将穿越所述第一局部网格的边界。

在其中一个实施例中，还包括：当前跟踪步结束时，对未穿越边界的待追踪粒子，根据待追踪粒子的第一曲线坐标、追踪时间和待追踪粒子在曲线坐标系下的逆变速度，得到所述待追踪粒子在第一局部网格中的第二曲线坐标。

在其中一个实施例中，还包括：获取所述第一局部网格和所述第二局部网格对应的曲线坐标系的坐标转换矩阵；

根据所述跟踪时间和所述穿越时间，得到剩余跟踪时间，并确定所述待追踪粒子在所述第一局部网格边界的第一边界坐标和第一穿越速度；

根据所述边界坐标、所述穿越速度和所述坐标转换矩阵，得到所述待追踪粒子在所述第二局部网格边界的第二边界坐标和第二穿越速度；

根据所述剩余跟踪时间、所述第二边界坐标和所述第二穿越速度，在所述第二局部网格中进行粒子跟踪。

在其中一个实施例中，还包括：建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为局部网格结构、全局网格非结构的多块网格，所述局部网格为结构化的曲线网格，每个所述局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同，所述全局网格由所述局部网格相互链接而成。

在其中一个实施例中，还包括：在根据所述边界坐标、所述穿越速度和所述坐标转换矩阵，得到所述待追踪粒子在所述第二局部网格边界的第二边界坐标和第二穿越速度之后，通过MPI非阻塞式的点对点通信将所述第二边界坐标和所述第二穿越速度发送到所述第二局部网格对应的目标计算分区；

通过所述目标计算分区接收所述第二边界坐标和所述第二穿越速度。

在其中一个实施例中，还包括：根据所述第二曲线坐标，得到待追踪粒子所在的网格索引为：

其中，

为所述网格索引，

为所述第二曲线坐标，int表示向下取整操作。

一种基于曲线坐标的并行粒子追踪装置，所述装置包括：

网格初始化模块，用于建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为局部网格结构、全局网格非结构的多块网格，所述局部网格为结构化的曲线网格，每个所述局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同；

第二曲线坐标确定模块，用于获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据所述第一曲线坐标得到待追踪粒子在第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据所述第一曲线坐标，以及所述第一局部网格和所述第二局部网格对应的曲线坐标系的转换关系，得到待追踪粒子在所述第二局部网格中的第二曲线坐标；

网格索引确定模块，用于根据所述第二曲线坐标，得到待追踪粒子所在的网格索引。

一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：

上述基于曲线坐标的并行粒子追踪方法、装置、计算机设备和存储介质，通过建立进行粒子追踪的流体模型，初始化流体模型的网格为局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同；获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据第一曲线坐标得到待追踪粒子在第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据第一曲线坐标，以及第一局部网格和第二局部网格对应的曲线坐标系的转换关系，得到待追踪粒子在第二局部网格中的第二曲线坐标；根据第二曲线坐标，得到待追踪粒子所在的网格索引。本发明数值实现保持了块结构网格下相对简洁的数据结构和粒子追踪实现过程，同时也很好的处理了分区间复杂拓扑问题，数值实现简洁，计算效率高，适合拉格朗日方法的大规模并行计算。

附图说明

图1为一个实施例中基于曲线坐标的并行粒子追踪方法的应用场景图；

图2为一个实施例中基于曲线坐标的并行粒子追踪方法的流程示意图；

图3为一个实施例中局部网格结构、全局网格非结构的多块网格的示意图；

图4为一个实施例中曲线坐标下粒子位移的示意图，其中a为一个实施例中曲线坐标下粒子位移的示意图，b为一个实施例中物理坐标下粒子位移的示意图；

图5为一个实施例中局部网格间粒子追踪方法的示意图；

图6为另一个实施例中基于曲线坐标的并行粒子追踪方法的流程示意图；

图7为时间发展混合层实施例中流动构型示意图和采用的网格示意图，其中a为时间发展混合层实施例中流动构型示意图，b为时间发展混合层实施例中采用的网格示意图；

图8为时间发展混合层实施例中粒子分布示意图，其中a为时间发展混合层实施例中非正交曲线网格下粒子方法计算粒子分布结果示意图，b为时间发展混合层实施例中直接数值模拟计算粒子分布结果示意图；

图9为圆柱绕流实施例中的网格示意图；

图10为圆柱绕流实施例中粒子分布图；

图11为一个实施例中基于曲线坐标的并行粒子追踪装置的结构框图；

图12为一个实施例中计算机设备的内部结构图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本申请提供的基于曲线坐标的并行粒子追踪方法，可以应用于如图1所示的应用环境中。其中，数字标号1表示第一局部网格，数字标号2表示第二局部网格，点A和点B的连线构成一部分网格边界。第一局部网格和第二局部网格为结构化的曲线网格，第一局部网格和第二局部网格对应的曲线坐标系不同，多个第一局部网格和第二局部网格构成全局网格，全局网格为非结构的多块网格，局部网格之间的分布是无序的。

在一个实施例中，如图2所示，提供了一种基于曲线坐标的并行粒子追踪方法，以该方法应用于图1中的第一局部网格和第二局部网格为例进行说明，包括以下步骤：

步骤202，建立进行粒子追踪的流体模型，初始化流体模型的网格；网格为局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同。

在计算流体力学中，拉格朗日粒子方法是研究流体力学问题的一种重要手段，在欧拉网格中，通过拉格朗日粒子处理固（液）相是计算流体力学中常见的问题，大量的粒子会导致拉格朗日分析需要大量内存，同时，对每一个粒子的连续追踪和输运会增加计算要求，特别是对于非结构化网格，繁琐耗时的粒子追踪导致拉格朗日方法的并行效率较低，大规模计算应用受到明显的限制。

本发明针对拉格朗日方法中的粒子追踪问题，如图3所示，构建的局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同。通过将每个局部网格的粒子追踪对应到一个计算分区进行处理，可以实现利用大规模并行计算的优势提升全局粒子追踪的效率。

曲线坐标下的追踪方法，是通过坐标变换，将传统在物理坐标下进行粒子追踪的过程，转换到曲线坐标下进行。传统方法基于粒子物理坐标，即直接通过物理坐标以及物理位移

进行粒子追踪，

其中

为粒子基于物理坐标的物理速度。

本发明则在曲线坐标下，根据曲线坐标下的位移

来进行粒子的追踪与定位，

其中

，

，

，

等为度量系数，其确定了坐标系间的映射关系，给定网格后，度量系数即被确定，不需要额外重复计算。

称为逆变速度，是基于曲线坐标的曲线速度。图4a为曲线坐标下粒子位移的示意图，图4b为物理坐标下粒子位移的示意图，从图4可以看出，粒子P从A输运到B，不管是从物理坐标，还是曲线坐标，粒子的位置是一致的，或者说就是重合的，只是从两个不同的坐标系下去看粒子的位移而已。由于粒子直接在计算空间进行追踪，因此在当前分区网格中，由于粒子的计算坐标

与网格单元索引

直接相关，计算坐标向下取整即为网格索引，即可直接确定粒子所在的网格，不需要额外计算，自动实现分区内粒子的高效定位与追踪。

步骤204，获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据第一曲线坐标得到待追踪粒子在第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据第一曲线坐标，以及第一局部网格和第二局部网格对应的曲线坐标系的转换关系，得到待追踪粒子在第二局部网格中的第二曲线坐标。

如图5所示，假设粒子初始位于给定第一局部网格中的A点，并向新的位置B点运动。如果B点跟A点位于同一分区，则B未穿越当前分区，如果B点在第二局部网格，则B穿越出了当前分区。

粒子是否穿越边界的判断，可以通过给定时间步长，按第一局部网格下的第一曲线坐标系计算B点坐标，如图5中，若粒子在B点的

和

坐标大于5或者小于1，即表明该粒子飞出了第一局部网格。也可以比较粒子到达边界的时间与给定的跟踪时间的大小，若在耗尽跟踪时间之前，粒子已经到达边界，即粒子到达边界的时间小于跟踪时间，那么，当跟踪时间耗尽时，粒子会飞出边界。

当前跟踪步结束时，对未穿越边界的待追踪粒子，粒子的网格位置索引可以由曲线坐标系下粒子坐标值直接得到；对穿越至第二局部网格的待追踪粒子，采用类似非结构网格追踪的方法，首先将粒子退回到边界处，即先确定粒子穿越网格的位置O，并在O位置处中断粒子追踪，计算剩余的输运时间

。在曲线坐标下，剩余的输运时间计算也十分简单，因为曲线坐标是正交的，粒子初始位置A到边界的最短距离由边界的坐标减去A的坐标即可直接得到，粒子逆变速度已知，从A到O的时间即可得到，总的跟踪时间减去

即可得到剩余的时间步长

。然后在边界O处，根据边界两侧的第一局部网格和第二局部网格对应计算分区的拓扑关系，进行曲线坐标间的坐标转换，即将当前分区的曲线坐标转换为目标分区的曲线坐标，如图5所示O点在分区1的曲线坐标为(5,2.6)，坐标转换的任务就是将其转变为分区2中的(1,2.6)。除了位置坐标需要转换外，粒子的逆变速度也要进行相应的转换。位置矢量和逆变速度矢量的坐标转换在数值上只需要乘以坐标转换矩阵即可完成。第二局部网格对应的计算分区接收该粒子，并从O处开始继续对该粒子进行追踪，根据剩余跟踪时间，转换后的粒子坐标和粒子速度，计算得到粒子在第二局部网格中的第二曲线坐标，重复前面的过程，直到剩余的时间步长变为0。

步骤206，根据第二曲线坐标，得到待追踪粒子所在的网格索引。

根据粒子的计算坐标

，计算坐标向下取整即为网格索引，可直接确定粒子所在的网格。

上述基于曲线坐标的并行粒子追踪方法中，通过建立进行粒子追踪的流体模型，初始化流体模型的网格为局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同；获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据第一曲线坐标得到待追踪粒子在第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据第一曲线坐标，以及第一局部网格和第二局部网格对应的曲线坐标系的转换关系，得到待追踪粒子在第二局部网格中的第二曲线坐标；根据第二曲线坐标，得到待追踪粒子所在的网格索引。本发明数值实现保持了块结构网格下相对简洁的数据结构和粒子追踪实现过程，同时也很好的处理了分区间复杂拓扑问题，数值实现简洁，计算效率高，适合拉格朗日方法的大规模并行计算。

在其中一个实施例中，还包括：在获取待追踪粒子在第一局部网格中的第一曲线坐标之后，获取追踪时间；计算待追踪粒子到第一局部网格和第二局部网格的边界的距离，以及到达边界的穿越时间；当穿越时间大于等于跟踪时间，待追踪粒子不会穿越第一局部网格的边界；当穿越时间小于跟踪时间，待追踪粒子将穿越第一局部网格的边界。

本发明判断粒子是否飞出当前分区的方法十分简单，不额外增加计算量。

在其中一个实施例中，还包括：当前跟踪步结束时，对未穿越边界的待追踪粒子，根据待追踪粒子的第一曲线坐标、追踪时间和待追踪粒子在曲线坐标系下的逆变速度，得到待追踪粒子在第一局部网格中的第二曲线坐标。

逆变速度，是基于曲线坐标的曲线速度。在单块网格中进行粒子追踪十分简单，根据待追踪粒子的第一曲线坐标、追踪时间和待追踪粒子在曲线坐标系下的逆变速度，可以直接得到待追踪粒子在第一局部网格中的第二曲线坐标。在并行计算中，绝大多数的粒子并不穿越内部边界，基于曲线坐标的方法降低了相应的计算代价。

在其中一个实施例中，还包括：获取第一局部网格和第二局部网格对应的曲线坐标系的坐标转换矩阵；根据跟踪时间和穿越时间，得到剩余跟踪时间，并确定待追踪粒子在第一局部网格边界的第一边界坐标和第一穿越速度；根据边界坐标、穿越速度和坐标转换矩阵，得到待追踪粒子在第二局部网格边界的第二边界坐标和第二穿越速度；根据剩余跟踪时间、第二边界坐标和第二穿越速度，在第二局部网格中进行粒子跟踪。

通过多区网格间基于块与块之间的拓扑信息实现粒子的坐标转换与通信，可以实现对复杂网格拓扑间的粒子连续追踪，避免了边界处丢失粒子的问题。局部网格对应的分区间粒子局部坐标转换矩阵在粒子初始化前即可生成，粒子追踪计算中，不需要更新。坐标转换的计算量相较传统方法在边界处的追踪和处理来说，也没有增加额外的计算量。

在其中一个实施例中，还包括：建立进行粒子追踪的流体模型，初始化流体模型的网格；网格为局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同，所述全局网格由所述局部网格相互链接而成。

在其中一个实施例中，还包括：在根据边界坐标、穿越速度和坐标转换矩阵，得到待追踪粒子在第二局部网格边界的第二边界坐标和第二穿越速度之后，通过MPI非阻塞式的点对点通信将第二边界坐标和第二穿越速度发送到第二局部网格对应的目标计算分区；通过目标计算分区接收第二边界坐标和第二穿越速度。

在其中一个实施例中，还包括：根据第二曲线坐标，得到待追踪粒子所在的网格索引为：

其中，

为网格索引，

为第二曲线坐标，int表示向下取整操作。

在一个具体实施例中，如图6所示，提供了一种基于曲线坐标的并行粒子追踪方法，包括以下步骤：

1）初始化网格，计算曲线坐标下的网格导数，以及网格分区间的坐标转换矩阵；

2）初始化粒子位置和粒子追踪时间，粒子位置为曲线坐标空间的位置，当前分区内所有粒子的追踪时间为

，

为粒子输运的总时间步长；

3）计算粒子从当前位置到分区各边界的距离

，并根据粒子逆变速度找出粒子将最可能穿越的边界以及到该边界的时间

；

4）如果

，则粒子的最终位置将仍在当前网分区，直接更新粒子坐标并结束当前步的追踪；

5）如果

，说明粒子飞出了当前分区边界。首先将粒子移动到粒子穿越点，并将相应缩减粒子追踪时间为

，此时

为剩余的粒子追踪时间；

6）将该粒子移动到该穿越边界的粒子链表中，并对所有该边界上的粒子进行分区间的坐标转换，将当前粒子的曲线坐标转换为相邻目标分区的曲线坐标，粒子的逆变速度转换为目标分区下的速度；

7）对应发送链表中的粒子，向目标分区进行MPI发送；

8）接收目标分区发送的粒子，并对所有粒子追踪时间大于零的粒子，重复上述步骤3到步骤8的操作，直到粒子发送链表为空，且所有粒子的追踪时间都变为0。

9）完成当前步所有粒子的追踪过程。

10）完成第n个时间步后，继续进行第n+1个时间步的粒子跟踪，直到完成预先设定的

次迭代。

在一个时间混合层实施例中，随时间发展的混合层流动示意图如图7a所示，由两股相反的平行流发展而成。分别采用粒子方法以及直接数值模拟方法对时间发展可压缩平面混合流进行模拟，对比验证曲线坐标下粒子追踪方法的准确性。采用的计算网格如图7b所示，为非正交的曲线网格。图8a和图8b分别为粒子方法计算粒子分布和直接数值模拟计算粒子分布的结果示意图，从图8可以看出，计算的粒子分布与直接数值模拟符合的很好，验证了本发明提出的曲线坐标下粒子追踪方法的准确性。

在另一个圆柱绕流实施例中，综合体现了多块网格间实现粒子追踪的主要问题。图9给出了圆柱绕流的网格示意图。该算例总共包含12个网格分区，每个分区网格数目相同。从图中可以看出分区1和分区2等几个块之间存在明显的局部坐标方向变化的情况。入口给定超声速入口条件，来流马赫数选择为2.5，因为在该马赫数下，激波正好能穿越多个网格分区，能更好的验证曲线坐标下粒子是否被正确输运。来流压力为101325Pa，来流温度为288.15K。上下以及出口边界均采用超声速出口条件。此外，用标量来对粒子进行标记，以显示化粒子的位置。从图10可以看出粒子能够正确穿越不同的网格分区，表明了曲线坐标下粒子追踪方法的正确性，以及其处理复杂网格拓扑的能力。

在另一个具体实施例中，粒子在边界处还可能被反射或者直接逃逸。具体的边界条件有：（1）周期条件：粒子穿越边界，并在另一周期处穿入计算域；（2）出口条件：粒子穿越边界，逃离计算域，释放粒子；（3）入口条件：生成粒子，并进入计算域；（4）反射条件：粒子在边界处反射。在曲线坐标下，由于边界与坐标曲面一致，边界条件的处理时十分简单的。对壁面边界，采用镜面反射，直接速度取反，重新进行追踪，直到消耗完所有时间步长；对于周期边界，粒子速度保持不变，发送到新的周期边界，消耗完剩余时间；对于出口边界，直接释放粒子。

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

在一个实施例中，如图11所示，提供了一种基于曲线坐标的并行粒子追踪装置，包括：网格初始化模块1102、第二曲线坐标确定模块1104和网格索引确定模块1106，其中：

网格初始化模块1102，用于建立进行粒子追踪的流体模型，初始化流体模型的网格；网格为局部网格结构、全局网格非结构的多块网格，局部网格为结构化的曲线网格，每个局部网格对应一个计算分区，每个曲线网格对应的曲线坐标系不同；

第二曲线坐标确定模块1104，用于获取待追踪粒子在第一局部网格中的第一曲线坐标，当前跟踪步结束时，对未穿越边界的待追踪粒子，根据第一曲线坐标得到待追踪粒子在第一局部网格中的第二曲线坐标，对穿越至第二局部网格的待追踪粒子，根据第一曲线坐标，以及第一局部网格和第二局部网格对应的曲线坐标系的转换关系，得到待追踪粒子在第二局部网格中的第二曲线坐标；

网格索引确定模块1106，用于根据第二曲线坐标，得到待追踪粒子所在的网格索引。

第二曲线坐标确定模块1104还用于在获取待追踪粒子在第一局部网格中的第一曲线坐标之后，获取追踪时间；计算待追踪粒子到第一局部网格和第二局部网格的边界的距离，以及到达边界的穿越时间；当穿越时间大于等于跟踪时间，待追踪粒子不会穿越第一局部网格的边界；当穿越时间小于跟踪时间，待追踪粒子将穿越第一局部网格的边界。

第二曲线坐标确定模块1104还用于在当前跟踪步结束时，对未穿越边界的待追踪粒子，根据待追踪粒子的第一曲线坐标、追踪时间和待追踪粒子在曲线坐标系下的逆变速度，得到待追踪粒子在第一局部网格中的第二曲线坐标。

第二曲线坐标确定模块1104还用于获取第一局部网格和第二局部网格对应的曲线坐标系的坐标转换矩阵；根据跟踪时间和穿越时间，得到剩余跟踪时间，并确定待追踪粒子在第一局部网格边界的第一边界坐标和第一穿越速度；根据边界坐标、穿越速度和坐标转换矩阵，得到待追踪粒子在第二局部网格边界的第二边界坐标和第二穿越速度；根据剩余跟踪时间、第二边界坐标和第二穿越速度，在第二局部网格中进行粒子跟踪。

第二曲线坐标确定模块1104还用于在根据边界坐标、穿越速度和坐标转换矩阵，得到待追踪粒子在第二局部网格边界的第二边界坐标和第二穿越速度之后，通过MPI非阻塞式的点对点通信将第二边界坐标和第二穿越速度发送到第二局部网格对应的目标计算分区；通过目标计算分区接收第二边界坐标和第二穿越速度。

网格索引确定模块1106还用于根据第二曲线坐标，得到待追踪粒子所在的网格索引为：

其中，

为网格索引，

为第二曲线坐标，int表示向下取整操作。

关于基于曲线坐标的并行粒子追踪装置的具体限定可以参见上文中对于基于曲线坐标的并行粒子追踪方法的限定，在此不再赘述。上述基于曲线坐标的并行粒子追踪装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图12所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种基于曲线坐标的并行粒子追踪方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本领域技术人员可以理解，图12中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims

1.一种基于曲线坐标的并行粒子追踪方法，其特征在于，所述方法包括：

建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为多块结构网格，所述多块结构网格中的每块网格为一个局部网格，所述局部网格为结构化的曲线网格，每个所述局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同；

获取待追踪粒子在第一局部网格中的第一曲线坐标，获取跟踪时间；计算所述待追踪粒子到所述第一局部网格和第二局部网格的边界的距离，以及到达所述边界的穿越时间；当所述穿越时间大于等于所述跟踪时间，待追踪粒子不会穿越所述第一局部网格的边界；当所述穿越时间小于所述跟踪时间，待追踪粒子将穿越所述第一局部网格的边界；

当前跟踪步结束时，对未穿越边界的待追踪粒子，根据待追踪粒子的第一曲线坐标、跟踪时间和待追踪粒子在曲线坐标系下的逆变速度，得到所述待追踪粒子在第一局部网格中的第二曲线坐标；

对穿越至第二局部网格的待追踪粒子，获取所述第一局部网格和所述第二局部网格对应的曲线坐标系的坐标转换矩阵；根据所述跟踪时间和所述穿越时间，得到剩余跟踪时间，并确定所述待追踪粒子在所述第一局部网格边界的第一边界坐标和第一穿越速度；根据所述边界坐标、所述穿越速度和所述坐标转换矩阵，得到所述待追踪粒子在所述第二局部网格边界的第二边界坐标和第二穿越速度；根据所述剩余跟踪时间、所述第二边界坐标和所述第二穿越速度，在所述第二局部网格中进行粒子跟踪；

2.根据权利要求1所述的方法，其特征在于，所述局部网格之间的拓扑为非结构的拓扑。

3.根据权利要求2所述的方法，其特征在于，在根据所述边界坐标、所述穿越速度和所述坐标转换矩阵，得到所述待追踪粒子在所述第二局部网格边界的第二边界坐标和第二穿越速度之后，包括：

通过MPI非阻塞式的点对点通信将所述第二边界坐标和所述第二穿越速度发送到所述第二局部网格对应的目标计算分区；

4.根据权利要求3所述的方法，其特征在于，根据所述第二曲线坐标，得到待追踪粒子所在的网格索引，包括：

根据所述第二曲线坐标，得到待追踪粒子所在的网格索引为：

其中，

为所述网格索引，

为所述第二曲线坐标，

表示向下取整操作。

5.一种基于曲线坐标的并行粒子追踪装置，其特征在于，所述装置包括：

网格初始化模块，用于建立进行粒子追踪的流体模型，初始化所述流体模型的网格；所述网格为多块结构网格，所述多块结构网格中的每块网格为一个局部网格，所述局部网格为结构化的曲线网格，每个所述局部网格对应一个计算分区，每个所述曲线网格对应的曲线坐标系不同；

第二曲线坐标确定模块，用于获取待追踪粒子在第一局部网格中的第一曲线坐标，获取跟踪时间；计算所述待追踪粒子到所述第一局部网格和第二局部网格的边界的距离，以及到达所述边界的穿越时间；当所述穿越时间大于等于所述跟踪时间，待追踪粒子不会穿越所述第一局部网格的边界；当所述穿越时间小于所述跟踪时间，待追踪粒子将穿越所述第一局部网格的边界；当前跟踪步结束时，对未穿越边界的待追踪粒子，根据待追踪粒子的第一曲线坐标、跟踪时间和待追踪粒子在曲线坐标系下的逆变速度，得到所述待追踪粒子在第一局部网格中的第二曲线坐标；对穿越至第二局部网格的待追踪粒子，获取所述第一局部网格和所述第二局部网格对应的曲线坐标系的坐标转换矩阵；根据所述跟踪时间和所述穿越时间，得到剩余跟踪时间，并确定所述待追踪粒子在所述第一局部网格边界的第一边界坐标和第一穿越速度；根据所述边界坐标、所述穿越速度和所述坐标转换矩阵，得到所述待追踪粒子在所述第二局部网格边界的第二边界坐标和第二穿越速度；根据所述剩余跟踪时间、所述第二边界坐标和所述第二穿越速度，在所述第二局部网格中进行粒子跟踪；

6.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至4中任一项所述方法的步骤。

7.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至4中任一项所述的方法的步骤。