粒子法のプログラムを実装しようとしたのですが、論文や本だと数式レベルでしか説明がなく、それをプログラムにして動くようにするまで苦労しました。
数式をプログラムに落とすこと自体は簡単なのですが、各定数の値が分からなかったりして、粒子が上手く動いてくれるまで大変でした。
なので、同じところでつまづく人も多いのではと思い、自分が書いた粒子法のプログラムを紹介しようと思います。下の動画を作ったときのプログラムです。SPH法です。
こちらがソースコードです。
※2009.11.4追記 こちらに整理したC++版のほか、Haskell版、OCaml版のコードがあります。
まずはシンプルにと考え、近傍粒子の探索は工夫せず、n(O^2)ですがわかりやすいアルゴリズムを採用しています。まず動くコードを書き、そこから高速化していくアプローチです。
また、計算量を減らした方が試しやすいのと、粒子法は３次元への拡張が簡単なことから、まずは２次元で実装しています。
プログラムの中身を別の記事で少しずつ説明していこうと思います。
粒子法のプログラム
粒子法のプログラム第１回（概要）粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
その他の解説エントリ
SPHによる巻き波のシミュレーション１SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと
固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation
剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies
動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
18 件のコメント:
大変貴重なソース・コードをありがとうございます。 今後を楽しみにしています。
コメントありがとうございます！
DEMに関して何かいいページがないか探しているところに、このブログに当たりました。
当方で使用しているのはFortranですが, C++コードを参考にさせていただきます。ありがとうございます。
コメントありがとうございます！DEMであれば、コードの大枠はそのままで力の計算の部分だけを差し替えれば良さそうですね！
できればＶＣ版も
このプログラムで使用しているパラメータ(SPH_ が頭につく定数)の意味を知りたいです。
> hiroshiさん
定数の意味ですが、以下のようになっています。
SPH_PMASS
→粒子１つの質量
SPH_INTSTIFF
→圧力係数
大きくすることで、非圧縮性を近似します
SPH_EXTSTIFF
SPH_EXTDAMP
→境界条件の係数
SPH_PDIST
→粒子の初期間隔
SPH_SIMSCALE
→シミューレションのスケール
（この係数で位置ベクトルを拡大しています）
SPH_RESTDENSITY
→定常密度
流体の密度がこの値に近づくように、粒子が動きます
SPH_VISC
→粘性係数
SPH_LIMIT
→境界条件における、速度の上限
初めまして、いつも楽しくホームページを見させていただいています。ソースコードを拝見させていただきました。壁との境界で、ペナルティー法のところが良くわからなかったのですが、詳しい説明か参考文献を教えていただけないでしょうか？
よろしくお願いします。
> Sakkyanさん
どうもありがとうございます。
ペナルティ法ですが、境界へ粒子がめり込んだ場合に、めり込みの大きさに応じた反発力を与えています。
詳細には、反発力は２の要素からなり、１つはめり込みの大きさに比例した項と、もう１つは境界面の法線方向に対する粒子速度に比例した項です。この２つの和がadjという変数に入ります。
このような、位置と速度に応じた反発力を与えるモデリングは「ばねダンパモデル」といわれます。
分からない点についての回答になっていますか？
回答ありがとうございます。よくわかりました。
定数の決め方や単位なども教えてもらっていいですか？
ソースコードを拝見させていただきました．変数DTは何を表しているのか教えていただけますでしょうか？
SPH法の具体的なコードを載せてくださっているのはここだけで、とても助かっています。
このプログラムを実装するにあたって、利用した数式が知りたいのですが、教えていただけないでしょうか。
こちらの論文(PDF)が、簡略化されて理解しやすいと思います。
kanonamaさん
非常に貴重なソースコードありがとうございます。大変参考になります。
１つお聞きしたいのですが、近傍粒子の検索について、c++コードの258行目の計算での
x+y*mx+z*mx*my
がなぜ検索の基準の値になるのですか？
よろしければ回答よろしくお願いします。
３次元の近傍マップを内部的には１次元の配列として持っており、３次元のインデックス (x, y, z) を１次元のインデックスに変換しているためです。
kanonamaさん
返信遅くなってすみません。
回答ありがとうございました。
コードを拝見させていただきました。
SPH法のサンプルコードがほしかったので、見つけたときは嬉しかったです。
10年以上も経っておりますが、返信いただけると良いのですが、、、
以前、掲載していた参考論文の著者名やタイトルを教えていただけないでしょうか？
コメントを投稿