[Paper]Real-Time Fluid Dynamics for Games

[Shimo888]

Jos Stam. Real-Time Fluid Dynamics for Games, 2003

Stab

• ゲームエンジン向けの流体ソルバー(格子法)
• 正確さより見た目重視のソルバー
• 安定性と計算速度が優れている(絶対発散しないらしい)
• StableFluidsと同じ著者

関連

#9 (Stable Fluids)

[Shimo888]

数理モデル

[密度]

$$\frac{\partial \rho}{\partial t} = -(\mathbf{u} \cdot \nabla)\rho + \kappa \nabla^{2} \rho + S$$

[流速]

$$\frac{\partial \mathbf{u}}{\partial t} = P\Big( - (\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} \Big)$$

Note

流速の式は論文の式から少し変更を加えた。
Projectionの式が抜けていたので、現論文Stable Fluid同様記述している。
密度の式は、StableFluidの時から減衰項が消えているが、シミュレーションしたい状況に依存するとおもうので気にしないことにした。

[Shimo888]

数値計算 -全体-

データ構造

#define IX(i,j) ((i)+(N+2)*(j))
static float u[size], v[size], u_prev[size], v_prev[size];
static float dens[size], dens_prev[size];

size=(N+2)*(N+2)であり、2次元グリッドを1次配列で表現している(速度の水平成分のセル(i,j)は、 u[IX(i,j)] という要素で与えられる)。

全体計算ステップ

while ( simulating )
{
    // UIから入力（ソースや力）を取得し、_prev配列に格納
    get_from_UI ( dens_prev, u_prev, v_prev );
    
    // 速度場を1ステップ進める
    vel_step ( N, u, v, u_prev, v_prev, visc, dt );
    
    // 密度場を1ステップ進める
    dens_step ( N, dens, dens_prev, u, v, diff, dt );
    
    // 計算結果の密度を描画
    draw_dens ( N, dens );
}

時間間隔 dt

[Shimo888]

数値計算 -密度-

0. 全体

void dens_step ( int N, float * x, float * x0, float * u, float * v, float diff, float dt )
{
    add_source ( N, x, x0, dt );
    SWAP ( x0, x ); 
    diffuse ( N, 0, x, x0, diff, dt );
    SWAP ( x0, x ); 
    advect ( N, 0, x, x0, u, v, dt );
}

1. ソースの追加

ユーザーの入力した値などを大体追加する

void add_source ( int N, float * x, float * s, float dt )
{
    int i, size=(N+2)*(N+2);
    for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i];
}

2. 拡散

[陰解法 + ガウスザイデル法]

void diffuse ( int N, int b, float * x, float * x0, float diff, float dt )
{
    int i, j, k;
    float a = dt*diff*N*N;
    for ( k=0 ; k<20 ; k++ ) { // 反復回数
        for ( i=1 ; i<=N ; i++ ) {
            for ( j=1 ; j<=N ; j++ ) {
                x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+
                                               x[IX(i,j-1)]+x[IX(i,j+1)])) / (1+4*a);
            }
        }
        set_bnd ( N, b, x );
    }
}

Note

上の方法で計算する
ラプラシアンを離散化して、次ステップの項とすると連立方程式が導ける
解が収束すると仮定して漸化式をたて反復しているだけ(それが上のガウスザイデル法)
具体的なステップは以下である。

1. ラプラシアンの離散化は中心差分を行うことで以下のように示せる

$$\nabla ^2 {\rho} \approx \frac{\rho(x- dr,y) + \rho(x+ dr,y) + \rho(x,y -dr) + \rho(x, y +dr) - 4 \rho(x, y)}{dr^2}$$

2. 上の式を用いて拡散方程式を陰解法で離散化すると、

$$\frac{\rho(t + dt) - \rho(t)}{dt} = a' \frac{\rho(x- dr,y) + \rho(x+ dr,y) + \rho(x,y -dr) + \rho(x, y +dr) - 4 \rho(t +dt)}{dr^2}$$

3. 上の式を $\rho(t + dt)$でくくって式変形すると

$$\rho(t+dt) = \frac{ \rho(t) + \frac{a' \, dt}{dr^2} \left( \rho(x- dr,y) + \rho(x+ dr,y) + \rho(x,y -dr) + \rho(x, y +dr) \right)}{1 + \frac{4a'\,dt}{dr^2}}$$

4. $a=\frac{4a',dt}{dr^2}$とおけば上の式が求まる
5. あとは、x_ijが20回の反復で収束すると仮定して計算を行うだけ

その他、格子全体の大きさが1であると規格化している ($dr N = 1$)

Note

ガウスザイデル法の詳細については別Issueでまとめた
並列計算する際は、ガウスザイデル法が使用できないので、ヤコビ法を用いる(k回目反復計算中にk回目の結果を使用するかどうかだけの違い)
係数行列をまとめずに、二重ループ内で係数行列に対応する部分を求めているので少しわかりにくいが、本質的な計算方法は同じ
近接の4点が影響をあたえるので、係数行列は疎行列であることはわかる
でも計算途中で境界条件設定したりしてるので厳密なガウスザイデル法とは異なる
でも「反復して計算すると解が近づくよ」程度に捉えることにした

[陽解法]

void diffuse ( int N, int b, float * x, float * x0, float diff, float dt )
{
    int i, j, k;
    float a = dt*diff*N*N;
    for ( k=0 ; k<20 ; k++ ) { // 反復回数
        for ( i=1 ; i<=N ; i++ ) {
            for ( j=1 ; j<=N ; j++ ) {
                x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+
                                               x[IX(i,j-1)]+x[IX(i,j+1)])) / (1+4*a);
            }
        }
        set_bnd ( N, b, x );
    }
}

Note

普通に陽解法で離散化すると上のようになる(拡散方程式とかシミュレーションするときにするやつ)。
連立方程式解く必要ないので計算量とか少ないが、数値的に不安定になっちゃうので通常は使用しない

3. 移流

void advect ( int N, int b, float * d, float * d0, float * u, float * v, float dt )
{
    int i, j, i0, j0, i1, j1;
    float x, y, s0, t0, s1, t1, dt0;

    dt0 = dt*N; // グリッド空間での時間ステップ
    for ( i=1 ; i<=N ; i++ ) {
        for ( j=1 ; j<=N ; j++ ) {
            // 時間を逆向きにたどり、過去の位置(x,y)を求める
            x = i - dt0*u[IX(i,j)]; 
            y = j - dt0*v[IX(i,j)];
            
            // 座標がグリッド範囲内に収まるように制限 (クランプ)
            if (x<0.5) x=0.5; if (x>N+0.5) x=N+0.5; 
            i0=(int)x; i1=i0+1;
            if (y<0.5) y=0.5; if (y>N+0.5) y=N+0.5; 
            j0=(int)y; j1=j0+1;
            
            // 補間のための重みを計算
            s1 = x-i0; s0 = 1-s1; 
            t1 = y-j0; t0 = 1-t1;
            
            // 双線形補間で過去の密度を求め、現在のセルに設定
            d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)] + t1*d0[IX(i0,j1)]) +
                         s1*(t0*d0[IX(i1,j0)] + t1*d0[IX(i1,j1)]);
        }
    }
    set_bnd ( N, b, d );
}

Note

移流項をバックトレースで表現するアイディアが、Stable Fluidの核の一つ。これで安定する。
数値的な詳細は、原論文（これじゃない）方の補足にまとめれてる。
双線形補間は、四角形の4点の値から補完する方法
まず上下の辺の値をx方向で補完して求める
次に求めた値を用いてy方向で補完する(図にかいたら簡単)

[Shimo888]

数値計算 -流速-

0: 全体

void vel_step ( int N, float * u, float * v, float * u0, float * v0, float visc, float dt )
{
    // 力の加算
    add_source ( N, u, u0, dt ); 
    add_source ( N, v, v0, dt );
    
    // 拡散
    SWAP ( u0, u ); 
    diffuse ( N, 1, u, u0, visc, dt );
    SWAP ( v0, v ); 
    diffuse ( N, 2, v, v0, visc, dt );
    
    // 射影
    project ( N, u, v, u0, v0 );
    
    // 自己移流
    SWAP ( u0, u ); 
    SWAP ( v0, v );
    advect ( N, 1, u, u0, u0, v0, dt ); 
    advect ( N, 2, v, v0, u0, v0, dt );
    
    // 再び射影
    project ( N, u, v, u0, v0 );
}

1. 射影

void project ( int N, float * u, float * v, float * p, float * div )
{
    int i, j, k;
    float h = 1.0/N;

    // 1. 速度場の発散(divergence)を計算する
    for ( i=1 ; i<=N ; i++ ) {
        for ( j=1 ; j<=N ; j++ ) {
            div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+
                                   v[IX(i,j+1)]-v[IX(i,j-1)]);
            p[IX(i,j)] = 0;
        }
    }
    set_bnd ( N, 0, div ); 
    set_bnd ( N, 0, p );
    
    // 2. ポアソン方程式を解き、圧力場pを求める (ガウス＝ザイデル法)
    for ( k=0 ; k<20 ; k++ ) {
        for ( i=1 ; i<=N ; i++ ) {
            for ( j=1 ; j<=N ; j++ ) {
                p[IX(i,j)] = (div[IX(i,j)] + p[IX(i-1,j)] + p[IX(i+1,j)] +
                                            p[IX(i,j-1)] + p[IX(i,j+1)]) / 4;
            }
        }
        set_bnd ( N, 0, p );
    }
    
    // 3. 圧力pの勾配(gradient)を速度場から差し引く
    for ( i=1 ; i<=N ; i++ ) {
        for ( j=1 ; j<=N ; j++ ) {
            u[IX(i,j)] -= 0.5 * (p[IX(i+1,j)] - p[IX(i-1,j)]) / h;
            v[IX(i,j)] -= 0.5 * (p[IX(i,j+1)] - p[IX(i,j-1)]) / h;
        }
    }
    set_bnd ( N, 1, u ); 
    set_bnd ( N, 2, v );
}

Note

一見意味不明だが、現論文のStableFluidを見れば納得できる。
まずは数式の定義から見る

1. ヘルムホルツホッジ分解の定義は以下である
   任意のベクトルuは、発散が0のベクトルとあるスカラー量pの勾配の和で表せる。

$$u = u_{div0} + \nabla p$$

2. このとき演算子Pは、uからu_{div0}を取り出す演算として定義できる

$$Pu = u_{div0}$$

3. 1の定義式の両辺で発散を求めると、以下の式が求める

$$\nabla \cdot u = \nabla ^2 p$$

つまりP演算は、3のポアソン方程式の解pを求め、uからgrad pを差し引くことに相当する

Note

上のコードでは
1と2でヘルムホルツホッジ分解から導出されるポアソン方程式をとく
３でヘルムホルツホッジ分解の、u_{div0}を求める
だけのシンプルコード

[Shimo888]

数値計算 -境界条件-

[反射境界条件]

void set_bnd ( int N, int b, float * x )
{
    int i;
    for ( i=1 ; i<=N ; i++ ) {
        // 左右の壁 (x[0, i] と x[N+1, i])
        x[IX(0, i)]   = b==1 ? -x[IX(1,i)] : x[IX(1,i)];
        x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)];
        
        // 上下の壁 (x[i, 0] と x[i, N+1])
        x[IX(i, 0)]   = b==2 ? -x[IX(i,1)] : x[IX(i,1)];
        x[IX(i, N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)];
    }
    // 四隅のセルは隣接するセルの平均値をとる
    x[IX(0, 0)]       = 0.5*(x[IX(1,0)]   + x[IX(0,1)]);
    x[IX(0, N+1)]     = 0.5*(x[IX(1,N+1)] + x[IX(0,N)]);
    x[IX(N+1, 0)]     = 0.5*(x[IX(N,0)]   + x[IX(N+1,1)]);
    x[IX(N+1, N+1)] = 0.5*(x[IX(N,N+1)] + x[IX(N+1,N)]);
}

Note

反射させているだけなのでめっちゃ簡単
bはxかyどっちの成分なのかを指定する
b=1のときは左右の壁で値が逆転するので、x成分
b=2のときは上下の壁で値が逆転するので、y成分
b=0のときは、境界値が固定になるだけ（Neuman境界条件に対応）