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

工程科学学报，第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 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 WOＲDS search algorithm; neighbor; smooth particle hydrodynamics; graphic processing unit; parallel computing; fluid sim￾ulation 收稿日期: 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 architec￾ture) 技术进行了 GPU 上的分子动力学( MD) 模拟． 左 颢睿等［5］提出的基于 GPU 的并行优化技术，优化后累 加和算法的运算速度相比标准并行算法提高了约 34 倍，相比 CPU 串行实现提高了约 70 倍． Goswami 等［6］ 使用 GPU 完成将 SPH 并行化，并且使用了 Z 排序的并 行化算法来降低 GPU 内存开销及访存开销． Domingu￾ez 和 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 和 Te￾schner［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 Flu￾id 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 ar￾chitecture，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) back￾ground 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 num￾bers 粒子数/ 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 ·