基于CUDA的SMAC算法并行化

常立博，杜慧敏，韓俊剛

(1.西安郵電大學 電子工程學院,陜西 西安 710121; 2.西安郵電大學 計算機學院,陜西 西安 710121)

?

基于CUDA的SMAC算法并行化

常立博1，杜慧敏1，韓俊剛2

(1.西安郵電大學 電子工程學院,陜西 西安 710121; 2.西安郵電大學 計算機學院,陜西 西安 710121)

改進SMAC(Simplified Marker and Cell)算法，增強其對流體模擬的實時處理能力。采用點差分格式對求解壓力場和速度更新的偏微分方程進行離散化；引入消除數(shù)據(jù)相關性的存儲算法以減少數(shù)據(jù)傳輸，并借助分層存儲機制提高訪存比，采用并行歸約增加線程并行度；在統(tǒng)一計算設備架構平臺下，對離散化的SMAC算法進行改進、優(yōu)化及并行化實現(xiàn)。純粹計算及多次迭代模擬實驗結果顯示，改進算法提速明顯，可實現(xiàn)對一般場景的實時模擬。

計算流體力學；統(tǒng)一計算設備架構；并行算法

流體流動現(xiàn)象的模擬不僅是計算機動畫、游戲、影視和飛行模擬等領域不可或缺的元素，也是科學計算可視化領域的研究熱點之一[1]。常用的流體模擬算法有基于Stable Fluids[2]和基于SMAC[3]兩類。Stable Fluids算法雖然可以較好地完成一般場景的模擬工作，但其本身存在嚴重的數(shù)值耗散，從而失去了模擬的物理真實感[1]。SMAC是MAC(Maker and Cell)方法的一種改進，具有良好的數(shù)值特性、求解精度和計算效率，已廣泛應用于計算機圖形學領域[4]。但基于SMAC的流體模擬算法存在復雜度較高、計算耗時過長等問題，無法用于實時性要求較高的場合。以提高計算密集型算法速度而發(fā)展起來的圖形處理器(Graphic Processing Unit, GPU)，為提高計算機動畫實時模擬速度提供了新的解決方法，尤其是統(tǒng)一計算設備架構(Compute Unified Device Architecture, CUDA)平臺的推出，進一步發(fā)揮了GPU的數(shù)據(jù)級并行計算能力，其強大的并行計算能力和浮點運算能力，以及可編程能力使得流體模擬加速成為可能。

文獻[5]利用CUDA提供的庫函數(shù)實現(xiàn)流體模擬的方法，模擬結果質量較好，但具體的實現(xiàn)原理并沒有說明；針對求解偏微分方程的迭代算法的并行化設計方法,采用GPU實現(xiàn)了算法加速，雖然最高已達到9倍的加速，但依然不能滿足實時流體模擬的需求[6]。文獻[7]通過對SMAC算法離散化過程進行分析，為該算法并行化提供了理論依據(jù)，基于OpenGL和Cg語言在NV40顯卡上得以實現(xiàn)，但限于當時的軟硬件條件，很多步驟并沒有在GPU上執(zhí)行，導致加速比都不是很高，無法達到實時處理的要求。

本文擬采用點差分方法完成SMAC算法的離散化操作，討論該算法在CUDA平臺上實現(xiàn)時的數(shù)據(jù)存儲算法、數(shù)據(jù)訪存機制及線程分配，采用并行SMAC算法實現(xiàn)流體實時模擬。

1　SMAC算法原理及CUDA平臺概述

1.1SMAC算法原理

流體模擬的基本方程N-S方程[8]采用交錯網(wǎng)格離散化整個求解域。二維N-S方程的無量綱表達式為[7]

(1)

(2)

(3)

其中，x和y分別表示流場的橫坐標和縱坐標，p表示壓力值，fx和fy分別表示流場在x和y方向受到的外力，u和v分別表示在x和y方向的流場速度，Re表示雷諾系數(shù)。

針對每個模擬節(jié)點，令中間變量F和G表示為

(4)

(5)

其中δt表示時間步進長度，n為模擬的迭代步數(shù)，將式(4)和式(5)分別代入式(1)和式(2)，可得到

(6)

(7)

為保證速度場滿足無源條件，引入壓力投影映像來修正速度場，將式(6)和式(7)代入式(3)，得到壓力項的求解方程為

(8)

顯式求解SMAC算法必須滿足式柯朗-弗里德里希斯-列維(Courant-Friedrichs-Lewy Condition，CFL)條件[9]，即最大速度值滿足

(9)

|umax|δt<δx,|vmax|δt<δy。

(10)

其中δx和δy分別表示x和y的步進長度，umax和vmax分別表示速度場在x和y方向的極大值。所以，在流場中計算更新仿真步長，其核心是找到速度場中的極大值。

1.2CUDA并行設計架構

CUDA將GPU視為一個并行數(shù)據(jù)計算的設備，對所進行的計算進行分配和管理。CUDA主要有線程執(zhí)行模式和存儲組織形式兩個顯著特點[10-11]。

(1)線程執(zhí)行模式

CUDA采用單指令多線程(Single Instruction Multiple Threads, SIMT)的執(zhí)行模式，即每個處理單元(Streaming Processor, SP)執(zhí)行1個線程、多個SP組成1個流處理器簇(Streaming Multiprocessor, SM)。每個SM的SP單元以連續(xù)的32個線程(稱為1個warp)為單位來創(chuàng)建、管理、切換和執(zhí)行指令，調試器不檢查指令流內(nèi)的依賴性，所以1個線程塊內(nèi)的線程數(shù)最好是32的整數(shù)倍。

(2)存儲組織形式

CUDA結構中存儲單元可分為全局存儲器、常數(shù)存儲器、共享存儲器和寄存器，如圖1所示。所有線程都可讀寫1塊相同的全局存儲器和1個只供讀取的常數(shù)存儲器；每個線程塊都擁有1個共享存儲器，同一線程塊內(nèi)的所有線程可以讀寫本塊內(nèi)的共享存儲；每個線程都有1個容量較小但速度快的私有寄存器。由于訪問共享存儲比訪問全局存儲快數(shù)百倍，所以針對同一個線程塊內(nèi)的線程之間的通信，應當盡可能地使用共享存儲以減少數(shù)據(jù)訪存延時。

圖1　CUDA存儲結構

2　SMAC算法改進

2.1SMAC算法離散化求解

采用點差分格式離散純二次項及混合二次項[12],則速度場在x和y方向二階偏導的具體差分格式可表示為

(11)

(12)

式中χ為人工粘性系數(shù)(χ=0.9)。

(13)

速度場主離散求解格式為

(14)

(15)

2.2SMAC算法并行實現(xiàn)步驟

根據(jù)SMAC離散化的過程，實現(xiàn)SMAC算法具體步驟如下。

步驟1初始化流場、邊界條件和流場中障礙物分布。

步驟2初始化速度場和壓力場。

步驟3計算每一個節(jié)點的F,G。

步驟4根據(jù)式(13)計算每一個節(jié)點的p。

步驟5根據(jù)式(14)和式(15)更新速度場。

步驟6根據(jù)式(9)和式(10)重新仿真步長。

步驟7根據(jù)壓力場計算流場的密度場。

步驟8繪制密度。

如果速度場變化大于閾值或超過仿真時間則結束模擬過程，否則從步驟3開始下一次迭代。

可以看出，計算各節(jié)點的F，G和p參數(shù)只與相鄰的4個節(jié)點和本節(jié)點的上一次迭代狀態(tài)相關，而與以前的狀態(tài)無關，具有時空局部性，可進行并行計算，從而減少計算時間，提高算法的運行速度。

3　基于CUDA加速的SMAC流體模擬

采用CPU和GPU協(xié)同的方式完成流體模擬。在CPU中設定流體的邊界條件，將初始化場景數(shù)據(jù)傳輸?shù)紾PU存儲中，并將每一格點上的計算任務分配給一個計算核上，當?shù)竭_時間閾值時完成流體模擬計算；最后將模擬完成的場景數(shù)據(jù)傳回CPU端，并利用圖形渲染混合接口實現(xiàn)流體模擬的可視化操作，具體過程如圖2所示。

圖2　CUDA加速的SMAC流體模擬流程

3.1數(shù)據(jù)預處理

在建立場景時，對塊間的ghost數(shù)據(jù)(邊緣填充數(shù)據(jù))進行預處理：在各塊的四周用相鄰塊的數(shù)據(jù)進行填充，各塊在計算之前，將填充好的ghost數(shù)據(jù)和本分塊數(shù)據(jù)一同合并讀入共享存儲。中央塊的ghost數(shù)據(jù)將由它周圍的8個數(shù)據(jù)塊中對應的數(shù)據(jù)組成。從而減少在模擬迭代過程中進行全局存儲訪問的操作次數(shù)。數(shù)據(jù)填充如圖3所示。

圖3　數(shù)據(jù)填充

假設流體模擬場景m×n，將其分解為T×T大小的塊。原始場景用I(i,j)表示，變換后的場景用O(i,j) 表示，則流體模擬場景數(shù)據(jù)預處理過程可劃分為兩個主要部分，偽代碼如下。

(1)對邊界及ghost元素的處理

for (i=0;i<(m+2(m/T-1)+2);i++)

for (j=0;j<(n+2(n/T-1)+2);j++)

if(i==0)

O(i,j) =左邊界；

else if(j==0)

O(i,j) =上邊界；

else if(i==(m+2(m/T-1)+1))

O(i,j) =右邊界；

else if(j==(n+2(n/T-1)+1))

O(i,j) =下邊界；

else if(m%((i-1)*T)==0) {

if(n%((j-2)*T)==0)

O(i,j) = I(i-1,j+1);//右上角ghost元素

else if(n%(j-1)*T)==0)

O(i,j) = I(i-1,j-1); //左上角ghost元素

else

O(i,j)= I(i-1,j); //上邊ghost元素

}

else if(m%((i-2)*T)==0) {

if(n%((j-2)*T)==0)

O(i,j) = I(i+1,j+1);//右下角ghost元素

else if(n%(j-1)*T)==0)

O(i,j) = I(i-1,j+1); //左下角ghost元素

else

O(i,j) = I(i+1,j); //下邊ghost元素

}

else if(n%((j-1)*T)==0)

O(i,j) = I(i,j-1);//左邊ghost元素

else if(n%((j-2)*T)==0)

O(i,j) = I(i,j+1);//右邊ghost元素

(2)對真實場景數(shù)據(jù)的處理

for(k=0;k

for(l=0;l

for(j=0;j

for(i=0;i

O(k*T+j+1,l*T+i+1) =I(k*T+j,l*T+i);

}

3.2數(shù)據(jù)傳輸單元設計

在場景初始化時將速度場、壓力場和受力場由CPU傳輸?shù)紾PU，完成場景模擬之后再將速度場和壓力場由GPU傳輸?shù)紺PU；而數(shù)據(jù)填充、對各節(jié)點的速度場和壓力場計算及更新均在GPU端完成，從而減少數(shù)據(jù)在顯存和主存中傳輸?shù)拇螖?shù)。

3.3核心計算單元設計

為了避免同1個warp中出現(xiàn)分支，使用1個線程計算1個節(jié)點的F、G、p信息及速度更新操作，將16×16或32×32個線程組成1個線程塊，并根據(jù)模擬場景自動確定線程塊數(shù)量。如果要模擬的場景不規(guī)則，則在進行模擬之前先對場景進行無用數(shù)據(jù)填充，再進行場景模擬。

利用共享存儲作為數(shù)據(jù)“中轉站”以降低非合并訪問。同時經(jīng)過數(shù)據(jù)填充處理過的狀態(tài)信息，塊與塊之間的信息是完全不相關的。利用CUDA提供的訪存延時隱藏技術，即在處理過程中，先啟動所有塊對全局存儲訪問，再進行計算，這樣幾乎可以隱藏所有的訪存延時。

3.4求解迭代步長單元設計

方程解隨著時間步的遞進而更新，直到計算時間溢出。設定空間步長為1，由式(10)確定各時間步之間的時間推進量。采用并行歸約方法計算每一個時間步的速度最大值[12]，具體如圖4所示。假如對N個數(shù)據(jù)進行求極大值操作，開啟N/2個線程分別進行兩兩比較大小，如0號線程求第0和第N/2+1個數(shù)據(jù)中的較大值，結果寫入第0個數(shù)據(jù)的存儲空間，同理可推其他線程計算的數(shù)據(jù)；每個循環(huán)參與計算的線程數(shù)減半，最后，0號單元的數(shù)據(jù)和1號單元的數(shù)據(jù)比較以計算出極大值，此算法將求極大值運算復雜度降為log2(n)。在實現(xiàn)過程中，首先在同一個線程塊內(nèi)完成速度場每行極大值的計算，然后將各線程塊的極大值(即各行的極大值)放在0號線程塊中計算出速度極大值。

圖4　并行歸約求極大值

4　實驗方法與結果分析

為了驗證改進后的SMAC算法在流體模擬中的有效性和效率，采用基于Intel i7 5500U的CPU、基于nVIDIA GTX560處理器的GPU和分別基于gcc 和nvcc3.0平臺的編譯器，分別完成了基于SMAC算法流體模擬的串行及并行實現(xiàn)。并且，對不同計算規(guī)模下的單次迭代及相同計算規(guī)模下多次迭代的性能進行了分析。

4.1單次迭代的性能分析

表1為模擬1個時間步長的運行時間分布，在不同計算規(guī)模下，CPU耗費的時間同GPU耗費的時間對比。其中GPU計算時間只包括GPU上由原始速度場u、v，力場fx、fy，初始壓力場p計算得到1個時間步的速度場和壓力場的時間，不包括GPU和CPU之間數(shù)據(jù)傳輸時間，即純粹計算時間。

表1　單步長的加速時間分布

從表1可看出，當模擬的流場較大(1 024×1 024)時，采用CPU實現(xiàn)的串行計算方法每秒僅可完成約6次流體場景計算，無法實現(xiàn)實時流體場景模擬。然而，采用并行SMAC算法每秒僅可完成數(shù)百次對流體場景模擬計算，所以完全可實現(xiàn)實時流體場景模擬。

4.2多次迭代性能分析

表2為計算規(guī)模(512×512)一定時，不同迭代次數(shù)CPU的時間同GPU的時間對比。其中GPU計算時間表示GPU+CPU的計算總時間，包括將原始速度場u、v，力場fx、fy，初始壓力場p由CPU傳輸?shù)紾PU，和將模擬后的速度場u、v、壓力場p傳輸回CPU的傳輸時間，和在GPU上多次迭代的運算時間。CPU計算時間表示CPU串行模擬時間。

表2　多步長的加速時間分布

從表2可以看出，針對中等規(guī)模(512×512)的流體模擬場景，在場景不是十分復雜(即在迭代次數(shù)小于100時求解方程收斂)的情況下，采用GPU并行實現(xiàn)方法模擬速度可達到20～30 fps，基本可滿足實時流體模擬的要求。

從實驗結果可以看出，采用點差分離散方法對SMAC算法進行離散化，并采用數(shù)據(jù)填充算法及各種優(yōu)化計算方面可以完成對一般流體場景進行實時模擬的工作。

5　結語

采用點差分離散方法實現(xiàn)SMAC算法中求解壓力場、速度更新等操作的并行化處理。然后，從減少數(shù)據(jù)傳輸、提高訪存比及增加線程并行度等方面，對離散化的SMAC算法在CUDA平臺下進行改進和優(yōu)化。最后，將改進后的算法在CUDA架構下實現(xiàn)。實驗結果表明，基于CUDA架構的SMAC并行算法提速效果明顯，可實現(xiàn)對一般場景的實時模擬工作。

[1]佟志忠, 姜洪洲, 韓俊偉. GPU加速的二維流體實時流動仿真[J/OL]. 哈爾濱工程大學學報, 2008, 29(3):278-284[2015-12-22]. http://d.wanfangdata.com.cn/Periodical/hebgcdxxb200803015. DOI:1006-7043( 2008) 03-0278-07.

[2]STAM J. Stable fluids[J/OL]. Acm Transactions on Graphics, 1999:121-128[2015-12-20]. http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/ns.pdf.

[3]AMSDEN A A, HARLOW F H. A Simplified MAC Technique for Incompressible Fluid Flow Calculations Physics[J/OL]. Journal of Computational , 1970, 6(2): 322-325[2015-12-20]. http://dx.doi.org/10.1016/0021-9991(70)90029-X.

[4]張永學, 曹樹良, 祝寶山. 求解不可壓三維湍流的隱式SMAC方法[J/OL]. 清華大學學報(自然科學版), 2005, 45(11):1561-1564[2015-12-30]. http://dx.chinadoi.cn/10.3321/j.issn%3a1000-0054.2005.11.031.

[5]KRAUS J, SCHLOTTKE M, ADINETZ A, et al. Accelerating a C++ CFD code with OpenACC[C/OL]// 2014 First Workshop on Accelerator Programming using Directives (WACCPD), New Orleans, LA:IEEE, 2014:47-54[2015-12-30]. http://dx.doi.org/10.1109/WACCPD.2014.11.

[6]AMADOR G, GOMES A. CUDA-Based Linear Solvers for Stable Fluids[C/OL]// 2010 International Conference on Information Science and Applications (ICISA), Texas , LA:IEEE, 2010:1-8[2015-12-30]. http://dx.doi.org/10.1109/ICISA.2010.5480268.

[7]SCHEIDEGGER C E, COMBA J L D, Da C R D. Practical CFD Simulations on Programmable Graphics Hardware using SMAC[J/OL]. Computer Graphics Forum, 2005, 24(24):715-728[2015-12-30]. http://dx.doi.org/10.1111/j.1467-8659.2005.00897.x.

[8]CHORIN A J. Numerical Solution of the Navier-Stokes Equations*[J/OL]. Mathematics of Computation, 1968, 22(104):17-34 [2016-1-1]. http://dx.doi.org/10.1016/B978-0-12-174070-2.50005-8.

[9]LAX P D, HERSH R, JELTSCHR, et al. The Courant-Friedrichs-Lewy (CFL) Condition[J/OL]. Communications on Pure & Applied Mathematics, 1962, 15(4):363-371[2016-1-1]. http://link.springer.com/978-0-8176-8394-8.DOI: 10.1007/978-0-8176-8394-8.

[10]KIRK D B. Programming Massively Parallel Processors[M/OL].北京：清華大學出版社, 2010：24-46[2016-1-10]. http://www.doc88.com/p-2939593408781.html.

[11]COOK S. CUDA Programming: A Developer’s Guide to Parallel Computing with GPUs[M/OL]. Morgan Kaufmann Publishers Inc., 2012：155[2016-1-10].http://download.csdn.net/detail/niexiao2008/5780159.

[12]HIRSCH C.Numerical computation of internal and external flows [M/OL]. Butterworth-Heinemann, 2007:491-539[2016-1-10].http://dx.doi.org/10.1016/B978-075066594-0/50053-9.

[13]李大力，張理論，徐傳福,等. 雅可比迭代的CPU/GPU并行計算及在CFD中的應用[C/OL]// 2012全國高性能計算學術年文集.張家界:中國計算機學會,中國軟件行業(yè)協(xié)會，2012：1-8[2016-1-10].http://www.ccf.org.cn/sites/ccf/gxnjs.jsp#8.

[責任編輯：祝劍]

SMAC algorithm parallelization based on CUDA

CHANG Libo1，DU Huimin1，HAN Jungang2

(1. School of Electronic Engineering, Xi’an University of Posts and Telecommunications, Xi’an 710121, China;2. School of Computer Science and Technology, Xi’an University of Posts and Telecommunications, Xi’an 710121, China)

SMAC (Simplified Marker and Cell) algorithm is improved in order to enhance the ability of real time processing for fluid simulation. The algorithm is discretized for solving pressure field and update speed used by central difference scheme. A data storage algorithm that can eliminate data correlation is proposed to reduce data transmission. Usage of memory is optimized based on CUDA and the parallel computing threads into appropriate blocks is organized to make better use of stream multiprocessor in GPU. The discretization of the SMAC algorithm is improved, optimized and parallel implemented under the CUDA (unified computing device architecture) platform. Result show that GPU-based implementation of SMAC can achieve a significant speedup compared with CPU based counterpart both in the SMAC algorithm and in real simulation. Therefore the improved algorithm can achieve real-time simulation of the general scene.

computational fluid dynamics, compute unified device architecture, parallel algorithm

10.13682/j.issn.2095-6533.2016.05.007

2016-03-14

國家自然科學基金資助項目(61136002)；西安市科技發(fā)展計劃資助項目(CXY1440(10))

常立博(1985-)，男，助教，碩士，從事計算機體系結構、高性能計算研究。E-mail: changlibo@xupt.edu.cn

杜慧敏(1967-)，女，教授，博士，CCF會員。從事計算機體系結構、計算機圖形學研究。E-mail: duhuimin0529@126.com

TP338.6

A

2095-6533(2016)05-0033-06