張鵬

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

1 緒論

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

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

流體模擬是21世紀流體力學的重要組成部分，在電影，游戲以及其它領域有廣泛的應用。但是，伴隨著流體計算量的增加，流體模擬的速度會變得非常慢。我們的目標就是研究在iOS系統里實現流體模擬，并且使用ARM的命令集進行并行加速。

2 流體模擬的含義

2.1 先前研究

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

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

2.2 數學方程式

圖1

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

這里P為投射算子。

2003年，Jos Stam發表了他的研究“Real-Time Fluid Dynamics for Games”[3]，使用密度替換速度，將上面方程式（1）變為：

3 在IOS系統的實時流體模擬的實現

OpenGL ES是OpenGL開發的基于嵌入式版本的3維計算機圖形API。我們使用OpenGL ES在IOS系統進行繪制，但是在IOS使用OpenGL ES進行繪制之前，我們需要先生成繪制環境。

根據Apple提供的在IOS使用OpenGL ES編程向導[4]，我們可以根據幾個步驟進行IOS繪圖環境設置：

設置EAGLContext→建立OpenGL繪圖內容→建立渲染緩沖→建立幀緩沖→展現

3.1 設計流體對象

以前面的方程式為基礎，我們設計流體對象包括部分如圖1所示。

3.2 流體對象的實現

為了解開原有的流體方程式，需要將連續方程變為離散方程。在歐拉方法的基礎上，我們利用有限差分法將原有流體方程式的各個項拆分為以下形式：

在這里，i為水平坐標，j表示為垂直坐標。為了判定邊界條件，我們追加了新的邊界，將網格大小設定為N+2。

3.2.1 外部力

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

3.2.2 擴散

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

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

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

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

3.2.3 對流

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

3.2.4 投射

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

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

最終可以求解得出U4

3.2.5 邊界條件

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

3.3 在iPhone的繪制

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

3.4 問題

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

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

4.1 以前算法的問題點

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

4.2 并行處理

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

因為neon指令集是按照線性順序存儲數據的，所以我們可以并行處理我們的數據。

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

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

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

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

Q2=從U加載數據(從索引i,j+1到i+3,j+1);Q3=從U加載數據(從索引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)相加的結果存儲在寄存器R里，然后將R里存儲的數據再寫入內存(圖4）。

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

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

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

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

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

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

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

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

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

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

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

使用neon指令集求解：

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

寄存器 R=vmlaq_f32(X,S,scalar); 存儲數據(內存，R);

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