Real-Time Fluid Dynamics for Gamesを読んでみる

前回、前々回に続き、Jos StamさんがGDC’03で発表した、Real-Time Fluid Dynamics for Gamesを読んでみたと思います。

Introduction

Fuildの流れはどこにでもあります、煙の立ち上る様や、雲、霧、川の流れ、そして海など。
プレイヤーを説得力のあるvisualの世界に没頭させるという、ゲームにおいて一つの重要な要素のために、fuildをゲームエンジンに実装するのは好ましいことです。
すでに多くの、ad-hocな手法による、スプライトや、パーティクルを用いた、fuildのようなエフェクトの実装の試みはたくさんあります。
しかしながら、これ等を正しい物理法則に従わせてアニメーションさせることは簡単ではありません。

我々は、Euler, Navier, Stokesの時代(1750~1850)に発展した、流体力学を利用するほうが、もっと良い方法だと思います。
これらの発展は、自然界に存在するfluidの流れの詳細な数学モデルである、Navier-Stokes方程式と呼ばれるものを導き出しました。
しかしながら、これ等の方程式は、非常な単純なケースの、解析ソリューションに利用されるのみでした。
1950年代に入り、研究者たちがコンピューターを用い、この方程式を解くために数値解析的手法を用いるようになるまで、発展はありませんでした。
一般的にこれ等の手法は、精度を重視し、非常に複雑で、計算時間の掛かるものです。これは、これ等の応用分野が、計算結果が物理的に正確である必要があったためです。
飛行機や橋にかかる負荷が正確に計算される必要があるのは、当たり前のことです。

CGやゲームにおいては、まったく逆で、それらしく見えて、はやく計算できることが重要です。
加えて、ソルバーが複雑すぎず、普通のPCやゲームコンソールやPDAなどで実装できることが重要です。
この論文では、これ等の要求にマッチするアルゴリズムを紹介します。
これ等のゴールを達成するため、まずは従来の計算物理学の知識よりスタートし、ビジュアルエフェクトを作るための、専用のアルゴリズムに発展させます。
物理的に正確な、(シミューレーションの)時間間隔にシビアなソルバーと異なり、我々のソルバーはstableで、決して”blow up”しません。

Fig.1 (訳注:オリジナルの論文を参照して下さい)

数学的に、fluidの状態は、ある時間における、velocity vector field(速度ベクトルが空間上の全ての点で定義された関数)でモデル化されます。
部屋に空気が充満している状態を想像して下さい。空気の速度は、熱源や、空気の流れなどで変化します。
たとえば、暖房機のそばの速度は、立ち上る熱のせいで、上方向に支配的なものでしょう。
タバコより立ち上る煙や、空中を舞うホコリを見れば分かるように、部屋の中の速度の分布は、とても複雑なものです。

Navier-Stokes方程式は、velocity fieldの時間的進展の詳細な定義です。
与えられた現在の速度と、外力の状態より、この方程式は、我々に、微小時間において、速度がどのように変化するかを、詳細に教えてくれます。
Fig.1(上)はこれ等の方程式を、簡易ベクトル表記であらわしたものです。
大まかにいって、この式は、右辺の三つの項による速度の変化を表したものです。
velocity field自身は、煙のパーティクルや、塵、葉っぱなどが動き始めるまで、それほどビジュアル的に面白くありません。
これらのオブジェクトの動きは、オブジェクト周辺の速度を、オブジェクトに対する力としてコンバートすることで計算されます。

塵の様に軽いオブジェクトは、通常、単にvelocity filedに従って運ばれます。
煙の場合は、粒子ごとにモデル化することは、非常にコストが掛かります。
そのため、このケースでは、煙の粒子は、煙の濃度(density)(全ての空間の点において、パーティクルの存在量を示す連続関数)で置換されます。
濃度は、通常0~1の範囲の値をとります。(煙が無ければ0。そうでなければ、粒子の存在量を示す)
このdensity fieldが、fluidのvelocity filed内で、進展する様は、Fig.1(下)に示される、詳細な数学的方程式で説明できます。

本論文の読者は、これ等の方程式を完全に理解することは期待されていません。
しかし、皆様にもこれ等の2つの式が非常に似通っていることは明白だと思います。
実は、この類似は、我々の手法の開発において、助けとなっています。
density(下)の方程式は、velocity(上)のものよりシンプルです。専門的に言えば、前者が線形で、後者が非線形です。
まず我々は、densityが、固定されたvector field内を動く様を計算するアルゴリズムを開発しました。
そしてこれが、velocity fieldの進展の計算にも同様に適用できるかも知れないことに気づきました。
本論文では、この歴史的な経緯に沿って進めていきます。まず始めに、どのようにdensityに関する方程式を解くかを示します。
これは、読者の、他の要素のソルバーに対する理解を深めるでしょう。使用されたコンセプトは、簡単に説明、映像化できます。
続けて、これ等のアイデアをもっと難しい問題である、velocity fieldのシミュレーションにもっていきます。

A Fuild in a Box

訳注:省略します。大まかな内容は、本論文では、説明の明確さを重視し、2次元のグリッドを用いてシミュレーションを行い、計算グリッド周辺に1要素余分に確保し、境界条件を適用します。
また、論文中のプログラムでは一次元配列を用いてメモリを確保するので、アクセス用のマクロを定義しています。またシミュレーションの時間間隔はdtとします。

Moving Densities

Fig.3 (訳注:オリジナルの論文を参照して下さい)

先に説明したとおり、まず、density filedの、時間的に変化しない、固定されたvelocity fieldに対するソルバーを説明します。
もう一度、Fig.1(下)で示された、densityに関する方程式について考察します。この式は、densityが3つの要因によって、時間経過に応じて変化する様を定義します。
これらの3つの要因とは、式の右辺の3つの項です。最初の項はdensityはvelocity fieldに従って動くことを示します。2番目の項は、densityは一定の確率で拡散することを示し、3番目の項は、外からの要因でdensityは上昇することを示します。
我々のソルバーはこれ等の項を、Fig.3で示したとおり、逆順で解いていきます。
初期状態のdensityより開始し、これ等の3つの項を、時間間隔ごとに、周期的にといていきます。

最初の項は、実装が簡単です。該当のフレームの外的要因の、sourceが配列s[]で提供されるとします。
この配列は、ゲームエンジンのsouceを検出する部分によって提供されます。我々のプロトタイプでは、この配列は、ユーザーのマウスの動きによって供されます。
sourceをdensityに加えるルーチンは単に以下のようになります。

Fig.4 (訳注:オリジナルの論文を参照して下さい)

2番目のステップでは、確率diffで示される拡散を計算します。diff>0の時、densityはグリッド上で広がっていきます。
まず始めに、一つのセル上で何が起きているかを考えます。このケースで、我々は、Fig.4で示す様に、直接隣接する4つのセルとdensityを交換するものとします。
対象セルのdensityは近隣のセルに移動することで減少しますが、同時に、近隣のセルから流れ込むことによって上昇もします。
結果的には、これらの差分である、
x0[IX(i-1, j)] + x0[IX(i+1, j)] + x0[IX(i, j-1)] + x0[IX(i, j+1)] -4*x0[IX(i, j)]
になります。

拡散のソルバーの考えられる実装は、単にこれ等の交換を各セルごとに計算し、元々の値に加算することです。これは結果的に以下の簡単な実装になります。

void diffuse_bad() {} (訳注:オリジナルの論文を参照して下さい)

set_bnd()ルーチンは、境界セルの値を設定するもので、後ほど説明します。
この拡散のルーチンは、素直で魅力的に見えるのですが、残念ながらうまく動作しません。
高い拡散率aを用いると、densityの値は振動し始め、負の値になり、最後には発散してしまい、ついにはシミュレーションが破綻してしまいます。
この挙動は、よく知られる問題で、unstableな手法では厄介な問題です。これらの理由から、我々はstableな手法を拡散のステップのために考えます。
我々の手法の背景となる基礎的なアイデアは、求めるべき拡散されたときのdensityを時間的に遡り、開始時点のdensity計算することです。
コードでは以下の通りです。
x0[IX(i, j)] = x[IX(i, j)] – a * (x[IX(i-1, j)] + x[IX(i+1, j)] + x[IX(i, j-1)] + x[IX(i, j+1)] – 4*x[IX(i, j)])

これは、未知数x[IX(i, j)]の線形システムになります。
我々は、この線形システムに対してマトリクスを作成し、通常の逆マトリクスを求めるルーチンを呼び出すことが出来るかもしれません。
しかしこれは、この問題に対しては過剰です。なぜならば、このマトリクスは非常に疎で、ごく一部の要素のみが、非ゼロだからです。
従って、我々はこの逆マトリクスに対して、単純で、反復的手法を用いることが出来ます。
実際にうまく動作し、最もシンプルな解法は、Gauss-Seidel法です。以下がその実装になります。

void diffuse() {} (訳注:オリジナルの論文を参照して下さい)

このバージョンの拡散ソルバーの美点は、unstableなバージョンと、ほとんどシンプルさは変わらないにもかかわらず、いかなるdiff,dt,Nを扱うことができ、どれだけ値を大きくても、シミュレーションが”blow up”しません。

Fig.5 (訳注:オリジナルの論文を参照して下さい)

densityソルバーの最後のステップである、与えられたvelocity fieldに従ってdensityが動く項について説明します。
Fig.5を参照して下さい。我々は、stableで”blow up”しない技法を使いたいと思います。
拡散の計算をした際と同様に、線形システムを構築し、Gauss-Seidel法を用います。
しかしながら、この線形方程式はvelocityに依存し、このことにより、解くのが困難になっています。
幸運なことに、もっと効果的な代替方法が存在します。
この新しい手法の肝となるのは、densityの動きを、粒子にモデル化することができれば、簡単に解く事が出来るということです。
この場合は単純に、velocity filed内で粒子をトレース必要があります。
たとえば、計算グリッドの各セルの中心を、粒子に見立てて、Fig.6(b)で示したように、velocity field内をトレースします。
問題は、これ等の粒子の値を、グリッドのセルの値へと変換し直さなければならないということです。
これをどのように行うかということは、必要無いことが明らかです。
もっと良い方法は、Fig.6(c)に示したとおり、1回のシミュレーションステップを経て、丁度グリッドのセルの中心に来る粒子を見つけることです。
これらの粒子が運ぶdensityの量は、粒子の開始点の近傍4セルの線形補完で簡単に求まります。
以上は、以下のdensityの更新手順につながります。
2つの計算グリッドを用意します:一つは、1計算ステップ前の状態のdenistyを保持し、もう片方に、更新した値を格納します。
更新した値を計算するために、各セルの中心位置から、逆方向にvelocity fieldをトレースします。
その後、前の状態のdensityを保持するグリッドの値を線形補完し、新しいグリッドのセルに値を格納します。

Fig.6 (訳注:オリジナルの論文を参照して下さい)

以下が、このアイデアの実装コードになります。単純な線形トレースを使用しています。

void advect() {} (訳注:オリジナルの論文を参照して下さい)

以上でdensityソルバーに関する説明を完了します。これ等の手順は、簡単に一つのルーチンへとまとめることが出来ます。
x0にはsourceが格納されているとします。

void dens_step() {} (訳注:オリジナルの論文を参照して下さい)

SWAPマクロは配列のポインタを交換します。

#define (訳注:オリジナルの論文を参照して下さい)

[ここまでで。感想]

氏のSiggraph’99のStable Fluidと異なり、比較的平易に書いてあるので読みやすいです。加えて実装のコードが理解を助けてくれます。
一方で、この発表は、背景になる部分をかなり省略したものだとも思います。しかし、Stable Fluidの論文とあわせて読むことで、その背景部分を補うことが出来ると思います。
densityの計算部分はStable Fluidの論文では余り触れられていませんが、最後のほうに少し記述があります。また、densityの計算自体は、流体シミュレーションの計算に必須というわけではなく、
velocity filedの計算さえ出来れば、その中を動く物質のモデル化は、濃度モデルでも粒子モデルでも好きなほうを用いればよいと思います。説明の都合上、velocity fieldに先立って説明されているに過ぎません。
文中で使用されているGauss Seidel法は、線形の連立方程式を解く際には有用で、対角成分が優位な場合に良好に収束するため、今回のケースでは、うまく動作することが想像できます。
移流項の計算で使用される、パーティクルのトレーサーは、Stable Fluidの論文ですと、2次のRunge Kuttaを用いたと記されていますが、本論文では、単純な線形補完で記述されています。
しかし、元となっているアイデアは同じです。

次回はvelocity filedの計算部分を読んでみたいと思います。

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中