面向Ghost光滑粒子动力学模拟的图形处理器快速邻居搜索算法.pdf_大学文库

工程科学学报，第37卷，第5期：661667,2015年5月 Chinese Journal of Engineering,Vol.37,No.5:661-667,May 2015 DOI:10.13374/j.issn2095-9389.2015.05.019:http://journals.ustb.edu.cn 面向Ghost光滑粒子动力学模拟的图形处理器快速邻居搜索算法刘旭，班晓娟四，陈昙吴，张雅斓北京科技大学计算机与通信工程学院，北京100083 ☒通信作者，E-mail:Banxj(@usth.cd.cn 摘要提出了一种全新的快速邻居搜索方法，该方法可提高基于光滑粒子动力学的流体模拟在图形处理器上的运行效率. 此外，这种新的邻居表建立方法可以对两种或者两种以上的粒子进行邻居搜索，使所有粒子能在同一背景网格下拥有独立的粒子属性.在此基础上，引入了Gost边界粒子以加强光滑粒子动力学方法在边界模拟上的准确性，从而使流体模拟更加真实.实验证明，与传统的基于图形处理器的光滑粒子动力学模拟相比，本文方法效率更高. 关键词搜索算法：邻居：光滑粒子动力学：图形处理器：并行计算：流体模拟分类号TP391.41 Fast neighbor search on GPU for Ghost SPH simulation LIU Xu,BAN Xiao-juan,CHEN Tan-hao,ZHANG Ya-lan School of Computer Communication Engineering,University of Science&Technology Beijing,Beijing 100083,China Corresponding author,E-mail:Banxj@ustb.edu.cn ABSTRACT This paper presents a novel fast neighbor searching method.By using this method,fluid simulation based on smooth particle hydrodynamics (SPH)can be parallelized easily and run on graphic processing unit (GPU)with high efficiency.The neigh- bor searching method can search two or more kinds of particles,while saving their information in the same background grid.Ghost boundary particles are introduced to improve the accuracy of boundaries,which can enhance the fidelity of the fluid simulation.Exper- iments show that the proposed method is more efficient than the traditional SPH method based on GPU. KEY WORDS search algorithm:neighbor;smooth particle hydrodynamics:graphic processing unit:parallel computing:fluid sim- ulation 在计算机图形学中，基于物理的流体模拟始终是 SPH流体模拟是一种基于拉格朗日观点的模拟方一项重要的研究问题，与传统建模的模拟方法相比其法，其将流体看成大量流体微元，通过微元间的相互作具有效果逼真、节省人力等特点.基于物理的流体模用来模拟流体.与其他的流体模拟方法相比，SPH利拟可分网格法、无网格法和半网格法.相比之下，无网用光滑核函数对流体的支配方程Navier-Stokes方程格法和半网格法在处理复杂的液体表面时效果更加真进行近似处理，因此其计算效率由粒子及其邻居粒子实，且不存在网格细分等问题.随着光滑粒子动力学的组织关系直接决定.此外，在传统图形处理器 (smooth particle hydrodynamics,SPH)方法的提出，基 (graphic processing unit,.GPU)并行算法的限制下，只于该方法的无网格流体模拟逐渐成为了流体模拟领域能使用单一性质的粒子进行运算，这对边界检测带来的研究热点，在许多领域有着重要应用0 了一定影响收稿日期：20150105 基金项目：国家自然科学基金资助项目(61272357,61300074)

工程科学学报，第 37 卷，第 5 期: 661--667，2015 年 5 月 Chinese Journal of Engineering，Vol． 37，No． 5: 661--667，May 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 05． 019; http: / /journals． ustb． edu． cn 面向 Ghost 光滑粒子动力学模拟的图形处理器快速邻居搜索算法刘旭，班晓娟，陈昙昊，张雅斓北京科技大学计算机与通信工程学院，北京 100083  通信作者，E-mail: Banxj@ ustb． edu． cn 摘要提出了一种全新的快速邻居搜索方法，该方法可提高基于光滑粒子动力学的流体模拟在图形处理器上的运行效率．此外，这种新的邻居表建立方法可以对两种或者两种以上的粒子进行邻居搜索，使所有粒子能在同一背景网格下拥有独立的粒子属性．在此基础上，引入了 Ghost 边界粒子以加强光滑粒子动力学方法在边界模拟上的准确性，从而使流体模拟更加真实．实验证明，与传统的基于图形处理器的光滑粒子动力学模拟相比，本文方法效率更高．关键词搜索算法; 邻居; 光滑粒子动力学; 图形处理器; 并行计算; 流体模拟分类号 TP391. 41 Fast neighbor search on GPU for Ghost SPH simulation LIU Xu，BAN Xiao-juan ，CHEN Tan-hao，ZHANG Ya-lan School of Computer ＆ Communication Engineering，University of Science ＆ Technology Beijing，Beijing 100083，China  Corresponding author，E-mail: Banxj@ ustb． edu． cn ABSTＲACT This paper presents a novel fast neighbor searching method． By using this method，fluid simulation based on smooth particle hydrodynamics ( SPH) can be parallelized easily and run on graphic processing unit ( GPU) with high efficiency． The neighbor searching method can search two or more kinds of particles，while saving their information in the same background grid． Ghost boundary particles are introduced to improve the accuracy of boundaries，which can enhance the fidelity of the fluid simulation． Experiments show that the proposed method is more efficient than the traditional SPH method based on GPU． KEY WOＲDS search algorithm; neighbor; smooth particle hydrodynamics; graphic processing unit; parallel computing; fluid simulation 收稿日期: 2015--01--05 基金项目: 国家自然科学基金资助项目( 61272357，61300074) 在计算机图形学中，基于物理的流体模拟始终是一项重要的研究问题，与传统建模的模拟方法相比其具有效果逼真、节省人力等特点．基于物理的流体模拟可分网格法、无网格法和半网格法．相比之下，无网格法和半网格法在处理复杂的液体表面时效果更加真实，且不存在网格细分等问题．随着光滑粒子动力学 ( smooth particle hydrodynamics，SPH) 方法的提出，基于该方法的无网格流体模拟逐渐成为了流体模拟领域的研究热点，在许多领域有着重要应用［1］． SPH 流体模拟是一种基于拉格朗日观点的模拟方法，其将流体看成大量流体微元，通过微元间的相互作用来模拟流体．与其他的流体模拟方法相比，SPH 利用光滑核函数对流体的支配方程 Navier--Stokes 方程进行近似处理，因此其计算效率由粒子及其邻居粒子的组织关系直接决定．此外，在传统图形处理器 ( graphic processing unit，GPU) 并行算法的限制下，只能使用单一性质的粒子进行运算，这对边界检测带来了一定影响．

·662 工程科学学报，第37卷，第5期在模拟流体过程中，粒子的数量有可能会达到几万、几十万甚至更多，因此在寻找周围的邻居粒子的过程中，遍历所有粒子找出满足条件的邻居粒子会带来巨大的计算量.Crespo等因在场景中加入了背景网格，使之在搜索邻居粒子时只需找到周围网格的粒子即可，无需大规模的遍历所有粒子.Dominguez等在文献B]中归纳了几种主流的、基本的邻居粒子搜寻方法. 除了在邻居搜索算法上进行优化之外，使用GPU 进行并行化来提高效率也是最近几年流行的方法.陈飞国等利用CUDA(compute unified device architec- ture)技术进行了GPU上的分子动力学(MD)模拟.左图1邻居粒子分布图颢容等同提出的基于GPU的并行优化技术，优化后累 Fig.1 Distribution of neighbor particles 加和算法的运算速度相比标准并行算法提高了约34 倍，相比CPU串行实现提高了约70倍.Goswami等率高且技术成熟的SPH模拟方法.SPH是一种无网格使用GPU完成将SPH并行化，并且使用了Z排序的并的拉格朗日方法，其描述粒子的基础就是Navier-一行化算法来降低GPU内存开销及访存开销.Domingu- Stokes方程.Navier--Stokes的表示形式如下： ez和Crespo提出了在GPU上的计算模型的优化，将 ou =-(mV)u-上Vp+nw+pg (1) 原本的邻居粒子搜索和系统更新也在GPU上实现，达 p (2) 到提高效率的目的.Valdez---Balderasa和Dominguez Vu=0. 在文献中使用了多图像处理器(Muli一GPUs)来应对式中，“为流体速度，p为压强，V为梯度算子，为拉粒子数量级达到GPU内存无法负荷的情况，此时可以普拉斯算子，g为外力场，p为粒子密度，v为黏性系联合多个GPU,存储空间以及计算效率将得到明显数，1为时间. 提升. (1)连续性方程.原本的SPH中密度求和公式为本文提出一种搜索邻居粒子的方法，该方法不仅 p=∑m,W (3) 加快了邻居粒子的索引速度，而且在同一背景网格中其中W=W(x。-x,),W为光滑核函数，a代表粒子可以索引多种类型的粒子，因此引入Ghost边界粒子 a,b同a.x。代表粒子a的位置，m,为粒子b的质量. 来进一步解决边界检测问题.其中，Ghost边界粒子是由于每个粒子的质量不变，所以由连续性方程可在固体和流体粒子之间加入的一层虚拟粒子层以得出密度的求解公式： 1相关知识告= (4) 1.1基本概念其中uh=u。-。本文在邻居搜索算法的基础上对一些相关术语作 (2)压强.本文采用泰勒方程来计算压强团，因以下定义为泰勒方程在计算密度变化较小的压强时效率较高. 定义1邻居粒子：在三维空间中，定义位于流体泰勒方程的表示如下：粒子周围的26个网格（如图1）内的粒子为该粒子的 (5) 邻居粒子，以此简化搜寻邻居粒子的方法，从而提升计算效率。其中y=7,在文献7]中，B≈1119kPa,Po=1000kg· 定义2网格序号：网格序号与每个粒子的属性 m 一一对应.将其存储在名为网格索引的数组里，对网 (3)动量方程.通过动量方程，利用压强和密度格索引进行排序，然后用首索引数组记录排序后的粒可以推演出粒子的加速度子在每个网格的起始位置和结束位置.例如，5号粒子 (6) 的网格坐标是(0,2)，则其网格序号就为3 dt L.2SPH理论根据加速度即可以更新粒子位置本文中所用到的SPH模拟方法是Becker和Te- (4)边界条件.Fedkiw等o所提出的“Ghost Flu- schner在2007年提出的弱可压SPH方法(weakly id Method'”,已经被广泛的应用于欧拉流体，包括耦合 compressible SPH,W/CSPH)）,它是一种模拟效果好、效的拉格朗日固体m.Schechter和Bridson☒提出了一

工程科学学报，第 37 卷，第 5 期在模拟流体过程中，粒子的数量有可能会达到几万、几十万甚至更多，因此在寻找周围的邻居粒子的过程中，遍历所有粒子找出满足条件的邻居粒子会带来巨大的计算量． Crespo 等［2］在场景中加入了背景网格，使之在搜索邻居粒子时只需找到周围网格的粒子即可，无需大规模的遍历所有粒子． Dominguez 等在文献［3］中归纳了几种主流的、基本的邻居粒子搜寻方法．除了在邻居搜索算法上进行优化之外，使用 GPU 进行并行化来提高效率也是最近几年流行的方法．陈飞国等［4］利用 CUDA( compute unified device architecture) 技术进行了 GPU 上的分子动力学( MD) 模拟．左颢睿等［5］提出的基于 GPU 的并行优化技术，优化后累加和算法的运算速度相比标准并行算法提高了约 34 倍，相比 CPU 串行实现提高了约 70 倍． Goswami 等［6］使用 GPU 完成将 SPH 并行化，并且使用了 Z 排序的并行化算法来降低 GPU 内存开销及访存开销． Dominguez 和 Crespo［7］提出了在 GPU 上的计算模型的优化，将原本的邻居粒子搜索和系统更新也在 GPU 上实现，达到提高效率的目的． Valdez--Balderasa 和 Dominguez［8］在文献中使用了多图像处理器( Multi--GPUs) 来应对粒子数量级达到 GPU 内存无法负荷的情况，此时可以联合多个 GPU，存储空间以及计算效率将得到明显提升．本文提出一种搜索邻居粒子的方法，该方法不仅加快了邻居粒子的索引速度，而且在同一背景网格中可以索引多种类型的粒子，因此引入 Ghost 边界粒子来进一步解决边界检测问题．其中，Ghost 边界粒子是在固体和流体粒子之间加入的一层虚拟粒子层． 1 相关知识 1. 1 基本概念本文在邻居搜索算法的基础上对一些相关术语作以下定义．定义 1 邻居粒子: 在三维空间中，定义位于流体粒子周围的 26 个网格( 如图 1) 内的粒子为该粒子的邻居粒子，以此简化搜寻邻居粒子的方法，从而提升计算效率．定义 2 网格序号: 网格序号与每个粒子的属性一一对应．将其存储在名为网格索引的数组里，对网格索引进行排序，然后用首索引数组记录排序后的粒子在每个网格的起始位置和结束位置．例如，5 号粒子的网格坐标是( 0，2) ，则其网格序号就为 3． 1. 2 SPH 理论本文中所用到的 SPH 模拟方法是 Becker 和 Teschner［9］在 2007 年提出的弱可压 SPH 方法 ( weakly compressible SPH，WCSPH) ，它是一种模拟效果好、效图 1 邻居粒子分布图 Fig． 1 Distribution of neighbor particles 率高且技术成熟的 SPH 模拟方法． SPH 是一种无网格的拉格朗日方法，其描述粒子的基础就是 Navier-- Stokes 方程． Navier--Stokes 的表示形式如下: u t = － ( u Δ ) u － 1 ρ Δ p + v 2 Δ u + ρg． ( 1) Δ u = 0． ( 2) 式中，u 为流体速度，p 为压强， Δ 为梯度算子， 2 Δ 为拉普拉斯算子，g 为外力场，ρ 为粒子密度，v 为黏性系数，t 为时间． ( 1) 连续性方程．原本的 SPH 中密度求和公式为 ρa = ∑b mbWab ． ( 3) 其中 Wab = W( xa － xb ) ，W 为光滑核函数，a 代表粒子 a，b 同 a． xa 代表粒子 a 的位置，mb 为粒子 b 的质量．由于每个粒子的质量不变，所以由连续性方程可以得出密度的求解公式: dρa dt = ∑b mbuab Δ Wab ． ( 4) 其中 uab = ua － ub ． ( 2) 压强．本文采用泰勒方程来计算压强［7］，因为泰勒方程在计算密度变化较小的压强时效率较高．泰勒方程的表示如下: P = B [ ( ρ ρ ) 0 γ ] － 1 ． ( 5) 其中 γ = 7，在文献［7］中，B≈1119 kPa，ρ0 = 1000 kg· m － 3 ． ( 3) 动量方程．通过动量方程，利用压强和密度可以推演出粒子的加速度． dua dt = － ∑b mb ( pa ρ 2 a + pb ρ 2 ) b Δ Wab + g． ( 6) 根据加速度即可以更新粒子位置． ( 4) 边界条件． Fedkiw 等［10］所提出的“Ghost Fluid Method”，已经被广泛的应用于欧拉流体，包括耦合的拉格朗日固体［11］． Schechter 和 Bridson［12］提出了一 · 266 ·

刘旭等：面向Ghost光滑粒子动力学模拟的图形处理器快速邻居搜索算法 ·663 种基于自由表面和固体边界条件的新的Ghost Fluid的 SPH模拟方法，即将周围的固体和气体，或者边界采样参数初始化、朝始化成一系列狭窄的粒子层，然后和场景中的流体粒子进传递参数行交互，产生边界条件. 场景粒子的邻居表建立 2面向Ghost SPH模拟的GPU快速邻居搜索算法边界虚拟粒子的邻居表建立在CPU上运行的算法分为两种，并行化算法和非并行化算法.并行化算法，即计算的对象在逻辑上要搜寻邻居粒子，计算密度，压强执行相同的操作，且对象相互之间没有数据依赖，在这里指的是对粒子的计算.本文需要将CPU上的非并行计算加速度化算法转化成GPU上的并行化算法. CPU端GP端和作用力本文将基于SPH方法的模拟过程分为五个部分 (如图2)，即系统初始化、邻居表建立、计算作用力、系更新粒子坐标统坐标更新和数据与图像显示.邻居表建立、计算作用力和系统坐标更新这三个部分是可以并行化处理俞历完所有粒子否的.对于每个粒子，邻居表建立只需要每一帧的粒子、&《输出帧初始位置，这对于每一帧来说是一个定值，而非计算 CPL绘制 GPU绘制值.计算作用力需要周围粒子的密度和压强，这看似粒子和 GP渲染否粒子和场景场存在着数据依赖，不能进行并行化，但是密度和压强只与自己以及周围粒子的初始位置有关，这与邻居表建图2 Ghost SPH方法的模拟过程 Fig.2 Processes of simulation based on Ghost SPH method 立是一致的，因此将密度和压强的计算并入邻居表建立中.计算作用力只需完成对合力和加速度的计算. 1-2 系统坐标更新依赖于每一帧的初始速度、位置以及加时钟周期块内共享速度线程内存当场景粒子的数量增大的时候，由于粒子的分布网格比较聚集，导致GPU的块内的共享存储器大小不足以存储块内所有的粒子的属性，而将粒子完全存储在全块0.0) 块(1.0 局内存中又会增加访问时长，访问时间相差数百倍，如 400-600 时钟周期图3所示.此外，GPU的并行特性可能会受到一些特全局内存殊情况的影响，比如流体粒子分布不均匀.在未达到块0.1) 块(1.1) 稳定状态时，计算密集部分可能会较慢，而计算稀疏部分可能较快，导致并行程度并不理想根据通用并行计算架构(compute unified device ar- 图3CUDA内存结构及访问时长 chitecture,CUDA)的内存特性，本方法将粒子存储在全 Fig.3 Memory structure and access time 局内存中，在计算邻居表时只将粒子的索引值拷贝到块的共享内存中，并建立邻居表.这样整个排序过程中绝大部分的内存访问将在共享内存完成，提高了内 5 6 存访问速度.邻居表的存在，使得计算部分中最为耗 090 自时的寻找周围粒子的步骤的时耗大大降低. 2 本节采用滑动动态向量法对图4中的粒子进行排图4粒子分布示意图序，将每个粒子对应的背景网格序号存储在名为网格 Fig.4 Distribution of particles 索引的数组里，将背景网格序号与每个粒子的属性一对于本文中加入的边界Ghost粒子，由于滑动动一对应，对网格索引进行排序，然后用首索引数组记录态向量法只适用于计算过程相同的粒子，而Ghost粒网格索引中每个背景网格的开始位置，如图5所示，后子只提供边界粒子对场景粒子的作用，本身不需要计续计算可以根据起始位置搜索到每个背景网格内的算，因此需要对滑动动态向量法进行拓展，存储Ghost 粒子粒子的属性但计算与流体粒子分开，对场景粒子的计

刘旭等: 面向 Ghost 光滑粒子动力学模拟的图形处理器快速邻居搜索算法种基于自由表面和固体边界条件的新的 Ghost Fluid 的 SPH 模拟方法，即将周围的固体和气体，或者边界采样成一系列狭窄的粒子层，然后和场景中的流体粒子进行交互，产生边界条件． 2 面向 Ghost SPH 模拟的 GPU 快速邻居搜索算法在 CPU 上运行的算法分为两种，并行化算法和非并行化算法．并行化算法，即计算的对象在逻辑上要执行相同的操作，且对象相互之间没有数据依赖，在这里指的是对粒子的计算．本文需要将 CPU 上的非并行化算法转化成 GPU 上的并行化算法．本文将基于 SPH 方法的模拟过程分为五个部分 ( 如图 2) ，即系统初始化、邻居表建立、计算作用力、系统坐标更新和数据与图像显示．邻居表建立、计算作用力和系统坐标更新这三个部分是可以并行化处理的．对于每个粒子，邻居表建立只需要每一帧的粒子初始位置，这对于每一帧来说是一个定值，而非计算值．计算作用力需要周围粒子的密度和压强，这看似存在着数据依赖，不能进行并行化，但是密度和压强只与自己以及周围粒子的初始位置有关，这与邻居表建立是一致的，因此将密度和压强的计算并入邻居表建立中．计算作用力只需完成对合力和加速度的计算．系统坐标更新依赖于每一帧的初始速度、位置以及加速度．当场景粒子的数量增大的时候，由于粒子的分布比较聚集，导致 GPU 的块内的共享存储器大小不足以存储块内所有的粒子的属性，而将粒子完全存储在全局内存中又会增加访问时长，访问时间相差数百倍，如图 3 所示．此外，GPU 的并行特性可能会受到一些特殊情况的影响，比如流体粒子分布不均匀．在未达到稳定状态时，计算密集部分可能会较慢，而计算稀疏部分可能较快，导致并行程度并不理想．根据通用并行计算架构( compute unified device architecture，CUDA) 的内存特性，本方法将粒子存储在全局内存中，在计算邻居表时只将粒子的索引值拷贝到块的共享内存中，并建立邻居表．这样整个排序过程中绝大部分的内存访问将在共享内存完成，提高了内存访问速度．邻居表的存在，使得计算部分中最为耗时的寻找周围粒子的步骤的时耗大大降低．本节采用滑动动态向量法对图 4 中的粒子进行排序，将每个粒子对应的背景网格序号存储在名为网格索引的数组里，将背景网格序号与每个粒子的属性一一对应，对网格索引进行排序，然后用首索引数组记录网格索引中每个背景网格的开始位置，如图 5 所示，后续计算可以根据起始位置搜索到每个背景网格内的粒子．图 2 Ghost SPH 方法的模拟过程 Fig． 2 Processes of simulation based on Ghost SPH method 图 3 CUDA 内存结构及访问时长 Fig． 3 Memory structure and access time 图 4 粒子分布示意图 Fig． 4 Distribution of particles 对于本文中加入的边界 Ghost 粒子，由于滑动动态向量法只适用于计算过程相同的粒子，而 Ghost 粒子只提供边界粒子对场景粒子的作用，本身不需要计算，因此需要对滑动动态向量法进行拓展，存储 Ghost 粒子的属性但计算与流体粒子分开，对场景粒子的计 · 366 ·

·664· 工程科学学报，第37卷，第5期算过程中也能找到周围的Ghost粒子. 因此，用另一张邻居表来存储Ghost粒子，进行相应的操作，让这张邻居表也按场景粒子的邻居表对应的背景网格来划分.这样，两张邻居表等同于一张扩展的邻居表.也就是说，在场景粒子搜寻邻居时，对两张邻居表的操作是基本相同的. 邻居表建立的步骤（如图5）： 6 2 2 5 (1)将所有的场景粒子对应到背景网格中（如图 4),并给场景粒子建立一个粒子索引数组. 粒子素引网格索引粒子索引网格索引 (2)计算场景粒子对应的网格序号，这里的序号 (b) 是将三维坐标转成线性排列的序号，存储在网格索引数组中（如图5(a). 网格1 (3)用基数排序将网格索引以网格序号的升序进 2 网格2 行排序，粒子索引数组随之改动，这样排序后，同一网网格3 格内的场景粒子在数组内是连续的，方便后面的邻居 6 网格4 表的建立（如图5(b)). 网格5 (4)最后，根据网格索引和粒子索引数组两个数 7 网格6 组，创建一个首索引数组，将排序后网格索引中每个背结桌网格景网格的首地址插入到首索引数组中.这样，根据排序后的网格序号网格索引和首索引数组就能定位到需网格索引首索引背景网格要的粒子首索引数组定义如下，首索引数组的长度比背景网格数量大1.即如共有n个背景网格，首索引数组的图5邻居表建立的具体过程.(a)建立索引，计算背景网格：长度为n+1,首索引数组中每一个位置代表一个背景 (h)进行基数排序：(c)插入首索引数组网格中的每一个网格，首索引数组中的内容为网格 Fig.5 Specific process of neighbor table establishment:(a)back- ground table and index initialization:(b)radix sort:(c)setting up 索引数组中网格编号的起始位置.此处可将首索引 begin array 数组视为图5(b)中网格索引数组的索引.图5(c) 描述了首索引数组的填写过程.通过首索引数组我居表来存储不同粒子的信息，这样在计算场景粒子的们可以准确计算每一个背景网格中所包含的粒子合作用力时，可以分别搜索相应邻居表来调用对应的数，如cell_3ell2为3说明背景网格2包含3个粒不同的周围粒子信息，计算多种流体粒子或者刚体粒子，再通过图5(b)中网格索引数组及粒子索引数组子、边界Ghost粒子等对其的影响.此外，与传统背景两个数组即可查找到相应的粒子.由于背景网格4 网格方法相比，该方法虽然在粒子排序等预处理上消为空，因此首索引数组中ce_4与cell5的值相同，耗了计算时间，但是在进行邻居搜索时，由于可以准确两值做差可知背景网格4包含0个粒子，首索引数组获得每网格中的粒子总数，且由于粒子已排序，所以可的最后一个网格cell_end旨在表明最后一个背景网以快速获取一个网格中的粒子信息，而不需要遍历所格中包含的粒子数. 有粒子或像传统的背景网格方法利用指针进行访问，同样地，Ghost粒子的邻居表也如上述步骤一样. 并且由于共享内存与全局内存的访问速度上的差异，虽然场景粒子和Ghost粒子自带的属性不全相同，但系统的运行效率得到了提高粒子的邻居关系只与粒子的当前帧位置有关.在这一作用力及加速度等计算，按照第1节所述方法进点上，场景粒子和Ghost粒子是等同的，因此可以进行行计算即可.过程中会频繁遇到邻居搜索的问题，由相同的邻居表建立的操作.然后计算作用力、加速度于邻居表已建立，并且粒子信息在计算密度时已从全并找出邻居后，再根据对应的邻居表，调用对应的粒子局内存中读取，此处访问开销不大.从图2看出，GPU 属性即可. 中的粒子数据可以按固定步长传回到CPU中渲染，也根据邻居表建立时的特性可以看出，本方法使不可以通过通用并行计算架构(CUDA)和开放图形库同粒子之间不互相干扰，但拥有相同的背景网格，且能 (Open Graphics Library,OPENGL)交互直接渲染.此处独立调用任意一种粒子的属性.对于多种流体或者不即使是传回到CPU渲染，也只传回粒子的位置信息同刚体采样出来的表面粒子，可以通过建立不同的邻即可

工程科学学报，第 37 卷，第 5 期算过程中也能找到周围的 Ghost 粒子．因此，用另一张邻居表来存储 Ghost 粒子，进行相应的操作，让这张邻居表也按场景粒子的邻居表对应的背景网格来划分．这样，两张邻居表等同于一张扩展的邻居表．也就是说，在场景粒子搜寻邻居时，对两张邻居表的操作是基本相同的．邻居表建立的步骤( 如图 5) : ( 1) 将所有的场景粒子对应到背景网格中( 如图 4) ，并给场景粒子建立一个粒子索引数组． ( 2) 计算场景粒子对应的网格序号，这里的序号是将三维坐标转成线性排列的序号，存储在网格索引数组中 ( 如图 5( a) ) ． ( 3) 用基数排序将网格索引以网格序号的升序进行排序，粒子索引数组随之改动，这样排序后，同一网格内的场景粒子在数组内是连续的，方便后面的邻居表的建立( 如图 5( b) ) ． ( 4) 最后，根据网格索引和粒子索引数组两个数组，创建一个首索引数组，将排序后网格索引中每个背景网格的首地址插入到首索引数组中．这样，根据排序后的网格序号网格索引和首索引数组就能定位到需要的粒子．首索引数组定义如下，首索引数组的长度比背景网格数量大 1．即如共有 n 个背景网格，首索引数组的长度为 n + 1，首索引数组中每一个位置代表一个背景网格中的每一个网格，首索引数组中的内容为网格索引数组中网格编号的起始位置．此处可将首索引数组视为图 5 ( b) 中网格索引数组的索引．图 5 ( c) 描述了首索引数组的填写过程．通过首索引数组我们可以准确计算每一个背景网格中所包含的粒子数，如 cell_3-cell_2 为 3 说明背景网格 2 包含 3 个粒子，再通过图 5( b) 中网格索引数组及粒子索引数组两个数组即可查找到相应的粒子．由于背景网格 4 为空，因此首索引数组中 cell_4 与 cell_5 的值相同，两值做差可知背景网格 4 包含 0 个粒子，首索引数组的最后一个网格 cell_end 旨在表明最后一个背景网格中包含的粒子数．同样地，Ghost 粒子的邻居表也如上述步骤一样．虽然场景粒子和 Ghost 粒子自带的属性不全相同，但粒子的邻居关系只与粒子的当前帧位置有关．在这一点上，场景粒子和 Ghost 粒子是等同的，因此可以进行相同的邻居表建立的操作．然后计算作用力、加速度并找出邻居后，再根据对应的邻居表，调用对应的粒子属性即可．根据邻居表建立时的特性可以看出，本方法使不同粒子之间不互相干扰，但拥有相同的背景网格，且能独立调用任意一种粒子的属性．对于多种流体或者不同刚体采样出来的表面粒子，可以通过建立不同的邻图 5 邻居表建立的具体过程 . ( a) 建立索引，计算背景网格; ( b) 进行基数排序; ( c) 插入首索引数组 Fig． 5 Specific process of neighbor table establishment: ( a) background table and index initialization; ( b) radix sort; ( c) setting up begin array 居表来存储不同粒子的信息，这样在计算场景粒子的合作用力时，可以分别搜索相应邻居表来调用对应的不同的周围粒子信息，计算多种流体粒子或者刚体粒子、边界 Ghost 粒子等对其的影响．此外，与传统背景网格方法相比，该方法虽然在粒子排序等预处理上消耗了计算时间，但是在进行邻居搜索时，由于可以准确获得每网格中的粒子总数，且由于粒子已排序，所以可以快速获取一个网格中的粒子信息，而不需要遍历所有粒子或像传统的背景网格方法利用指针进行访问，并且由于共享内存与全局内存的访问速度上的差异，系统的运行效率得到了提高．作用力及加速度等计算，按照第 1 节所述方法进行计算即可．过程中会频繁遇到邻居搜索的问题，由于邻居表已建立，并且粒子信息在计算密度时已从全局内存中读取，此处访问开销不大．从图 2 看出，GPU 中的粒子数据可以按固定步长传回到 CPU 中渲染，也可以通过通用并行计算架构( CUDA) 和开放图形库 ( Open Graphics Library，OPENGL) 交互直接渲染．此处即使是传回到 CPU 渲染，也只传回粒子的位置信息即可． · 466 ·

刘旭等：面向Ghost光滑粒子动力学模拟的图形处理器快速邻居搜索算法 ·665 本文的实验是经典的Dam Break的模拟，图6是 3算法仿真与性能评估 10个流体粒子和50O0个边界Ghost粒子的Dam Break 3.1实验环境实验.图中绿色为流体粒子，红色为Ghost边界粒子. 在Windows7操作系统中，Visual Studio2010开发图6(a)为系统初始状态.图6(b)为有Ghost粒子参平台下，用通用并行计算架构(CUDA)实现了GPU上与的模拟在0.15s时的结果.为了对比需要，在图6 基于Ghost边界粒子的邻居搜索算法的SPH流体模 (b)中隐藏了Ghost粒子的显示.图6(c)为无Ghost 拟.硬件环境：CPU是Intel至强E5V2637,内存80 粒子时系统运行到0.15s的模拟效果.从图6(b)和 GB.显卡为NVIDIA Quadro K4000. 图6(c)两幅图中可以看出，图6(c)的边界处粒子出 3.2实验结果现了不齐整、不自然的现象，这是由于流体边界处的密本文中提出的基于Ghost边界粒子的邻居搜索算度过小导致的，Ghost边界粒子的加入弥补了密度上的法不仅能较大地利用GPU的并行特性以提高模拟的缺失，使流体表面更加光滑，更加真实.图6是对模拟效率，而且在边界上的真实性得到了提高效果的展示，表1和表2对运行效率进行了分析 a 图6同一时刻有无Chost边界粒子的效果对比（图中为粒子运行到0.15s时遗坝的模拟截图）.(a)初始状态：(b)有Ghost边界粒子：(c) 没有Ghost边界粒子 Fig.6 Comparison of results with and without Ghost particles (the state of the simulation of dam break at 0.15s).(a)initial state (b)running with Ghost particles,(c)running without Ghost particles 以下是用本文的算法和CPU链式搜索方法国、多表215×103粒子数下不同算法在每个计算步骤的耗费时间线程CPU链式搜索方法以及GPU的式搜索方法进行 Table 2 Cost time of each step for different algorithms with 15 x 10 particles 对比，见表1. 邻居表建立计算力耗时/ 更新坐标表1不同粒子数下不同算法的耗费时间算法耗时/ms ms 耗时/ms Table 1 Running time of different algorithms at different particles num- 本文算法 2.1 3.9 0.3 bers GPU算法 0.7 6.4 0.3 粒子数/ 本文算法/ GPU CPU多线 CPU CPU算法 1 49.2 0.5 103 ms 算法/ms 程/ms 算法/ms 15 6.3 7.4 23 51 图7是在系统运行到同一时刻0.15s时，对是否 30 12 15 50 145 使用GPU加速，以及是否使用Ghost边界粒子进行实 70 31 哈 122 320 验效果对比.实验证明，无论是否有Ghost边界粒子参与计算，本文所述的GPU加速算法对模拟结果均不产表1中的耗费时间是指邻居表建立，计算压强、密生影响.图7的第一列为该系统在CPU上运行的结度、作用力以及更新坐标所需的耗费时长.由表1可果，第二列为该系统用本文所述方法运行的结果。图7 以看出，本文的算法较CPU的链表搜索算法有非常明第一行为无Ghost运行结果，第二行为有Ghost运行结显的提升，即便是与四线程的CPU算法也有接近6的果.图8则完整展示了溃坝的运行结果.为了展示需加速比.由于邻居搜索的部分所占的耗时比重并不要，此处对流体表面进行了重构高，本算法对GPU链表搜索方式没有显著的效果，但是也有一定的效率的提升 4结论从表2中可以看出，邻居搜索在整体计算中所占本文中提出了一种面向Ghost粒子的GPU邻居搜比重并不高，这导致在同为GPU上的算法得到的加速索方法.这种方法可以确保SPH能在GPU上正确运比不是非常明显，更大的加速比是在GPU和CPU计算行.实验表明，本算法相比于传统的CPU算法和CPU 性能上的提升并行算法，在效率上有明显提升：与GPU链式搜索方

刘旭等: 面向 Ghost 光滑粒子动力学模拟的图形处理器快速邻居搜索算法 3 算法仿真与性能评估 3. 1 实验环境在 Windows 7 操作系统中，Visual Studio 2010 开发平台下，用通用并行计算架构( CUDA) 实现了 GPU 上基于 Ghost 边界粒子的邻居搜索算法的 SPH 流体模拟．硬件环境: CPU 是 Intel 至强 E5 V2637，内存 80 GB．显卡为 NVIDIA Quadro K4000． 3. 2 实验结果本文中提出的基于 Ghost 边界粒子的邻居搜索算法不仅能较大地利用 GPU 的并行特性以提高模拟的效率，而且在边界上的真实性得到了提高．本文的实验是经典的 Dam Break 的模拟，图 6 是 104 个流体粒子和 5000 个边界 Ghost 粒子的 Dam Break 实验．图中绿色为流体粒子，红色为 Ghost 边界粒子．图 6( a) 为系统初始状态．图 6( b) 为有 Ghost 粒子参与的模拟在 0. 15 s 时的结果．为了对比需要，在图 6 ( b) 中隐藏了 Ghost 粒子的显示．图 6 ( c) 为无 Ghost 粒子时系统运行到 0. 15 s 的模拟效果．从图 6( b) 和图 6( c) 两幅图中可以看出，图 6( c) 的边界处粒子出现了不齐整、不自然的现象，这是由于流体边界处的密度过小导致的，Ghost 边界粒子的加入弥补了密度上的缺失，使流体表面更加光滑，更加真实．图 6 是对模拟效果的展示，表 1 和表 2 对运行效率进行了分析．图 6 同一时刻有无 Ghost 边界粒子的效果对比( 图中为粒子运行到 0. 15 s 时溃坝的模拟截图) ． ( a) 初始状态; ( b) 有 Ghost 边界粒子; ( c) 没有 Ghost 边界粒子 Fig． 6 Comparison of results with and without Ghost particles ( the state of the simulation of dam break at 0. 15 s) . ( a) initial state ( b) running with Ghost particles，( c) running without Ghost particles 以下是用本文的算法和 CPU 链式搜索方法［3］、多线程 CPU 链式搜索方法以及 GPU 的式搜索方法进行对比，见表 1．表 1 不同粒子数下不同算法的耗费时间 Table 1 Ｒunning time of different algorithms at different particles numbers 粒子数/ 103 本文算法/ ms GPU 算法/ms CPU 多线程/ms CPU 算法/ms 15 6. 3 7. 4 25 51 30 12 15 50 145 70 31 36 122 320 表 1 中的耗费时间是指邻居表建立，计算压强、密度、作用力以及更新坐标所需的耗费时长．由表 1 可以看出，本文的算法较 CPU 的链表搜索算法有非常明显的提升，即便是与四线程的 CPU 算法也有接近 6 的加速比．由于邻居搜索的部分所占的耗时比重并不高，本算法对 GPU 链表搜索方式没有显著的效果，但是也有一定的效率的提升．从表 2 中可以看出，邻居搜索在整体计算中所占比重并不高，这导致在同为 GPU 上的算法得到的加速比不是非常明显，更大的加速比是在 GPU 和 CPU 计算性能上的提升．表 2 15 × 103粒子数下不同算法在每个计算步骤的耗费时间 Table 2 Cost time of each step for different algorithms with 15 × 103 particles 算法邻居表建立耗时/ms 计算力耗时/ ms 更新坐标耗时/ms 本文算法 2. 1 3. 9 0. 3 GPU 算法 0. 7 6. 4 0. 3 CPU 算法 1 49. 2 0. 5 图 7 是在系统运行到同一时刻 0. 15 s 时，对是否使用 GPU 加速，以及是否使用 Ghost 边界粒子进行实验效果对比．实验证明，无论是否有 Ghost 边界粒子参与计算，本文所述的 GPU 加速算法对模拟结果均不产生影响．图 7 的第一列为该系统在 CPU 上运行的结果，第二列为该系统用本文所述方法运行的结果．图 7 第一行为无 Ghost 运行结果，第二行为有 Ghost 运行结果．图 8 则完整展示了溃坝的运行结果．为了展示需要，此处对流体表面进行了重构． 4 结论本文中提出了一种面向 Ghost 粒子的 GPU 邻居搜索方法．这种方法可以确保 SPH 能在 GPU 上正确运行．实验表明，本算法相比于传统的 CPU 算法和 CPU 并行算法，在效率上有明显提升; 与 GPU 链式搜索方 · 566 ·

工程科学学报，第 37 卷，第 5 期图 7 同一时刻有无 Ghost 边界粒子在 CPU 和 GPU 上运行的效果对比( 图中为粒子运行到 0. 15 s 时溃坝的模拟截图) . ( a) 无 Ghost 粒子在 CPU 上的运行结果; ( b) 无 Ghost 粒子使用本文算法在 GPU 的运行结果; ( c) 有 Ghost 粒子在 CPU 上的运行结果; ( d) 有 Ghost 粒子使用本文算法在 GPU 的运行结果 Fig． 7 Comparison of results with and without Ghost particles running on CPU and GPU ( the state of the simulation is dam breaking at 0. 15 s． ) : ( a) running on CPU without Ghost particles; ( b) running on GPU without Ghost particles; ( c) running on CPU with Ghost particles; ( d) running on GPU with Ghost particles 图 8 完整的溃坝模拟: 第一排从左至右分别为第 0 秒、1. 2 秒和 2. 4 秒时的运行结果，第二排从左至右分别为第 4. 8 秒、6 秒和 7. 2 秒的模拟结果 Fig． 8 Simulation of dam breaking: on first row from left to right are the simulation results at 0，1. 2 and 2. 4 s; on the second row from left to right are the simulation results at 4. 8，6 and 7. 2 s 法相比，建表时间稍长，但是搜索效率高，整体性能有一定的提升．此外，由于算法引入了 Ghost 粒子，在边界处的粒子更为真实、稳定．实验证明该方法不会对模拟结果产生影响．参考文献［1］ Gingold Ｒ A，Monaghan J J． Smoothed particle hydrodynamics: theory and application to non-spherical stars． Mon Not Ｒ Astron Soc，1977，181( 3) : 375 ［2］ Crespo A J C，Dominguez J M，Barreiro A，et al． GPUs a new tool of acceleration in CFD: efficiency and reliability on smoothed particle hydrodynamics methods． PLoS ONE，2011，6( 6) : 1 ［3］ Dominguez J M，Crespo A J C，Gómez-Gesteira M． Neighbour lists in smoothed particle hydrodynamics． Int J Numer Methods Fluids，2010，67( 12) : 2026 ［4］ Chen F G，Ge W，Li H J． The simulation of complex multiphase flow of molecular dynamics implemented on GPU． Sci China Ser · 666 ·