一種基于gpu集群的并行流體仿真加速方法
【專利摘要】本發(fā)明是一種基于GPU集群的并行流體仿真加速方法。針對流體仿真計算量大、并行度高的特點(diǎn)，設(shè)計了自動化的并行加速方法，進(jìn)行流體的模擬與算法研究。提出并實(shí)現(xiàn)了一種在同一節(jié)點(diǎn)內(nèi)多個GPU之間、多個節(jié)點(diǎn)之間的負(fù)載均衡算法。在給定空間中使用基于位置的流體仿真算法模擬流體的物理行為，通過減少分支以及縮小臨域搜索范圍對算法進(jìn)行加速。對流體所在的空間進(jìn)行分割，每個節(jié)點(diǎn)處理一個子空間，在節(jié)點(diǎn)內(nèi)根據(jù)GPU數(shù)量進(jìn)行進(jìn)一步分割，從而完成流體仿真在GPU集群上的并行化。
【專利說明】一種基于GPU集群的并行流體仿真加速方法

【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于計算機(jī)應(yīng)用領(lǐng)域，具體地說是一種基于GPU集群的并行流體仿真加速 方法，該方法可用于影視、游戲、醫(yī)學(xué)等各個領(lǐng)域中，模擬流體的行為。

【背景技術(shù)】
[0002] 隨著虛擬現(xiàn)實(shí)技術(shù)的發(fā)展，流體動畫在電影特技、娛樂游戲、軍事仿真、醫(yī)學(xué)仿真、 數(shù)字媒體中廣泛應(yīng)用，人們對真實(shí)感有了越來越高的需求。為了仿真這些復(fù)雜的場景效果， 一幀一幀手工渲染動畫自然滿足不了需求，近年來，計算機(jī)硬件性能飛速發(fā)展，極大的促進(jìn) 了計算機(jī)圖形學(xué)的深入研究。流體仿真更是多種多樣，從早期簡單的大面積水體仿真，如平 靜的水面、靜止的煙霧，到現(xiàn)在渲染更為豐富細(xì)致的動態(tài)流體，如飛濺、泡沫、沸騰、多流體 融合、多相流，流體交互仿真等，流體和其物理屬性已經(jīng)緊密的結(jié)合在一起。因此，基于物理 的流體動畫仿真成為計算機(jī)圖形學(xué)和虛擬現(xiàn)實(shí)領(lǐng)域中的研究熱點(diǎn)之一。
[0003] 近年來被廣泛應(yīng)用的光滑粒子流體動力學(xué)（Smoothed Particle Hydrodynamics) 算法是一種基于物理的流體仿真算法。該類方法使用粒子離散化問題域，通過加權(quán)求和鄰 居粒子的貢獻(xiàn)值來計算粒子的物理量。該算法能夠很好的處理自由表面、變形邊界、運(yùn)動交 界面以及大變形等問題。但是SPH算法使用了剛性方程以維持其不可壓縮性，這使得力的 計算結(jié)果較大，只能使用足夠小的時間步長。
[0004] 基于位置的流體仿真（Position Based Fluid)算法在基于位置的動力學(xué) (Position Based Dynamics7F1DB)框架中添加了密度的迭代求解。通過一組位置限制公式 來維持一個固定的密度。該方法能夠得到同SPH算法類似的不可壓縮性以及收斂性，同時 繼承了 PBD的幾何穩(wěn)定性，從而可以在實(shí)時仿真中使用更大的時間步。同時改善了 SPH算 法對時間步長的限制，允許使用更大的時間步長，在相同時間內(nèi)可以獲得更大的粒子位置 更新，同時保留了基于粒子的流體仿真所固有的優(yōu)勢，具有廣闊的應(yīng)用前景。
[0005] 為了使得仿真效果更加逼真、細(xì)節(jié)更豐富，粒子規(guī)模被不斷提高。PBF算法的計算 量非常大，CPU已經(jīng)難以滿足PBF方法的計算需求，浮點(diǎn)性能超過CPU幾十倍的GPU開始被 用于加速PBF方法。最近十余年，GPU在浮點(diǎn)性能和可編程性方面都取得了長足的進(jìn)步。浮 點(diǎn)性能方面，目前GPU無論在單精度浮點(diǎn)運(yùn)算還是雙精度浮點(diǎn)運(yùn)算都達(dá)到了數(shù)Tflops?？?編程性方面，隨著CUDA的問世，GPU從傳統(tǒng)的圖形領(lǐng)域跨入了通用計算領(lǐng)域。越來越多的 應(yīng)用通過GPU加速獲得了大幅度的性能提升。GPU是一種并行硬件架構(gòu)，特別適合于加速計 算密度大的并行算法。


【發(fā)明內(nèi)容】

[0006] 本發(fā)明的目的是在GPU集群上實(shí)現(xiàn)PBF的自動化加速算法，使得可以用同一種負(fù) 載均衡策略、通信管理策略、存儲管理策略來實(shí)現(xiàn)基于GPU集群的流體仿真算法。讓開發(fā)人 員能夠快速在GPU集群上實(shí)現(xiàn)仿真算法，將研發(fā)精力從算法實(shí)現(xiàn)轉(zhuǎn)移至對仿真算法本身的 研究，只投入少量精力即可獲得GPU集群的加速。
[0007] 為了實(shí)現(xiàn)上述目的，對流體物理行為的模擬采用基于位置的流體仿真算法（PBF) 實(shí)現(xiàn)，并通過減少分支，縮小臨域搜索空間對算法進(jìn)行優(yōu)化。按照節(jié)點(diǎn)數(shù)量、節(jié)點(diǎn)內(nèi)GPU數(shù) 量和粒子分布情況對全局空間進(jìn)行劃分，將計算任務(wù)均衡的劃分至各個GPU中，在一些聯(lián) 系性較強(qiáng)的步驟執(zhí)行后對粒子的一些關(guān)鍵屬性進(jìn)行全局同步，使其能夠正確的并行的快速 的在GPU集群上完成整個計算任務(wù)。在計算過程中，根據(jù)各個GPU的計算時間，對全局的 GPU進(jìn)行負(fù)載均衡。

【專利附圖】

【附圖說明】
[0008] 圖1示出本發(fā)明中并行加速方法所使用平臺的整體邏輯結(jié)構(gòu)；
[0009] 圖2示出本發(fā)明基于粒子分布情況與服務(wù)器節(jié)點(diǎn)數(shù)量劃分空間；
[0010] 圖3示出本發(fā)明每個子空間周圍的ghost粒子層；
[0011] 圖4示出本發(fā)明縮小臨域搜索范圍。

【具體實(shí)施方式】
[0012] 下面結(jié)合附圖對本發(fā)明作進(jìn)一步說明。
[0013] 并行加速方法所使用平臺的整體邏輯結(jié)構(gòu)如圖1所示。具體步驟如下：
[0014] (1)按照服務(wù)器節(jié)點(diǎn)數(shù)量和粒子的空間分布情況對空間進(jìn)行二維ORB劃分，如圖2 所示。統(tǒng)計粒子在三個方向上坐標(biāo)分布的方差，取方差最大的兩個方向作為切割方向。首 先沿方差最大的方向?qū)鼍扒懈顬閮蓚€長方體，然后在兩個長方體中繼續(xù)選擇其余一個方 向切割為更小的兩個長方體，循環(huán)在兩個方向上切割直到子空間的數(shù)量和GPU數(shù)量相同；
[0015] (2)統(tǒng)計每個節(jié)點(diǎn)的所屬子空間中粒子在ORB劃分的兩個方向上的方差，選定方 差大的方向?qū)υ撟涌臻g進(jìn)行一維劃分，將整個仿真空間沿該方向切割為多個長方體子空 間。將粒子數(shù)據(jù)發(fā)送至所屬節(jié)點(diǎn)的相應(yīng)GPU中，并包含該粒子群周圍的一些冗余粒子信息 用于后續(xù)計算，稱之為ghost粒子層，如圖3所示；
[0016] (3)根據(jù)粒子上一個時間步受力情況（模擬開始時初始化為重力）和速度，計算每 個粒子的速度變化；
[0017] (4)將GPU的粒子空間劃分為多個網(wǎng)格，如圖4所示，網(wǎng)格尺寸為光滑半徑的一半， 以縮小搜索范圍，統(tǒng)計粒子在空間網(wǎng)格中的分布情況，統(tǒng)計每一個網(wǎng)格的起始終止粒子，用 于后續(xù)步驟的臨域搜索，這里按（2)步驟中的劃分方向合并多個網(wǎng)格，以減少分支；
[0018] (5)計算每個粒子的密度P。粒子密度由標(biāo)準(zhǔn)SPH密度估算公式得出，如粒子i 的密度P i是關(guān)于其所有臨域粒子位置坐標(biāo)的函數(shù)：

【權(quán)利要求】
1. 一種基于GPU集群的并行流體仿真加速方法，其特征在于包括以下步驟： (1) 在全局空間中按照節(jié)點(diǎn)數(shù)量和粒子分布情況對其在某一二維平面進(jìn)行劃分，而非 三維空間上進(jìn)行ORB劃分，每個節(jié)點(diǎn)對應(yīng)一塊子空間； (2) 在每個節(jié)點(diǎn)的所屬子空間中根據(jù)GPU數(shù)量和該子空間中粒子分布差異最大化的坐 標(biāo)軸方向?qū)υ撟涌臻g進(jìn)行一維劃分，并沿該方向在對其子空間內(nèi)粒子進(jìn)行矩陣變換，將粒 子數(shù)據(jù)發(fā)送至所屬節(jié)點(diǎn)的相應(yīng)GPU中； (3) 各GPU使用基于位置的流體仿真算法（PBF)進(jìn)行計算，多次迭代修正粒子位置，再 計算限制因子項(xiàng)，更新粒子位置之后分別對限制因子項(xiàng)和粒子位置進(jìn)行兩次全局同步，同 步后本時間步結(jié)束； (4) 根據(jù)各GPU的計算時間，對各個節(jié)點(diǎn)各個GPU進(jìn)行負(fù)載均衡處理，調(diào)整全局空間的 劃分，根據(jù)新的劃分在GPU之間交換粒子，而后重復(fù)步驟（3)，處理下一個時間步。
2. 如權(quán)利要求1所述的基于GPU集群的并行流體仿真加速方法，其特征在于：在所述 步驟（1)中在全局空間中按照節(jié)點(diǎn)數(shù)量和粒子分布情況對其在某一二維平面進(jìn)行劃分，具 體方法是： 按照服務(wù)器節(jié)點(diǎn)數(shù)量和粒子的空間分布情況對空間進(jìn)行二維ORB劃分，ORB方式采用 多層次的方式切割，選定兩個方向進(jìn)行劃分；統(tǒng)計粒子在三個方向上坐標(biāo)分布的方差，取方 差最大的兩個方向作為切割方向；首先沿方差最大的方向?qū)鼍扒懈顬閮蓚€長方體，然后 在兩個長方體中繼續(xù)選擇其余一個方向切割為更小的兩個長方體，循環(huán)在兩個方向上切割 直到子空間的數(shù)量和GPU數(shù)量相同；此切割方法目的在于減小節(jié)點(diǎn)間交界面面積，同時不 至于使傳輸模型過于復(fù)雜，最小化各節(jié)點(diǎn)傳輸數(shù)據(jù)量。
3. 如權(quán)利要求1所述的基于GPU集群的并行流體仿真加速方法，其特征在于：在所屬 步驟（2)中對每個節(jié)點(diǎn)的全局子空間進(jìn)行劃分，具體方法是： 計算粒子坐標(biāo)集在步驟（2)選取兩個方向上的方差，沿方差最大的坐標(biāo)軸粒子的分布 差異相對最大，按照粒子在該方向的數(shù)量分布，將整個仿真空間沿該方向切割為多個長方 體子空間；此切割方法目的在于在使節(jié)點(diǎn)間傳輸事務(wù)次數(shù)最小化的前提下，節(jié)點(diǎn)內(nèi)GPU間 的傳輸數(shù)據(jù)量和傳輸事務(wù)次數(shù)最小化。
4. 如權(quán)利要求1所述的基于GPU集群的并行流體仿真加速方法，其特征在于：在所屬 步驟（3)中各GPU使用基于位置的流體仿真算法（PBF)進(jìn)行計算并在計算限制因子項(xiàng)和校 正粒子位置兩個步驟中進(jìn)行全局同步兩個關(guān)鍵步驟，包括步驟如下： (3. 1)根據(jù)粒子上一個時間步受力情況和速度，計算每個粒子的速度變化；所述受力 情況模擬開始時初始化為重力； (3. 2)根據(jù)粒子速度，預(yù)測每個粒子的下一步位置并更新； (3. 3)將GPU的粒子空間劃分為多個網(wǎng)格，網(wǎng)格尺寸為光滑半徑的一半，以縮小搜索范 圍，統(tǒng)計粒子在空間網(wǎng)格中的分布情況，統(tǒng)計每一個網(wǎng)格的起始終止粒子，用于后續(xù)步驟的 臨域搜索，這里按（2)步驟中的劃分方向合并多個網(wǎng)格，以減少分支； (3.4)計算每個粒子的限制因子項(xiàng)，并全局同步粒子的限制因子項(xiàng)； (3. 5)根據(jù)粒子的臨域粒子限制因子項(xiàng)，計算每個粒子的位置偏移，更新粒子位置并全 局同步； (3. 6)對每個粒子進(jìn)行邊界檢測； (3. 7)對步驟（3. 4)、（3. 5)、（3. 6)進(jìn)行多次重復(fù)迭代，以校正（3. 2)中的粒子的預(yù)測 位置，當(dāng)?shù)竭_(dá)最大迭代次數(shù)或位置偏移小于某一閾值后，進(jìn)行下一步； (3. 8)根據(jù)校正后的粒子位置，與上一時間步中粒子的位置，計算粒子的速度，計算粒 子的受力情況，包括渦流和粘滯力； (3. 9)在分布粒子或重新分布粒子N個時間步后，若全局?jǐn)?shù)據(jù)傳輸量大于某一閾值M， 則表明粒子目前的分布狀態(tài)不適用于當(dāng)前的空間切割方式，此時重新按照節(jié)點(diǎn)和GPU的數(shù) 量劃分整個空間，并計算重新分布后第一次全局同步的數(shù)據(jù)傳輸總量m，設(shè)置閾值M為a*m， 其中a是一個大于1的常量，在實(shí)現(xiàn)中根據(jù)問題的具體情況調(diào)整大??； (3. 10)重復(fù)步驟（3. 1)至（3. 9)，直至整個模擬過程被結(jié)束。
5.如權(quán)利要求1所述的基于GPU集群的并行流體仿真加速方法，其特征在于：在所屬 步驟（4)中各節(jié)點(diǎn)各GPU需要進(jìn)行負(fù)載均衡，具體方法是： 將整個集群組織成二叉樹結(jié)構(gòu)，葉節(jié)點(diǎn)為服務(wù)器節(jié)點(diǎn)，非葉節(jié)點(diǎn)代表左右子樹上的節(jié) 點(diǎn)間的交界面；對于每一個非葉節(jié)點(diǎn)，我們根據(jù)左右子樹節(jié)點(diǎn)在上一個時間步的平均耗時 來調(diào)整它們的交界面。當(dāng)左子樹的節(jié)點(diǎn)的計算時間高于右側(cè)節(jié)點(diǎn)時，交界面向左側(cè)移動一 段距離，一個網(wǎng)格層的粒子被交換到右側(cè)節(jié)點(diǎn)上，反之則將右側(cè)節(jié)點(diǎn)上的一個網(wǎng)格層中的 粒子交換給左側(cè)節(jié)點(diǎn)；逐層向上迭代重復(fù)此過程，直至根節(jié)點(diǎn)； 節(jié)點(diǎn)內(nèi)GPU的負(fù)載均衡方式同節(jié)點(diǎn)間相似，但因?yàn)槭且痪S劃分，故在不需要沿葉節(jié)點(diǎn) 向上重復(fù)，只計算葉節(jié)點(diǎn)一層即可。
【文檔編號】G06F9/46GK104360896SQ201410737607
【公開日】2015年2月18日 申請日期:2014年12月4日 優(yōu)先權(quán)日:2014年12月4日
【發(fā)明者】沈旭昆, 黃冠喆 申請人:北京航空航天大學(xué)