CN106484532B - 面向sph流体模拟的gpgpu并行计算方法 - Google Patents

Abstract

本发明公开了一种面向SPH流体模拟的GPGPU并行计算方法。该方法将SPH方法中粒子的计算任务进行了分组来迎合GPGPU中的线程层次结构，同一个协作线程队列（CTA）中的线程能够处理逻辑上相似的工作。这个分组方法针对了GPGPU中的CTA调度规则开发并在GPGPU上实现，从而能快速分配任务的同时增加系统的cache命中率。以此分配方法为基础，片内共享内存能够被用来缓存来自全局内存中读入的数据。并且本发明能够根据全局内存的存储访问特性充分存储器带宽。同时，通过对GPGPU的SIMD特性的利用，由于使用共享内存所带来的线程同步开销也能被有效规避。

Description

面向SPH流体模拟的GPGPU并行计算方法

技术领域

本发明属于计算机图形学与高性能计算领域，具体地说是一个新的基于通用计算图形处理器(GPGPU)并行方法，用于以SPH算法为基础的流体模拟计算。涉及到计算机图形学以及仿真，计算流体动力学，单指令流多数据流(SIMD)架构，GPGPU等。

背景技术

SPH方法是一种无网格拉格朗日数值方法。该方法在近年来被广泛运用于流体模拟、天体动力学、分子动力学以及生物模拟中。作为一个粒子方法，SPH需要大量的粒子来模拟高精度大规模的场景。而大量的高密度粒子会导致海量的计算需求。所以对于用户来说，他们总是被迫于在模拟效果和系统性能之间做出权衡。因此，提升SPH方法的性能总是非常有意义的。

如今，非常多的研究致力于提升SPH方法的性能。其中一个主题是从系统体系结构的角度来做性能的提升。由于SPH方法本身非常适合并行化，非常多的研究得益于近年来并行硬件以及开发环境的发展并成功地将SPH方法和它们结合起来。然而，大多数这些研究，尤其是那些利用GPU的研究工作，没有能够充分地利用这些新硬件和开发环境的能力。比如说，主流的传统SPH的GPGPU实现中，尽管这种实现非常易于实现，并且相比起CPU代码来说，有1到2个数量级的性能提升，但是GPU在这种实现中只是被单纯地看作是一个SIMT处理器，而没有能够运用到很多GPU中的高级特性。而这种做法导致了传统的实现让SPH成为了内存受限型的应用。于是，也有一些研究试图去利用GPU中的各种内存空间来改进其在内存访问方面的高额代价，特别是共享内存。然而，已有的共享内存解决方案并没有带来理想中的性能提升，甚至要比传统方法慢上20％～70％。这个解决方案使用了显式的CTA中线程同步的方法来保证数据的正确性，但是线程同步带来的非常高昂的性能开销；同时该方法过多地占用共享内存使得其在大多数设备上运行时有着较低的设备占用率；其次，其低效的全局内存访问模式没有很好地利用全局内存的访问特性；而其随机的任务分配策略会导致较低的cache命中率。

总而言之，尽管利用GPGPU设备的高并行处理能力，SPH的性能得到了显著的提升，但是其潜力仍然远远没有被充分挖掘。

发明内容

本发明的目的是提供一个新的GPU并行计算架构来加速以SPH方法为基础的流体模拟在GPGPU设备上的性能。该方法包含了一个新的面向SPH流体模拟的GPGPU并行计算方法。

本发明的目的是这样实现的：

一种面向SPH流体模拟的GPGPU并行计算方法，该方法包括以下具体步骤：

步骤1：设SPH方法中的光滑核半径为h，将整个模拟空间划分为以h为边长的互不相交的正方体单元格；

步骤2：使用基于GPGPU的计数排序算法，将空间内的粒子根据它们所在的单元格来进行排序，在计数排序算法中，首先计算每个单元格中所包含的粒子数量，然后通过并行前缀求和算法计算每个单元格中的粒子在内存中的偏移量，之后根据每个粒子所在单元格在内存中的偏移量以及其在该单元格中的位置进行复制；

步骤3：根据每个单元格中粒子的数量，为接下来的粒子交互阶段中每一个CTA分配任务，此过程完全在GPGPU上实现：

i)计算每个单元格所包含的粒子在内存中的范围，其结果保存在两个数组中，数组中每个元素描述一个单元格的相关信息：(1)单元格中的粒子在内存中的偏移量；(2)单元格所包含的粒子数量B；这步由个GPU线程并行完成，为空间中所包含的单元格的数量；

ii)计算每一个单元格将在步骤4中所需要用来处理其内部粒子的交互的CTA数量，并将结果保存在数组B’中；所需的CTA数量为其中，S_cta表示一个CTA中包含了线程数量，表示单元格c_j中所包含的粒子的数量，这些数据保存在数组B中；这步由个GPU线程并行完成；

iii)对数组B’中的元素做前缀求和并将结果保存在数组O’中；O’中的元素表示每个单元格所需的任务在任务队列T中的起始位置；这步由个GPU线程并行完成；

iv)T中每个元素t_i＝{c_i，k_i}，其中c_i表示当前这个任务需要处理哪个单元格中的粒子，而k_i代表自己是第几个处理该单元格的任务；这步由个GPU线程并行完成；

v)获取总共需要的CTA数量S_grid；该值通过将O’中最后一个值与最后一个非空单元格所需的CTA数量相加得到；

步骤4：进行粒子交互；粒子交互包括两个独立的遍历邻居粒子阶段：读取邻居粒子的位置与质量来计算每个离散粒子表示的液体密度；通过邻居粒子的密度信息来计算每个离散粒子的受力；其中，遍历邻居粒子的方法如下(两个阶段相似，因此不重复阐述)：

i)每个线程将自己所负责的粒子的必要数据从全局内存中读取到寄存器中；如果该线程不负责任何粒子，则跳过此步骤；

ii)一个CTA中前9个线程计算它们所负责的那个单元格周围的单元格的粒子在内存中的保存范围，结果保存在共享内存中；周围的单元格的粒子在内存中的范围是9个连续区域即3D场景下；其余线程跳过此步骤；

iii)读取一个邻居粒子区域中最多S_cta个粒子，保存到共享内存中；CTA中每个线程最多负责读取一个；

iv)凡是有粒子负责的线程通过在寄存器中自己负责的粒子数据以及共享内存中最多S_cta个邻居粒子的数据来计算自己负责的粒子的密度或是受力；

v)S_cta设为32，也就是一个线程束中所包含的线程的数量；因此内存访问与计算之间不需要显式的线程同步；

步骤5：根据步骤4中计算出的每个粒子的受力情况，来移动粒子；移动粒子的时候要考虑外力影响(如重力)与边界情况；当粒子在边界上时给予一个反向的力将粒子弹回模拟空间内；

步骤6：一个时间步中所需要的计算就完成了；之后，输出粒子的位置信息，然后计算下一个时间步中粒子的变化，回到步骤1。

本发明的有益效果

本发明所提出的方法相对于经典的SPH方法的GPGPU实现有60％的性能提升。这个提升是建立在没有对数值方法进行任何修改，并且没有任何精度损失的情况下的。本发明的方法能够运用于其他基于SPH方法的流体模拟算法中并得到显著的性能提升。同时，本发明也适用于那些与SPH方法有着相同的底层算法(固定网格法)的其他算法。

附图说明

图1为本发明基于固定网格法的SPH流体模拟的流程图；

图2为本发明计算任务如何分配给CTA的示意图；

图3为本发明任务分配算法的流程图；

图4为本发明27个正方体邻接单元格分成9个长方体单元格示意图；

图5为本发明中，每次邻居粒子数据访问时内存读取的模式示意图；图中，每个方块代表一个邻居粒子，灰色方块为一次能读取的粒子。

具体实施方式

以下结合附图对本发明作详细说明。

(1)基于SPH方法的流体模型

SPH方法是一个无网格粒子方法。作为一个基于拉格朗日的方法，持续的模拟空间被视为由一些离散的粒子所组成的。每个粒子的属性(比如密度，受力)是通过光滑核函数W(r，h)带权插值计算得到的：

其中，m和ρ分别表示粒子所代表的质量和密度。W(r，h)中的h表示光滑核半径。

在计算每个粒子的受力之前，需要首先计算每个粒子所代表的密度ρ和产生压力p：

ρ_i＝∑_jm_jW(r_i-r_j，h)， (2)

p_i＝κ(ρ_i-ρ₀)， (3)

其中，κ是气体常数，ρ₀是静止密度。接着，通过解N-S方程来计算每个粒子的受力f：

其中，g是外力密度场，v是速度，μ是黏性系数，σ是张力系数。

于是，通过每个粒子的受力情况以及时间步长，可以计算出在当前时间步结束后粒子的位移了。

(2)查找固定半径内的邻居粒子

由上述公式中所用到的光滑核函数中知道：一个在r_i的粒子只会被那些位于r_j位置的粒子影响，如果|r_i-r_j|≤h的话。这是因为只有这种情况下光滑核函数W(r，h)才不为0。因此，SPH方法的实现就能映射为固定半径内邻居(FNN)问题。

FNN问题是处理在多维度欧几里德空间中，所有在给定距离内的点与点之间的交互的问题。而固定网格法是目前最为流行的用来解决这个问题的方法之一。这是由于这个方法非常适合并行执行并且很容易地运用到各种数值方法中去。在这种方法中，整个空间被分割为边长为h的互不相交的正方形单元格，设这些单元格的集合为空间中的每个粒子必然被放到其中一个单元格中。在邻居粒子遍历阶段，每个粒子只需要遍历以下它所在单元格周围的单元格中的粒子就能找到所有离它在h以内的邻居粒子了。使用固定网格法需要两个步骤：建立邻居信息和遍历邻居数据。前者是预处理阶段，用排序算法来完成，如果使用基于非比较的排序算法的话，时间复杂度为O(wn)，而邻居遍历阶段的时间复杂度为O(3^kwnN)，其中w是排序用关键字的长度，n是空间中粒子的数量，k是空间纬度(在三维SPH流体模拟中是3)，N是每个单元格中的粒子数上限。

(3)基于CTA调度特性的任务分配算法

固定网格法已经被证明是在低维空间中最适合用来处理FNN问题的方法。本发明也是基于这个方法的。基于固定网格法的SPH流体模拟的流程图如图1所示，图中，黑色的块是本发明的重点。灰色块是关于建立邻居信息的模块。本发明是通过共享内存的使用来解决原先SPH的GPGPU实现中内存受限的问题。而共享内存在使用上由非常多的限制，它是被一个CTA共享的。因此，一个CTA中的线程所处理的粒子必须要有近似的内存访问需求，不然共享内存将很难顾及所有线程的访问需求，并且出现大量的程序分支导致的串行执行。另一方面，一个CTA中所包含的线程在整个GPU内核程序执行的过程中都是固定的，而每个单元格中的粒子数量可以是任意的。所以，本发明中需要一个任务分配的过程，让每个CTA知道它们所负责的单元格以及自己是第几个负责该单元格的CTA。图2为计算任务如何分配给CTA的示意图。图中，虚线表示那些没有计算需求的空闲线程；左边第二行的单元格所分配的任务被省略了。图2展示了一个被固定网格法所划分的空间，这个空间中有若干的粒子被插入在不同的单元格中。在这张图中，假设CTA的大小S_cta的大小是5，表示单元格c_j中所包含的粒子的集合。那么任务分配的输出需要满足如下规则：

单元格c_j被分配到的CTA数量是

如果那么必然有一些没有任何粒子需要计算的闲线程存在于被分配给单元格c_j的CTA中。

在任务分配的最后，需要计算出之后粒子数值计算的GPGPU内核程序所使用的整个CTA数量S_grid，而每个CTA所被赋予的任务则存储在数组T中。T可以表示为其中是表示该任务是第几个处理单元格c的任务。以图3为例，S_grid为12，t₉＝{θ，1}如果是使用行优先顺序的话。

在GPGPU平台中，最简单的实现上述功能的方法是使用原子操作。每个线程负责一个单元格，计算该单元格需要多少CTA来计算其内部的粒子，然后调用原子操作获取个随机但是全局唯一的值来保存这个单元格的任务在整个任务数组T中位置并把任务写进去。但是这种做法会导致非常低的cache命中率。这是由于以CUDA为代表的GPGPU平台的CTA调度器是根据巡回制来给每个多处理器安排任务的。换句话说，在线程网格中相邻的CTA更有可能同时被执行。如果同时被执行的CTA有近似的内存访问需求，那么可以得到更高的cache命中率。而基于原子操作的做法无法保证这一点。因此，本发明不但要并行快速地在GPU中实现这个任务分配的过程，同时T中相邻任务也要有关联性，也就是说，同一个单元格的任务要在T中相邻，而相邻单元格的任务也要尽可能相邻。

详细的流程如图3所示，此图使用了图2中场景为例子。首先，在建立好邻居信息后，很容易地得到了每个单元格中所包含的粒子的数量以及它们在内存中的起始索引。假设每个单元格所包含的粒子数量保存在数组B中，空间上相邻的单元格在B中也是邻近的(优先顺序为X轴＞Y轴＞Z轴)。然后个线程被用来计算每个单元格需要多少个CTA并把结果存放到数组B’中。之后对B’进行前缀求和得到数组O’。O’中的数据就是每个单元格的任务在T中的起始位置。然后个线程将任务写进T中。而总CTA数量S_cta可以很容易的计算到，只需要将O’中最后一个元素和B’中最后一个非0值相加即可。

(3)全新的粒子交互方法

这里主要是说明粒子如何来遍历其邻居粒子的过程。这部分充分利用了GPU的线程和内存层次结构的特性。

首先，每个线程从T中读取它们所在的CTA所负责的任务。然后，将每个线程所负责的粒子的数据从全局内存读入到寄存器中。接着，CTA中前9个线程负责读取它们所负责的单元格周围的单元格中的粒子在内存中的存储范围。在3D空间中，每个单元格有27个相邻单元格(包括它自己)，由于粒子在内存中的存储规则是按照其所在的单元格在空间中的位置来存放的，并且按照X轴优先的规则。因此，周围的27个正方体邻近单元格可以视为9个长方体的连续单元格，每个长方体单元格包含了Y轴和Z轴坐标相同而X轴不同的3个连续的正方体单元格，如图4所示。然后初始化以下读取参数，包括现在正在读取哪个邻接长方体单元格以及该单元格中已经读取了多少粒子了。

之后进入读取邻居数据并进行计算的循环。和传统的方法不同，本发明并不是每次读取一个邻居粒子然后做计算，而是读取最多S_cta个。首先，S_cta个线程将当前正在读取的长方体单元格中没有被读取的前S_cta个邻居粒子的数据从全局内存读到共享内存中。由于粒子在内存中存放的顺序是按照其所在的单元格的，因此这些邻居粒子的数据在内存中都是连续的，仅仅需要很少的内存访问次数就能得到大量的数据。这种做法还能使得每次内存访问的带宽有效占用率达到最高。因此，这种访问模式中，每个线程束所产生内存读取指令所导致的访存次数是最低，并且得到的有效数据是最高的，如图5所示。这种模式下，如果当前长方体单元格中的粒子数量达不到S_cta，将读取少于S_cta个数量的粒子，并且在读取完毕后将当前读取的长方体单元格指向下一个。如果全部的长方体单元格已经都读完了，那么就跳出循环，当前CTA的邻居访问算是完成了。读取邻居数据之后，就可以根据存放在共享内中最多S_cta个邻居粒子的数据与自己所负责的存放在寄存器中的粒子数据进行计算了。共享内存的性能和L1cache是近似的，因此访问速度是非常快的。

这里值得注意的是，如果一个CTA中包含的线程S_cta超过32个，也就是一个线程束的数量。一般来说，为了保证计算的正确性，必须在数据读取前后显式地调用CTA中线程同步指令来保证数据操作的正确性。这是由于不同线程束在执行过程中时异步执行的。由于线程同步的开销非常非常的大，因此这是需要想办法要避免的。另一方面，在同一个线程束中的线程在执行时必然是完全同步的。如果一个CTA中的线程只有32个，那么这种天然的线程同步特性就可以被运用到保证数据正确性上。所以在本发明中，将S_cta设为32。

至此，全新的粒子交互方式已经说明完毕。现以图2中的单元格ε为例，说明这一部分具体的实现步骤：

1)由于单元格ε中只包含了5个粒子，因此，分配了一个包含了32个线程的CTA用来负责该单元格的计算。CTA中前5个线程与单元格ε中的粒子一一对应。

2)首先，该CTA中前5个线程将自己所负责的粒子的数据从全局内存读入寄存器中。

3)然后，该CTA中前9个线程将单元格ε周围9个长方体单元格中的粒子在内存中的范围计算出来并存储在共享内存中。由于图2是一个二维实例图，因此周围的长方体单元格只显式了3个，分别是单元格αβγ，单元格δεζ以及单元格ηθι。

4)之后，试图读取第一个长方体单元格αβγ中的粒子，由于该长方体单元格中的粒子有8个，因此CTA中前8个线程负责将这8个粒子的数据从全局内存读入到共享内存中。由于单元格αβγ已读完，下一次将读取单元格δεζ中的粒子数据。

5)然后根据在寄存器中的粒子数据以及在共享内存中的例子数据，CTA中前5个线程为他们所负责的粒子计算新的属性。

6)计算完后，读取新的邻居粒子的信息。这次读取的是单元格δεζ中的粒子数据。由于单元格δεζ中包含20个粒子，所以CTA中前20个线程负责读取这些粒子的信息并保存到共享内存中，读取完后依旧执行计算。

7)直到所有邻居长方体单元格中的数据都被读取完毕并完成计算，该CTA所负责的粒子的计算任务就完成了。

Claims (1)

1.一种面向SPH流体模拟的GPGPU并行计算方法，其特征在于该方法包括以下具体步骤：

步骤1：设SPH方法中的光滑核半径为h，将整个模拟空间划分为以h为边长的互不相交的正方体单元格；

步骤2：使用基于GPGPU的计数排序算法，将空间内的粒子根据它们所在的单元格来进行排序，在计数排序算法中，首先计算每个单元格中所包含的粒子数量，然后通过并行前缀求和算法计算每个单元格中的粒子在内存中的偏移量，之后根据每个粒子所在单元格在内存中的偏移量以及其在该单元格中的位置进行复制；

步骤3：根据每个单元格中粒子的数量，为接下来的粒子交互阶段中每一个CTA分配任务，此过程完全在GPGPU上实现：

i)计算每个单元格所包含的粒子在内存中的范围，其结果保存在两个数组中，数组中每个元素描述一个单元格的相关信息：(1)单元格中的粒子在内存中的偏移量；(2)单元格所包含的粒子数量B；这步由个GPU线程并行完成，为空间中所包含的单元格的数量；

ⅱ)计算每一个单元格将在步骤4中所需要用来处理其内部粒子的交互的CTA数量，并将结果保存在数组B’中；所需的CTA数量为其中，S_cta表示一个CTA中包含的线程数量，表示单元格c_j中所包含的粒子的数量，这些数据保存在数组B中；这步由个GPU线程并行完成；

ⅲ)对数组B’中的元素做前缀求和并将结果保存在数组O’中；O’中的元素表示每个单元格所需的任务在任务队列T中的起始位置；这步由个GPU线程并行完成；

ⅳ)T中每个元素t_i＝{c_i,k_i}，其中c_i表示当前这个任务需要处理哪个单元格中的粒子，而k_i代表自己是第几个处理该单元格的任务；这步由个GPU线程并行完成；

ⅴ)获取总共需要的CTA数量S_grid；该值通过将O’中最后一个值与最后一个非空单元格所需的CTA数量相加得到；

步骤4：进行粒子交互；粒子交互包括两个独立的遍历邻居粒子阶段：读取邻居粒子的位置与质量来计算每个离散粒子表示的液体密度；通过邻居粒子的密度信息来计算每个离散粒子的受力；这两个阶段遍历邻居粒子的实现流程是相同的，具体步骤如下：

i)每个线程将自己所负责的粒子的必要数据从全局内存中读取到寄存器中；如果该线程不负责任何粒子，则跳过此步骤；

ii)一个CTA中前9个线程计算它们所负责的那个单元格周围的单元格的粒子在内存中的保存范围，结果保存在共享内存中；周围的单元格的粒子在内存中的范围是9个连续区域；其余线程跳过此步骤；

iii)读取一个邻居粒子区域中最多S_cta个粒子，保存到共享内存中；CTA中每个线程最多负责读取一个；

iv)凡是有粒子负责的线程通过在寄存器中自己负责的粒子数据以及共享内存中最多S_cta个邻居粒子的数据来计算自己负责的粒子的密度或是受力；

v)S_cta设为32，也就是一个线程束中所包含的线程的数量；因此内存访问与计算之间不需要显式的线程同步；

步骤5：根据步骤4中计算出的每个粒子的受力情况，来移动粒子；移动粒子的时候要考虑外力影响与边界情况；当粒子在边界上时给予一个反向的力将粒子弹回模拟空间内；

步骤6：一个时间步中所需要的计算就完成了；之后，输出粒子的位置信息，然后计算下一个时间步中粒子的变化，回到步骤1。