基于iOS的实时流体模拟的并行处理

张鹏

（郑州科技学院信息工程学院 河南省郑州市 450064）

1 绪论

随着移动智能手机的最新发展，越来越多的人使用智能手机生活。其中的一些代表智能手机，比如：携带Apple iOS系统的iPhone以及iPad，三星携带Android系统的galaxy，以及华为的手机等。

这些智能手机的处理器都具有体积非常紧凑，且功耗低的特点。现在大部分智能手机的系统都是由ARM处理器构成的，由于设计以及存储等其它缺点的存在，由ARM构成的设备性能要远低于PC电脑的配置性能。随着智能手机性能的多样化，人们使用它观看电影，玩游戏或者做更多的事情，在这种情况下，我们期待智能手机能够提高现在的处理器性能从而使我们做更多的工作。

流体模拟是21世纪流体力学的重要组成部分，在电影，游戏以及其它领域有广泛的应用。但是，伴随着流体计算量的增加，流体模拟的速度会变得非常慢。我们的目标就是研究在iOS系统里实现流体模拟，并且使用ARM的命令集进行并行加速。

2 流体模拟的含义

2.1 先前研究

在物理学中，纳维斯托克斯方程(Navier-Stokes equations)用来描述流体的运动。流体模拟，一般是利用计算机求解流体方程式后，用数字模拟的形式在计算机里模拟流体的运动。现在流体模拟的技术分为几种，最常用的是欧拉方法。这种方法将流体领域划分为特定的空间，将空间分为网格组成，利用有限差分法求解流体方程式的速度，压力等。

第二种是拉格朗日方法，将流体看做是一个个粒子，每个粒子代表流动的流体，具有流体的特征。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+Δ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 以前算法的问题点

为了解开流体方程式，传统的方法是利用有限差分方法进行求解。基本思路是将连续空间划为离散网格构成的空间，然后使用有限差分法进行求解。但是，当网格变大的时候，不仅方程式求解的速度会变得很慢，而且速度和密度的流动也会变得很慢。将方程式进行并行处理可以解决速度变慢的问题。

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指令集是线性顺序处理，因此半拉格朗日算法不适用，所以对对流不使用并行处理。