2009年2月23日月曜日

粒子法のプログラム第1回(概要)


粒子法のプログラムを実装しようとしたのですが、論文や本だと数式レベルでしか説明がなく、それをプログラムにして動くようにするまで苦労しました。

数式をプログラムに落とすこと自体は簡単なのですが、各定数の値が分からなかったりして、粒子が上手く動いてくれるまで大変でした。

なので、同じところでつまづく人も多いのではと思い、自分が書いた粒子法のプログラムを紹介しようと思います。下の動画を作ったときのプログラムです。SPH法です。



こちらがソースコードです。

※2009.11.4追記 こちらに整理したC++版のほか、Haskell版、OCaml版のコードがあります。

まずはシンプルにと考え、近傍粒子の探索は工夫せず、n(O^2)ですがわかりやすいアルゴリズムを採用しています。まず動くコードを書き、そこから高速化していくアプローチです。

また、計算量を減らした方が試しやすいのと、粒子法は3次元への拡張が簡単なことから、まずは2次元で実装しています。

プログラムの中身を別の記事で少しずつ説明していこうと思います。

粒子法のプログラム
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)

その他の解説エントリ
SPHによる巻き波のシミュレーション1
SPHによる巻き波のシミュレーション2
SPHによる巻き波のシミュレーション3
このあとやりたいこと

固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション(プレビュー)
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧


--

17 件のコメント:

匿名 さんのコメント...

大変貴重なソース・コードをありがとうございます。 今後を楽しみにしています。

kamonama さんのコメント...

コメントありがとうございます!

Masahiko さんのコメント...

DEMに関して何かいいページがないか探しているところに、このブログに当たりました。
当方で使用しているのはFortranですが, C++コードを参考にさせていただきます。ありがとうございます。

kamonama さんのコメント...

コメントありがとうございます!DEMであれば、コードの大枠はそのままで力の計算の部分だけを差し替えれば良さそうですね!

a さんのコメント...

できればVC版も

hiroshi さんのコメント...

このプログラムで使用しているパラメータ(SPH_ が頭につく定数)の意味を知りたいです。

kamonama さんのコメント...

> hiroshiさん

定数の意味ですが、以下のようになっています。

SPH_PMASS
→粒子1つの質量

SPH_INTSTIFF
→圧力係数
 大きくすることで、非圧縮性を近似します

SPH_EXTSTIFF
SPH_EXTDAMP
→境界条件の係数

SPH_PDIST
→粒子の初期間隔

SPH_SIMSCALE
→シミューレションのスケール
 (この係数で位置ベクトルを拡大しています)

SPH_RESTDENSITY
→定常密度
 流体の密度がこの値に近づくように、粒子が動きます

SPH_VISC
→粘性係数

SPH_LIMIT
→境界条件における、速度の上限

Sakkyan さんのコメント...

初めまして、いつも楽しくホームページを見させていただいています。ソースコードを拝見させていただきました。壁との境界で、ペナルティー法のところが良くわからなかったのですが、詳しい説明か参考文献を教えていただけないでしょうか?
よろしくお願いします。

kamonama さんのコメント...

> Sakkyanさん
どうもありがとうございます。

ペナルティ法ですが、境界へ粒子がめり込んだ場合に、めり込みの大きさに応じた反発力を与えています。

詳細には、反発力は2の要素からなり、1つはめり込みの大きさに比例した項と、もう1つは境界面の法線方向に対する粒子速度に比例した項です。この2つの和がadjという変数に入ります。

このような、位置と速度に応じた反発力を与えるモデリングは「ばねダンパモデル」といわれます。

分からない点についての回答になっていますか?

大介 さんのコメント...
このコメントは投稿者によって削除されました。
Sakkyan さんのコメント...

回答ありがとうございます。よくわかりました。
定数の決め方や単位なども教えてもらっていいですか?

chatnoir さんのコメント...

ソースコードを拝見させていただきました.変数DTは何を表しているのか教えていただけますでしょうか?

匿名 さんのコメント...

SPH法の具体的なコードを載せてくださっているのはここだけで、とても助かっています。
このプログラムを実装するにあたって、利用した数式が知りたいのですが、教えていただけないでしょうか。

kamonama さんのコメント...

こちらの論文(PDF)が、簡略化されて理解しやすいと思います。



匿名 さんのコメント...

kanonamaさん
非常に貴重なソースコードありがとうございます。大変参考になります。
1つお聞きしたいのですが、近傍粒子の検索について、c++コードの258行目の計算での
x+y*mx+z*mx*my
がなぜ検索の基準の値になるのですか?
よろしければ回答よろしくお願いします。

kamonama さんのコメント...

3次元の近傍マップを内部的には1次元の配列として持っており、3次元のインデックス (x, y, z) を1次元のインデックスに変換しているためです。

匿名 さんのコメント...

kanonamaさん
返信遅くなってすみません。
回答ありがとうございました。