粒子法のプログラムを実装しようとしたのですが、論文や本だと数式レベルでしか説明がなく、それをプログラムにして動くようにするまで苦労しました。
数式をプログラムに落とすこと自体は簡単なのですが、各定数の値が分からなかったりして、粒子が上手く動いてくれるまで大変でした。
なので、同じところでつまづく人も多いのではと思い、自分が書いた粒子法のプログラムを紹介しようと思います。下の動画を作ったときのプログラムです。SPH法です。
こちらがソースコードです。
※2009.11.4追記 こちらに整理したC++版のほか、Haskell版、OCaml版のコードがあります。
まずはシンプルにと考え、近傍粒子の探索は工夫せず、n(O^2)ですがわかりやすいアルゴリズムを採用しています。まず動くコードを書き、そこから高速化していくアプローチです。
また、計算量を減らした方が試しやすいのと、粒子法は３次元への拡張が簡単なことから、まずは２次元で実装しています。
プログラムの中身を別の記事で少しずつ説明していこうと思います。
粒子法のプログラム
粒子法のプログラム第１回（概要）粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
その他の解説エントリ
SPHによる巻き波のシミュレーション１SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと
固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation
剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies
動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
大変貴重なソース・コードをありがとうございます。 今後を楽しみにしています。
返信削除コメントありがとうございます！
返信削除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年以上も経っておりますが、返信いただけると良いのですが、、、
以前、掲載していた参考論文の著者名やタイトルを教えていただけないでしょうか？