Semi-Lagrangian法

提供: Akionux-wiki
Share/Save/Bookmark
移動: 案内検索

概要

Robert Bridson著Fluid Simulation for Computer Graphics[1]を参考にSemi-Lagrangian法をまとめる。 コンピュータグラフィックスにおけるSemi-Lagrangian法は視覚効果と計算パフォーマンス優先なので精度に関しては参考にすべき出ないところがあるはずなので、高精度化に関しては他の文献等を参考されたい。

アルゴリズム概観

  • 適切なの値を決定する。
  • 移流項(advection)を適用する。
 (1)
  • 重力項を適用する。
  • Sweeping
  • 境界条件を適用。
  • 非圧縮状態にする(projection)。#非圧縮化を参照。
 (2)

実装

1ステップ進行の例:

void advance_one_step(Grid &grid, Particles &particles, double dt)
{
	for(int i=0; i<5; ++i)
		particles.move_particles_in_grid(0.2*dt);
	particles.transfer_to_grid();
	grid.save_velocities();
	grid.add_gravity(dt);
	grid.compute_distance_to_fluid();
	grid.extend_velocity();
	grid.apply_boundary_conditions();
	grid.make_incompressible();
	grid.extend_velocity();
	grid.get_velocity_update();
	particles.update_from_grid();
}

Time stepの決定

semi-Lagrangian法ではどんなに大きなをとっても安定である。

  • CFL条件

CFL条件は3人の応用数学者R. Courant, K. Friedrichs, H. Lewyの頭文字から名前が付けられた、収束条件である。 時間依存偏微分方程式の解におけるある点には初期時間のあるドメインに依存する。例えば、等速運動の場合は速度ベクトルで経過時間分戻った点に依存するし、拡散方程式は初期状態のあらゆる点に依存する。このように正確な解を得るにはより正確なドメイン依存性を指定してやる必要がある。これがCFL条件である。


すなわちタイムステップの条件


情報が伝播する速度、言い換えれば音速をと置けば、


はCFL数と呼ばれる。

Grid

MACグリッド

MAC(Marker-And-Cell)グリッド[2] を構成する。MACグリッドはスタッガードグリッドの一種である。 1/2インデックス刻みのグリッド(これをスタッガードグリッドと呼ぶ)上において、におけるの微分を以下のように近似する。


スタッガードでないグリッド上の微分における前方か後方かのバイアスを無くす。精度はである。
MACグリッドは圧力と非圧縮性を正確に扱うために用いられる。具体的には以下のような配列になる:





圧力のグリッドのマス目の数に対して、となる。

実装

まずは自作の3次元配列の構造体を定義しておくと便利。スタッガードグリッドを構成するため、内挿(interpolation)を頻繁に使うので、3次元配列のメソッドで内挿が求められるのが望ましい。

//内挿メソッド
T bilerp(int i, int j, int k, T fx, T fy, T fz)
{ return (1-fx)*((1-fy)*((1-fz)*(*this)(i,j,k)
+fz*(*this)(i,j,k+1))+fy*((1-fz)*(*this)(i,j+1,k)+fz*(*this)(i,j+1,k+1)))
+fx*((1-fy)*((1-fz)*(*this)(i+1,j,k)+fz*(*this)(i+1,j,k+1))+fy*((1-fz)*(*this)(i+1,j+1,k)
+fz*(*this)(i+1,j+1,k+1))); }

Gridクラスのプロパティに速度場、圧力などのArray3構造体を組み込む。値の取り出しにはセルのセンター(インデックスがであるもの)とセルの端(インデックスがであるもの)の2種類を用意しておくと良い。

移流項(Advection)

移流方程式(1):


を解くことで次の時間における場を求めたい。ルーチンは:


一次元では、


これをMACグリッド上でオイラー近似する:


について解く:

しかしながら、この差分方程式は非常に不安定である。仮想粒子(Particle)を移動させ、移動後のグリッドへ値を内挿(Interpolation)する手法がStam(1999)[3]によって導入された[4]

仮想粒子による内挿

仮想粒子の移動の微分方程式は、


終点から時間を巻き戻して始点に戻る移動をオイラー近似すると、


前方オイラー近似が適切な場合もあるが、このような後方高次のRunge-Kutta法を使ったほうが圧倒的に良い結果が得られる[1]。2次のRunge-Kutta法が推奨される:



までやってきた仮想粒子から以下のような式で移流後の状態を得る:

仮想粒子はグリッド上に載らないので、内挿が必要である。例えば1次元では、の仮想粒子はに戻ってきて、に存在する場合、とすると、なので、アップデートは

境界条件

仮想粒子がドメインを越境した場合、異なるドメイン間の速度や圧力の調節が必要となる。詳しくはRobert Bridson(2008)[1]の6章を参照。

数値粘性による散逸

例えば、一次元の移流方程式:


が近似の誤差によって実際には


になり、右辺に粘性に類似した項が現れる。にすれば数値粘性項をなくすことができるが、そのためには計算コストがかさむし、粘性のある流体を扱うならば実際の粘性が大きいので数値粘性は無視できるようになる。もしくは、数値粘性があることによって実際の粘性があるように振る舞うとも言える。

非圧縮化

流体シミュレーションの中心部に当たるプロセスである。ルーチンとしては


の部分に当たる。
元の数式は式(2)であり、

 (2-a)

かつ

 (2-b)

である。

速度場の更新式

式(2-a)を変形して:


離散化:


変形:

 (3)

3次元におけるアップデートは次のようになる:

 (3-a)
 (3-b)
 (3-c)

MACグリッドのおかげでの計算が簡単かつ正確に求めることができる。
ただし、これらの式は流体を含んだセルで、かつ速度と圧力が有効なセルのみに適用される。まずはセルに流体または固体または空気のみが含まれたセルを扱うこととし、セル中に複数の相が混ざったセルの扱いは後回しにする。

ダイバージェンスを0にする条件式

更新後の速度場は式(2-b)を満たす必要がある:


境界条件

境界条件で最も簡単なのは自由境界の場合で、液体と空気の境界の場合に相当し、Dirichlet境界条件と呼ばれる。境界が固体に接している場合はもう少し難しいが、Neumann境界条件と呼ばれ、 次の境界条件も満たす必要がある:


Neumann境界条件のとき、セルが流体でセルが固体の時、


において、に等しいので


という圧力の更新になる。
projectルーチンの境界のセルでは次のようなアップデートを行う:

u(i,j,k)=usolid(i,j,k);
u(i+1,j,k)=usolid(i+1,j,k);

ダイバージェンスの離散化

は3次元では


であり、MACグリッド上では

 (4)

と表せる。MACグリッドでない一般的な連結グリッドだとダイバージェンスの計算はnullな空間で計算されるため、うまく求めることができない。
ダイバージェンスは圧力と同様にスタッガードグリッドの整数インデックスに格納され、圧力の式の右辺(rhs, right-hand side)は次のように求める:

scale = 1/dx;
loop over i,j,k where cell(i,j,k)==FLUID;
   rhs(i,j,k) = -scale * (u(i+1,j,k)-u(i,j,k) +v(i,j+1,k)-v(i,j,k) +w(i,j,k+1)-w(i,j,k));

圧力方程式

式(2-b)へ式(4)を代入:

これに式(3-a), (3-b), (3-c)を代入:

変形して、

 (5)

実は、式(5)はポアソン問題の近似式である。

流体セル(i,j,k)の隣(i,j,k+1)が空気セルならば、とすればよい。一方、(i,j,k)の隣(i,j,k+1)が固体セルならばを代入するが、処理は少々面倒になり、

  • の係数を1減らす。言い換えれば、の係数は隣合う非固体セルの個数である。
  • 右辺(rhs)にを加える変更が必要となる。

流体セルの周囲の固体セルによる変更は次のようなコードになる:

scale = 1 / dx;
loop over i,j,k where cell(i,j,k)==FLUID:
    if cell(i-1,j,k)==SOLID:
        rhs(i,j,k) -= scale * (u(i,j,k) - usolid(i,j,k));
    if cell(i+1,j,k)==SOLID:
        rhs(i,j,k) += scale * (u(i,j,k) - usolid(i,j,k));
    if cell(i,j-1,k)==SOLID:
        rhs(i,j,k) -= scale * (v(i,j,k) - vsolid(i,j,k));
    if cell(i,j+1,k)==SOLID:
        rhs(i,j,k) += scale * (v(i,j,k) - vsolid(i,j,k));
    if cell(i,j,k-1)==SOLID:
        rhs(i,j,k) -= scale * (w(i,j,k) - wsolid(i,j,k));
    if cell(i,j,k+1)==SOLID:
        rhs(i,j,k) += scale * (w(i,j,k) - wsolid(i,j,k));

不完全コレスキー分解による前処理付き共役勾配法

ここからは式(5)を連立方程式として具体的に解く方法となる。
式(5)を次のような行列表現にする:


は係数の行列、は未知の圧力の行列、は式(5)の右辺に境界条件を施した行列である。

  • Modified Incomplete Cholesky conjugate gradient(不完全コレスキー分解&共役勾配法)
  • 共役勾配法(conjugate gradient method, CG)

解に近づくまでイテレートする方法の1つ。グリッドの幅に比例して望ましい精度の解を得るために必要な計算量が増える。より計算量が少ない簡易な方法としてGauss-Seidel法などがある。Bridsonの本では不完全コレスキー分解を前処理として行った(preconditioned)共役勾配法(Preconditioned Conjugate gradient, PCG)を用いる。
前処理を行う理由はAが単位行列に近いほどイテレーションの収束は楽なので、


として前処理でを単位行列に近づけるのである。

  • 収束

収束するかどうかを知るために番目の推定値に対する残余(residual)を定義する:


が0に近づくほど誤差が小さい。
がどの程度小さくなれば十分かを知るため、物理的考察に戻ろう。であり、式(2-b)にある通り物理的にも0にしたい値であった。つまり、残余がどれだけ0に近いかという問題はどれだけダイバージェンスが0に近いかという問題と同じである。十分に小さな値をとすると、tolの次元は1/時間である。そしてはイテレーション過程の速度場をいくらか圧縮するために必要な時間の下限である。ゆえに、tolはシミュレーション時間に反比例するはずである。しかし、シミュレーション時間は場合によって数秒であったりインタラクティブな場合は無限だったりするので、実用的には程度にする[1]。より小さなtolにするほどよりダイバージェンスの誤差が小さな結果が求められるが、同時に計算に時間がかかるトレードオフな問題である。流体シミュレーションの計算時間のほとんどは圧力ソルバに使われるので、このtolの値は最も多くのチューニングが求められる部分である。
また、浮動小数点の精度が原因で収束しない場合もあるので、イテレーション回数に制限をかける必要がある。100回程度が推奨される[1]
もうひとつ、イテレーションにおいて圧力の初期推定値をいくらにするかという問題がある。良い初期値を設定できればイテレーションが早く終わらせることができる。流体ドメインがほとんど変形しない場合には前回のステップの値を初期値とするのが良いが、一般的には流体シミュレーションは流体ドメインが大きく変形するので、圧力の値が前のステップと次のステップの間で大きく変わる場合が多いので、初期値を前回のステップの値にしようと常に0にしようと、どちらも変わらない場合が多い。そのため、初期値は全て0にするのが推奨される[1]
前処理付き共役勾配法のアルゴリズムは次のようになる:

  • 初期推定値をとし、残余ベクトルをとする。(この時点でならばを返す。)
  • 補助ベクトルをとし、検索ベクトルとする。
  • 条件(r≦tolまたはイテレーション回数上限を超える)までループ:
    • 補助ベクトル

  • を返す。

PCGは最も計算時間がかかる部分なので、ベクトル処理が最適化されているBLASやLAPACKを頼るのが良いだろう[1]

コレスキー分解

行列は行列の積に分解することができる:


これはコレスキー分解と呼ばれる。もとのとなるのでまたはを解く問題になる。

不完全コレスキー分解(IC)

前処理には多くの選択肢があるが、Bridsonの本では不完全コレスキー分解(Incomplete Cholesky, IC)を用いている。
コレスキー分解を採用しない理由は、が0成分がとても少なくても、に大量に無駄な0成分が現れ、3次元空間で解を求めようとするとメモリ不足になるためである。
不完全コレスキー分解はで0の場所から0成分を作るのを避けるようになっている。
が次のように下三角行列と対角行列に分解できると仮定する:


すると、


ここでは対角行列である。
3次元では、

修正不完全コレスキー分解(MIC)


フーリエ解析をすると、ICは高周波のノイズ除去に優れ、MICは低周波のノイズ除去に優れていることがわかる。そのため実際にはICおよびMICに重み付けをして結果を求めるのがよい。重み付けは普通は0.97程度で、大きなグリッドならば1に近づけると良い[1]

この節を書こうとした人は途中で投げ出してしまいました。今後場合によっては加筆されるかもしれません。


参考

  1. 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Robert Bridson, Fluid Simulation for Computer Graphics, A K Peters, Ltd., 2008.
  2. F. Harlow and J. Welch,"Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface.", Phys. Fluids 8 2182-2189, 1965.
  3. J. Stam, "Stable Fluids", In SIGGRAPH 99 Conference Proceedings, Annual Conference Series, pp.121-128, New York: ACM, 1999.
  4. この方法が流体力学の分野において通用しているのかどうかは要調査である。