CN104360896B - 一种基于gpu集群的并行流体仿真加速方法 - Google Patents

Abstract

本发明是一种基于GPU集群的并行流体仿真加速方法。针对流体仿真计算量大、并行度高的特点，设计了自动化的并行加速方法，进行流体的模拟与算法研究。提出并实现了一种在同一节点内多个GPU之间、多个节点之间的负载均衡算法。在给定空间中使用基于位置的流体仿真算法模拟流体的物理行为，通过减少分支以及缩小临域搜索范围对算法进行加速。对流体所在的空间进行分割，每个节点处理一个子空间，在节点内根据GPU数量进行进一步分割，从而完成流体仿真在GPU集群上的并行化。

Description

一种基于GPU集群的并行流体仿真加速方法

技术领域

本发明属于计算机应用领域，具体地说是一种基于GPU集群的并行流体仿真加速方法，该方法可用于影视、游戏、医学等各个领域中，模拟流体的行为。

背景技术

随着虚拟现实技术的发展，流体动画在电影特技、娱乐游戏、军事仿真、医学仿真、数字媒体中广泛应用，人们对真实感有了越来越高的需求。为了仿真这些复杂的场景效果，一帧一帧手工渲染动画自然满足不了需求，近年来，计算机硬件性能飞速发展，极大的促进了计算机图形学的深入研究。流体仿真更是多种多样，从早期简单的大面积水体仿真，如平静的水面、静止的烟雾，到现在渲染更为丰富细致的动态流体，如飞溅、泡沫、沸腾、多流体融合、多相流，流体交互仿真等，流体和其物理属性已经紧密的结合在一起。因此，基于物理的流体动画仿真成为计算机图形学和虚拟现实领域中的研究热点之一。

近年来被广泛应用的光滑粒子流体动力学(Smoothed Particle Hydrodynamics)算法是一种基于物理的流体仿真算法。该类方法使用粒子离散化问题域，通过加权求和邻居粒子的贡献值来计算粒子的物理量。该算法能够很好的处理自由表面、变形边界、运动交界面以及大变形等问题。但是SPH算法使用了刚性方程以维持其不可压缩性，这使得力的计算结果较大，只能使用足够小的时间步长。

基于位置的流体仿真(Position Based Fluid)算法在基于位置的动力学(Position Based Dynamics，PDB)框架中添加了密度的迭代求解。通过一组位置限制公式来维持一个固定的密度。该方法能够得到同SPH算法类似的不可压缩性以及收敛性，同时继承了PBD的几何稳定性，从而可以在实时仿真中使用更大的时间步。同时改善了SPH算法对时间步长的限制，允许使用更大的时间步长，在相同时间内可以获得更大的粒子位置更新，同时保留了基于粒子的流体仿真所固有的优势，具有广阔的应用前景。

为了使得仿真效果更加逼真、细节更丰富，粒子规模被不断提高。PBF算法的计算量非常大，CPU已经难以满足PBF方法的计算需求，浮点性能超过CPU几十倍的GPU开始被用于加速PBF方法。最近十余年，GPU在浮点性能和可编程性方面都取得了长足的进步。浮点性能方面，目前GPU无论在单精度浮点运算还是双精度浮点运算都达到了数Tflops。可编程性方面，随着CUDA的问世，GPU从传统的图形领域跨入了通用计算领域。越来越多的应用通过GPU加速获得了大幅度的性能提升。GPU是一种并行硬件架构，特别适合于加速计算密度大的并行算法。

发明内容

本发明的目的是在GPU集群上实现PBF的自动化加速算法，使得可以用同一种负载均衡策略、通信管理策略、存储管理策略来实现基于GPU集群的流体仿真算法。让开发人员能够快速在GPU集群上实现仿真算法，将研发精力从算法实现转移至对仿真算法本身的研究，只投入少量精力即可获得GPU集群的加速。

为了实现上述目的，对流体物理行为的模拟采用基于位置的流体仿真算法(PBF)实现，并通过减少分支，缩小临域搜索空间对算法进行优化。按照节点数量、节点内GPU数量和粒子分布情况对全局空间进行划分，将计算任务均衡的划分至各个GPU中，在一些联系性较强的步骤执行后对粒子的一些关键属性进行全局同步，使其能够正确的并行的快速的在GPU集群上完成整个计算任务。在计算过程中，根据各个GPU的计算时间，对全局的GPU进行负载均衡。

附图说明

图1示出本发明中并行加速方法所使用平台的整体逻辑结构；

图2示出本发明基于粒子分布情况与服务器节点数量划分空间；

图3示出本发明每个子空间周围的ghost粒子层；

图4示出本发明缩小临域搜索范围。

具体实施方式

下面结合附图对本发明作进一步说明。

并行加速方法所使用平台的整体逻辑结构如图1所示。具体步骤如下：

(1)按照服务器节点数量和粒子的空间分布情况对空间进行二维ORB划分，如图2所示。统计粒子在三个方向上坐标分布的方差，取方差最大的两个方向作为切割方向。首先沿方差最大的方向将场景切割为两个长方体，然后在两个长方体中继续选择其余一个方向切割为更小的两个长方体，循环在两个方向上切割直到子空间的数量和GPU数量相同；

(2)统计每个节点的所属子空间中粒子在ORB划分的两个方向上的方差，选定方差大的方向对该子空间进行一维划分，将整个仿真空间沿该方向切割为多个长方体子空间。将粒子数据发送至所属节点的相应GPU中，并包含该粒子群周围的一些冗余粒子信息用于后续计算，称之为ghost粒子层，如图3所示；

(3)根据粒子上一个时间步受力情况(模拟开始时初始化为重力)和速度，计算每个粒子的速度变化；

(4)将GPU的粒子空间划分为多个网格，如图4所示，网格尺寸为光滑半径的一半，以缩小搜索范围，统计粒子在空间网格中的分布情况，统计每一个网格的起始终止粒子，用于后续步骤的临域搜索，这里按(2)步骤中的划分方向合并多个网格，以减少分支；

(5)计算每个粒子的密度ρ。粒子密度由标准SPH密度估算公式得出，如粒子i的密度ρ_i是关于其所有临域粒子位置坐标的函数：

其中，m_j是临域粒子的质量，h是光滑半径，W是光滑核函数。

(6)计算每个粒子的限制公式C，粒子i的限制公式如下：

其中，ρ_i是粒子i的密度，ρ₀是静态密度。

(7)计算每个粒子的限制因子项，求粒子i的限制因子项λ_i的公式如下：

(8)全局同步限制因子项λ，每个GPU同相邻GPU交换相邻区域内的粒子的λ；

(9)计算位置偏移Δp，当一个粒子的临域粒子少于一定数量导致其密度无法达到静态密度时，会造成负向的压强，将出现粒子结块和聚集的现象，可以通过在光滑核函数中增加一个附加压力项S_corr来解决：

其中这里Δq是一个小于光滑半径的长度，n是一个常数，一般取值为4；

(10)根据粒子位置偏移更新粒子位置，并全局同步新的粒子位置，进行边界检测；

(11)对步骤(5)、(6)、(7)、(8)、(9)、(10)进行多次重复迭代，以校正步骤(3)中的粒子的预测位置，当到达最大迭代次数或位置偏移小于某一阈值后，进行下一步；

(12)根据校正后的粒子位置，与上一时间步中粒子的位置，计算粒子的速度，计算粒子的受力情况，包括涡流和粘滞力；

(13)对各节点各GPU进行负载均衡，将整个集群组织成二叉树结构，叶节点为服务器节点，非叶节点代表左右子树上的节点间的交界面。对于每一个非叶节点，我们根据左右子树节点在上一个时间步的平均耗时来调整它们的交界面。当左子树的节点的计算时间高于右侧节点时，交界面向左侧移动一段距离，一个网格层的粒子被交换到右侧节点上，反之则将右侧节点上的一个网格层中的粒子交换给左侧节点。逐层向上迭代重复此过程，直至根节点。节点内GPU的负载均衡方式同节点间相似，但因为是一维划分，故在不需要沿叶节点向上重复，只计算叶节点一层即可。

(14)在分布粒子或重新分布粒子N个时间步后，若全局数据传输量大于某一阈值M，则表明粒子目前的分布状态不适用于当前的空间切割方式，重新分布粒子，执行步骤(1)(2)，并计算重新分布后第一次全局同步的数据传输总量m，重新设置阈值M为a*m，其中a是一个大于1的常数。

(15)重复步骤(3)至(14)，直至整个模拟过程被结束。

Claims (1)

1.一种基于GPU集群的并行流体仿真加速方法，其特征在于包括以下步骤：

(1) 在全局空间中按照节点数量和粒子分布情况对其在某一二维平面进行划分，而非三维空间上进行ORB划分，每个节点对应一块子空间；

(2) 在每个节点的所属子空间中根据GPU数量和该子空间中粒子分布差异最大化的坐标轴方向对该子空间进行一维划分，并沿该方向对其子空间内粒子进行矩阵变换，将粒子数据发送至所属节点的相应GPU中；

(3) 各GPU使用基于位置的流体仿真算法(PBF) 进行计算，多次迭代修正粒子位置，再计算限制因子项，更新粒子位置之后分别对限制因子项和粒子位置进行两次全局同步，同步后本时间步结束；

(4) 根据各GPU的计算时间，对各个节点各个GPU进行负载均衡处理，调整全局空间的划分，根据新的划分在GPU之间交换粒子，而后重复步骤(3)，处理下一个时间步，

所述步骤(1)中在全局空间中按照节点数量和粒子分布情况对其在某一二维平面进行划分，具体方法是：

按照服务器节点数量和粒子的空间分布情况对空间进行二维ORB划分，ORB方式采用多层次的方式切割，选定两个方向进行划分；统计粒子在三个方向上坐标分布的方差，取方差最大的两个方向作为切割方向；首先沿方差最大的方向将场景切割为两个长方体，然后在两个长方体中继续选择其余一个方向切割为更小的两个长方体，循环在两个方向上切割直到子空间的数量和GPU数量相同；此切割方法目的在于减小节点间交界面面积，同时不至于使传输模型过于复杂，最小化各节点传输数据量；

所述步骤(2)中在每个节点的所属子空间中根据GPU数量和该子空间中粒子分布差异最大化的坐标轴方向对该子空间进行一维划分具体方法是：

计算粒子坐标集在步骤(2)选取两个方向上的方差，沿方差最大的坐标轴粒子的分布差异相对最大，按照粒子在该方向的数量分布，将整个仿真空间沿该方向切割为多个长方体子空间；此切割方法目的在于在使节点间传输事务次数最小化的前提下，节点内GPU间的传输数据量和传输事务次数最小化；

所述步骤(3)中各GPU使用基于位置的流体仿真算法(PBF) 进行计算并在计算限制因子项和校正粒子位置两个步骤中进行全局同步两个关键步骤，包括步骤如下：

(3.1) 根据粒子上一个时间步的受力情况和速度，计算每个粒子的速度变化；在模拟开始时，将受力情况初始化为重力；

(3.2) 根据粒子速度，预测每个粒子的下一步位置并更新；

(3.3) 将GPU的粒子空间划分为多个网格，网格尺寸为光滑半径的一半，以缩小搜索范围，统计粒子在空间网格中的分布情况，统计每一个网格的起始终止粒子，用于后续步骤的临域搜索，这里按(2)步骤中的划分方向合并多个网格，以减少分支；

(3.4) 计算每个粒子的限制因子项，并全局同步粒子的限制因子项；

(3.5) 根据粒子的临域粒子限制因子项，计算每个粒子的位置偏移，更新粒子位置并全局同步；

(3.6) 对每个粒子进行边界检测；

(3.7) 对步骤(3.4)、(3.5)、(3.6)进行多次重复迭代，以校正(3.2)中的粒子的预测位置，当到达最大迭代次数或位置偏移小于某一阈值后，进行下一步；

(3.8) 根据校正后的粒子位置，与上一时间步中粒子的位置，计算粒子的速度，计算粒子的受力情况，包括涡流和粘滞力；

(3.9) 在分布粒子或重新分布粒子N个时间步后，若全局数据传输量大于某一阈值M，则表明粒子目前的分布状态不适用于当前的空间切割方式，此时重新按照节点和GPU的数量划分整个空间，并计算重新分布后第一次全局同步的数据传输总量m，设置阈值M为a*m，其中a是一个大于1的常量，在实现中根据问题的具体情况调整大小；

(3.10) 重复步骤(3.1)至(3.9)，直至整个模拟过程被结束；

所述步骤(4)中各节点各GPU需要进行负载均衡，具体方法是：

将整个集群组织成二叉树结构，叶节点为服务器节点，非叶节点代表左右子树上的节点间的交界面；对于每一个非叶节点，根据左右子树节点在上一个时间步的平均耗时来调整它们的交界面；当左子树的节点的计算时间高于右侧节点时，交界面向左侧移动一段距离，一个网格层的粒子被交换到右侧节点上，反之则将右侧节点上的一个网格层中的粒子交换给左侧节点；逐层向上迭代重复此过程，直至根节点；

节点内GPU的负载均衡方式同节点间相似，但因为是一维划分，故不再需要沿叶节点向上重复，只计算叶节点一层即可。