（计算机应用技术专业论文）gpu在基于sph方法的混合流模拟中的应用.pdf

g p u 在基于s p i t 方法的混合流模拟中的应用 摘要 6 p u 在基于s p h 方法的混合流模拟中的应用 专业：计算机应用技术 硕士生：熊晨晖 指导老n - 纪庆革副教授 摘要 流体现象在我们的日常生活中扮演着重要的角色，是我们周围环境的一个重 要组成部分，对其进行真实感的模拟，是影视、游戏、计算机仿真等行业非常迫切 的需要。此类研究早期，主要通过构造的方法对一些简单的流体现象进行仿真， 近年来随着图形硬件性能的提高以及人们对图形真实感的追求，该领域学者开始 陆续提出基于物理的流体仿真算法，从而可以用计算机生成更多逼真多样的流体 效果。由于流体的丰富特性，这些算法并不能通用于各种流体现象，随着应用背 景的不同以及仿真环境的改变，这些算法被不断的改进和完善。 流体中复杂而又有趣的现象主要源自于交互过程，然而在计算机图形学界对 流体与静态或动态固体交互较多研究的同时，却忽略了如空气与水或水与油等不 同类型流体之间的交互，而这种混合流体之间的交互却是研究很多重要流体现象 的基础。 本文在总结现有混合流仿真研究方法的基础上，采用一种新颖的基于s p h 的方法对混合流交互现象进行模拟。s p h 方法在流体模拟方面有着许多其它仿真 方法所无法比拟的优点，但这种方法时间复杂度较高，很难保证模拟效果的实时 性。本文所采用方法的新颖之处在于实验中利用了可编程g p u 对模拟过程加速， 通过结合g p u 编程框架与s p h 方法，显著提高了模拟速度。实验结果表明，在 相同粒子数情况下，采用g p u 加速的模拟方法的平均帧率是传统方法的1 0 倍以 上。同时本文还引入流体运动固有的表面张力和混合流体之间的界面张力，较逼 真的模拟出不相溶液体之间的交互情况，从而使得实验效果在达到实时性的同 时，也月 - 匕e , 1 z , k e t 好地满足真实感的要求。 关键词：流体仿真，混合流，s p h ，g p u g p u 在基于s p h 方法的混合流模拟中的应用目录 t h ea p p l i c a t i o no fg p ui ns p h - - b a s e dm i x e df l o w ，1 j 3 l mu l a n o n m a j o r ：c o m p u t e ra p p l i c a t i o nt e c h n o l o g y n a m e ：c h e n h u ix i o n g s u p e r v i s o r ：a s s o c i a t ep r e f e s s o rq i n g g ej i a b s t r a c t f l u i d s ，w h i c ha r ea ni m p o r t a n tc o m p o n e n to fo u rs u r r o u n d i n ge n v i r o n m e n t ， p l a ya ni m p o r t a n tr o l ei no u re v e r y d a yl i f e s ot h e r ei sp r e s s i n gn e e df o rs i m u l a t i n g f l u i d sw i t hh i 曲r e a l i t yi ni n d u s t r yj u s tl i k em o v i e ，c o m p u t e rg a m i n ga n dc o m p u t e r s i m u l a t i o na n ds oo n e a r l i e rr e s e a r c h e sm a i n l ys i m u l a t e ds i m p l ef l u i dp h e n o m e n o n s b yp a r a m e t e rm o d e l i n g i nr e c e n ty e a r s ，谢t l li m p r o v e m e n to ft h eg r a p h i ch a r d w a r e p e r f o r m a n c ea n dm a n sp u r s u i tf o rr e a l i s t i cs e n s eo fg r a p h i c s ，m o r ea n dm o r e r e s e a r c h e r si nt h ef i e l dp r o p o s e dp h y s i c a l l yb a s e ds i m u l a t i o na l g o r i t h m so nf l u i d s ，s o a st og e n e r a t em o r el i f e l i k ea n dd i v e r s ef l u i de f f e c t sb yc o m p u t e r b e c a u s eo ft h e c o l o r f u lc h a r a c t e r i s t i co ff l u i dp h e n o m e n a , t h e s ea l g o r i t h m sc a n n o tb ea p p l i e dt oa l l k i n d so ff l u i dp h e n o m e n o n w i t hd i f f e r e n ta p p l i c a t i o nb a c k g r o u n d sa n dc h a n g e so f t h es i m u l a t i n ge n v i r o n m e n t s ，t h e s ea l g o r i t h m sa r ec o n s t a n t l yi m p r o v e d f l u i dc o m p l e xa n di n t e r e s t i n gp h e n o m e n o ni sm a i n l yd e r i v e df r o mi n t e r a c t i v ep r o c e s s w h i l ei n t e r a c t i o n so ff l u i d s 、 ，i t hs t a t i co rd y n a m i cs o l i d sh a sc a u g h ts o m e a t t e n t i o ni nc o m p u t e rg r a p h i c sl a t e l y , t h em u t u a li n t e r a c t i o no fd i f f e r e n tt y p e so ff l u i d s s u c ha sa i ra n dw a t e ro rw a t e ra n dw a xh a sr e c e i v e dm u c hl e s sa t t e n t i o na l t h o u g h t h e s et y p e so fi n t e r a c t i o na r et h eb a s i sf o ra v a r i e t yo fi m p o r t a n tp h e n o m e n a o nt h eb a s i so fs u m m i n gu pt h ee x i s t i n gm i x e df l o ws i m u l a t i o nr e s e a r c h e s ，w e a d o p tan o v e ls p h - b a s e dm e t h o dt os i m u l a t ef l u i d - f l u i di n t e r a c t i o n s p hm e t h o dh a sa l o to fa d v a n t a g e so t h e rf l u i ds i m u l a t i o nm e t h o d sc a n tm a t c hw i t h ，b u ti t st i m e i l l g p u 在基于s p h 方法的混合流模拟中的应用 目录 c o m p l e x i t yi ss oh i g ht h a tt h es i m u l a t i o nr e a l t i m ec a n tb es a t i s f i e d t h en o v e l t yo f t h ew a ya d a p t e di nt h i sp a p e rl i e si nt h a ti tu s e sp r o g r a m m a b l eg p ut os p e e du p s i m u l a t i o n t h e r ei sa p p a r e n ti n c r e a s ei ns i m u l a t i o ns p e e da f t e rc o m b i n i n gg p u p r o g r a m m i n gf r a m e w o r ka n ds p h m e t h o d t h er e s u l t so fe x p e f i m e n t ss h o w st h a tt h e f r a m er a t ew i t l lg p ui sm o r et h a n10t i m e st h a to fn o r m a ls p hw a y b ym o d e l i n gt h e n a t u r a ls u r f a c et e n s i o no ff l u i dm o v e m e n t sa n dt h ei n t e r f a c et e n s i o nb e t w e e nm i x e d f l o w s ，w ec a na l s os i m u l a t ev i v i df l u i d f l u i di n t e r a c t i o n ，w h i c hw i l lm a k et h ee f f e c t m o r er e a l i s t i cw h i l er e a l t i m e k e y w o r d s ：f l u i ds i m u l a t i o n ，m i x e df l o w , s p h ，g p u i v 论文原创性声明内容： 本人郑重声明：所呈交的学位论文，是本人在导师的指导下， 独立进行研究工作所取得的成果。除文中已经注明引用的内 容外，本论文不包含任何其他个人或集体已经发表或撰写过 的作品成果。对本文的研究做出重要贡献的个人和集体，均 已在文中以明确方式标明。本人完全意识到本声明的法律结 果由本人承担。 学位论文作者签名：礁葭晖 1 日期：2 0 10 年5 月17 日 学位论文使用授权声明 本人完全了解中山大学有关保留、使用学位论文的规定，即： 学校有权保留学位论文并向国家主管部门或其指定机构送 交论文的电子版和纸质版，有权将学位论文用于非赢利目的 的少量复制并允许论文进入学校图书馆、院系资料室被查 阅，有权将学位论文的内容编入有关数据库进行检索，可以 采用复印、缩印或其他方法保存学位论文。 保密论文保密期满后，适用本声明。 学位论文作者签名：髭奸 日期：20 10 年5 月17 日 导师签名：磁 j 日期：20 10 年5 月17 日 g p i j 在基于s p t 方法的混合流模拟中的应用第一章绪论 1 1 问题提出 第一章绪论帚一早殖下匕 目前，在计算机图形学中，物理现象模拟越来越多的在动画制作里得到应用。 即便是复杂的流体运动，也因其千变万化的形态，惟妙惟肖的视觉效果，而得以 在游戏场景、电影及动画制作中被广泛地模拟。 流体模拟可以根据系统所包含的流体种类数目分为几类，但目前大部分的工 作都是针对只有一种流体成分的单元流进行的模拟，而对包含两种流体成分的二 元流研究工作则少之又少，尤其是二元流中的液液交互现象，其作为很多重要现 象衍生的基础，却很少得到关注。 本文采用一种基于s p h 方法的算法对二元流液液交互现象进行模拟，因为 用基于粒子的方法来模拟混合流必须用到大量的粒子，这将严重影响模拟效果的 实时性，因此本文引入了基于物理显卡的g p u 加速方法解决这一难题，使得模 拟出的效果更加美观，更加逼真。 1 2 研究背景 流体现象在我们的日常生活中扮演着重要的角色，是我们周围环境的一个重 要组成部分。在影视制作、计算机游戏以及虚拟仿真等行业，对流体仿真进行真 实感的模拟有着非常迫切的需要。定义一个完整的流体系统，在空间上不仅包括 流体覆盖的求解区域，还包括界定流体求解区域的边界区域。如果排除边界因素 而仅从流体求解系统所含的不同种类流体的数目来看，目前的研究工作可以大致 地分为以下几类： ( 1 ) 单元流( s i n g l e c o m p o n e n tf l u i d ) ，即所要模拟的流体系统中只包含有一 种流体成分。比较常见的现象有水f l2 1 、烟【3 1 和火焰4 1 等等。 ( 2 ) 二元流( b i n a r yf l u i d s ，b i n a r ym i x t u r e s ，t w o f l u i d s ) ，即流体系统中有两种 流体成分，这其中又包括单相二元流与多相二元流【5 1 。多卡l ( m u l t i p h a s e ) 的是指 两种流体成分有两种不同的状态，如气液混和系统，典型的是气泡现象【6 】。单相 g p u 在基于s p h 方法的混合流模拟中的应用 第一章绪论 ( s i n g l e p h a s e ) 的二元流主要就是液- 液混和系统，如油和水的混合现象，奶酪与咖 啡的交融现象等。本文不考虑气液混和的多相流系统，而仅仅以单相的液液混合 系统作为研究对象。 ( 3 ) 多元流( m u l t i c o m p o n e n t sf l u i d s ) ，即系统中包含有三种以上不同流体成 分，如水中的气泡冲击水面的浮油，会导致浮油破碎【5 j 。然而到目前为止，在图形 学领域，作者还没有见到针对这类现象的模拟。当然如果要在流体系统中加入固 体相态，就演变成了另外的一大类问题，即流固耦合t n f 口- j 题。本文只关注纯粹的 流体系统，流固耦合的问题不在本文的考虑范围内。 当前很多流体类的模拟工作都是针对单元流的现象进行模拟，并且已经取得 非常逼真的效果，s t a m 引入半拉格朗日方法( s e m i l a g r a n g i a nm e t h o d ) 来求解 纳维斯多克斯公式( n a v i e r s t o k e se q u a t i o n ( n s e ) ) 中的对流项【8 】，使得n s e 的 稳定求解摆脱了c f l ( c o u r a n tf r i e d r i c h sl e w y ) 条件的限制【9 】，在此基础上j o s s t a m 等人引入“涡流限制( v o r t i c i t yc o n f i n e m e n t ) ，来弥补由于数值耗散而丢 失的细节【l 】，如烟雾中小尺度的翻滚现象，大大增加了烟雾模拟的真实性。f o s t e r 等人采用水平集( l e v e ls e t ) 方法来追踪水的自由表面【3 】，e n r i e 出t 等人则进一步 提出粒子水平集( p a r t i c l el e v e ls e t ) 方法，为水的自由表面增加了更多的细节， 大大改善了自由表面模拟的真实性【4 】。文献【1 0 】对纳维斯多克斯公式( n s e ) 的 质量守恒方程进行完善，引入质量损失项，以此描述底面对水的吸收，模拟了水 在吸收性表面上的流动。更为重要的是，随着半拉格朗日( s e m i l a g r a n g i a n ) 方 法的引入【引，模拟的时间成本大幅下降，于是流体模拟的成果就被广泛用到了电 影等实际领域当中，这进一步激发了实际应用领域对流体模拟的需求。 然而，在流体系统中多加入一种流体，其模拟复杂程度并不是简单地线性增 加，因为加入新的流体之后，不仅仅要考虑到各种流体自身内部的相互作用，还 要考虑到各种流体之间的相互作用。通常流体之间的相互作用是十分复杂的，这 也是为什么到目前为止，针对二元流体现象研究工作相对较少的原因，同时，这 些研究工作也有较多值得改进和完善的地方。 文献 1 1 通过修正质量守恒方程，解决了多相流表面密度不连续的问题，但 是对于长时间的模拟会出现不稳定的结果。文献【1 2 】和文献 1 3 】在他们的研究中 使用了经过修正后的密度求和法。前者在修正的s p h 流体方程中结合对流扩散 2 g p i j 在基于s p i t 方法的混合流模拟中的应用 第一章绪论 方程，集中于复杂环境中可混合流体研究，但是在他们的模型中并未考虑张力的 影响因素，并且只允许在封闭系统的模拟或者具有周期性边界条件的系统模拟中 使用到压力计算。后者的工作主要集中在液滴振动现象和二维剪切流的变形之类 的研究工作中。文献 1 4 1 对这一类的工作在不可压缩性方面进行了扩展。 除了这些文献之外，早期有些文献也针对多相流模拟的不连续问题进行过研 究1 1 5 1 6 ，1 7 1 ，而文献 1 8 ，1 9 ，2 0 ，2 1 贝j j 主要关注气泡泡沫的现象模拟等。 1 3 研究意义 流体现象近年来已经成为基于物理动画制作的一个热门领域。如壁泉、海浪、 火焰、烟云等等，这些流体因其形状千变万化，从而能展现各种形态不一的复杂 效果，也正因如此，流体极大地丰富了虚拟世界带给人们的视觉感观。 流体的定义为：在静态平衡中不能支撑剪应力的物质，或者更直接的说，流 体是因为抵抗变形性而不断流动的物质。为了将流体保持在固定的地方，必须将 它们盛放于水池或其它容器中，这也是为什么在流体模拟中液体与固体边界的交 互研究工作占主导地位的原因了。液固交互现象在很多文献中都有叙述，这里 不再重复。 另一个重要的现象就是不同种类的流体间的交互现象，这类现象在计算机图 形学研究里还没有得到足够的重视，然而，我们可以从液液交互过程中观察到很 多有趣的现象。比如在水沸腾的时候，水不断地从液态变成气态，这是气液交 互现象；当把水倒入杯中的时候，杯中的空气受到流体作用形成气泡；在熔岩灯 中，两种不同极性与密度的流体混合在一起，产生有趣的现象等等。 本文采用一种基于s p h 的方法来模拟液液交互现象。该方法的基本思想是 将连续的流体( 或固体) 用相互作用的质点组来描述，各个物质点上承载各种物理 量，包括质量、速度，位置，密度，颜色等，通过求解质点组的动力学方程和跟 踪每个质点的运动轨迹，求得整个系统的力学行为。从原理上说，只要质点的数 目足够多，就能精确地描述力学过程。该方法是一种纯l a g r a n g e 方法，对于 e u l e r i a n 的基于网格方法地言，模拟混合流或者多相流会相对比较棘手。在这样 的情况下，基于粒子的模拟算法很自然地就成了较理想的选择。在一个粒子系统 中，每个单独的粒子都拥有各自的属性值，具有不同属性值的粒子能够任意地混 g p i j 在基于s p h 方法的混合流模拟中的应用第一章绪论 合，在需要时，粒子可以动态地生成或者从系统当中删除。 然而基于粒子算法一个比较明显的不足在于，当要计算单个粒子受到的力或 其它属性值时，需要对其所有邻近粒子的属性值进行插值求取，因此必须要搜索 与该粒子相邻的其它粒子，这将是一件时间复杂度很高的工作，而且要高仿真的 去模拟混合流液液交互现象，必须用到大量的粒子，这将严重地影响模拟效果的 实时性。本文在之前的研究基础上，引入g p u 物理加速显卡，解决了传统混合 流模拟帧数小，效果迟滞的问题。在流体模拟研究中，利用g p u 结合s p h 方法 来模拟混合流液液交互现象还比较少见，这对于以后影视动画制作与工业生产都 具有十分重要的意义。 1 4 本文的结构 本文的组织如下：第二章主要介绍s p h 方法的原理，以及g p u 编程相关知 识背景的介绍；第三章详细介绍如何应用s p h 方法结合g p u 编程实现对液液 混合流交互现象的模拟；第四章对实验模拟的结果进行综合分析，指出文中所用 方法优势所在及相应的不足；第五章为总结与展望。 4 g p u 在基于s p h 方法的混合流模拟中的应用第二章s p h 仿真原理及g p u 编程 第二章s p h 仿真原理及g p u 编程概述 本文采用基于s p h 的方法来模拟混合流交互现象，并通过g p u 对模拟进行 加速，因此在实验之前有必要对s p h 方法及g p u 编程的知识背景进行相关说明。 在本章前半部分会讲述s p h 方法模拟流体的仿真原理，后半部分会介绍g p u 的 发展及c p u + g p u 异构编程的相关知识。 2 1s p h 仿真原理 光滑粒子动力学( s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ，s p h ) 方法是一种用于求 解连续介质动力学的无网格方法，它通过引入空间场函数和核函数的概念，将基本 方程离散。l u c y l 2 2 1 和g i n g o l d 等人【2 3 1 在1 9 7 7 年首先提出了这种纯拉格朗日流体 动力学方法。起初，该方法用于天体物理学，可处理三维流体自引力问题。后来，l i b e r s k y 等人【2 4 1 将材料强度及完整的应力、应变张量引入到线弹性s p h 方法中。 将该方法用于求解固体动力学问题。对于用一般有限元方法不能处理的问题如边 界不连续问题、大变形问题和冲击破坏问题等，采用s p h 方法可以比较容易求解。 但s p h 方法主要有两个缺陷 2 5 1 ：( 1 ) 自由边界精度难以保证；( 2 ) 在计算中会出现 张力不稳定振荡【2 6 1 。 尽管s p h 方法出现较早，但它在解决一些复杂问题如流体力学、爆炸力学时具 有其他数值方法如有限元方法无法比拟的优点。近来的一些研究，如修正光滑粒 子动力学( c o r r e c t e ds m o o t h e dp a r t i c l eh y d r o d y n a m i c s ) 方法【2 7 1 ，规则光滑粒 子动力学( r e g u l a r i z e ds m o o t h e dp a r t i c l eh y d r o d y n a m i c s ) 方法【2 8 , 2 9 1 ，b 一样条 有限元插值的s p h 方法1 3 0 】等的提出改善了该方法的计算精度和稳定性，并解决了 磁流体动力学( m a g n e t oh y d r o d y n a m i c s ，m h d ) 振动问题【2 8 2 9 ，拉伸和压缩的材料 动力学稳定性问题【3 0 1 ，使s p h 方法的潜力得到进一步发挥。m o n a g h a n 对s p h 方法 作了较全面的介绍，指出了该方法的一些新进展【3 。在国内，s p h 方法的研究也正 逐渐趋于成熟【2 5 , 3 2 , 3 3 1 。 g p l i 在基于s p h 方法的混合流模拟中的应用 第二章s p h 仿真原理及g p u 编程 2 1 1s p h 基本思想 光滑粒子流体动力学，本质上是一种基于粒子系统的插值方法。其基本思想 是用粒子来离散化表示流体，每个粒子代表了一定的流体质量，其携带的各种数 值量将光滑地分布于以该粒子中心为球心，半径为h 的区域内( 该区域称为粒子 的影响域) 。当需要对问题域内某个位置x 的数值量么( x ) 进行采样时，先确定对 该z 处有影响的所有粒子f ( x 位于这些粒子的影响域内，也称这些粒子位于x 支 持域上) ，这些粒子所携带的对应数值量为彳，再对所有4 分布于石处的贡献量 进行积分，就得到了该待采样的数值量么( x ) 。其中很关键的一点是，如何将其 所携带的量光滑地分布于其影响域内，也即如何设计出光滑核函数( s m o o t h i n g k e r n e lf u n c t i o n ) w ( r ，乃) 用以解决了这一分布问题，其中，为粒子中心到影响域上 某点的向量，h 为分布半径。w ( r ，厅) 其应满足这些性质： 1 ) 、对称性，也即其对所携带量分布应以粒子中心为对称中心，呈中心对称 分布。用公式来表示则为：w ( r ，办) = ( 一，办) ； 2 ) 、归一性，也即对所有分布于影响域q 内的数值量进行积分时，其和为 该粒子所携带的数值量总量。用公式表示为：少( ，h ) d r = 1 ； q 3 ) 、当h 趋向于零时具有狄拉克函数性质。 ! i m w ( x - x ，h ) = a ( x - x ) h - - o 4 ) 、紧支性，即对粒子影响域以外的分布量为零，用公式表示为： w ( r ，办) = o 当h 办时。 流体域 l 一一一一一一一一一一一一一一一一一一一一一一一一一一 矿十 影响域 + 一卜。一。1 图2 1 光滑核函数曲线示意图f 3 4 1 6 g p i 在基于s p h 方法的* 城模拟中 月第二章s p h 仿真月8 厦c 叫犏程 州瑚= 击 o - 2 5 ( 2 圳3 1-1 15 d7 5 d0 d d s 2 1 i 一2 + o 3 1 5 0 d 2 式中 一：h 此。 当光滑核函数确定以后，在连续的情况下位于流体域内位置z 的数值量 a ( x ) 的值可以积分得到公式( 21 ) ： a ( x ) * i a ( x 7 0 一， ) 出( 2 1 ) o 其梯度及拉普拉斯运算形式分别为公式( 22 ) 和( 2 - 3 ) ： v a ( x ) ； a ( x ) v w ( x z ， ) 出( 2 2 ) 二 v 2 一( x ) zh ( ，2 w ( xz 7 ，h ) d x ( 2 _ 3 ) o 图2 - 2 光滑粒子示意图 对该连续形式转化为支持域内所有粒子叠加求和的离散化形式，粒子及其 影响域如图2 - 2 所示。若使用粒子的体积a r 来取代连续积分中粒子i 处无穷小 体元d r ，则粒子的质量啊可以表示为啊= n ，其中n 为粒子i 的密度， ( i = 1 , 2 ，n ) ，n 为在z 支持域内的粒子总量，可阻得到一( x ) 的离散化粒子近 似公式( 24 ) ： 爿：p ( x ，似一， ) 出* 兰爿( ) 矿( x - x , , h ) a v , ( 2 _ 4 ) = 姜一( ) ( z 一王，柳去( 一k ) ( 2 5 ) 即 似) = 姜詈训x - x , h ) ( 2 卸 g p u 在基于s p h 方法的混合流模拟中的应用 第二章s p h 仿真原理及g p u 编程 相应的梯度及拉普拉斯运算的粒子近似式为( 2 - 7 ) 和( 2 - 8 ) ： v a ( 班喜罢4 v 形( x - x ，, h ) ( 2 - 7 ) v 2 砸) = 喜罢矿形( x - x i , h ) ( 2 - 8 ) 对于每个粒子，其质量是预先定义的，因此每个粒子的密度可由公式( 2 - 9 ) 得到： p ( x ) = 喜鲁以肜( x - x j , h ) = 善n 形( x - x j , h ) ( 2 _ 9 ) 2 1 2s p h 基本方程 连续介质流体动力学的守恒方程为公式( 2 - 1 0 ) ： 竺：刃u 俄 掣：土v 盯 ( 2 1 0 ) a t p 塑：一1 盯：v p u 公式( 2 1 0 ) 中：p 是密度；u 是速度向量；万是应力张量；e 是能量。 对公式( 2 1 0 ) 中的三个方程进行变换，则连续介质的守恒方程变成s p h 的粒 子方程，有公式( 2 11 卜( 2 1 3 ) ： 警一z ，m j ( u j 删 ( 2 1 1 ) 百d u , 一手( 考一v ( 2 - 1 2 ) 生d t 手哆考( 删叫 ( 2 - 1 3 ) 式中，= v 形( 薯一t ，厅) 。在上面三组公式中，对于一个特定的粒子f ，它的密 度随时间变化由公式( 2 1 1 ) 决定，加速度由公式( 2 1 2 ) 决定，内部能量变化 由公式( 2 1 3 ) 决定。 8 g p u 在基于s p i t 方法的混合流模拟中的应用第二章s p h 仿真原理及g p u 编程 2 1 3 关于s p h 粒子支持域与影响域 在前面论述中本文用到了支持域及影响域的概念，因为这两个概念极易混 淆，故在此有必要做进一步区分。 在无网格方法的广义范围内，支持域是这样定义的：场点x = ( x ，y ，z ) 的支持 域是这样的一个领域，即域内所有的点的信息都被用于决定在x 处的点的信息。 影响域则是一个点对周围产生影响的范围。因此，在无网格方法中影响域只与节 点相关，而支持域则是与任意位置点x 相关，其域内可能只有一个节点，但也可 能有许多个节点。 由以上定义可得： 当应用支持域的概念时，我们所考虑的问题是基于场点x 上的；当应用影响 域的概念时，我们所考虑的问题只基于节点上。 若节点在位置点x 的支持域内，且位置点x 受到节点f 的影响，则位置点x 也 在节点2 的影响域内。 若场点x 刚好也是节点j ，则此点将同时具有支持域和影响域。在这种情况下， 节点z 的影响域可能与场点x 的支持域相同。 作为一种粒子法，s p h 方法只在粒子上对场变量进行近似计算。换句话说， 即位置x 永远位于粒子( 节点) 上。因此，在s p h 方法中，粒子同时具有支持 域和影响域。 在s p h 方法中，任一粒子的支持域和影响域都与粒子的光滑长度h 紧密相 关。光滑粒子h 与因子k 的乘积即确定了光滑函数作用的支持域或影响域。光滑 长度可以随时间或空间的变化而变化，在一些情况下，光滑长度也可以随维数变 化，即在1 维中是变量，2 维中是矢量，3 维中是张量。因此，在s p h 中粒子的 支持域和影响域直接与粒子的领域砌相关。 本文之所以会对这两种不同领域进行讨论是因为在s p h 方法的历史中，这 两种领域的概念产生了两种不同的粒子近似模型。在s p h 方法中使用支持域和 影响域的概念可分别产生对粒子f 处的场变量进行近似的分散模型和集中模型。 在分散模型内，通过应用那些影响域覆盖到粒子f 的粒子来对粒子f 进行近似计 9 g p u 在基于s p h 方法的混合流模拟中的应用第二二章s p h 仿真原理及g p u 编程 算。在集中模型内，通过应用粒子f 的支持域内的粒子来对粒子f 进行近似计算。 由于两个粒子的光滑长度不一定相同，因此有可能发生如下情况，即粒子f 位 于粒子，的影响域内，但却不影响粒子，的求解计算。这种情况常发生在不涉及 粒子性质，只使用这些区域对场变量进行插值计算的无网格方法中。然而，对于 无网格粒子法，如s p h ，此不平衡或非对称影响可能会导致求解失真【3 5 】。在s p h 方法中，若矗比h ，小，则粒子，的影响域可能覆盖粒子f ，但粒子f 的影响域不一 定覆盖粒子，。因此，有可能粒子，对粒子f 施加了一个作用力，而粒子f 没有对 应的施加力于粒子上。这明显地与牛顿第三定律相违背。为了解决此问题，常应 用一些均值法求解两个相互作用粒子的光滑长度，并由此来获取两粒子的支持域 和影响域。在这种情况下，上述所提的分散模型将与集中模型一样。因此，s p h 中粒子的支持域和影响域实际上是一样的。 2 2g p u 编程概述 目前，主流计算机中的处理器主要是中央处理器c p u 和图形处理器g p u 。 传统上，g p u 只负责图形渲染，而大部分的处理都交给了c p u 。 然而，2 l 世纪人类所面临的重要科技问题，如卫星成像数据的处理、基因 工程、全球气候准确预报、核爆炸模拟等，数据规模已经达到t b 甚至p b 量级， 没有万亿次以上的计算能力是无法解决的。与此同时，我们在日常应用中( 如游 戏、高清视频播放) 面临的图形和数据计算也越来越复杂，对计算速度提出了严 峻挑战。 g p u 在处理能力和存储器带宽上相对c p u 有明显优势，在成本和功耗上也 不需要付出太大代价，从而为这些问题提供了新的解决方案。可编程的g p u 已 发展成为一种高度并行化、多线程、多核的处理器，具有杰出的计算功率和极高 的存储器带宽，如图2 3 所示。 g p u 其实是由硬件实现的一组图形函数的集合，这些函数主要用于绘制各 种图形所需要的运算。这些和像素，光影处理，3 d 坐标变换等相关的运算由 g p u 硬件加速来实现。图形运算的特点是大量同类型数据的密集运算如图 形数据的矩阵运算，g p u 的微架构就是面向适合于矩阵类型的数值计算而设计 1 0 g p u 枉基于s p h j 法的* 台流模拟中的应用第= 章s p h 仿真原理肚g p u 编程 的，大量重复设计的计算单元，这类计算可以分成众多独立的数值计算大量 数值运算的线程而且数据之间没有像程序执行的那种逻辑关联性。 1 一 ：，一 g 7 1 匹妄兰之兰苎：堂 7 t 一 7 。p 7 ；。一“ 一 图2 - 3c p u 和g p u 的每秒浮点运算次数和存储器带宽 c p u 和g p u 之间浮点功能之所以存在这样的差异，原因就在于g p u 专 为计算密集型、高度并行化的计算而设计图2 3 显示的正是这种情况，因而 g p u 的设计能使更多晶体管用于数据处理，而非数据缓存和流控制，如图2 - 4 所示。 cpu g p u 圈2 4g p u 中的更多晶体管用于数据处理 c3|。_d g p u 在基于s p h 方法的混合流模拟中的应用 第二章s p h 仿真原理及g p u 编程 更具体地说，g p u 专用于解决可表示为数据并行计算的问题在许多数 据元素上并行执行的程序，具有极高的计算密度( 数学运算与存储器运算的比 率) 。由于所有数据元素都执行相同的程序，因此对精密流控制的要求不高；由 于在许多数据元素上运行，且具有较高的计算密度，因而可通过计算隐藏存储器 访问延迟，而不必使用较大的数据缓存【3 6 1 。 数据并行处理会将数据元素映射到并行处理线程。许多处理大型数据集的应 用程序都可使用数据并行编程模型来加速计算。在3 d 渲染中，大量的像素和顶 点集将映射到并行线程。类似地，图像和媒体护理应用程序( 如渲染图像的后期 处理、视频编码和解码、图像缩放、立体视觉和模式识别等) 可将图像块和像素 映射到并行处理线程。实际上，在图像渲染和处理领域之外的许多算法也都是通 过数据并行处理加速的从普通信号处理或物理仿真一直到数理金融或数理 生物学。 2 2 1c p u + g p u 异构并行 传统的c p u + g p u 异构并行处理的典型任务是图形实时渲染。在这类应用 中，c p u 负责根据用户的输入和一定的规则( 如游戏的a i ) 确定在下一帧需要 显示哪些物体，以及这些物体的位置，再将这些信息传递给g p u ，由g p u 绘制 这些物体并进行显示。两者的计算是并行的：在g p u 绘制当前帧的时候，c p u 可以计算下一帧需要绘制的内容。 在这些处理中，c p u 负责的是逻辑性较强的事务计算，g p u 则负责计算密 集度高的图形渲染。为了满足事务计算的需要，c p u 的设计目标是使执行单元 能够以很低的延迟获得数据和指令，因此采用了复杂的控制逻辑和分支逻辑，以 及大量的缓存来提高执行效率；而g p u 必须在有限的面积上实现很强的计算能 力和很高的存储器带宽，因此需要大量地执行单元来运行更多相对简单的线程， 在当前线程等待数据时就切换到另一个处于就绪状态等待计算的线程【3 7 1 。简而言 之，c p u 对延迟将更敏感，而g p u 则侧重于提高整体的数据吞吐量。 2 2 2 从g p g p u 到c u d a 传统上的g p u 应用被局限于处理图形渲染计算任务，无疑是对计算资源的 1 2 g p u 在基于s p h 方法的混合流模拟中的应用第二章s p h 仿真原理及g p u 编程 极大浪费。随着g p u 可编程性的不断提高，利用g p u 完成通用计算的研究逐渐 活跃起来。将g p u 用于图形渲染以外领域的计算称为g p g p u ( g e n e r a l p u r p o s e c o m p m i n go ng r a p h i c sp r o c e s s i n gu n i t s ，基于g p u 的通用计算) 。g p g p u 计算通 常采用c p u + g p u 异构模式，由c p u 负责执行复杂逻辑处理和事务管理等不适 合数据并行的计算，由g p u 负责计算密集型的大规模数据并行计算。这种利用 g p u 强大处理能力和高带宽弥补c p u 性能不足的计算方式在发掘计算机潜在的 性能，在成本和性价比方面有显著优势。但是，传统的g p g p u 受硬件可编程性 和开发方式的制约，应用领域受到了限制，开发难度也很大。 2 0 0 7 年，n v i d i a 推出了c u d a ( c o m p m e u n i f i e dd e v i c ea r c h i t e c t u r e ，统 一计算设备架构) 。c u d a 是一种将g p u 作为数据并行计算设备的软硬件体系， 到目前为止，c u d a 的版本已经进行了三次版本提升，功能不断完善，能够很好 地支持新硬件的特性。 c u d a 不需要借助于图形学a p i ，并采用了比较容易掌握的类c 语言进行开 发。其核心有三个重要抽象概念：线程组层次结构、共享存储器、屏蔽同步( b a r r i e r s y n c h r o n i z a t i o n ) ，可轻松将其作为c 语言的最小扩展级公开给程序员，保证熟 悉c 语言等标准编程语言的程序员能够迅速掌握c u d a ，从而比较平稳地从 c p u 过渡到g p u ，而不必重新学习语法。当然，要开发高性能的g p u 通用计算 程序，开发人员仍然需要掌握并行算法和g p u 架构方面的知识。 与以往的g p u 相比，支持c u d a 的g p u 在架构上有了显著的改进，这两 项改进使c u d a 架构更加适用于g p u 通用计算。一是采用了统一处理架构，可 以更加有效地利用过去分布在顶点渲染器和像素渲染器的计算资源；二是引入了 片内共享存储器，支持随机写入和线程间通信【3 7 1 。 c u d a 为开发人员有效利用g p u 的强大性能提供了条件。自推出后，c u d a 被广泛应用于石油勘探、天文计算、流体力学模拟、分子动力学仿真、生物计算、 图像处理、音视频编解码等领域，在很多应用中获得了几倍、几十倍、乃至上百 倍的加速比【3 刀。 2 2 3c u d a 软件体系 c u d a 的软件堆栈由三层构成：c u d al i b r a r y 、c u d ar u n t i m ea p i 、c u d a g p u 在基于s p h 方法的混合流模拟中的应用第二章s p h 仿真原理及g p u 编程 d r i v e ra p i ，如图2 5 所示。c u d a 的核心是c u d ac 语言，它包含对c 语言的 最小扩展集和一个运行时库，使用这些扩展和运行时库的源文件必须通过n v c c 编译器进行编译。 匣二二二二二二 图2 5c u d a 软件体系 c u d ac 语言编译得到的只是g p u 端代码，而要管理g p u 资源，在g p u 上分配显存并启动内核函数，就必须借助c u d a 运行时a p i ( r u n t i m ea p i ) 或者 c u d a 驱动a p i ( d r i v e ra p i ) 来实现。在一个程序中只能使用c u d a 运行时a p i 与c u d a 驱动a p i 中的一种，不能混合使用。 2 2 4c u d ac 语言 c u d ac 语言为程序员提供了一种用c 语言编写设备端代码的编程方式， 包括对c 的一些必要扩展和一个运行时哭库。本节主要介绍语言上的扩展， c u d a 对c 的扩展主要包括一下几个方面： 引入了函数类型限定符。用来规定函数是在h o s t 还是在d e v i c e 上执行， 以及这个函数是从h o s t 调用还是从d e v i c e 调用。这些限定符有： 1 4 g p u 在摹于s p h 方法的混合流模拟中的应用第二章s p i - ! 仿真缴理及g p u 编程 _ d e v i c e 、_ h o s t _ 并g l _ g l o b a l _ 。 引入了变量类型限定符。用来规定变量被存储在哪一类存储器上。传统 的在c p u 上运行的程序，编译器能自动决定将变量存储在c p u 的寄存 器还是内存中。在c u d a 编程模型中，一共抽象出多达8 种不同的存储 器。为了区分各种存储器，必须引入一些限定符，包括：d e v i c e 、s h a r e d 和_ c o n s t a n t 。注意，此处的_ d e v i c e _ 与之前说的一d e v i c l 限定符含义不 同。 引入了内置矢量类型。如c h a r 4 、u s h o r t 3 、d o u b l e 2 、d i m 3 等，它们是由 基本的整型或浮点型构成的矢量类型，通过x 、y 、z 、w 访问每一个分 量，在设备端代码中各矢量类型有不同