流体计算的主要问题就是解NS方程，至于怎么解就是各家之长了。Alias的大牛Stam在上世纪末提出了的稳定解算流体的做法，10年之后来瞅瞅……

把流体解算过程看作一个积分过程，但积分一不小心就会整成不稳定的，原因很多，如步长大于某种方法的稳定步长啦，没把质量守恒、能量守恒考虑进去啦，等等。

Stam用了个隐式积分的方法解决稳定积分的问题，以前是直接显式积分（加dt），现在是倒退回dt的时间，拿t-dt时刻的量，但这额外需要计算一个线性方程组。由于这个线性系统非常稀疏，可以使用Gauss-Seidel迭代法解这个系统（试下来性能令人满意）。

好，稳定的积分方法解决了，现在看质量守恒。通常根据NS方程几个项（平流、扩散、外力作用）积分下来，通常不会是质量守恒的。如果把符合质量守恒的状态可视化成所有一个平面的点集，那么每次积分的结果都可能把一个点（状态）带离这个平面，也就是破坏了质量守恒的条件。那么我们要在积分末尾进行一个projection的动作，把状态投影到这个平面上，也算满足了质量守恒的先决条件。可能这个投影的动作不那么科学，但至少满足了“稳定”的需求。投影的依据是Helmholtz-Hodge分解法，这个方法指出：任何向量场（我们所积的速度场）都可以唯一分解为一个收敛的向量场（即不发散）和一个标量的梯度场。
w = u + del(q) ==> u = w - del(q)
（del是del operator，向量微分算子，是一个记号，和向量一起使用才有微分的意义）
那么我们只要取这个收敛的向量场作为速度场即可，这个操作就是投影。至于求这个q，又需要解一个线性系统，
del(w) = del^2(q)
这是个泊松方程，也存在一些数值方法可以解它。总之得出q后也得到了最后的投影操作，也完成了整个解算的过程。

下面是在PS3上试验的结果，二维、用CPU更新的纹理，还没捣腾出截屏功能，凑合看……