基于GPU的原油液體仿真效率優(yōu)化研究

劉賢梅閆沖

摘要：針對有桿抽油系統(tǒng)的采油過程，原油液體在油管內(nèi)受力而不斷運(yùn)動。為了將此過程形象、逼真的展現(xiàn)出來以便于相關(guān)工作人員觀察和分析，需建立一個高精度、高效率、易使用的可視化環(huán)境。本文在VC++6.0環(huán)境下，運(yùn)用粒子系統(tǒng)技術(shù)模擬原油液體的運(yùn)動過程，并采用通過GPU加速和非GPU加速對比的方式優(yōu)化了系統(tǒng)效率。通過數(shù)據(jù)對比表明：隨著粒子數(shù)的增加，通過采用GPU加速運(yùn)算方式的系統(tǒng)效率大大提高。

關(guān)鍵詞：可視化；粒子系統(tǒng)；GPU

中圖分類號：TP393文獻(xiàn)標(biāo)識碼：A

1引言

隨著油田鉆井技術(shù)的發(fā)展，有桿采油系統(tǒng)中的可視化問題已經(jīng)成為了油田行業(yè)的一個重大問題，特別是在采油作業(yè)的過程中，需要時刻定性、定量的了解井筒中設(shè)備的運(yùn)動狀況以及原油液體狀態(tài)變化，方便相關(guān)技術(shù)人員觀察與分析，以確保采油過程順利進(jìn)行。

GPU的并行處理技術(shù)已經(jīng)在各個領(lǐng)域廣泛應(yīng)用。1999年，文獻(xiàn)1[1]在求解輻射度方程時利用顏色來標(biāo)示物體序號，利用深度緩沖來輔助計算物體可見性。2002年，文獻(xiàn)22]等利用圖形硬件進(jìn)行實時的運(yùn)動規(guī)劃，即用顏色來標(biāo)記多邊形。2004年，文獻(xiàn)3[3]等利用圖形硬件的多邊形光柵化過程和深度緩沖檢測來計算二維和三維的通用Voronoi圖。2006年，文獻(xiàn)4[4]等采用均勻網(wǎng)格來計算場景中的紋理，并建立三角形鏈表紋理，整個計算過程都放在GPU上進(jìn)行，提高了效率。自此，GPU的并行計算技術(shù)在碰撞檢測、數(shù)學(xué)運(yùn)算、像素模擬等方面得到了廣泛的應(yīng)用。

本文充分考慮圖形處理技術(shù)的并行性，以VC++6.0為開發(fā)平臺，采用GPU加速的方法，運(yùn)用粒子系統(tǒng)技術(shù)模擬了原油液體在油管內(nèi)的運(yùn)動狀態(tài)，建立一個高精度、高效率、易使用的可視化環(huán)境。本系統(tǒng)對液體仿真過程進(jìn)行優(yōu)化，是GPU技術(shù)在液體仿真和油田領(lǐng)域的一次新的突破。最終通過數(shù)據(jù)對比表明：隨著粒子數(shù)的增加，通過采用GPU加速運(yùn)算的方法使系統(tǒng)效率大大提高。

2粒子系統(tǒng)

2.1粒子系統(tǒng)簡介

迄今為止，粒子系統(tǒng)作為一種有效的圖形生成算法通常被用來模擬不規(guī)則的模糊物體，其基本思想就是將無數(shù)微小粒子作為最基本的元素來表示不規(guī)則物體[5-6]。根據(jù)粒子系統(tǒng)的應(yīng)用，大致分為四種類型：

1）隨機(jī)粒子系統(tǒng)。在隨機(jī)粒子系統(tǒng)中，粒子的屬性隨機(jī)變化，例如煙花的綻放過程、流水過程等。

2）流體粒子系統(tǒng)。流體粒子系統(tǒng)中的粒子運(yùn)動線路和運(yùn)動軌跡受流體力學(xué)的影響。例如，工廠里的煙霧擴(kuò)散情況、瀑布、以及高溫高壓下的巖漿流動情況。

3）方向粒子系統(tǒng)。方向粒子系統(tǒng)中的粒子在活動過程中相互作用、相互影響。在方向粒子系統(tǒng)中，所有粒子的屬性具有關(guān)聯(lián)性，主要用于綿軟不規(guī)則的物體，例如紡織物等。

4）結(jié)構(gòu)化粒子系統(tǒng)。結(jié)構(gòu)化粒子系統(tǒng)是由具有一定組織結(jié)構(gòu)的較小粒子組成的系統(tǒng)。其基本單元在一定的時間內(nèi)變化量較小，例如微風(fēng)下的花草樹木。

2.2粒子的基本性質(zhì)

1）粒子的物質(zhì)性。在粒子系統(tǒng)中，在任何時間和空間上，粒子與粒子、粒子與其他物體之間都不會發(fā)生重疊，具有絕對的物質(zhì)性。

計算技術(shù)與自動化2015年12月

第34卷第4期劉賢梅等：基于GPU的原油液體仿真效率優(yōu)化研究

2）粒子的基本屬性。粒子系統(tǒng)中的每一個粒子都有其基本屬性，包括運(yùn)動屬性和靜態(tài)屬性兩種。靜態(tài)屬性是指屬性物體的屬性不會隨時間的變化而變化，例如大小、顏色、生命周期、質(zhì)量等。運(yùn)動屬性是指物體的屬性會隨時間的變化而產(chǎn)生變化，這個變化可以是線性的，例如速度、位移等。

3）粒子的生命機(jī)制。每個粒子都有其生命周期，一般粒子要經(jīng)歷產(chǎn)生、活動、消亡基本過程。粒子的生命閾值到達(dá)上限是會被自動刪除，在此過程中會有新的粒子產(chǎn)生。

4）粒子的活動機(jī)制。粒子在其生命周期內(nèi)會按照預(yù)先設(shè)定好的運(yùn)動規(guī)律在場景中進(jìn)行運(yùn)動，直到生命周期結(jié)束。

5）粒子的繪制。根據(jù)被模擬物體的基本特性設(shè)定粒子的靜態(tài)屬性，再分析被模擬物體的運(yùn)動特性設(shè)定粒子的動態(tài)屬性。系統(tǒng)中的粒子按照產(chǎn)生、活動、消亡、更新的順序進(jìn)行全生命周期的活動。

2.3粒子系統(tǒng)模型

粒子系統(tǒng)是按照時間順序進(jìn)行全生命周期活動的系統(tǒng)[7]。粒子系統(tǒng)中的每一個粒子都要經(jīng)歷產(chǎn)生、活動、消亡這三個過程。在進(jìn)行每一次粒子繪制的時候，整個粒子系統(tǒng)要維持一個穩(wěn)定的循環(huán)體系，對于活動中的粒子，我們要根據(jù)被模擬對象的屬性來更新粒子的屬性，其運(yùn)行基本步驟主要包括五部分：

1）在指定區(qū)域隨機(jī)生成粒子，并賦予粒子隨機(jī)屬性；

2）根據(jù)被模擬物體的基本特性來歸納總結(jié)粒子的基本活動規(guī)律，設(shè)置好粒子的靜態(tài)屬性或者動態(tài)屬性，例如初始位置、顏色、大小、生命周期等。；

3）根據(jù)粒子生命周期，將粒子的當(dāng)前存在時間與生命上線進(jìn)行對比，判斷粒子的生命狀態(tài)，刪除生命值為零的粒子并根據(jù)需要更新粒子；

4）根據(jù)粒子的動態(tài)屬性或靜態(tài)屬性，模擬粒子的位移變換；

5）將活動粒子進(jìn)行渲染，最終描繪出被模擬物體。

3可編程GPU技術(shù)

3.1GPU工作流程

GPU處理復(fù)雜和大量圖形具有很高的效率[8]，主要是由于：

1）減少了GPU與CPU之間的數(shù)據(jù)通信

當(dāng)整個應(yīng)用針對圖形生成的時候，不再需要在CPU與GPU之間進(jìn)行多次數(shù)據(jù)交換，這樣就可以將CPU解放出來做其他的任務(wù)。這些優(yōu)勢使得GPU比CPU更適用于流處理計算，因此GPU也被認(rèn)為是一個SIMD的并行機(jī)或者流處理器，可以用于處理大規(guī)模數(shù)據(jù)集，使應(yīng)用得到加速[9]。

2）向量運(yùn)算架構(gòu)

GPU本質(zhì)上是一個向量計算模型架構(gòu)，計算單元豐富，它的這種架構(gòu)使它在處理大規(guī)模向量運(yùn)算時體現(xiàn)出較高的性能[10]，GPU運(yùn)算過程如圖1所示。

3.2GPU數(shù)據(jù)存儲與訪問

在GPU運(yùn)算中，大量數(shù)據(jù)需要存儲在紋理中，通過紋理坐標(biāo)來訪問紋理中的數(shù)據(jù)[11]。紋理相當(dāng)于一個存儲顏色值的二維數(shù)組，紋理中每個元素都存儲著一個四維向量的顏色值。在進(jìn)行GPU通用計算時，將大量的數(shù)據(jù)存儲到紋理中[12]。紋理有多種像素格式，可以根據(jù)每種數(shù)據(jù)的精度和值域來選擇合適的像素格式進(jìn)行存儲。如果存儲的數(shù)據(jù)是布爾型或者是精度很低的整型，可以用低精度的8bitsARGB2222格式來存儲；如果是高精度的浮點(diǎn)型數(shù)據(jù)，可以選用相應(yīng)的浮點(diǎn)型紋理格式來存儲[13]。在GPU通用計算中，浮點(diǎn)像素格式是比較通用的像素格式，其每個顏色通道都可以存儲一個32位的浮點(diǎn)數(shù)據(jù)[14-15]。

在紋理中通過紋理坐標(biāo)對數(shù)據(jù)進(jìn)行索引。紋理坐標(biāo)是一個二維的實數(shù)向量，一般紋理坐標(biāo)的值域在O到1之間，索引屬性值需要精確定位到紋理中的像素。例如在一個m×n尺寸的紋理中索引第x列第y行的像素則紋理坐標(biāo)的計算方法是：

U=x/m+1/（2×m）；V=y/n+1/（2×n）。

其中U和V分別為UV兩個方向的紋理坐標(biāo)值。

4系統(tǒng)實現(xiàn)與分析

4.1GPU下油管內(nèi)液體粒子的數(shù)據(jù)組織

實驗描述：模擬在上沖程過程中，原油液體在泵筒中由下往上的運(yùn)動狀態(tài)，通過機(jī)組實驗數(shù)據(jù)對比在粒子數(shù)變化幅度相同的情況下，CPU和GPU的效率變化。

在上沖程過程中，抽油桿和柱塞上行，原油液體進(jìn)入泵筒內(nèi)部，此時液柱受到重力G（豎直向下）、推力ft（豎直向上）、摩擦阻力fm（沿著泵筒向下）。在泵筒內(nèi)液柱向上運(yùn)動，根據(jù)原油粒子特點(diǎn)，需要將以下粒子屬性儲存到紋理中：

Speed：用三維向量表示液體粒子的流動速度；

Live：粒子的生命閾值，存儲為布爾型的標(biāo)量，True表示粒子處于活動中，F(xiàn)alse表示粒子消亡；

Position：用三維向量存儲粒子在世界坐標(biāo)系下的位置。

由于液體粒子數(shù)量很多、大小不一，因此存放粒子時應(yīng)采取紋理多樣化的方式，根據(jù)不同精度選擇合適的像素格式進(jìn)行有效存儲。例如：當(dāng)存儲低精度數(shù)據(jù)時，采用低精度的8 bits ARGB 2222格式，反之當(dāng)存儲高精度數(shù)據(jù)時，采用高精度的32 bits ARGB 8888格式。在GPU的通用計算中，為了滿足計算的精度要求，通常浮點(diǎn)像素格式的每一個顏色通道都將存儲一個8位或者32位的浮點(diǎn)數(shù)據(jù)。在存儲數(shù)據(jù)的過程中既要考慮效率問題又要考慮存儲空間的合理利用問題，通常的像素格式都是具有RGBA四個顏色通道，粒子的速度和位置都是三維向量，這樣把粒子的速度和位置存放到紋理像素中的時候就會剩余一個通道，所以將粒子的Live屬性儲存到最后一個通道里，這樣就充分合理的利用了存儲空間。

格式選取完畢后開始創(chuàng)建存儲粒子數(shù)據(jù)的紋理，確保每一個粒子對應(yīng)一個像素點(diǎn)，這樣像素點(diǎn)的最大值等于粒子的最大數(shù)目。每一個粒子的屬性值存儲在通過紋理中并通過各自的紋理坐標(biāo)進(jìn)行檢索。

4.2液體粒子屬性計算

在液體粒子產(chǎn)生、活動、消亡的過程中，把握并控制好粒子的屬性，準(zhǔn)確實時計算出粒子的速度、生命值、位置是繪制液體粒子的核心。首先編寫像素著色器程序，通過GPU的通用計算將粒子屬性的計算問題具體化，計算出每一個液體粒子的屬性值并存儲到紋理中以備使用。

利用GPU通用計算來計算每一個粒子屬性的時候，粒子要經(jīng)歷產(chǎn)生、活動。消亡三個過程。由于在GPU下不能動態(tài)為粒子分配存儲空間，所以每一個活動中的粒子的生命值由最大值降為零時，將此粒子的生命值、位置、速度進(jìn)行初始化，作為新粒子重新加入到系統(tǒng)中。粒子系統(tǒng)開始運(yùn)行后，粒子的數(shù)量隨著時間的推移不斷增加，粒子的屬性以二維索引的形式存儲在紋理中，在進(jìn)行GPU通用計算時需要將二維索引轉(zhuǎn)化為一維索引。對于二維索引來說，需要依靠紋理坐標(biāo)來將其轉(zhuǎn)換為一維索引，轉(zhuǎn)換公式為：M=（U+V×H）×W，其中M是二維索引轉(zhuǎn)換為一維索引的序列號，U和V為紋理坐標(biāo)的兩個分量，W和H為兩個紋理坐標(biāo)分量的大小。液體粒子屬性值計算步驟如下：

1）液體在受力中的體積和數(shù)量的計算

原油液體在油管內(nèi)流動的過程中，粒子系統(tǒng)在每一幀都要保持產(chǎn)生新的粒子以維持液體的流動，這樣就要對新生成的粒子的速度、位置、生命值進(jìn)行初始化。在上沖程過程中，油管內(nèi)的液柱隨著柱塞的上升而升高，所有液體粒子由下部初始位置升到最高點(diǎn)，初始點(diǎn)的所有粒子位置要添加一個隨機(jī)值：P=P0+Prant，其中P是粒子的初始位置，P0是油管內(nèi)液柱上升的初始位置，Prand是一個三維隨機(jī)向量，此向量決定液柱的粗細(xì)。如圖2為上沖程過程中原油液體進(jìn)入泵筒內(nèi)的示意圖。此時泵筒內(nèi)原油液體受重力G（豎直向下）、推力ft（豎直向上）、摩擦阻力fm（沿著泵筒向下）、泵筒內(nèi)壓力fp。如圖2中x表示柱塞在時刻t的位移，表示柱塞位移為x時的液面高度，Lw是尾管長度。設(shè)任意時刻活塞t運(yùn)動到x的位置，進(jìn)入泵筒內(nèi)原油液體體積為Vt。根據(jù)受力分析，將整個泵筒內(nèi)空間分為兩部分，Va為原油液體未充滿泵筒前真空部分的體積，Vb為下半部分體積，其中Va=（x-y）A，Vt=Va+Vb，Vb=（Rs（x）-Rs（0））（1-Lw）（Ls+y）ABo，其中Rs是溶解油氣比、Bo為原油體積系數(shù)與水的體積系數(shù)、A為壓力系數(shù)。假設(shè)單位粒子體積為1，則有：1×Vnum=Vt，其中Vnum為粒子隨體積變化的數(shù)量。

2）活動粒子屬性值的計算

上沖程過程中，液體粒子受力向上運(yùn)動，其位置和速度在不斷地變化著。下面分析粒子的受力狀況，根據(jù)里的平衡原理得出：fm+fp+G=ft，vm=vm-1+a×Δt，其中Vm是粒子的速度，Vm-1是上一幀粒子的速度，a是粒子的加速度，Δt是一幀持續(xù)的時間。在Y方向上，粒子的加速度由受力變化確定，在X和Z方向上粒子的速度分量不變。Pm=Pm-1+（a×Δt）×Δt，其中Pm是粒子的位置，Pm-1是粒子在上一幀的位置，Δt是一幀持續(xù)的時間。

3）消亡粒子屬性值的計算

當(dāng)液體粒子在油管內(nèi)達(dá)到最高點(diǎn)時，粒子的生命值降為零，此時的粒子屬性不在發(fā)生變化。在此粒子重新生成并再次加入粒子系統(tǒng)之前，既不需要計算其屬性值，也不用被渲染。

4.3液體粒子繪制

計算完液體粒子的位置和速度，下一步需要讀取粒子的屬性值并進(jìn)行繪制。運(yùn)用紋理獲取技術(shù)（Vertex Texture Fetch）來獲取粒子的屬性值，用tex2Dlod函數(shù)在頂點(diǎn)著色器讀取紋理像素中的粒子屬性值進(jìn)行繪制，這樣就不需要通過CPU傳輸，大大降低了CPU與GPU之間的通信，提高了渲染速度。下面進(jìn)行粒子的繪制步驟：

1）建立頂點(diǎn)緩沖區(qū)

選擇圓形作為繪制液體粒子的基本元素，組成每個粒子的最基本單位設(shè)定為四邊形。首先建立頂點(diǎn)緩沖區(qū)，頂點(diǎn)的個數(shù)為四倍的粒子數(shù)，然后建立頂點(diǎn)索引，將每四個頂點(diǎn)組合在一起組成四邊形，粒子以四邊形的方式被繪制出來。創(chuàng)建出來的頂點(diǎn)緩沖區(qū)需要預(yù)先設(shè)計其格式。首先給液體粒子貼上紋理（此處選擇普通水滴紋理），利用Texcoord0（紋理通道0的紋理坐標(biāo)）給粒子的頂點(diǎn)的紋理坐標(biāo)分別賦值為（0，1）（1，0）（0，1）（1，1），這樣就可以通過紋理坐標(biāo)0將水滴的紋理映射在粒子上，并判斷出水滴紋理在粒子上的位置。其次要對粒子的屬性值進(jìn)行索引，將粒子的紋理屬性存儲在Texcoord1（紋理通道1的紋理坐標(biāo)）中，這樣就可以索引到粒子的屬性值了。

2）頂點(diǎn)繪制

首先讀取粒子的紋理屬性并獲取粒子在世界坐標(biāo)系下的位置，由于粒子必須要正對著屏幕，所以要將粒子變換到相機(jī)坐標(biāo)下，利用紋理坐標(biāo)將粒子的四個頂點(diǎn)分別繪制在粒子的四個角上。由于GPU會將所有頂點(diǎn)繪制出來，即無論是活動中的頂點(diǎn)還是已經(jīng)消亡的頂點(diǎn)，GPU都會無選擇性的將其全部繪制出來，因此把消亡的粒子的頂點(diǎn)位置全部設(shè)定到一起，這樣粒子就不會顯示出來。頂點(diǎn)位置的計算公式如下：

Pt=（P×Mv+（U0，V0，0）×Size×Live）×Mp，其中Pt是最終計算的頂點(diǎn)位置，P是屬性紋理中的粒子位置，U0和V0分別是頂點(diǎn)紋理坐標(biāo)0的兩個分量，Size是粒子的大小，Live是粒子的生命值，Mv和Mp分別是相機(jī)矩陣和投影矩陣。當(dāng)粒子消亡時，粒子的Live值為零，粒子不需要被繪制出來；反之當(dāng)粒子活動時，粒子的Live值為1，粒子頂點(diǎn)的位置會被正確的計算出來。

5效果與結(jié)論

本實驗硬件環(huán)境為：Intel2.8GHz處理器，4G內(nèi)存，NVIDIAGeForceGT640M顯卡，操作系統(tǒng)是64位Window7系統(tǒng)，開發(fā)工具為VC++6.0和Directx9.0。圖4.1為原油液體的流動效果圖。表1為獨(dú)立環(huán)境下CPU模擬的原油液體和GPU模擬的原油液體幀速率的比較。如圖4.2為表一數(shù)據(jù)的對比圖。

結(jié)論：針對有桿抽油系統(tǒng)采油過程中原油液體在油管內(nèi)的受力運(yùn)動，本文設(shè)計并建立了一個能滿足高精確度、高效率、易使用的可視化環(huán)境。從表中我們可以看出，當(dāng)原油液體粒子數(shù)較少時，CPU實現(xiàn)的速度比GPU還要快，但隨著粒子數(shù)的增大，CPU實現(xiàn)的幀速率快速下降，而GPU實現(xiàn)的幀速率下降速度緩慢，當(dāng)粒子數(shù)達(dá)到300萬時，GPU仍然能夠?qū)崟r的繪制，而CPU實現(xiàn)每秒只有6.8幀，這說明GPU對大量粒子的處理能力遠(yuǎn)遠(yuǎn)高于CPU。

參考文獻(xiàn)

[1]HopfMatthias，ErtlThomas.Accelerating 3D convolutionusinggraphicshardware.[J] Proceedings of IEEE Visualization.1999，28（05）：56-70，62.

[2]FRANCISCO J Seron，Rafael Rodrigues，Eva Cerezo，AlfredoPina.AddingSupport for Highlevel Skeletal Animation[J]. IEEE Transactions onvisualization and computer graphics.2002，28（05）：120-134，122.

[3]Peter Kipfer，Mark Segal，Rüdiger Westermann.UberFlow：a GPUbased particle engine[J]. Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on Graphics HardwareHWWS04.2004， 28（05）：128-132，130.

[4]ARON E Lefohn，Shubhabrata Sengupta，Joe Kniss，et al.Glift：Generic，efficient，randomaccess GPU data structures[J]. ACM Transactions on Graphics.2006， 28（05）：120-134，123.

[5]田蘇昕.虛擬現(xiàn)實中大規(guī)模粒子系統(tǒng)的研究[J].沈陽航空工業(yè)學(xué)院學(xué)報，2009，30（3）：149-150，160.

[6]徐陽東.基于粒子系統(tǒng)的不規(guī)則景物建模研究[J].山東師范大學(xué)學(xué)報，2009， 28（5）：128-132，138.

[7]白玉.基于粒子系統(tǒng)的煤礦火災(zāi)模擬系統(tǒng)[J].吉林大學(xué)學(xué)報，2010，8（2）：47-49.

[8]王濤. 基于GPU的程序分析與并行化研究[J].解放軍信息工程大學(xué)學(xué)報，2010，8（2）：47-51.

[9]胡杰. CPU-GPU異構(gòu)平臺計算模型的研究與應(yīng)用[J].大連理工大學(xué)學(xué)報，2011，29（9）：44-48.

[10]郭忠明. 基于CUDA的并行圖像處理性能優(yōu)化[J].大連理工大學(xué)學(xué)報，2012，99（4）：49-58.

[11]方旭東. 面向大規(guī)?？茖W(xué)計算的CPUGPU異構(gòu)并行技術(shù)研究[J].國防科學(xué)技術(shù)大學(xué)學(xué)報，2009，22（9）：165-169.

[12]BUCK I，F(xiàn)OLEY T，HORN D，SUGERMAN J，HANRAHAN P.Brook for GPUs： Stream Computing on Graphics Hardware[J].Proceedings of SIGGRAPH，2004，28（4）：55-59.

[13]Jizhu Lu，et.al.HMMerCell： High Performance Protein Profile Searching on the Cell/B.E.Processor[J].Proceeding of IEEE International Symposium on Performance Analysis of Systems and SoftwareISPASS，2008，33（7）：94-98.

[14]S.García，A.Fernández，J.Luengo et al.A study of statistical techniques and performance measures for geneticsbased machine learning： accuracy and interpretability[J].Soft Computing，2009，83（7）：104-108.

[15]Zhe Wei， Yixiong Feng， Jianrong Tan et al.. Multi-objective performance optimal design of largescale injection molding machine[J]. The International Journal of Advanced Manufacturing Technology 2009， 30（3）：166-170，160.