粒子法シミュレーション記述言語によるn-Body(多体問題)シミュレーション

以前のエントリでは、粒子法シミュレーションを記述するための言語として、まず１つの質点についての自由落下を記述できるようにしてシミュレーションする例を示しました。

今回は、それを進めて、複数質点についてのn-Bodyシミュレーションを記述できるようにしたので紹介します。

n-Bodyシミュレーションの記述

n-Bodyシミュレーションを、数式で以下の様に表現した場合を考えます。



ここで、xi, vi, ai, fiは、それぞれ、ある質点iの位置、速度、加速度、力です。mは質点の質量、gは重力定数です。質点の質量はすべての質点について同じとします。

提案する言語を使うと、このシミュレーションを、以下のようにHaskellの組み込みDSLとして記述できます。

nBody =
  do define         "x" (integral (var "v" <*> dt))
     define         "v" (integral (var "a" <*> dt))
     defineWithType "a" (var "f" </> m) ExpTypeVector
     define         "f"
                    (let r = norm (var' "x" <-> var "x")
                         k = m <*> m <*> g </> r </> r
                         n = (var' "x" <-> var "x") </> r
                     in sigma (k <*> n))
     initialConditionV "x" x0
  where dt = constantS 0.01
        m  = constantS 1
        g  = constantS 9.8
        x0 = V.fromList [ (0,0,0), (2,0,0), (0,2,0) ]

記述の説明

"define"によって、ある質点iの物理量について方程式を定義しています。

"defineWithType"は、"define"と同じように物理量を定義しますが、併せて型を指定する点が異なります。"define"を使った場合には、その式の型は他の式から推論されますが、他の式から型を推論できない場合には、"defineWithType"を使って型も併せて指定します。現在使える型は、スカラーとベクトルの２種類です。

"initialConditionV"は、integralを使って定義された物理量の初期条件を与えるためのものです。最後の"V"は、その物理量の型がベクトルであることを示しています。引数として、ベクトルの配列を渡しています。速度vについては初期条件が省略されており、その場合は0ベクトルが使われます。

"integral"は、積分を表します。内部では陽解法で計算しています。

"var", "var'"は、引数で与えたシンボルで表される変数です。var'はsigmaの中でのみ使うことができます。varとvar'の違いについては、以下のsigmaの中で説明します。

"sigma"は、自身を除いた質点jについてオペランドの式を計算した値の合計です。オペランド中の var "x" がxiを、var' "x" がxjを表します。iが一定で、jをi以外の各質点について変えた場合の値を合計したものを返します。

dt、m、gといった定数は、constantS関数でリフトさせて使います。最後の"S"は、その定数がスカラーであることを示しています。

実行

上記で定義したnBodyを、以下のようにrunSimMachine関数に渡すことで計算します。

main = mapM_ (printRow . flip machineRegisterValueV "x")
     $ take 100 $ runSimMachine nBody

※コード全体についてはこちらを見てください

実行速度について

現時点では、この言語はインタプリタとして実装されているため、実行速度はあまり速くありません。

同じシミュレーションをHaskellで直接書いた場合(ソースコード)と速度を比較すると、以下のようになりました。SimLangがこの言語、NBodyがHaskellで直接書いた場合です。

$ time ./SimLang > /dev/null

real    0m2.900s
user    0m2.869s
sys     0m0.020s

$ time ./NBody > /dev/null 

real    0m0.846s
user    0m0.829s
sys     0m0.011s

Haskellで直接書いた場合に比べて、1/3〜1/4の速度になっています。

実行速度については、Temlate Haskellを使って、Haskellで直接書いた場合と同等のコードに（広義の）コンパイルすることで効率化できます。

課題

現時点での実装には、以下のような課題があります。

• 0割算などの例外処理
• デバッグのしやすさ、記述にエラーがあったときのメッセージの改善

ロードマップ

今後のロードマップとして以下のようなことを考えています。

• より強力な型推論
• Template Haskellによる（広義の）コンパイル
• テスト
• 設計の解説
• 組み込みDSLや言語抽象についての解説
• SPHの記述（分岐や近傍探索の導入）
• 剛体シミュレーションの記述（階層表現の導入）
• 固液連成シミュレーションの記述（ドメインの導入）
• ベンチマーク

関連エントリ

このエントリに関連するページです。

粒子法シミュレーションを記述するための言語を設計中
【まとめ】埋め込みDSL関係の論文(1)
【まとめ】埋め込みDSL関係の論文(2)

まとめページ

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

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

動画の一覧

【まとめ】CG、シミュレーション関係の論文
【論文】Particle-Based Fluid Simulation ...
【論文】Two-way Coupling of Fluids to ...

そのほかの論文まとめ（準備中）

--