今回は、粒子法のプログラムでどのように密度と圧力を計算しているかを説明します。
粒子の密度と圧力の計算は次のコードになります。
void calculate_density_and_pressure( Particles* ps )
{
double r2, c, sum, H2;
vec dr;
H2 = H*H;
FOR_EACH_PARTICLES( ps, p )
{
sum = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = (p->r - pj->r) * SPH_SIMSCALE;
r2 = norm2( dr );
if ( H2 > r2 )
{
c = H2 - r2;
sum += c * c * c;
}
}
p->rho = sum * SPH_PMASS * Poly6Kern;
p->p = ( p->rho - SPH_RESTDENSITY ) * SPH_INTSTIFF;
p->rho = 1.0 / p->rho; // take inverse for later calculation
}
}
それぞれの粒子について、影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。そのあと、足し合わせた重みに、粒子質量についての係数とカーネル係数をかけています。密度が計算された後、圧力を密度から計算します。
上のやり方だと、計算量のオーダーは、各粒子について他の粒子それぞれとの計算になるので、粒子の数をnとするとO(n^2)となります。しかし、影響半径外の粒子との計算は不要なので、他粒子すべてと計算する必要は本来ありません。なので、近傍粒子をあらかじめマップに入れておき、それを参照するようにすれば、計算量をO(n)にすることができます。
また、粒子pと粒子pjについて対称なので、その対称性を利用すれば、計算量をさらに半分にすることが可能です。
密度と圧力は、このように計算されています。いったん素直な実装になっていますが、上の２つの点で高速化することができます。
次回は、粒子にかかる力をどのように計算しているかという説明です。
粒子法のプログラム
粒子法のプログラム第１回（概要）粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
その他の解説エントリ
SPHによる巻き波のシミュレーション１SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと
固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation
剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies
動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。
