使用FFT的WebGL墨滴模拟

Evgeny Demidov写了很多精彩的WebGL例子。其中我最感兴趣的是墨滴模拟。我把计算方法做了点改动,自己实现了一个:(点击图片)

FFT based droplet

原来的页面用的方法是Jos Stam的Real-Time Fluid Dynamics for GamesStable 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)$$

求逆变换时,反过来利用公式就行了。这样向量基的FFT算法也可应用了。本仿真中就用了向量基的FFT,比起做两个方向的FFT要快不少。不过我还没懂究竟怎样的代码在WebGL中较快,因我发现这和运行在CPU上的代码有不同。

其实之中还有不少细节我还不完全清楚,还望大家多多指点。

posted @ 2012-04-30 21:29 枣木夹子 阅读(...) 评论(...) 编辑 收藏