使用FFT的WebGL墨滴模拟
Evgeny Demidov写了很多精彩的WebGL例子。其中我最感兴趣的是墨滴模拟。我把计算方法做了点改动，自己实现了一个：（点击图片）
原来的页面用的方法是Jos Stam的Real-Time Fluid Dynamics for Games和Stable Fluid，使用的是交错网格。不过在输运（不知道这是不是Advect的正确翻译）变量时，有些小问题，作者在页面上也做了注明。我总觉得其中以雅可比迭代解泊松方程显得有点“原始”，自己以前接触过用快速傅立叶解泊松方程，所以想试试。实际上如果使用FFT，不需要构造泊松方程，见Jos Stam的A Simple Fluid Solver based on the FFT，本例就是用的此方法。
关于流体力学的基本知识，推荐Robert Bridson的SIGGRAPH 2007 Course。前列诸文献写得很好，这里只简要介绍一下基本知识和我使用的方法。对于这个仿真，可认为流体运动遵循不可压缩纳维叶－斯托克斯方程：
$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \nabla \mathbf{u} + \frac1\rho \nabla p = \mathbf{f} + \nu \nabla\nabla \mathbf{u}$$$$\nabla \mathbf{u} = 0$$前一个是动量方程，由牛顿第二定律得来，后一个为质量守恒方程。其中$\mathbf{f}$为外力产生的加速度，在本例中$\mathbf{f}=\mathbf{g}c$，$\mathbf{g}$为重力加速度，$c$为墨的浓度。而墨的浓度又会随液体流动变化：
$$\frac{\partial c}{\partial t} + \mathbf{u} \nabla c = 0$$综合起来，共是三个偏微分方程。一种较简单的方法“分裂法”，即是将复杂方程分成几个简单方程分别求解，再将结果叠加起来。本例中，先施加重力，接着输运速度，再令速度满足质量守恒方程，最后用得到的速度输运墨的浓度。其中最复杂的是令速度满足质量守恒方程，即使速度场为无源场。这过程又被称为亥姆霍兹分解，指的是将矢量场分解为无源场和无旋场之和。本例中使用快速傅立叶来完成此分解，不需要计算压力值。利用FFT，可将速度表示为：
$$u_x = \sum_{k = 0}^{M-1}\sum_{l = 0}^{N-1}\mathcal{F}(u_x)_{kl}e^{2\pi i(\frac kM x+\frac lN y)}$$$$u_y = \sum_{k = 0}^{M-1}\sum_{l = 0}^{N-1}\mathcal{F}(u_y)_{kl}e^{2\pi i(\frac kM x+\frac lN y)}$$使速度场为无源场，即$\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y} = 0$，将上两式代入，可得对每一个$k,l$：$\frac kM\mathcal{F}(u_x)_{kl} + \frac lN\mathcal{F}(u_y)_{kl} = 0$。按照这个关系，将FFT的结果投影到和$(k,l)$垂直的方向上，就得到了一个无源场。“A Simple Fluid Solver based on the FFT”文章中的图形形象地表达了这个过程。
在实现方面，如果是用C之类的语言的话，实现起来应比较容易，因为有现成的快速傅立叶库可用。不过我没找到之前有人写过基于WebGL的快速傅立叶实现，只好自己写一个。由于这方面我之前没接触过，实现起来费了不少力气。我采用的方法是The FFT on a GPU。尽管我也只是半懂，感觉其中的“Index Magic”还是很有意思的。文中介绍的是基为2时的情形，但要扩展到基为其它值以及向量基的情况是很容易的。因为计算的是两组实数的FFT，还采用了文中的将两实数组成一复数计算FFT的策略。这一步骤在文中被称为“Tangle/Untangle”。不过这一过程可以在完成二维FFT后再进行，不必在每一方向的FFT后都加上一步，因为二维实数FFT也有和一维相应的共轭对称关系。采用和文中相似的记法，$M\times N$的二维实数$f(x,y),g(x,y)$对应的FFT结果为$F(u,v),G(u,v)$，令$h(x,y)=f(x,y)+jg(x,y)$，计算FFT得$H(u,v)$。可由$H(u,v)$得到$F(u,v),G(u,v)$：
$$F(u,v)_R = \frac12 \left(H(u,v)_R+H(M-u,N-v)_R\right)$$$$F(u,v)_I = \frac12 \left(H(u,v)_I-H(M-u,N-v)_I\right)$$$$G(u,v)_R = \frac12 \left(H(u,v)_I+H(M-u,N-v)_I\right)$$$$G(u,v)_I = -\frac12 \left(H(u,v)_R-H(M-u,N-v)_R\right)$$求逆变换时，反过来利用公式就行了。
其实之中还有不少细节我还没完全清楚，还望大家多多指点。
TAG:WebGL FFT
