基于iOS的實(shí)時流體模擬的并行處理

張鵬

（鄭州科技學(xué)院信息工程學(xué)院 河南省鄭州市 450064）

1 緒論

隨著移動智能手機(jī)的最新發(fā)展，越來越多的人使用智能手機(jī)生活。其中的一些代表智能手機(jī)，比如：攜帶Apple iOS系統(tǒng)的iPhone以及iPad，三星攜帶Android系統(tǒng)的galaxy，以及華為的手機(jī)等。

這些智能手機(jī)的處理器都具有體積非常緊湊，且功耗低的特點(diǎn)?，F(xiàn)在大部分智能手機(jī)的系統(tǒng)都是由ARM處理器構(gòu)成的，由于設(shè)計(jì)以及存儲等其它缺點(diǎn)的存在，由ARM構(gòu)成的設(shè)備性能要遠(yuǎn)低于PC電腦的配置性能。隨著智能手機(jī)性能的多樣化，人們使用它觀看電影，玩游戲或者做更多的事情，在這種情況下，我們期待智能手機(jī)能夠提高現(xiàn)在的處理器性能從而使我們做更多的工作。

流體模擬是21世紀(jì)流體力學(xué)的重要組成部分，在電影，游戲以及其它領(lǐng)域有廣泛的應(yīng)用。但是，伴隨著流體計(jì)算量的增加，流體模擬的速度會變得非常慢。我們的目標(biāo)就是研究在iOS系統(tǒng)里實(shí)現(xiàn)流體模擬，并且使用ARM的命令集進(jìn)行并行加速。

2 流體模擬的含義

2.1 先前研究

在物理學(xué)中，納維斯托克斯方程(Navier-Stokes equations)用來描述流體的運(yùn)動。流體模擬，一般是利用計(jì)算機(jī)求解流體方程式后，用數(shù)字模擬的形式在計(jì)算機(jī)里模擬流體的運(yùn)動?，F(xiàn)在流體模擬的技術(shù)分為幾種，最常用的是歐拉方法。這種方法將流體領(lǐng)域劃分為特定的空間，將空間分為網(wǎng)格組成，利用有限差分法求解流體方程式的速度，壓力等。

第二種是拉格朗日方法，將流體看做是一個個粒子，每個粒子代表流動的流體，具有流體的特征。Vortices方法和Lattice Boltzmann方法是最近常用的求解方法。這篇論文使用的流體模擬是根據(jù)歐拉方法進(jìn)行研究。

2.2 數(shù)學(xué)方程式

圖1

根據(jù)不可壓縮性流體的質(zhì)量守恒定律和動量守恒定律，我們得到不可壓縮牛頓流體的連續(xù)性的納維斯托克斯方程。1999年，Jos Stam發(fā)表了“Stable Fluids”[2]研究，為了解出對流式，使用了半拉格朗日方法(semi-lagragian)。此方法解決了大時間間隔數(shù)值的不穩(wěn)定問題，隨后，使用霍奇分解將納維斯托克斯方程變?yōu)椋?/p>

這里P為投射算子。

2003年，Jos Stam發(fā)表了他的研究“Real-Time Fluid Dynamics for Games”[3]，使用密度替換速度，將上面方程式（1）變?yōu)椋?/p>

3 在IOS系統(tǒng)的實(shí)時流體模擬的實(shí)現(xiàn)

OpenGL ES是OpenGL開發(fā)的基于嵌入式版本的3維計(jì)算機(jī)圖形API。我們使用OpenGL ES在IOS系統(tǒng)進(jìn)行繪制，但是在IOS使用OpenGL ES進(jìn)行繪制之前，我們需要先生成繪制環(huán)境。

根據(jù)Apple提供的在IOS使用OpenGL ES編程向?qū)4]，我們可以根據(jù)幾個步驟進(jìn)行IOS繪圖環(huán)境設(shè)置：

設(shè)置EAGLContext→建立OpenGL繪圖內(nèi)容→建立渲染緩沖→建立幀緩沖→展現(xiàn)

3.1 設(shè)計(jì)流體對象

以前面的方程式為基礎(chǔ)，我們設(shè)計(jì)流體對象包括部分如圖1所示。

3.2 流體對象的實(shí)現(xiàn)

為了解開原有的流體方程式，需要將連續(xù)方程變?yōu)殡x散方程。在歐拉方法的基礎(chǔ)上，我們利用有限差分法將原有流體方程式的各個項(xiàng)拆分為以下形式：

在這里，i為水平坐標(biāo)，j表示為垂直坐標(biāo)。為了判定邊界條件，我們追加了新的邊界，將網(wǎng)格大小設(shè)定為N+2。

3.2.1 外部力

因?yàn)閕Phone采用的是觸屏模式，所以將我們的觸屏做為外部的力源。通過我們手指對屏幕的控制做為我們的壓力源。假設(shè)之前時刻為T，那么現(xiàn)在時刻為T+ΔT,如果追加外部力的話，可以用下面方程式表示：

3.2.2 擴(kuò)散

每個網(wǎng)格只能和周圍相鄰的網(wǎng)格進(jìn)行交互，而且通過上面的方程式（4）和（6）得到在T時刻的速度和密度，為了得到T+ΔT時刻的速度使用以下方程式求解。在方程中，h表示網(wǎng)格的長度，k表示擴(kuò)散的強(qiáng)度由于擴(kuò)散比率大的話，會造成密度數(shù)值的不穩(wěn)定，從而最終造成發(fā)散現(xiàn)象。為了避免這種現(xiàn)象，我們考慮采用Jos Stam論文[3]的方法：

圖2：在索引(i,j)的數(shù)據(jù)的計(jì)算

圖3：在索引(i+1,j)數(shù)據(jù)的計(jì)算

下一步，我們采用高斯賽德爾迭代算法求解。

3.2.3 對流

關(guān)于對流的求解，1999年Jos Stam使用半拉格朗日算法(Semi-Lagrangian)求解?；舅枷胧峭蟮雇薚時刻，來計(jì)算現(xiàn)在時刻：T+ΔT 時刻的粒子位置。此方法解決了流體模擬時間步長過大時數(shù)值的不穩(wěn)定問題。

3.2.4 投射

計(jì)算完對流后，我們可以結(jié)束密度的計(jì)算。但是對于速度場，我們必須要投射速度。我們假定U3是已求解出擴(kuò)散后的速度，根據(jù)Stable Fluids[2]提供的方法，我們將投射算子P投射于速度U3，利用霍奇分解并且使用?算子，可以得到：

圖4：寄存器 Q4=Add(temp1,temp2);存儲數(shù)據(jù)(i,j)

最終可以求解得出U4

3.2.5 邊界條件

求解邊界條件問題的方法有許多種。這里，我們簡單的假定流體是放在一個類似箱子里。當(dāng)密度和速度遇到邊界的時候，取相反的方向。

3.3 在iPhone的繪制

我們求解出速度和密度的數(shù)值后，可以使用OpenGL ES的函數(shù)繪制出來。在渲染管道圖元中，我們選擇使用GL_LINE繪制速度，使用GL_TRIANGLE_FAN來繪制密度，顏色上使用混合色。

3.4 問題

在IOS測試流體模擬的時候，當(dāng)流體的網(wǎng)格增大的時候，流體模擬的速度場和密度場會變的非常慢，所以我們決定進(jìn)行并行處理流體模擬，從而加快流體模擬的速度和密度運(yùn)動。iPhone的處理器是ARM架構(gòu)，而ARM提供一套并行計(jì)算的命令集合:neon指令集,我們決定使用neon指令集來進(jìn)行流體模擬的并行處理。

4 利用Neon指令集的并行處理

4.1 以前算法的問題點(diǎn)

為了解開流體方程式，傳統(tǒng)的方法是利用有限差分方法進(jìn)行求解?；舅悸肥菍⑦B續(xù)空間劃為離散網(wǎng)格構(gòu)成的空間，然后使用有限差分法進(jìn)行求解。但是，當(dāng)網(wǎng)格變大的時候，不僅方程式求解的速度會變得很慢，而且速度和密度的流動也會變得很慢。將方程式進(jìn)行并行處理可以解決速度變慢的問題。

4.2 并行處理

在式（3）,（4）,（5）,（6）中，我們利用(i,j)數(shù)據(jù)的上，下，左，右來計(jì)算(i,j)位置的數(shù)據(jù)，如圖2所示。比如，求解下列數(shù)據(jù)類似的計(jì)算方法可以計(jì)算在索引(i+2,j)數(shù)據(jù)的計(jì)算。如圖3所示。

因?yàn)閚eon指令集是按照線性順序存儲數(shù)據(jù)的，所以我們可以并行處理我們的數(shù)據(jù)。

（1）比如，我們定義4個寄存器Q0，Q1，Q2，Q3：

Q0寄存器加載從索引(i-1,j)到(i+2,j)的數(shù)據(jù)，Q1寄存器加載從索引(i+1,j)到(i+4,j)的數(shù)據(jù)，Q2寄存器加載從索引(i,j+1)到(i+3,j+1)的數(shù)據(jù)，Q3寄存器加載從索引(i,j-1)到(i+3,j-1)的數(shù)據(jù)。

（2）我們使用寄存器使用計(jì)算。比如，使用neon指令集計(jì)算擴(kuò)散方程式（6）：

Q0=從U加載數(shù)據(jù)(從索引i-1,j到i+2,j);Q1=從U加載數(shù)據(jù)(從索引i+1,j到i+4,j)；

Q2=從U加載數(shù)據(jù)(從索引i,j+1到i+3,j+1);Q3=從U加載數(shù)據(jù)(從索引i,j-1到i+3,j-1)。

寄存器 Temp1=Add(Q1,Q0); Temp1=Add(Temp1,Q2); Temp1= Add(Temp1,Q3);

寄存器 Temp2=Multiply(Temp1,scalar); Temp2=Add(Q4,Temp2);

（3）將(i,j)到(i+3,j)相加的結(jié)果存儲在寄存器R里，然后將R里存儲的數(shù)據(jù)再寫入內(nèi)存(圖4）。

使用neon指令集求解投射方程式（5）和（7）：

寄存器 Q0=從U加載數(shù)據(jù)(從索引i-1,j到i+2,j);寄存器 Q1=從U加載數(shù)據(jù)(從索引i+1,j到i+4,j);

寄存器 Q2=從V加載數(shù)據(jù)(從索引i,j+1到i+3,j+1);寄存器 Q3=從V加載數(shù)據(jù)(從索引i,j-1到i+3,j-1)。

寄存器 temp1=Substract(Q1,Q0); 寄存器 temp2=Substract(Q2,Q3);

寄存器 R=Add(temp1,temp2); R=Multiply(scalar,R); 存儲數(shù)據(jù)(內(nèi)存(i,j),R));

使用neon指令集求解方程式（5）：

寄存器 Q0=從P加載數(shù)據(jù)(從索引i-1,j到i+2,j);寄存器 Q1=從P加載數(shù)據(jù)(從索引i+1,j到i+4,j)；

寄存器 Q2=從P加載數(shù)據(jù)(從索引i,j+1到i+3,j+1);寄存器 Q3=從P加載數(shù)據(jù)(從索引i，j-1到i+3,j-1);

寄存器 Temp=Add(Q0,Q1); Temp=Add(temp,Q2); Temp=Add (temp,Q3);

寄存器 R=Multiply(temp,scalar); 存儲數(shù)據(jù)(內(nèi)存(i,j),R);

外部力：U (x,T+ΔT)= U (x,T)+ ΔTSource(x,T)

使用neon指令集求解：

寄存器 X=從U加載數(shù)據(jù)；寄存器 S=從外部力加載數(shù)據(jù)；寄存器 scalar=存儲 ΔT.

寄存器 R=vmlaq_f32(X,S,scalar); 存儲數(shù)據(jù)(內(nèi)存，R);

由于neon指令集是線性順序處理，因此半拉格朗日算法不適用，所以對對流不使用并行處理。