粒子法シミュレーションの並列化 - www.iii.u-tokyo.ac.jp for III website.

Vol. 46
情報処理学会論文誌
No. 4
Apr. 2005
粒子ベースシミュレーションの並列化
原田 隆宏†
越塚 誠一††
田中 正幸††
河口 洋一郎††
現在のプロセッサは並列アーキテクチャに向かっている.そのため,その計算機性能を十分に発揮
するためにデータ並列処理を行なうアルゴリズムを開発する必要がある.本論文では並列化されたリ
アルタイム粒子ベースシミュレーション手法を提案する.粒子ベースシミュレーションを Graphics
Processing Units (GPUs) を並列計算機として用いて高速化する.並列化したことによって,リア
ルタイムシミュレーション今まで計算可能だった粒子数の十倍以上の粒子数を計算できるようになっ
た.本論文は粒子ベースシミュレーションにおいて流体と剛体を相互作用させる手法も提案する.こ
の手法の利点は計算コストが低いだけでなく,そのアルゴリズムを並列化できることにもある.流体,
剛体,そしてそれらの相互作用をリアルタイムで計算することによって,並列化された粒子ベースシ
ミュレーションの能力を示す.
Parallelization of Particle-based Simulations
Takahiro Harada,† Masayuki Tanaka ,†† Seiichi Koshizuka††
and Yoichiro Kawaguchi††
Since processors are shifting towards parallel architecture at the present time, it is important to develop data-parallel algorithms to exploit their computational power. This paper
presents a parallelized real-time particle-based simulation method. Particle-based simulation
is accelerated by using Graphics Processing Units (GPUs) as a parallel computation platform.
With parallelized particle-based simulation, the number of particles which can be simulated
in real-time is increased by more than an order of magnitude. This paper also presents a
particle-based method to interact fluids and rigid bodies. Rigid bodies are represented by a
set of particles. The benefits of this method are not only low computational cost but also parallelism of its algorithm. By simulating fluids, rigid bodies and their interactions in real-time,
we demonstrate the capability of the parallelized particle-based simulation.
1. 序
ション手法の開発はコンピュータグラフィックスにお
論
いて重要な研究課題の 1 つである.映像制作以外にお
我々の身近に存在するものの大半は物理法則に従っ
いても,ゲームなどのリアルタイムアプリケーション
ている.例えば流体の複雑な挙動は手で作成するのは
でも物理シミュレーションを用いることができ,これ
困難だが,支配方程式を解くことによって求めること
らに応用するためにはリアルタイムに計算しなければ
ができるため,コンピュータグラフィックスにおいて物
なければならない.
理シミュレーションは重要である.近年,流体をはじ
また現在のプロセッサテクノロジは周波数を向上さ
め剛体,弾性体などのシミュレーションについて様々
せるのではなく,効率性を向上させる方向に向かって
な研究が行なわれてきた.映像制作において求めら
いる.並列アーキテクチャを備えたプロセッサも多く
れる物理シミュレーションは正確さではなく,計算速
存在する.Cell Broadband Engine (Cell BE) Archi-
度や安定性に重点が置かれるため,高速なシミュレー
tecture は汎用計算用に開発されたマルチコアプロセッ
サであり,Graphics Processing Units (GPUs) はグ
ラフィックス処理に特化した並列プロセッサである.
† 東京大学
The University of Tokyo
東京都文京区本郷 7-3-1
t[email protected]
†† 東京大学
The University of Tokyo
さらに CPU もマルチコア化されてきている.これら
の計算能力を最大限に活用することができる並列計算
アルゴリズムを開発していかなければならない.
粒子ベースシミュレーションはメッシュを用いない
1
情報処理学会論文誌
2
Apr. 2005
計算手法である.粒子間で接続情報を持たないため,
分布の場合に適したデータ構造を用いて SPH を高速
近傍粒子を毎タイムステップで探索しなければならな
化した29) .自由表面流れのシミュレーションは計算コ
い.視覚的に十分な結果を得るためには,多くの粒子
ストが高く,リアルタイムアプリケーションに応用す
を用いなければならないが,多くの粒子の挙動をリア
るためにはまだ研究の余地がある.
ルタイムで計算するには計算コストが高かった.既存
剛体シミュレーションに関しても多く研究が行なわ
の GPU を並列計算機としてシミュレーションに用い
れてきた.Baraff は剛体間の接触の際に働く力を解
た研究は全て計算ノードの接続関係が計算中に変化し
析的に求める手法を開発し2) ,剛体の接触に関して研
ないという共通点を持っていた.例えば格子を用いた
究を行なった3)4)5)6) .剛体の衝突に力積ベースの手
流体計算では各計算点はその計算点を囲む点のみと相
法を用いた研究も行なわれた37)22)45) .Mirtich は衝
互作用し,布のシミュレーションでは質点同士の接続
突に関係している剛体のみを衝突時刻に戻す Time-
は計算中に変化しない24)13) .しかし粒子ベースシミュ
warp アルゴリズムを開発した37) .バネとダンパを用
レーションは,計算粒子は自由に動くことが可能であ
いるペナルティ法を衝突応答に用いた研究もある39)46) .
り,タイムステップごとに相互作用する計算粒子が変
Hasegawa らはこのペナルティ法を用いてリアルタイ
化する.そのため,計算粒子の近傍粒子を探索しなく
ムシミュレーションを行なっているが,剛体数は非常
てはならない.本研究では計算粒子の近傍粒子探索を
に少なく,凸形状のみしか扱えないという問題があ
並列化し,計算ノードの接続関係が動的に変化する粒
る25)26) .Guendelman らは形状表現にポリゴンと符
子ベースシミュレーションも並列化することが可能で
号付き距離関数を用いて非凸形状の剛体に関して衝
あり,並列計算機上で計算することができることを示
突計算を行なった20) .Kaufman らは多数の剛体間の
す.本研究では GPU を並列計算機として用いた.そ
接触を高速に解く手法を開発した.しかし剛体シミュ
の結果,スカラープロセッサ 1 個では得ることのでき
レーションではまず衝突検出を行なわなければ成らず,
ないパフォーマンスを得ることができた.本論文では
剛体シミュレーション自体を高速化するためには,高
粒子ベースシミュレーションで流体と剛体の相互作用
速な衝突検出手法が必要である28) .Dingliana らは剛
を計算する手法も提案する.この手法では剛体は粒子
体を階層構造を持つ球で表現し衝突計算を行なう手法
で表され,計算コストが低く,並列処理可能であると
を開発し 2 次元計算を行ない14) ,Suzuki らはこの手
いう利点がある.本手法を実装し,図 1 を初めとする
法を 3 次元に拡張した48) .Bell らは剛体を複数の大
様々なシーンの計算を行ない,手法の検証を行なった.
きさの異なる球で表現して衝突計算を行なった7) .こ
のように多く剛体計算の研究が行われてきたが,リア
2. 関 連 研 究
ルタイムで多数の剛体を用いたシミュレーションを行
Foster and Metaxas がコンピュータグラフィックス
の分野では初めて 3 次元ナビエストークス方程式を解
17)
いて流体計算を行ない
,Stam が Semi-Lagrangian
なった研究はない.
剛体と流体のカップリングの研究も行なわれてき
た.Yngve らは圧縮性流体を用いた爆発の計算を行な
法を CG 分野に導入した47) .Foster and Fedkiw は
い,圧縮性流体と剛体とのカップリングを行なった51) .
16)
自由表面追跡に Level Set 法を導入した.
.その後,
G´enevaux らは剛体を粒子をバネで繋いだ弾性体で表
界面追跡手法の改良15) ,粘弾性流体19) ,融解9) ,混相
し,流体のマーカー粒子と剛体粒子間の力を計算する
流36) ,2 次元と 3 次元のカップリング27) ,Octree35) ,
ことによって流体と剛体のカップリングを行なった18) .
四面体メッシュを用いた計算30) などの研究が行なわれ
Carlson らは剛体を流体として解き,剛体に変形の拘束
てきた.また流体計算手法はこれらの格子を用いた計
条件を課すことで剛体と流体のカップリングを行なっ
算手法の他に格子を使わず粒子を用いる Lagrangian
た10) .Takahashi らは流体と飛沫と泡の計算を行な
の手法がある.Premoze らは Moving Particle Semi-
い,さらに剛体と流体のカップリングも行なった49) .
Implicit (MPS) Method を
44)
,Cani and Desbrun
Takahashi らは Volume of Solid という各ボクセル
は Smoothed Particle Hydrodynamics (SPH) をコ
での剛体の占める割合を導入して剛体部分を特定し
8)
ンピュータグラフィックスに応用した .Lagrangian
た.そして剛体境界部分に圧力計算の境界条件を課し,
手法は移流によって質量損失が起こらないため,リア
計算された圧力を用いて剛体にかかる力を計算した.
ルタイムシミュレーションに向いていることを M¨
uller
M¨
uller らは SPH で流体を解き,Lagrangian メッシュ
らは示し40) ,他にも混相流42) ,相転移50) の研究も行
で弾性体の計算を行ない,弾性体メッシュ表面に仮想
なわれてきた.Kipfer and Westermann は疎な粒子
粒子群を毎ステップ生成して流体粒子とそれらの間
Vol. 46
粒子ベースシミュレーションの並列化
No. 4
3
図1
グラスと流体のリアルタイムシミュレーション.積み上げられたグラスを流体で満たし,
グラスを投げ込んだ.
Fig. 1 Real-time simulation of glasses and liquid. Glass tower is filled with liquid
and a glass is thrown onto the scene.
で相互作用を計算し,数千粒子でリアルタイムシミュ
レーションを行なった
41)
φ(x) =
∑
.Guendelman らは薄膜と
mj
j
φj
W (x − xj )
ρj
(3)
流体の相互作用を弱連成で解いた21) .Klingner らは
mj , ρj , xj は,それぞれ粒子 j の質量,密度,座標で
境界適合格子を用いて,流体の圧力のポワソン方程式
あり,W は重み関数である.
に剛体の圧力境界条件を組み込み,強連成で相互作用
の計算を行なった31) .Chentanez らは Klingner らの
手法を拡張し弾性体と流体をカップリングさせた11) .
自由表面流れ流体シミュレーションと剛体シミュレー
流体の密度は式 (3) を用いて近傍粒子の質量の重み
付け和で求められる.
∑ ρj
ρ(x) =
mj W (x − xj )
ρj
j
ションをカップリングする研究は行なわれているが,
=
∑
それらをリアルタイムで計算するという研究はほとん
ど行なわれていない.
mj W (x − xj )
(4)
j
また状態方程式より流体の圧力は
GPU は処理によっては CPU よりも高速な計算が可
p = k(ρ − ρ0 ) + p0
(5)
能であり,様々な研究が行なわれてきた43) .Amada ら
と計算され p0 , ρ0 はそれぞれ基準となる圧力と密度で
は SPH を GPU で高速化したが,全ての計算を GPU
ある.
で行なうことはできず,近傍粒子探索は CPU で行なっ
式 (2) を粒子で計算するためには圧力項と粘性項を
た1) .Kolb らの手法では粒子位置における物理量を
離散化する必要がある.対称性を考慮した粒子 i に働
求める際にまず格子上の値を求め,それらを補間して
く圧力 fipress と粘性力 fivis は近傍粒子の物理量を用
求めるため補間による数値拡散を引き起こしてしま
いて以下のように計算される.
∑ pi + pj
fipress = −
∇Wpress (rij )
mj
2ρj
う32) .
3. 流体シミュレーション
fivis = µ
3.1 支配方程式
mj
j
非圧縮流体の支配方程式は以下に示す連続の式と運
vj − vi
∇Wvis (rij )
ρj
(6)
(7)
また本研究では M¨
uller らの用いた重み関数を用い
た40) .
動量保存の式である.
Dρ
=0
Dt
(1)
DU
1
= − ∇P + ν∇2 U + g
Dt
ρ
(2)
ρ, U, P, ν, g はそれぞれ流体の密度,圧力,速度,動
粘性係数そして重力である.
3.2 Smoothed Particle Hydrodynamics
SPH では座標 x での物理量は近傍粒子の持つ物理
量の重み付け和で計算される.
∑
j
4. 剛体シミュレーション
剛体の運動は並進運動と回転運動に分けて計算する.
4.1 並 進 運 動
剛体に働く力 F が剛体の重心の運動量 P の時間微
分である.
dP
=F
dt
(8)
また剛体の速度 v は重心の運動量 P と剛体の質量 M
情報処理学会論文誌
4
Apr. 2005
図 2 粒子による形状表現.
Fig. 2 Particle representations.
を用いて以下のように求まる.
v=
1
P
M
(9)
図 3 4,096 個のチェスの駒のリアルタイムシミュレーション.
Fig. 3 Real-time simulation of 4,096 chess pieces.
剛体の速度 v を用いて剛体の重心位置 x の時間微分は
dx
=v
dt
(10)
と表される.
4.2 回 転 運 動
剛体に働く力 F は剛体の回転運動を生み,角運動
量を変化させる.角運動量 L の時間微分は
dL
=r×F
dt
(11)
と表される.角速度 w は角運動量 P と時刻 t での慣
性テンソル I(t) を用いて
−1
w = I(t)
P
(12)
図 4 2,731 個のトーラスのリアルタイムシミュレーション.
Fig. 4 Real-time simulation of 2,731 toruses.
と計算される.ここで時刻 t での慣性テンソルの逆行
列 I(t)−1 は,初期状態での慣性テンソル I(t) と時刻
る.ポリゴンモデルから粒子データの生成は Harada
t での回転行列 R(t) を用いて
and Koshizuka の手法を用いた23) .剛体粒子を表す
−1
I(t)
−1
= R(t)I(0)
−1
R(t)
(13)
粒子の径を流体粒子の径と同じにすることで流体との
と計算する.そして式 (12) で求めた角運動量 w を用
相互作用の計算も容易に行なえるようになる.この剛
いて,クオータニオンの変化量 dq を以下のように計
体を構成している粒子を剛体粒子と呼ぶ.図 2 にポリ
算する.
ゴンモデルと粒子で表現した形状の例を示す.粒子の
( )
dq = [cos
( )
θ
θ
, asin
]
2
2
(14)
ここで回転軸 a と回転角 θ は
w
a=
|w|
θ = |wdt|
解像度を高くするほど物体の形状の近似は良くなり,
解像度を低くするほど計算コストが下がるが,計算精
度も同時に下がる.つまり粒子による表現は計算精度
(15)
を変えることによって計算速度をコントロール可能な
のでコンピュータグラフィックスでは有用である.
(16)
と計算される.
5.2 剛体の衝突
剛体の衝突の検出は各剛体を構成している剛体粒子
計算されたクオータニオンの変化量と時刻 t でのク
の衝突を検出することによって行なう.そして衝突に
オータニオン q(t) を用いて時刻 t + dt でのクオータ
よって剛体に生じる力は剛体粒子に生じる力の和をと
ニオンを以下のように更新する.
ることで計算する.剛体粒子を導入することによって
q(t + dt) = dq × q(t)
(17)
5. 相 互 作 用
衝突検出を容易に行なうことができるようになる.衝
突検出が行なわれた後は,粉体の計算手法である個別
要素法を応用して計算する12) .
5.1 形 状 表 現
5.2.1 反 発 力
本研究では剛体を均一な径を持つ粒子 (ある大きさ
個別要素法の計算では様々な反発力の計算方法が研
を持った球) の集合体として近似して相互作用を計算す
究されているが,本研究では線形バネとダッシュポッ
Vol. 46
粒子ベースシミュレーションの並列化
No. 4
5
図 5 流路内にある流体と剛体.剛体の流体に対する比重はそれぞれ 0.5, 1, 2 である.
Fig. 5 A liquid and a rigid body are in a channel. Relative densities of rigid
bodies to the liquid are 0.5, 1,and 2, respectively.
トを用いて反発力を計算する.このモデルでは衝突し
まずは剛体を構成している粒子を流体として粒子にか
ている粒子間にめり込み量に比例した反発力と相対速
かる力を計算する10) .剛体の粒子 i が流体から受ける
度に比例した減衰力を働かせる.2 つの粒子 i と j が
力 fi,f luid は圧力勾配による力と粘性力なので式 (6)
あり,これらの粒子の直径を d とするとそれらの距離
と式 (7) を用いて
|rij | が d より小さいときに,以下のバネによる力 fs
とダンパによる減衰力 fd を働かせる.
fs = −k(d − |rij |)
rij
|rij |
fd = η(vj − vi )
fi,f luid = fipress + fivis
(23)
と計算される.
(18)
剛体を構成する粒子に働く力は剛体間の衝突による
力 fi,rigid と流体から受ける力 fi,f luid であるので,剛
(19)
ここで k, η はバネ定数と減衰定数であり,これらと反
発係数の関係を定式化することはできないため経験的
に値を決定した.また rij = rj − ri であり ri , rj は
体を構成する粒子にかかるこれらの力の和を計算する
ことで,重心に働く力を計算することができる.
∑
Fc =
(fi,rigid + fi,f luid )
(24)
i∈RigidBody
また粒子に働く力と剛体の重心からの座標 r′i を用い
粒子 i, j の位置ベクトルである.
5.2.2 摩 擦 力
て,その粒子に働く力が生む剛体のトルクを計算する
クーロンの法則によると動摩擦力は垂直応力に比例
ことができる.剛体全体に働く力は粒子にかかるトル
した力であり,相対運動方向に働く.よって粒子間の
動摩擦力は
vt
ft = −µ|fn |
|vt |
クの和を取り
Tc =
∑
r′i × (fi,rigid + fi,f luid ) (25)
i∈RigidBody
(20)
と計算できる.計算された力によって剛体の重心座標
と計算される.しかし式 (20) は剪断方向の相対速度
とクオータニオンを更新する.そして剛体を構成し
vt が 0 のときは計算することができなくなる.そこで
ている粒子の座標をその粒子の初期配置,剛体の重心
vt がある値よりも小さい場合は動摩擦力は剪断方向
座標とクオータニオンから計算し,粒子の座標も更新
の速度 vt に比例した力として以下のように計算する.
する.
ft = kt vt
(21)
動摩擦力は静止摩擦力が最大値を超えた場合に働く.
7. 近傍粒子探索
静止摩擦は Lee and Herrmann が行なっているよう
剛体間の衝突の計算も剛体と流体の相互作用の計算
に衝突した粒子のペアとその相対座標を保持しておく
でも,各粒子の周囲に存在する粒子との距離から力を
ことで計算することができる33) .しかし静止摩擦力を
計算するため,近傍粒子を探索しなければならない.
計算しなくても視覚的に十分な挙動を計算することが
総当たりで計算すると粒子数 n の 2 乗のオーダーの
できたので,本研究では静止摩擦力は導入しなかった.
計算コストがかかり,粒子数が増加すればするほど計
これらの反発力と摩擦力から球 i に働く衝突力
算コストは高くなってしまう.効率的に近傍粒子を探
fi,rigid は
fi,rigid = fs + fd + ft
索するために,計算領域を包括する計算格子 (グリッ
(22)
と計算される.
6. 流体と剛体の相互作用
流体と剛体の相互作用は Carlson らの手法のように
ド) を導入する38) .グリッドを構成する立方体要素が
ボクセルであり,各ボクセルにそれらが囲む計算領域
を割り当てる.そしてボクセルには割り当てられた計
算領域内に存在する粒子の番号を格納する.こうする
ことである粒子の周囲に存在する粒子はその粒子が格
情報処理学会論文誌
6
Apr. 2005
図 6 256 個のチェスの駒と液滴を水槽に落とした.
Fig. 6 256 Chess pieces and balls of liquid are dropped into a tank.
8.1 データ構造
GPU 上で物理計算を行なうため,物理量は全てテ
クスチャとしてビデオメモリ上に保持される.流体シ
ミュレーションで各粒子が持つ物理量は座標と速度と
密度であるため,それぞれテクスチャを用意する.計
算を行なうにあたり,流体の密度を計算しなければな
らないが,密度はそれぞれのタイムステップで粒子座
標から計算されるため,1 枚のテクスチャを用意する.
また剛体シミュレーションでは各剛体は座標,クオー
タニオン,並進運動量,角運動量であるため,それぞ
れにテクスチャを用意する.また慣性テンソルと初期
粒子配置は同じ形状の剛体ならば共通のデータなので
剛体ごとに用意する.流体と剛体シミュレーションで
必要なこれらの物理量は計算要素 1 つに対して 1 ピク
セルを割り当てる.つまり流体粒子 1 個と剛体 1 個に
図 7 1 個のグラスに流体を流入させた.
Fig. 7 Liquid is poured into a glass.
1 ピクセルを割り当て,その物理量をピクセルに保持
する.また本手法では近傍粒子探索を行なうためにグ
リッドを構築するため,グリッドをメモリ上に確保し
納されているボクセルの周りのボクセル内に格納され
なければならない.計算領域を囲むグリッドは 3 次元
ている粒子に限定され,計算コストは粒子数 n のオー
格子であるが,GPU は 3 次元格子に出力することは
ダーになり効率よく探索することができる.本研究で
できないため,Harris らが行なったように 2 次元格子
は一辺の長さが粒子の直径と同じ大きさのボクセルを
を複数枚並べた 1 枚のテクスチャとしてメモリ上に確
用いた.すると粒子が最密構造の配置を取ったとき,
保する24) .グリッドを構成するボクセル 1 個に 1 ピ
1 個のボクセルに入る粒子数は 4 個である.
クセルを割り当てる.
8. GPU を用いた並列計算
8.2 グリッド構築
グリッド構築はボクセル内に存在する粒子番号を,
今まで述べてきた流体と剛体シミュレーションとそ
対応するボクセルに格納する処理である.データを並
れらの相互作用は全ての計算を並列化することが可能
列処理するとそれぞれのボクセルには最終的に 1 個
である.本研究では GPU を並列計算機として用いて
の値しか保持されないため,この処理では 1 個のボク
シミュレーションを並列化させる方法を述べていく.
セルに数個粒子が存在する場合には正しくグリッドを
本章ではまず GPU 上でのデータ構造について説明し,
構築することができない.そこでこの処理を数段階に
実装の詳細を述べていく.粒子からグリッドを構築す
分けることによって,正しく処理することができるよ
る処理は流体と剛体のシミュレーションに共通の処理
うになる.これは GPU を用いた処理においても言え
であるため,まず GPU 上でのグリッド構築手法につ
ることであり,マルチパスレンダリングを行なうこと
いて述べる.それから流体と剛体シミュレーションの
によってグリッドを正しく構築することができる.こ
計算について順に述べていく.
のデータ分配処理は GPU ではそれぞれの粒子に対応
Vol. 46
粒子ベースシミュレーションの並列化
No. 4
7
図 8 流体を水槽に流入させた.パドルは 1 軸によって空間に固定されている.
Fig. 8 Liqid is poured into a tank. Axes of paddles are fixed in the space.
付けられた頂点を用意し,それらを対応する粒子のグ
4 パス目では基本的に 2 パス目で行なった処理と同じ
リッド内の座標に粒子番号を色として出力することで
処理を行なう.しかし 3 パス目では RG チャンネル,
行なうことができる.ここからマルチパスでグリッド
4 パス目では RGB チャンネルへの書き込みを防ぐた
を正しく構築する手法を述べていく.
め,カラーマスクを用いる.そして深度バッファは初
あるボクセル b に格納される小さい順に並んだ粒
子番号 i0 , i1 , i2 , i3 を考える.そしてこの粒子番号を
期化せずに使い続け,ステンシルバッファはそれぞれ
のパスごとに初期化する.
ボクセル b に割り当てられたピクセルの Red, Green,
8.3 流体シミュレーション
Blue, Alpha (RGBA) チャンネルにそれぞれ格納する
8.3.1 密度と力の計算
ことを目標として 4 パスで処理を行なう.頂点は粒子
流体の計算ではまず密度の計算を行なう.この計算
番号の小さい順にレンダリングする.用いるバーテッ
では近傍粒子を探索する必要がある.グリッドを参照
クスシェーダは頂点に対応する粒子座標に対応するグ
することで粒子 i が格納されているボクセルの周辺の
リッドテクスチャ内の位置に頂点を変位させ,フラグ
ボクセルに格納されている粒子番号を得ることができ
メントシェーダでは粒子番号を色と深度値として書き
る.これらの粒子番号から粒子座標を読み出し,式 (4)
出す.
を用いて密度を計算する.粒子に働く力を計算する際
1 パス目では i0 を R チャンネルにレンダリングす
も近傍粒子を探索する必要があるが,この計算も密度
る.i0 に対応する頂点は最も小さい深度値を持つ頂点
の計算と同様にグリッドを参照することで近傍粒子番
としてレンダリングされるため,デプスバッファを最
号を得ることができる.そしてその粒子番号から計算
大値で初期化し,小さい値がデプステストを合格する
した密度,座標,速度を読み出し,圧力項と粘性項の
ように設定することで i0 を書き入れることができる.
計算を式 (6) と式 (7) を用いて行なう.また剛体粒子
2 パス目では i1 を G チャンネルにレンダリングする
に関しては衝突による力を計算する.
が,このときにすでに値を書き入れた R チャンネルに
8.3.2 速度と座標の更新
上書きしないためにカラーマスクを行ない R チャン
計算された粒子に働く力を用いて速度と座標をオイ
ネルには値を書き入れないようにする.1 パス目で用
ラー陽解法で更新する.これらの処理は各粒子の持つ
いた深度バッファを用いるがデプステストは大きな値
物理量のみを用いて計算することができるので,それ
が合格するように設定することで i0 は書き込まれな
ぞれのテクスチャの各粒子に割り当てられたピクセル
くなる.しかしデプステストだけだと最も深度値が大
の値を参照するだけで行なうことができる.
きい i3 が最終的に書かれてしまう.i2 , i3 の書き込み
8.4 剛体シミュレーション
を防ぐため,ステンシルテストを用いる.ステンシル
8.4.1 並進運動量と角運動量の更新
ファンクションでは書き込まれるごとにステンシル値
剛体の並進運動量の更新は式 (8) の計算を行なう必
を増加させるように設定し,ステンシルテストは 1 よ
要があり,剛体を構成する粒子に働く力の総和を取る.
り大きい値を不合格にするように設定する.こうする
また角運動量の更新は式 (11) の計算を行なう必要が
ことで各ピクセルには 1 回値がレンダリングされたあ
あり,この計算では粒子に働く力以外に,剛体の重心
とは書き込まれなくなる.つまり i1 が書き込まれた
に対する粒子の相対位置ベクトル r′ が必要なので,剛
後は値を書き込むことができなくなるの G チャンネ
体のクオータニオンと重心座標と初期配置を用いて計
ルに選択的に i1 を書き込むことができる.3 パス目と
算する.
情報処理学会論文誌
8
8.4.2 座標とクオータニオンの更新
計算した速度を用いて式 (10) を計算し座標を更新
する.また計算された角速度を用いて式 (17) を計算
しクオータニオンを更新する.これらの計算は各剛体
Apr. 2005
表 1 剛体粒子数とフレームレート.
Table 1 Number of rigid particles and frame rates.
Figure 3
Figure 4
Number of particles
16,384
16,384
FPS
133.3
193.7
に割り当てられた物理量のみを用いた計算なので,そ
れぞれの剛体に割り当てられたピクセルの値のみの参
照で計算することができる.
9. 結
表 2 チェスの駒の数とフレームレート.
Table 2 Number of chess pieces and frame rates.
Number of particles
256
1,024
4,096
16,384
果
本手法を Core 2 X6800 CPU,2.0GB のメモリ,
GeForce 8800GTX を搭載した PC 上で実装した.プ
FPS
636.9
355.9
133.3
40.7
ログラムは C++, OpenGL, C for Graphics を用い
て開発した.本章で取り上げる計算結果は全てリア
個のグラスを積み重ねたシーンに 2 カ所から流体を流
ルタイムで計算した結果であり,それらの計算結果を
入させたシミュレーション結果である.流体がグラス
オフラインでレイトレーシングしたものである.流体
を上から満たしていった.そして流体流入後,1 個の
は計算粒子それぞれに濃度分布を与え,陰関数曲面を
グラスを投げ入れ,積み重なったグラスを崩れ落とし
34)
Marching Cubes を用いて抽出した
.
た.剛体の相互作用では静止摩擦力の計算をせず,動
流体の計算に比べ剛体の計算は 2 つの理由から計算
摩擦力の計算しか行なっていないが,このような複数
コストが低い.1 つ目の理由として,流体計算では近
の剛体が積み重なったシーンでも剛体をほぼ静止させ
傍粒子探索を密度の計算と圧力項の計算において 2 回
ることができた.次の計算例は密度のテストである.
行わなければならないが,剛体計算では衝突計算で 1
図 5 では流路が流体で満たされ,それぞれ異なった密
回行うだけで済む.2 つ目の理由は衝突計算を行なう
度を持つ剛体を落とした.流体は圧力の均衡を保つた
とき,粒子 i が入っているボクセルと,そのボクセル
め,高さを変えていることがわかる.図 6 では水槽を
に隣接しているボクセルに格納されている粒子番号を
少量の水で満たし,256 個のチェスの駒を落とした.
取り出さなければならない.しかし流体計算を行なう
その後液滴を水槽に落とした.チェスの駒の密度とし
ときには影響半径が粒子径よりも大きいときには,そ
て 2 種類の値を用いたため,重い駒は沈み,軽い駒は
れ以上のボクセル内に近傍粒子が存在している可能性
浮いている.図 7 には 1 個のグラスに水を注ぎ,グ
があるため,より多くのボクセルから粒子番号を取り
ラスが水によって倒されるシミュレーションの結果を
出さなければならない.
示す.図 8 では 1 軸で空中に固定された 5 個のギア
まず剛体計算のみの結果を示す.図 3 に 4,096 個の
が配置されたシーンにおいて 2 カ所から流体を流入
チェスの駒,図 4 に 2,731 個のトーラスを落としたシ
させた.流体によってギアが回転し,その抵抗により
ミュレーション結果を示す.フレームレートは表 1 に
流れが乱れている様子が計算されている.図 9 に示
示すように 133.3 と 193.7 である.このフレームレー
すシミュレーションでは水をインタラクティブに操作
トは 1 タイムステップ計算を行ない,レンダリングも
し波を作り,その後,水面に船を浮かべた.この波に
行なったものであるため,計算自体はより高速に行な
よって船が倒される様子が示されている.これらの計
うことができている.これらの計算では 16,384 個の
算に用いた粒子数と 1 タイムステップの計算時間を表
剛体粒子を用いている.両方の計算で用いた剛体粒子
4 に示す.約 6 万粒子を用いた流体と剛体の計算は約
数が同じだがフレームレートが異なっているのはレン
17fps で計算できている.
ダリングしなければならない剛体数が多いことに起因
している.剛体数を変えてフレームレートを測定した
10. 結
論
結果を表 2 に示す.また剛体計算を CPU で実装して
本研究では粒子を用いたシミュレーションを並列化
計算速度を計測し,本手法と比較したものを表 4 に示
しリアルタイムで計算する手法を提案した.GPU 上
す.この結果はレンダリングを行なわず計算のみを行
で粒子を用いたシミュレーションを並列化し,流体,
なったものである.剛体数 16,384 において CPU で
剛体そしてそれらの相互作用をリアルタイムで計算し
の計算より 24.5 倍高速に計算できている.
た.本手法を用いることによって流体シミュレーショ
次に流体と剛体の連成計算結果を示す.図 1 では 10
ンに関しては,従来リアルタイムシミュレーションで
Vol.46
粒子ベースシミュレーションの並列化
No.4
9
図 9 波と船.
Fig. 9 A wave and ships.
表 3 CPU での計算時間 (time (a)) と GPU での計算時間
(time (b)) (ミリ秒).
Table 3 Simulation time on CPU (time (a)) and on GPU
(time (b)) in milliseconds.
Number of chess pieces
1,024
4,096
16,384
Time(a)
12.5
51.8
195.1
Time(b)
1.41
2.81
7.96
表 4 流体粒子数とフレームレート.
Table 4 Number of fluid particles and frame rates.
Figure
Figure
Figure
Figure
Figure
Figure
1
5
6
7
8
9
Number of particles
49,153
8,193
64,512
12,289
64,512
245,760
FPS
17.1
91.3
16.2
75.2
15.2
3.85
用いられていた粒子数の数十倍の計算をリアルタイム
で行なうことができた.また剛体計算では多数の剛体
があるシーンをリアルタイムで計算した.またそれら
の相互作用もリアルタイムで計算した.本研究ではリ
アルタイムシミュレーション結果のみ示したが,本手
法を用いてオフラインシミュレーションの高速化も可
能である.高解像度のシミュレーションを用いた映像
制作において計算時間は重要なため,本手法は映像制
作においても有用である.
本研究では SPH を用いた流体と剛体のシミュレー
ションを行なったが,本研究で示した手法は他の粒子
法にも適用することが可能であるため,粒子を用いた
弾性体やシェルの計算にも応用することができる.ま
た本研究では流体のレンダリングはリアルタイムシ
ミュレーション結果から表面抽出をオフラインで行っ
たが,表面抽出をリアルタイムで行なう手法を開発し,
シミュレーションとレンダリングを全てリアルタイム
で行うことのできる手法も研究していく.最後に本研
究では並列計算のプラットフォームとして GPU を用
いたが,本手法は他のプロセッサ上でも同様に並列計
算することが可能であると考える.そこで今後マルチ
コアの CPU や PPU などを用いた並列化されたリア
ルタイムシミュレーションを研究していく.
参
考
文
献
1) T.Amada, M.Imura, Y.Yasumoto, Y.Yamabe,
and K.Chihara. Particle-based fluid simulation
on gpu. In 2004 ACM Workshop on GeneralPurpose Computing on Graphics Processors,
2004.
2) D. Baraff. Analytical methods for dynamic
simulation of non-penetrating rigid bodies.
Computer Graphics, 23(3):223–232, 1989.
3) D.Baraff. Curved surfaces and coherence for
non-penetrating rigid body simulation. Computer Graphics, 24(4):19–28, 1990.
4) D. Baraff.
Coping with friction for nonpenetrating rigid body simulation. Computer
Graphics, 25(4):31–40, 1991.
5) D. Baraff. Fast contact force computation
for nonpenetrating rigid bodies. In Computer Graphics Proc., Annual Conference Series, pages 23–34, 1994.
6) D.Baraff. Interactive simulation of solid rigid
bodies. IEEE Computer Graphics and Applications, 15:63–75, 1995.
7) N.Bell, Y.Yu, and P.J. Mucha. Particle-based
simulation of granular materials. In Proc. of the
2005 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pages 77–86,
2005.
8) M.P. Cani and M.Desbrun. Animation of deformable models using implicit surfaces. IEEE
Transactions on Visualization and Computer
Graphics, 3(1):39–50, 1997.
9) M. Carlson, P. Mucha, B.V. Horn III, and
G.Turk. Melting and flowing. In Proc. of ACM
SIGGRAPH Symposium on Computer Animation, pages 167–174, 2002.
10) M. Carlson, P.J. Mucha, and G. Turk. Rigid
fluid: Animating the interplay between rigid
bodies and fluid. ACM Transactions on Graphics, 23(3):377–384, 2004.
11) N.Chentanez, T.G. Goktekin, B.E. Feldman,
and J.F. O’Brien.
Simultaneous coupling
of fluids and deformable bodies. In Proc.
10
情報処理学会論文誌
of ACM SIGGRAPH/Eurographics Symposium
on Computer Animation, pages 83–89, 2006.
12) P.A. Cundall and O.D.L. Strack. A discrete numerical model for granular assemblies.
Geotechnique, 29:47–65, 1979.
13) Cyril.Z. Cloth simulation on the gpu. In ACM
SIGGRAPH Sketches, No.39, 2005.
14) J. Dingliana and C. O’Sullivan.
Graceful
degradation of collision handling in plysically
based animation. Computer Graphics Forum,
19(3):239–247, 2000.
15) D.Enright, S.Marschner, and R.Fedkiw. Animation and rendering of complex water surfaces. ACM Transactions on Graphics, 21:721–
728, 2002.
16) N.Foster and R.Fedkiw. Practical animation
of liquids. In Proc. of SIGGRAPH, pages 15–
22, 2001.
17) N.Foster and D.Metaxas. Realistic animation
of liquids. Graphical Models and Image Processing, 58(5):471–483, 1996.
18) O. G´enevaux, A. Habibi, and J.M. Dischler.
Simulating fluid-solid interaction. In Graphics
Interface, pages 31–38, 2003.
19) T.Goktekin, A.W. Bargteil, and J.F. O’Brien.
A method for animating viscoelastic fluids.
ACM Transactions on Graphics, 23:464–467,
2004.
20) E. Guendelman, R. Bridson, and R. Fedkiw.
Nonconvex rigid bodies with stacking. ACM
Transactions on Graphics, 22:871–878, 2003.
21) E. Guendelman, A. Selle, F. Losasso, and
R.Fedkiw. Coupling water and smoke to thin
deformable and rigid shells. ACM Transactions
on Graphics, 24:910–914, 2005.
22) J.K. Hahn. Realistic animation of rigid bodies. Computer Graphics, 22(4):299–308, 1988.
23) T. Harada and S. Koshizuka. Fast solid voxelization using graphics hardware. Transactions of Japan Society for Computational Engineering and Science, 2006.
24) M.J. Harris, W.V. Baxter, T. Scheuermann,
and A. Lastra. Simulation of cloud dynamics
on graphics hardware. In Proc. of the SIGGRAPH / Eurographics Workshop on Graphics Hardware, pages 92–101, 2003.
25) S. Hasegawa, N. Fujii, Y. Koike, and M. Sato.
Real-time rigid body simulation based on volumetric penalty method. In Proc. of the 11th
Symposium on Haptic Interfaces for VirtualEnvironment and Teleoperator Systems, pages
326–332, 2003.
26) S.Hasegawa and M.Sato. Real-time rigid body
simulation for haptic interactions based on con-
Apr. 2005
tact volume of polygonal objects. Computer
Graphics Forum, 24(3):946–956, 2004.
27) G. Irving, E. Guendelman, F. Losasso, and
R.Fedkiw. Efficient simulation of large bodies
of water by coupling two and three dimensional
techniques. ACM Transactions on Graphics,
25:812–819, 2006.
28) D.M. Kaufman, T. Edmunds, and D.K. Pai.
Fast frictional dynamics for rigid bodies.
ACM Transactions on Graphics, 24(3):946–
956, 2005.
29) P. Kipfer and R. Westermann. Realistic and
interactive simulation of rivers. In Proc. of the
2006 Conference on Graphics Interface, volume
137, pages 41–48, 2006.
30) B.M. Klinger, B.E. Feldman, N. Chentanez,
and J.F. O’Brien. Fluid animation with dynamic meshes. ACM Transactions on Graphics, 25:820–825, 2006.
31) B.M. Klingner, B.E. Feldman, N.Chentanez,
and J. O’Brien. Fluid amination with dynamic meshes. ACM Transactions on Graphics, 25(3):820–825, 2006.
32) A.Kolb and N.Cuntz. Dynamic particle coupling for gpu-based fluid simulation. In Proc.
of 18th Symposium on Simulation Technique,
pages 722–727, 2005.
33) J.Lee and H.J. Herrmann. Angle of repose and
angle of marginal stability: molecular dynamics of granular particles. Journal of Physics A,
26:136–146, 1993.
34) W.E. Lorensen and H.E. Cline. Marthing
cubes: A high resolution 3d surface construction algorithm. In Proc. of the 14th Annual
Conference on Computer Graphics and Interactive Techniques, pages 163–169, 1987.
35) F. Losasso, F. Gibou, and R. Fedkiw. Simulating water and smoke with and octree data
structure. ACM Transactions on Graphics,
23:457–462, 2004.
36) F.Losasso, T.Shinar, A.Selle, and R.Fedkiw.
Multiple interacting liquids. ACM Transactions on Graphics, 25:812–819, 2006.
37) B.Mirtich. Timewarp rigid body simulation.
In Proc. of SIGGRAPH 2000, pages 193–200,
2000.
38) B.K. Mishra. A review of computer simulation of tumbling mills by the discrete element
method: Parti-contact mechanics. International Journal of Mineral Processing, 71(1):73–
93, 2003.
39) M.Moore and J.Wilhelms. Collision detection
and response for computer animation. Computer Graphics, 22(4):289–298, 1988.
Vol. 46
No. 4
粒子ベースシミュレーションの並列化
40) M. M¨
uller, D. Charypar, and M. Gross.
Particle-based fluid simulation for interactive
applications. In Proc. of SIGGRAPH Symposium on Computer Animation, pages 154–159,
2003.
41) M. M¨
uller, S. Schirm, M. Teschner, B.
Heidelberger, and M. Gross. Interaction of
fluids with deformable solids.
Journal of
Computer Animation and Virtual Worlds,
15(3):159–171, 2004.
42) M. M¨
uller, B. Solenthaler, R. Keiser, and
M.Gross. Particle-based fluid-fluid interaction.
In Proc. of Siggraph Symposium on Computer
Animation, pages 237–244, 2005.
43) John D. Owens, David Luebke, Naga Govindaraju, Mark Harris, Jens Kr¨
uger, Aaron E.
Lefohn, and Timothy J. Purcell. A survey of
general-purpose computation on graphics hardware. Eurographics 2005, State of the Art Reports, pages 21–51, 2005.
44) S. Premoze, T. Tasdizen, J. Bigler, A. Lefohn,
and R.T. Whitaker. Particle-based simulation of fluids. Computer Graphics Forum,
22(3):401–410, 2003.
45) H. Schmidl and V.J. Milenkovic. A fast impulsive contact suite for rigid body simulation.
IEEE Transactions on Visualization and Computer Graphics, 10(2):189–197, 2004.
46) J.M. Snyder, A.R. Woodbury, K. Fleischer,
B.Currin, and A.H. Barr. Interval methods for
multi-point collisions between time-dependent
curved surfaces. In Proc. of the 20th Annual
Conference on Computer Graphics and Interactive Techniques, pages 321–334, 1993.
47) J.Stam. Stable fluids. In Proc. of ACM SIGGRAPH 99, pages 121–128, 1999.
48) K.Suzuki, J.Kubota, and H.Ohtsubo. Rigid
body simulation using volume based collision
detection. Structural Engineering / Earthquake
Engineering, 22:185–192, 2003.
49) T. Takahashi, H. Fujii, A. Kunimatsu, K.
Hiwada, T. Saito, K. Tanaka, and H. Ueki.
Realistic animation of fluid with splash and
foam. Computer Graphics Forum, 22(3):391–
400, 2003.
50) M. Wicke, P. Hatt, M. Pauly, M. M¨
uller, and
M.Gross. Versatile virtual materials using implicit connectivity. In Proc. of Point-Based
Graphics, pages 73–82, 2006.
51) G.D. Yngve, O’Brien J.F., and J.K. Hodgins.
Animating explosions. In Proc. of ACM SIGGRAPH 2000, pages 29–36, 2000.
11
(平成 16 年 11 月 28 日受付)
(平成 17 年 2 月 4 日採録)
原田 隆宏(正会員)
昭和 56 年生.平成 18 年東京大学
大学院工学系研究科システム量子工
学専攻修士課程修了.同年東京大学
大学院情報学環学際情報学府助手.
平成 19 年 4 月助教.計算力学とコ
ンピュータグラフィックスの研究に従事.日本機械学
会,画像電子学会,ACM 各会員.
田中 正幸
昭和 57 年生.平成 19 年東京大
学大学院工学系研究科システム量子
工学専攻修士課程修了.同年株式会
社東芝生産技術センター.粒子ベー
スのマルチフィジックスシミュレー
ションの研究に従事.日本機械学会会員.
越塚 誠一
昭和 37 年生.昭和 61 年東京大学
大学院工学系研究科原子力工学専攻
修士課程修了.同年東京大学工学部
助手.平成 3 年博士(工学).同年
講師.平成 5 年助教授.平成 16 年
教授.連続体の力学シミュレーションの研究に従事.
特に粒子法の開発をおこなう.平成 17 年に丸善より
「粒子法」を出版.
河口 洋一郎
昭和 27 年生.昭和 53 年東京教育
大学大学院修了(現,筑波大学大学
院).平成 4 年筑波大学芸術学系助
教授.平成 10 年東京大学大学院工
学系研究科・人工物工学センター教
授.平成 12 年東京大学大学院情報学環教授.昭和 52
年から CG に着手し,世界的 CG アーティストとし
て活躍中.第 100 回ベネチアビエンナーレ日本代表芸
術家に選ばれる.