WebGPU で実装したリアルタイム 3D 流体シミュレーションの紹介
ブラウザ上で動作するリアルタイム 3D 流体シミュレーションを,WebGPU を使って実装しました.
Demo(WebGPU に対応したブラウザが必要です.)
Repository
Demo Video
本記事では,主に以下の 3 つについて説明します.
- シミュレーションの実装に WebGPU を使った理由
- 流体の運動をシミュレーションするのに用いた SPH 法
- 特に,GPU 上での近傍探索の高速化について詳しく説明します.
- 水面のレンダリングに用いた Screen-Space Rendering
WebGPU を実装に使った理由
このシミュレーションでは,水面のレンダリングには WebGPU の頂点シェーダー・フラグメントシェーダーを,SPH 法の計算にはコンピュートシェーダーを用いています.今回,実装言語として WebGPU を用いたのは,なんといってもコンピュートシェーダーを使ってみたかったからです.
ブラウザ上で GPGPU を使った計算をする場合,基本的には,WebGL を使うか WebGPU を使うかの二択です.しかし,筆者の WebGL の経験は,以前書いた流体シミュレーションで円を大量に描画した程度であり,ほぼ皆無に等しいです.この状態から,いきなり WebGL のシェーダーを駆使して GPGPU で計算するコードを書くのは少々ハードルが高いように感じられました.そこで,新進気鋭の API である WebGPU のコンピュートシェーダーに思い切って入門し,シミュレーションのロジックを書いてみることにしました.
SPH 法について
このシミュレーションは,Smoothed-particle hydrodynamics(SPH, SPH 法)に基づいて流体の運動を再現しています.SPH 法とは,流体を有限個の粒子の集合体とみなして,それらの粒子の動きを通じて流体の物理的挙動をシミュレーションする方法です.シミュレーションの "Particle" のチェックボックスを押すと,以下のように描画が粒子モードに切り替わり,流体が粒子の集合体として近似されている様子がよくわかると思います.
粒子どうしの相互作用を定めるカーネル関数や,各粒子の密度・圧力の計算式は,筆者が以前書いた ブラウザ上でヌルヌル動作する流体シミュレーションを Rust + wasm-bindgen-rayon で実装する という記事の中で紹介してあるものと同じです(カーネル関数は 3 次元版のものを用いていますが).もし SPH 法自体の詳細を知りたい場合には,この記事を参照してください.
SPH 法のパフォーマンスのボトルネックになるのは,各粒子について,その近くにある粒子を探索する近傍探索のパートです.そのため,この近傍探索を高速化することが,大量の粒子をリアルタイムにシミュレーションする鍵となります.今回は,GPU 上で近傍探索を行ったため,CPU 上での近傍探索の高速化とは異なる方法を実装する必要がありました.以下では,はじめに CPU 上での近傍探索の高速化について説明し,その後に GPU 上での近傍探索をどう実装したのかを説明します.
なお,以下では簡単のために 2 次元平面上での近傍探索について説明しますが,実際の実装では 3 次元空間上での近傍探索を行っています.適宜読み替えてください.
CPU 上での近傍探索
SPH 法の近傍探索では,各粒子について,その粒子から一定の距離
すると,各粒子について,距離
以前の記事では,各セルについて,そのセルに属する粒子のインデックスをVec<u32>
として保持しました.格子全体は,Vec<Vec<u32>>
のような 2 次元ベクタとして保持されることになります.
GPU 上での近傍探索
上記の近傍探索のロジックは,GPU にそのまま移植することはできません.なぜなら,上記の方法は各セルのベクタのメモリを動的に確保しなければいけないのに対して,GPU 上では基本的に動的メモリ確保ができないからです.ゆえに,GPU 上で近傍探索を高速化するためには,固定長のメモリで済ませる方法を考える必要があります.
固定長のメモリで GPU 上の近傍探索を高速化するための有効な方法が,FAST FIXED-RADIUS NEAREST NEIGHBORS: INTERACTIVE MILLION-PARTICLE FLUIDS というスライドで説明されている,カウントソートに基づく手法です.この手法は,以下のステップからなります.なお,粒子の個数とセルの個数は固定されているとします.
- 各セルについて,そのセルに属する粒子の個数を記録するためのカウンタを 0 に初期化する
- i 番目のセルのカウンタを
count[i]
で表すことにします.
- i 番目のセルのカウンタを
- 各粒子について,その粒子が属するセルのカウンタをインクリメントする.
- 並列にインクリメントするので,データ競合が起こらないように
atomicAdd
命令を使う必要があります.具体的には,粒子のセルが属するセルの ID がcellID
のとき,atomicAdd(&count[cellID], 1)
のようにします.
- 並列にインクリメントするので,データ競合が起こらないように
-
count
の累積和をとる.- 累積和は前から順番に計算していくので,一見 GPU で並列化するのは難しそうですが,うまいこと並列計算するための手法がいろいろと研究されています.
- この累積和を
prefixSum
と表すことにします.
- 前ステップで計算した累積和
prefixSum
を用いて,セル番号の小さい順に粒子を並べ替える.- いわゆるカウントソートです.
このように粒子を並べ替えると,同じセルに属する粒子がメモリ上で連続するようになります.具体的には,cellID
番目のセルに属する粒子が,[prefixSum[cellID], prefixSum[cellID + 1])
番目に存在することになります.よって,各粒子の近傍探索の際には,周辺の 9 個のセルについて,[prefixSum[cellID], prefixSum[cellID + 1])
の範囲にある粒子を見ればよいです.
上記のアルゴリズムで必要になるのは,count
配列,prefixSum
配列,そして粒子の並び替え先の配列 3 つのみです.前の 2 つの配列のサイズはセルの個数に比例し,最後の配列のサイズは粒子の個数に比例します.セルの個数と粒子の個数はともに固定されていると仮定したので,必要になるメモリも固定長です.すなわち,これで格子分割による近傍探索が固定長のメモリでできたことになります.これが今回のシミュレーションで実装した,GPU 上の近傍探索の全容です.
格子分割による高速化の効果
百聞は一見に如かずということで,ここまで説明してきた格子分割による近傍探索が,近傍粒子を全探索するのに比べてどのくらい高速なのかを実際に見てみます.
格子分割で探索空間を限定した場合
近傍粒子を全探索した場合
違いは一目瞭然で、リアルタイムなシミュレーションには、格子分割で探索空間を限定することが事実上必須になることが分かると思います。
水面のレンダリング
SPH 法で計算された粒子の位置情報から,水面をレンダリングする方法として今回用いたのは,Screen-Space Fluid Rendering という手法です.この手法は,リアルタイムに美しい水面をレンダリングするための強力な手法として知られています.
Screen-Space Fluid Rendering は,以下のステップからなります.
- 粒子を球体とみなして深度をレンダリングする
- フィルタを用いて深度をぼかす
- 今回は Bilateral filter を使いました.
- ぼかした深度から screen-space で法線を計算する
- 上の 3 ステップとは別に,各ピクセルの厚み情報を計算しておく
- 加算合成を用いて,各粒子の厚みをピクセルごとに足し合わせていくことで実現できます.
- 法線と厚みの情報を用いて,水面をレンダリングする
ここでは,ポイントとなる法線の計算と,法線・厚みの情報を用いてどう水面をレンダリングするかについて解説します.
法線の計算
このステップでは,フィルタでぼかした深度情報を入力として,各ピクセルに対応する法線ベクトルを計算します.その計算方法は,以下のようなシンプルな方法です.いま注目しているピクセルのカメラ座標系での座標をviewPos
とします(これはテクスチャ座標と深度情報から計算できます).
- いま見ているピクセルの 1 つ右のピクセルについて,カメラ座標系での座標を深度から計算し,
viewPos
との差分をとる.この差分をddx
とする. - いま見ているピクセルの 1 つ上のピクセルについて,カメラ座標系での座標を深度から計算し,
viewPos
との差分をとる.この差分をddy
とする. -
ddx
とddy
の外積を計算して,それを法線とする.
ただ,上記の方法をそのまま実装すると,深度が急激に変化するエッジ付近で妙なノイズが現れることがあります.これを防ぐ方法としてこのスライドで紹介されているテクニックが,右(あるいは上)のピクセルだけではなく,左(あるいは下)のピクセルについても差分を計算しておき,奥行きの差分の絶対値が小さいほうをddx
(あるいはddy
)として採用する方法です.このテクニックにより,エッジ付近でノイズが現れるのを防ぐことができます.
法線・厚みの情報を用いて水面をレンダリングする
法線と厚みの情報が求められたら,いよいよ水面のレンダリングができます.今回は,以下の 3 つの物理法則を考慮してレンダリングしました.
- フレネルの式
- 鏡面反射
- Beer's law
それぞれの物理法則について,簡単に説明します.
光が異なる物質間の界面に入射すると,一部は反射し,一部は屈折します.この現象をフレネル反射といいます.また,フレネル反射において,反射する光の割合と入射角との関係を具体的に定めた式をフレネルの式といいます.コンピュータグラフィックスの世界では,偏光を考慮する正確なフレネルの式ではなく,以下の Schlick の近似式が用いられることが多いようです(今回のシミュレーションもこの式を用いています).
ここで,
フレネル反射とは別に,Blinn-Phong モデルに基づく鏡面反射も考慮しています.詳しい説明については他の資料にゆずるとして,鏡面反射を実装することにより,以下のように水面がキラキラと光るのを再現することができます.
Beer's lawは,光が物質中を通過するときに,その強度が距離に関して指数的に減衰するという法則です.この法則を用いると,以下のように「水の奥の部分ほど深い青に見える」という効果を再現することができます.
実装では,ステップ 4 で求めた流体の厚みに関して,光の強度が指数的に減衰するようにしています.また,RGB のうち,B が R, G よりもゆるやかに減衰するようにすることで,水中が青っぽく見えるようにしています.
Bonus: 10 万粒子をシミュレーションしてみる
冒頭に載せたデモでは,マシンへの負担を大きくしすぎないように,シミュレーションする粒子の最大個数を 4 万個に制限しています.ただ,そこそこの GPU を積んだマシンであれば,粒子の個数が 10 万個くらいになってもリアルタイムにシミュレーションできるようです.以下が 10 万個の粒子のシミュレーション結果になります(RTX 3060 Mobile という GPU を積んだノートパソコン上でのシミュレーションです).
壮観です.やはり 10 万個ともなるとだいぶ見ごたえがありますね!
(デモには,10 万個の粒子をシミュレーションするモードを「隠しモード」として実装しておきました.Number of Particles の下にあるボタンのうち,デフォルトで非表示にしている 5 つ目のボタンを,開発者ツールなどを使って表示することで利用可能になります.4 万粒子のシミュレーションがスムーズに動くのであれば,試してみるといいかもしれません.)
まとめ & 所感
以下がまとめです.
- ブラウザ上でリアルタイムに動作する 3D 流体シミュレーションを,WebGPU を用いて実装した.
- シミュレーションには SPH 法を用い,そのロジックは WebGPU のコンピュートシェーダーで実装した.
- GPU 上における近傍探索の高速化には,カウントソートに基づく手法を用いた.これにより,格子分割による近傍探索の高速化が,固定長のメモリで可能になる.
- 水面のレンダリングには,Screen-Space Rendering を用いた.
- レンダリングのための物理法則としては,フレネル反射,鏡面反射,Beer's law を考慮した.
今回のシミュレーションで,はじめて WebGPU のコンピュートシェーダーに入門してみたわけですが,総じて快適に実装することができたように感じます.特に,試作品として TypeScript で実装した近傍探索のロジックを,特段の困難なくコンピュートシェーダーに移植できたのは嬉しかったポイントです.ハイクオリティな流体シミュレーションを数多く実装している saharan さんの記事には,粒子法の近傍探索の GPGPU 化について「(WebGL の)シェーダだとかなり厳しいです。不可能ではないですが。」と書かれており,WebGL で実装した場合にはここまですんなりと行かなかったのではないかと予想しています.
次の目標は,さらに大規模な流体シミュレーションを実装することです.今回の流体シミュレーションでリアルタイムに扱える粒子の個数は,たかだか 3 万くらいが関の山といったところです(Bonus のところで示したように,GPU によっては 10 万粒子レベルのリアルタイムシミュレーションも可能ですが).一方で,たとえば David という方の Fluid Particles という流体シミュレーションは,粒子の数が 18 万くらいになってもそれなりの速さで動作します.この規模のシミュレーションをブラウザ上で実現するためには,近傍探索を必要としない流体シミュレーションの手法を用いる必要があると考えています.たとえば,PIC/FLIP 法や MLS-MPM といった,粒子法と格子法のハイブリッドな手法は,近傍探索を必要とせず高速なようです(ちなみに David の Fluid Particles は PIC/FLIP 法です).
というわけで,現在は PIC/FLIP 法や MLS-MPM の理論と実装についてちまちま勉強しているところです.最終的には,これらのアルゴリズムを WebGPU のコンピュートシェーダーで実装し,より大規模なシミュレーションをぜひ実現したいです.
参考文献
-
Screen-Space Fluid Rendering に関する GDC2010 のスライド
- 今回,水面のレンダリングに用いた,Screen-Space Fluid Rendering について詳しく説明している有名なスライドです.レンダリングの各ステップについて,対応するコードや画像を交えながら詳細に説明してあります.
-
解説: Just a Pool
- ハイクオリティな流体シミュレーションを数多く実装している saharan さんの記事です.Screen-Space Fluid Rendering を含めた,リアルタイムな水面のレンダリング手法について詳しく書かれており,大いに参考にしました.
-
解説: Jelly
- おなじく saharan さんの記事です.この記事では,流体シミュレーションの 3 つの手法(格子法・粒子法・ハイブリッド)が体系的に説明されており,流体シミュレーションの全体像を把握するのにとても参考になりました.
-
WebGPU Fundamentals
- WebGPU を体系的に学習できるサイトです.WebGPU の基本的事項はすべてここで学びました.頭から読んでいくのはしんどかったので,シミュレーションの実装の中で必要になった箇所をつまみ食いして読む感じで進めました.
- このサイトの作者である Greggman は,後述の wgsl offset computer や,wgsl 向けの行列演算ライブラリである wgpu-matrix の作者でもあり,WebGPU 界では知られた存在です.
-
wgsl offset computer
- WebGPU では,構造体のメモリレイアウトを整えたうえでシェーダにデータを渡す必要がありますが,この作業を手計算でやるのは苦痛です.wgsl offset computer はその作業を自動化してくれる便利ツールであり,今回のシミュレーションの実装を縁の下で支えてくれました.
-
Babylon.js の Screen-Space Rendering に関する説明
- Screen-Space Rendering の実装にあたり,手本となるコードがないとやっぱりきついです.このサイトでは,Babylon.js というグラフィックのライブラリ内で用いられている Screen-Space Rendering の原理とその実装について詳しく説明されています.Screen-Space Rendering の実装で一番参考にしたのはこのサイトです.
-
Coding Adventure: Simulating Fluids
- 流体シミュレーションを一から実装する動画です.この動画のシミュレーションは Unity 上で動くように実装されていますが、これをブラウザ上で動かせるようにしたいと思ったのが,今回のシミュレーションを実装したきっかけです(詳しくは 以前の記事).
Discussion