Stable Fluids(2)

前回に引き続き、Jos Stamさんが1999年のSiggraphで発表した論文を読んでみたいと思います。

2.2 Method of Solution

Eq.5を初期条件u0 = u(x, 0)から、⊿tごとに時間を進めて解いていきます。
まず、とある時間tにおいてはvelocity fieldは既知として、t+⊿t時間のfieldの状態を計算したいとします。
Eq.5を4つのステップを用いて、時間間隔⊿t分の進展を計算してきます。
直前の時間の解である、w0(x) = u(x, t)を開始点とし、順をおって、Eq.5の右辺の各項を計算してきます。
その後投影操作を行い、発散無しのfieldにします。手順の概略をFigure.1に示しました。(訳注:Fig.1はオリジナルの論文を参照してください)
Fig1_

時間t+⊿tにおける解は、最後のvelocity fieldの、u(x,t+⊿t) = w4(x)で与えられます。シミュレーションは、この手順を繰り返すことで計算することが出来ます。
ここから、各ステップでどのように計算を行うか、詳細を説明していきます。
まず、解くのが最も簡単な項は、外力fの加算(add force)です。外力fが時間間隔⊿tの間に、それほど変化しないと仮定すると、
Fig3_
この式は⊿tの間に受ける外力の影響の計算の良い近似となります。
インタラクティブなシステムでは、外力は各時間間隔の始めにのみ適用されるので、これは良い近似となります。

次の項では、流体自身の移流もしくは対流(advection or convection)を計算します。
流体内のどこかにある流れは、-(u・∇)uに従って、伝達されていきます。
この項は、Navier-Stokes方程式を非線形にしています。
ForsterとMetaxasはこれを有限差分法を用いて解いています。彼らの手法は、計算の時間間隔が十分に小さい場合、⊿t<⊿γ/|u|(⊿γは計算グリッドの大きさ)の場合は、安定した解が得られます。つまり、計算グリッドが細かい場合か、もしくはfieldの速度が早い場合は、非常に小さな時間間隔を用いて計算する必要があります。

これに対して、本手法ではまったく違ったアプローチの、無条件にstableなsolverを用います。
時間間隔が大きくても問題にはならず、シミュレーションは"blow up"しません。
本手法では、偏微分方程式の解法として知られる、特性曲線法(method of characteristics)を利用しています。
これは、stableなsolverの導出において非常に重要なので、詳細をAppendix Aに記載しておきます。

しかしながら、本手法は直感的に理解することが出来ます。
すべての流体のパーティクルは、時間間隔ごとに、流体自身の速度によって動かされます。
したがって、ある点xにおける、t+⊿tにおける速度を取得するために、点xより、velocity field:W1内を⊿t時間分遡ってトレースします。これは流れの一部を表す、パスp(x, s)を定義します。
点xにおける新しい速度は、現在xにいるパーティクルの速度で表され、⊿t時間前には、古い位置にいました。
Fig4_

Figure2が上式を説明したものです。(訳注:Fig2はオリジナルの論文を参照して下さい)この手法は幾つかの利点があります。
最も重要なことは、無条件にstableなことです。明らかに、上式を見れば分かるとおり、新しいfiledが持つ最大値は、絶対に以前のfieldが持つ値より大きくなりません。
次に、この手法は実装が簡単です。この手法が実際に必要とするのは、パーティクルの追跡と線形補完です。(次の章を参照)
この手法は、stableで実装が簡単という、コンピューターグラフィックスの流体のソルバーとして、望ましい2つの特性を保持しています。
同じような手法をユーザーが定義したvelocity fieldの中を動くdensityを動かす際にも利用します。[19]
特性曲線法は他のリサーチでも様々な形で使われております。flow fieldの可視化[13,18]や、気体シミュレーションのレンダリング[21, 5]などで利用されています。
本手法での利用は、これらとは根本的に異なります。今までの研究者が動的に動かさなかった、velocity fieldの更新に使います。

3番目のステップでは、粘性の影響を解く為の拡散方程式になります。
Fig5_
このスタンダードな方程式には、たくさんの数学的手法が確立されています。
最も素直にこの方程式を解く方法は、微分オペレーター∇2を離散化し、FosterとMetaxasが行ったように、ある時間間隔おきに計算する方法です。[7]
しかし、この手法は、粘性が大きいときに不安定になります。そのため、我々は陰的な手法(implicit method)を用います:
Fig6_

I は単位オペレーターになります。微分オペレーターを離散化すると、この式は、w3の疎な線形システムになります。このような線形システムは効率的に解く事が出来ます。

4番目のステップでは、ここまでの演算結果のfieldを、発散の無いものにするための投影を行います。
前の章で示した通り、Eq.4で定義された、Poisson問題を解く必要があります。
Fig7_

この投影演算では、良いPoissonソルバーが必要です。FosterとMetaxasは、これに似た方程式を、緩和法を用いて解いています。
しかし緩和法は、収束が悪く、一般的に多くの反復を必要とします。FosterとMetaxasは、少ない反復回数で良い結果を得られたと報告しています。
しかし、我々は異なった手法を移流(advection)の計算で使用したので、もっと精度の高い手法が必要です。

特性曲線法は、fieldが発散無しに近い状態のとき、より高精度になります。
ビジュアル的な観点で、もっと重要なのは、この投影演算がfieldを渦巻状に、もっとぐにゃぐにゃした動きにすることです。
これらの理由より、我々はもっと高精度のソルバーを投影演算に使用します。

Poisson方程式が、空間的に離散化されると、疎な線形システムになります。
つまり、投影演算と、粘性の影響のステップでは、どちらも大規模な線形方程式を解く必要があります。
たとえば、マルチグリッド法は、疎線形システムを線形時間で解くことができます。[10]
我々の、移流(advection)の計算もちょうど線形なので、我々の提案するアルゴリズムの複雑度は、O(N)になります。FosterとMetaxasのソルバーも同様の複雑度になります。
この性能は、複雑な流体に対しては、概念的には最適な状態です。どのようなアルゴリズムも、最低でも各計算グリッドを考慮して計算する必要があるからです。

Appendix A Method of Characteristics

特性曲線法は、移流方程式の類を解くために使うことができます。
Eq_AP1

aはscalar fieldで、vは定常なvector fieldになります。a0はt=0の時の、fieldを示します。
p(x0, t)は、field vにおける、位置x0,時間t=0からの流れの、特性曲線を表すとします。
Eq_AP2

ā(x0, t) = a(p(x0, t))を、位置x0,時間t=0から特性曲線に沿った値とします。
時間経過に伴う、この値の変位は、微分法の連鎖律に基づき、以下の式で表すことができます。
Eq_AP3
”
これは、このscalar値が、この流線(訳注:特性曲線)に沿っては、変化しないことを示します。
特筆すべきは、ā(x0, t) =ā(x0, 0) =a0(x0)が既知であるということです。
したがって、fieldの初期状態と特性曲線は、移流の問題の解を完全に定義します。
fieldの位置x、時間tにおける状態は、位置xより時間を遡り、特性曲線にそってトレースし、位置x0を取得し、
fieldの初期状態におけるその位置の値を評価することで、計算することが出来ます。
Eq_AP4

我々は、この手法を、移流方程式を、時間[t, t+⊿t]にわたって解くのに利用しています。
このケースでは、v=u(x,t)で、a0は時刻tにおける、velocity fieldの要素いずれかひとつに相当します。

[ここまでで。感想]

2.2章と、Appendix Aを見てみました。
またもや数値計算に疎いのでなんの事やらです。
いよいよ、前回の記事のEq.5を計算するわけですが、計算は各項を逐次的におこなっていくようです。
まず、外力の加算は問題ないと思います。
次に移流項(advect)の計算では、特性曲線法という手法を用いて計算されています。
Appendix Aで、移流方程式の特性曲線が、流れそのものであることと、特性曲線上でスカラー値aが変化しないことが示されています。
従って、scalar fieldの初期状態が分かれば、特性曲線をトレースして値をサンプリングすることで、移流項を計算したことになるというわけです。
⊿t時間分の進展を計算するには、特性曲線の⊿t時間前の位置が分かれば良いわけで、それはvector fieldを⊿t時間分遡ってトレースすることで得られるというわけです。
もっともこれは、論文中でも触れられているとおり、直感的にも理解することができます。流体の流れは、流れそのものによって流されるといってるに過ぎません。
3つ目の粘性項については、要は拡散方程式なので、時間を追うごとに拡散していくわけですが、これを離散化し、差分法(計算対象の周辺と平均化するオペレーター)で計算すると、条件によっては解が不安定になるので、これを陰的手法を用いて解くようにしているそうです。
陰的手法は、要は拡散後に得られる結果に相当するw3に対して、拡散と逆にあたるオペレーターを適用すると、拡散前の結果が手に入る、と立式します。
これは線形の連立方程式になるので、これを解く事で、w2からw3を計算するようです。
4つ目の投影操作では、ポワソン方程式を解く必要がありますが、詳細についてはあまり書いてありません。空間的に離散化して、疎線形システムになると書いてあるので、
差分法で式を立てて、連立方程式を解くものと思われます。右辺の∇・w3は、w3がこの時点で既知なので、離散化する際に計算してしまうのかもしれません。

結局、具体的には、シミュレーショングリッドの大きさの(疎)線形システムを、拡散項と投影操作の際に解く必要がありそうです。
しかし、どちらの場合でも疎な線形方程式になるのは自明です。なぜなら、個々の計算グリッドはその近辺のグリッドから最も強く影響を受けるからです。
したがって、対角成分とその近辺にのみ係数が存在し、他はほとんどが0となるはずです。ですから、グリッドの大きさの割には高速に解が求まると思われます。
移流項の計算では、vector fieldをトレースするための計算が必要ですが、トレースの具体的方法はまだ述べられていません。

[次回は]

次回は引き続きこの論文を読んでいこうかと思っていましたが、本論文はここまでにしたいと思います。
ちなみに
2.3章では周期的境界条件を適用した際のFFTの活用について
2.4章ではField内に存在する物質の動きについて
3章では実装について
4章では実際の計算時間や機材について
5章ではまとめ
が述べられています。
3章では実装について述べられていますが、氏のGDC03の発表の方がより実践的であるため、こちらを読むことで十分だと判断しました。
次回はこちらを読み進めていきたいと思います。


Reference
Stable Fluids [Jos Stam SIGGRAPH’99]
http://dl.acm.org/citation.cfm?id=311548


コメントを残す

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

WordPress.com ロゴ

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

Twitter 画像

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

Facebook の写真

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

Google+ フォト

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

%s と連携中