一种改进GPU加速策略在光滑粒子流体动力学方法中的应用.pdf

收稿日期:修回日期:录用日期:基金项目:国家重点研发计划(Y F E ),国家自然科学基金面上项目()作者简介:管延敏(),高级工程师,从事船舶水动力学研究通信作者:杨彩虹,博士,讲师,电话(T e l):;E m a i l:y c h n i h a o c o m文章编号:()D O I:/j c n k i j s j t u 一种改进G P U加速策略在光滑粒子流体动力学方法中的应用管延敏,杨彩虹,康庄,周利(江苏科技大学 船舶与海洋工程学院,江苏 镇江 ;哈尔滨工程大学 船舶工程学院,哈尔滨 )摘要:为解决粒子的无序化易引起的图形处理器(G P U)内存访问冲突问题和提高计算效率,通过建立粒子重排序技术提出了一种改进的G P U加速策略将该加速策略应用于光滑粒子流体动力学(S P H)方法中对三维带障碍物溃坝进行模拟,并与实验结果对比对算法进行验证,获得了较高的计算精度基于此算例,通过在不同硬件设施上进行模拟分别对粒子重新编号的效果和算法的求解效率比较研究结果表明,粒子重新编号技术可以保证稳定的单步运行时间,能够有效解决G P U S P H算法显存访问冲突问题;该G P U加速的S P H并行算法能够大幅提高S P H方法求解效率,随着粒子数量的增加,其大幅缩短计算时间的优势愈发明显,为扩大S P H方法解决三维数值模拟的适用性提供了可能关键词:光滑粒子流体动力学;并行计算;溃坝问题;计算效率中图分类号:U 文献标志码:AA p p l i c a t i o no f a n I m p r o v e dG P UA c c e l e r a t i o nS t r a t e g y f o rt h eS 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 sM e t h o dG U ANY a n m i n,Y ANGC a i h o n g,KANGZ h u a n g,ZHO UL i(S c h o o l o fN a v a lA r c h i t e c t u r ea n dO c e a nE n g i n e e r i n g,J i a n g s uU n i v e r s i t yo fS c i e n c ea n dT e c h n o l o g y,Z h e n j i a n g ,J i a n g s u,C h i n a;C o l l e g eo fS h i p b u i l d i n gE n g i n e e r i n g,H a r b i nE n g i n e e r i n gU n i v e r s i t y,H a r b i n ,C h i n a)A b s t r a c t:I no r d e rt os o l v et h ep r o b l e m o fg r a p h i c sp r o c e s s i n gu n i t(G P U)m e m o r ya c c e s sc o n f l i c t sp o s s i b l yc a u s e db yt h ed i s o r d e ro fp a r t i c l e sa n de n h a n c et h ec o m p u t a t i o ne f f i c i e n c y,a ni m p r o v e dG P Ua c c e l e r a t i o ns t r a t e g y i sp r o p o s e db ye s t a b l i s h i n gp a r t i c l er e o r d e rt e c h n o l o g y T h ea c c e l e r a t i o ns t r a t e g yi sa p p l i e dt ot h es 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)m e t h o dt os i m u l a t et h ed a m b r e a k i n g w i t ho b s t a c l e s i nt h r e ed i m e n s i o n s,a n dt h ea l g o r i t h mi sv e r i f i e db yc o m p a r i n gw i t ht h ee x p e r i m e n t a l r e s u l t s,w h i c ho b t a i n e dah i g hc a l c u l a t i o na c c u r a c y B a s e do nt h i sb e n c h m a r ke x a m p l eo f t h eS P H,t h es t u d i e so nt h ee f f e c to fp a r t i c l er e n u m b e r i n ga n dt h es o l u t i o ne f f i c i e n c yo f t h ea l g o r i t h ma r ec o n d u c t e db yc o m p a r i n gt h es i m u l a t i o n s o f d i f f e r e n t h a r d w a r e f a c i l i t i e s T h e r e s u l t s i n d i c a t e t h a t t h ep a r t i c l e r e o r d e r t e c h n o l o g yc a ne n s u r eas t a b l es i n g l e s t e pr u n n i n gt i m e,a n dc a ne f f e c t i v e l ys o l v et h ep r o b l e mo fg r a p h i cc a r dm e m o r ya c c e s s c o n f l i c t s t h a t c o mm o n l ye x i s t i n t h eG P U S P Ha l g o r i t h m F u r t h e r m o r e,t h eG P Up a r a l l e l a l g o r i t h m第 卷 第期 年月上 海 交 通 大 学 学 报J OUR NA LO FS HANGHA I J I A OT ON GUN I V E R S I T YV o l N o A u g 上海交通大学学报第 卷c a ng r e a t l y i m p r o v e t h es o l u t i o ne f f i c i e n c yo f t h eS P H m e t h o d,a n dw i t ht h e i n c r e a s eo fp a r t i c l en u m b e r,t h ea d v a n t a g eo f d r a s t i c a l l y r e d u c i n g t h e c o m p u t a t i o n t i m eb e c o m e sm o r eo b v i o u s T h em e t h o dp r o p o s e d i nt h i sp a p e rp r o v i d e st h ep o s s i b i l i t yt oe x p a n dt h ea p p l i c a t i o no ft h eS P H m e t h o dt os o l v e Dn u m e r i c a ls i m u l a t i o n s K e yw o r d s: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);p a r a l l e lc o m p u t i n g;d a mb r e a k i n g;c o m p u t a t i o n a le f f i c i e n c y光滑粒子流体动力学(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)对模拟自由表面流动问题具有先天优势,但随着粒子数量的增加,尤其是三维流动问题,计算效率急剧降低,极大限制了该方法在大规模计算中的适用性利用图形处理器(G P U)众核架构开展多核并行计算的计算机统一设备架构(C U D A)技术,其强大的并行特性非常适合解决大规模的高级计算问题,尤其是计算密集型问题 S P H计算能力密集的特性使其能容易地在G P U上实现并行运算,因此一些学者开始尝试将G P U技术应用于S P H方法 H a r a d a等首次描述了S P H方法在G P U上 的 实 现,C r e s p o等将G P U加 速 的S P H方法应用于模拟复杂的自由表面流动他们的研究均表明在单个G P U上的S P H模拟比在单核中央处理器(C P U)上进行的S P H模拟要快两个数量级目前S P H方法在G P U平台下的计算研究发展迅速,G P U加速的S P H模型已应用于颗粒流体流动的模拟和浅水方程的求解 ,并通过结合自适应粒子技术实现了高效计算国内,徐锋在实现了基于G P U众核架构的并行S P H算法的同时,结合所使用的G P U硬件上的特点,对并行算法进行了优化,将计算效率提高了 倍以上;金善勤等提出一种基于粒子对的并行计算方法并将其与改进的S P H方法结合,实现了超过 的加速比;杨志国等将G P U算法应用于二维溃坝模拟中也获得了数量级的加速效果但是,G P U的系统设计理念与C P U正好相反,G P U面对的是类型一致、互不关联的大量数据,其显存的读写对算法求解效率影响很大随着流场的不断演变,S P H粒子的无序化很容易导致多个线程同时读写同一地址而引起访问冲突,极大地影响了计算效率的稳定性针对此问题,并行索引排序方法,、Z索引排序方法 及其改进方法 相继被提出,这些方法通过优化S P H粒子索引存储方式使同单元粒子在G P U显存中尽可能相邻,一定程度上改善了G P U显存访问的不连续问题,但没有从实质上改变单元粒子的无序化为此,本文提出了一种粒子重新编号技术,实现了S P H粒子的有序排列和G P U显存的连续访问,并将该方法应用于三维带障碍物溃坝模拟,通过与实验数据比较以及不同硬件设施上不同算法求解效率的对比,验证了本文方法的精确性和高效性 S P H的控制方程S P H控制方程的离散形式可表示为didtjj(ujui)Wi jVj()duidtjmjpiipjji jWi jg()式中:i、pi、ui、Vj分别为流体粒子i的密度、压力、速度和体积;t为时间;mj和j为粒子i支持域中粒子j的质量和密度;g为重力加速度;Wi j为核函数,本文采用W e n d l a n d 的C 核函数:W(r,h)dq(q)()(q)式中:r为两粒子间距离;h为光滑长度;q为粒子间相对距离,qr/h;d为常数,对于三维问题,d/h人工黏性为i j ci ji ji j,ui jri j,ui jri j()式中:i jh ui jri jri j(h)()为控制速度耗散强度的系数,通常为 ,取 ;平均声速ci j (cicj);i j为平均密度;粒子速度差ui juiuj;粒子间距ri jrirj;(h)为添加项,避免分母为导致数值发散假定流体弱可压缩,采用M o n a g h a n等 提出的人工压缩法求解压强:pcpa t m()式中:为参考密度,取 k g/m;pa t m为大第期管延敏,等:一种改进G P U加速策略在光滑粒子流体动力学方法中的应用 气压力;为常数,;c为人工声速,为确保密度波动低于,本文取c g H,H为水深C U D A G P U架构下的S P H模型 并行算法的实现流程针对现代C UD A G P U架构设计了一种高效的S P H流体数据,以提高S P H流体模拟在G P U上的数据处理速度并行代码使用C语言和C UD A编程语言进行开发,大多数源代码都是C P U和G P U所共有的,可以在C P U或G P U上执行,也可以在没有启用G P U的工作站上运行,只使用C P U实现并行算法实现的流程如图所示,粒子信息的初始数据存储在C P U中,将数据传输到G P U中之后,后续的所有计算都在G P U中进行,即所有涉及粒子循环的任务执行都是由G P U的并行结构来完成的在完成控制方程求解和粒子物理量更新之后,将更新后的粒子物理量从G P U再次传输到C P U进行保存和输出图S P H并行算法的实现流程F i g P a r a l l e l a r i t h m e t i cm o d eo fS P H 粒子搜索S P H方法由于核函数的紧支性,只有粒子支持域内的相邻粒子才会相互作用,因粒子间的相对位置不断变化,每个时刻都需要进行最近相邻粒子的搜索,所以搜索算法极大地影响着S P H的计算效率本文采用M o n a h a n 的链表搜索算法,通过使用控制网格搜索粒子的相邻质点其主要思想为:在计算域建立临时网格,网格大小和粒子支持域大小一致,然后建立每个粒子与所处网格的关联形成关联链表对于三维问题,计算域被划分为立方单元,每个单元与其周围的单元(共 个)组成紧邻单元,如图所示粒子根据它们在区域中的位置被放置在单元格中,在搜索粒子时只需搜索粒子所在的紧邻单元即可,形成最近相邻粒子对图三维相邻粒子链表搜索法示意图F i g S k e t c hf o r t h r e e d i m e n s i o n a ln e i g h b o r l i s t s e a r c hm e t h o d 粒子重新编号技术邻近粒子搜索的G P U并行算法的构建与在C P U上使用的过程有所不同,在G P U上需要建立线程与粒子之间的关联,并为每个粒子单独分配一个计算线程所有粒子按编号排入一个连续的标签数组,并给固定数目的连续粒子分配一个线程块在流场模拟过程中,随着时间步不断推进,S P H粒子分布与初始粒子分布差异越来越大,链表搜索法中S P H粒子不再是顺序分布,如图(a)所示,这就导致线程结束并行访问内存时,极易导致多个线程同时对同一内存地址进行读写操作而产生访问冲突,降低并行算法的求解效率为保证算法效率稳定性,提出了一种粒子重新编号方法,图显示了用于生成按单元格重新编号法的粒子标签数组简单示意图该方法主要包含以下步骤:()遍历所有S P H粒子,计算每个粒子所在网格索引,统计该索引所在网格存储S P H粒子数量,存储于数组I D C o u n t中()创建数组I D B e g i n,记录新粒子编号中每个网格中首粒子编号,对于网格m有I D B e g i nm 上海交通大学学报第 卷图粒子重新编号算法示例F i g E x a m p l eo fn e i g h b o r l i s t r e o r d e rI D B e g i nmI D C o u n tm;()清空数组I D C o u n t;()遍历所有S P H粒子,重新统计该索引所在网格存储S P H粒子数量,根据粒子所在网格索引,对其重新编号,如对于某S P H粒子,其对于网格索引为n,其新的粒子编号为I D B e g i nnI D C o u n tn,其中I D C o u n tn 为统计至该粒子时网格n中存储的S P H粒子数量,最终得到的重新粒子编号效果如图(b)所示结果验证与分析以S P H法验证自由表面流动问题的基准测试案例 溃坝问题为研究对象,对三维带障碍物的溃坝进行数值模拟,通过与实验数据进行对比验证本文并行算法的可靠性和有效性 溃坝模型荷兰海洋研究所(M a r i n)的溃坝试验 被广泛认为是验证S P H自由表面流动的基准算例试验包含一个溃坝流与障碍物的碰撞,如图所示水箱长 m、宽m、高m,水柱被储存在水箱一端,长 m、宽m、高 m,并在试验开始瞬间释放障碍物设置在水流下游,随着挡水墙的拆除,由于重力作用,流体逐渐淹没水箱干床并与障碍物发生碰撞试验通过设置个垂直高度探头(H、H 和H)测量不同位置的水深,具体位置如图所示图试验配置和试验数据的测量位置(m)F i g E x p e r i m e n t a lc o n f i g u r a t i o na n dm e a s u r e m e n tp o s i t i o no f e x p e r i m e n t a l d a t a(m)结果对比数值模拟选取初始粒子间距为 m、流场实粒子总数为 ,时间步长为 m s 为验证本文粒子重新编号方法的准确性,对不同位置处的水深展开比较研究图所示为H、H 和H 位置处不同时刻的水深粒子重新编号前后计算值和实验值的比较图中:hw为水柱高度可见粒子重新编号前后计算值基本一致,不会对精度造成影响;在初始水柱中间(H)、箱 体 中 部(H)、障 碍 物 前 缘 附 近(H),本文计算值与实验值的趋势和大小都很接近图中H 探头清晰地再现了溃坝的整个过程,在最初的s内水柱坍塌,相应地,这段时间的水位也不断下降,而在其他探头中,水流依次到达H 和H 在 s后,反射的水波撞击左墙后反向移动,第次撞击障碍物以及右墙,同时,水深达到第个峰值(H 的 水 深 峰 值 时 刻 大 概 为 s,H 为 s,H 为 s)本文计算的水深与实验水深随时间的变化趋势大致相同,表明本文方法具有良好的计算精度第期管延敏,等:一种改进G P U加速策略在光滑粒子流体动力学方法中的应用 图实验测量和本文模拟的探测点处的垂直水柱高度比较F i g C o m p a r a t i o no f v e r t i c a lw a t e rh e i g h t s a t t h ed e t e c t i o np o i n tm e a s u r e de x p e r i m e n t a l l ya n ds i m u l a t e di nt h i sp a p e r效率测试结果与分析 测试平台通过 比 较C P U串 行、C P U并 行、G P U并 行S P H算法三维带障碍物溃坝数值模拟计算耗时,验证本文所采用G P U算法求解效率为保证计算结果的适用性,分别在不同的C P U、G P U硬件上运行S P H算法,具体配置如表所示表C P U和G P U配置T a b H a r d w a r ec o n f i g u r a t i o n so fC P Ua n dG P U设备类型型号核数主频/GH zC P UI n t e lC o r e i F X e o n G o l d C L G P UG T X s u p e r G e F o r c eR T X G e F o r c eR T X T i G P U算法粒子重新编号效果对比图为实粒子总数 时不同硬件条件下粒子重新编号和未重新编号G P U加速S P H算法单步运行时间(ts)对比图中:sn为计算步数由图可见,本文采用的粒子重新编号算法在不同硬件上都获得了稳定的单步运行时间,而粒子未重新编号时,随着流场中粒子的无序化导致G P U显存访问冲突,其单步运行时间呈对数增长,算例表明本文所采用的粒子重新编号方法可以保证稳定的单步运行时间,是有效的图粒子重新编号和未重新编号单步运行时间对比F i g C o m p a r i s o no fs t e pr u n n i n gt i m eb e t w e e nr e o r d e ra n dn o n r e o r d e rm e t h o d 算法求解效率对比图为 实 粒 子 数 为 时C P U并 行 和G P U并行算法迭代 步单步运行时间对比,可见G P U并行算法都有良好的计算效率,而求解效率稳定性弱于C P U并行算法受S P H方法部分算法、函数间存在串行关系影响,计算效率未能随核数增加而线性增加,以I n t e lC o r ei F为参照,各硬件核数、效率比如表所示,可见相对计算成本(核数效率比)随核数的增加而增大图为不同实粒子数下C P U并行、G P U并行S P H算法单步平均用时(tm)对比图中:np为实粒子数可见随着粒子数量的增加,C P U并行 算法运行时 间显著增加,G P U并行算法大幅缩短计算时间的优势愈发明显 上海交通大学学报第 卷图C P U并行和G P U算法单步运行时间对比F i g C o m p a r i s o no f s t e pr u n n i n gt i m eb e t w e e nC P Up a r a l l e l a n dG P Up a r a l l e l图不同实粒子总数下C P U与G P U运行时间的对比F i g C o m p a r i s o no f r u n n i n gt i m eb e t w e e nC P Ua n dG P Ua td i f f e r e n tn u m b e ro fp a r t i c l e s表各硬件核数和效率比T a b N u m b e ro fc o r e sa n de f f i c i e n c yr a t i oo fd i f f e r e n th a r d w a r e s硬件核数比效率比核数/效率I n t e lC o r e i FX e o n G o l d C L G T X s u p e r G e F o r c eR T X G e F o r c eR T X T i 为进一步验证本文粒子重新编号算法的有效性,对不同实粒子数下G P U算法并行效率与开源软件D u a l S P H y s i c s进行了对比,如图(b)所示,可见同实粒子总数、同硬件设备条件下,本文方法平均单步运行时间均小于D u a l S P H y s i c s软件,算例表明本文粒子重新编号方法具有良好的效率优势结论运用粒子重新编号技术开发了一套高效G P U S P H并行算法,将该算法应用于三维带障碍物溃坝问题,并对算法求解效率进行了比较研究,得到以下结论:()粒子重新编号前后计算值基本一致,不会对精度造成影响,与试验值的对比表明本文所采用的方法精确有效()粒子重新编号技术能够有效解决G P U S P H算法中的显存访问冲突问题()G P U并行算法能够大幅提高S P H方法求解效率,随着粒子数量的增加,其大幅缩短计算时间的优势愈发明显参考文献:刘肃肃,胡祎乐,余音基于G P U的近场动力学模拟的并行化方法J上海交通大学学报,():L I US u s u,HU Y i l e,YU Y i n P a r a l l e lc o m p u t i n gm e t h o do fp e r i d y n a m i c m o d e l sb a s e do n G P UJJ o u r n a lo fS h a n g h a iJ i a oT o n gU n i v e r s i t y,():HA R A D AT,KO S H I Z UKAS,KAWAGU CH IYS 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 so n G P U sCC o m p u t e rG r a p h i c sI n t e r n a t i o n a l P e t r o p o l i s,B r a z i l:C o m p u t e rG r a p h i c sS o c i e t y,:C R E S P O AC,D OM I NGU E ZJM,B A R R E I R O A,e t a l G P U s,an e wt o o l o f a c c e l e r a t i o n i nC F D:E f f i c i e n c ya n dr e l i a b i l i t yo ns 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 sm e t h o d sJ P l o SO n e,():第期管延敏,等:一种改进G P U加速策略在光滑粒子流体动力学方法中的应用 HEY,B AY L Y A E,HA S S AN P OUR A,e t a lAG P U b a s e dc o u p l e dS P H D EM m e t h o df o rp a r t i c l e f l u i df l o ww i t hf r e es u r f a c e sJ P o w d e rT e c h n o l o g y,:X I A X,L I AN G QA G P U a c c e l e r a 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(S P H)m o d e l f o rt h es h a l l o ww a t e re q u a t i o n sJ E n v i r o n m e n t a lM o d e l l i n g&S o f t w a r e,:X I ONG Q,L IB,X UJ G P U a c c e l e r a t e da d a p t i v ep a r t i c l es p l i t t i n ga n dm e r g i n gi nS P HJC o m p u t e rP h y s i c sC o mm u n i c a t i o n,():徐锋基于众核架构的并行S P H算法的研究与实现D上海:上海交通大学,X U F e n g R e s e a r c h a n d i m p l e m e n t a t i o n o f t h es 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 sa l g o r i t h mb a s e do nm u l t i c o r ea r c h i t e c t u r eD S h a n g h a i:S h a n g h a iJ i a oT o n gU n i v e r s i t y,金善勤,郑兴,段文洋基于G P U并行的改进S P H方法对黏性流场的模拟J哈尔滨工程大学学报,():J I NS h a n q i n,Z HE N GX i n g,D UAN W e n y a n g V i s c o s i t yf l o ws i m u l a t i o nu s i n gi m p r o v e dS P H m e t h o db a s e do n G P U p a r a l l e lc a l c u l a t i o nJJ o u r n a lo fH a r b i nE n g i n e e r i n gU n i v e r s i t y,():杨志国,黄兴,郑兴,等G P U在S P H方法模拟溃坝问题的应用研究J哈尔滨工程大学学报,():YANG Z h i g u o,HUANG X i n g,Z HE NG X i n g,e ta lT h ea p p l i c a t i o nr e s e a r c ho f G P U i nt h eS P Hm e t h o dt os i m u l a t et h ed a m b r e a k i n gp r o b l e mJJ o u r n a lo fH a r b i nE n g i n e e r i n gU n i v e r s i t y,():车庆首,李传文,张轶,等GA P I:G P U加速的移动对象并行索引方法J计算机科学与探索,():CHEQ i n g s h o u,L IC h u a n w e n,Z HANG Y i,e ta lGA P I:G P Ua c c e l e r a t e dp a r a l l e lm e t h o df o r i n d e x i n gm o v i n go b j e c t sJ J o u r n a lo fF r o n t i e r so fC o m p u t e rS c i e n c ea n dT e c h n o l o g y,():I HM S E N M,AK I N C IN,B E C K E R M,e ta lAp a r a l l e lS P Hi m p l e m e n t a t i o no nm u l t i c o r eC P U sJC o m p u t e rG r a p h i c sF o r u m,:聂霄不可压缩S P H流体的真实感模拟及其加速技术研究D成都:电子科技大学,N I EX i a o S t u d yo nr e a l i s t i cs i m u l a t i o na n da c c e l e r a t i o nt e c h n i q u e so fi n c o m p r e s s i b l eS P H f l u i d sDC h e n g d u:U n i v e r s i t yo fE l e c t r o n i cS c i e n c ea n dT e c h n o l o g yo fC h i n a,WE N D L AN D H P i e c e w i s ep o l y n o m i a l,p o s i t i v ed e f i n i t ea n dc o m p a c t l ys u p p o r t e dr a d i a l f u n c t i o n so fm i n i m a l d e g r e eJ A d v a n c e s i nc o m p u t a t i o n a lM a t h e m a t i c s,():MONA GHANJJ,G I NG O L D R AS h o c ks i m u l a t i o nb yt h ep a r t i c l em e t h o dS P HJ J o u r n a l o fC o m p u t a t i o n a lP h y s i c s,():MONA GHANJ J P a r t i c l em e t h o d s f o r h y d r o d y n a m i c sJC o m p u t e rP h y s i c sR e p o r t,():K L E E F S MANK M T,F E KK E NG,V E L DMANAEP,e t a l Av o l u m e o f f l u i db a s e ds i m u l a t i o nm e t h o d f o rw a v e i m p a c t p r o b l e m sJ J o u r n a l o fC o m p u t a t i o n a lP h y s i c s,:(本文编辑:王一凡)