だいぶ前の記事でCPUで大量の粒子を計算することをやりました。
このときは、あまり大量の粒子をけいさんすることはできませんでしたが、今回のGPUを使った計算では、以下のように万を超える粒子を”ほぼ”リアルタイムで計算することができました！
実行環境は、GPU: GeForce RTX 3060Ti、CPU: Intel Core i7-12700 2.10 GHzです。
GPUは、CPUが数個～数十個のコアを持つのに対し、GPUは数千個の小さなコアを持ち、多数の演算を同時並行で処理できます。
WebGPUはWeb上でGPUを使った大規模な並列計算を行える最新の技術です。従来のWebGLと比べて低レベルAPIへのアクセスが可能で、計算処理に特化した機能も備えています。
一方でDEM（個別要素法/Discrete Element Method）は、多数の粒子の挙動をシミュレーションし、それらの相互作用からマクロな現象を再現する解析手法です。土木工学や地盤工学、粉体工学の分野で広く用いられています。
DEMのような大量の粒子間の相互作用計算はGPUの並列処理能力と非常に相性がよく、比較的新しい技術であるWebGPUに適用すれば、ブラウザ上でリアルタイムのシミュレーションができるかなと思って、今回実装しました。
コードとデモページはこちら↓
WebGPUは現在、Chrome、Edge、Firefoxなどの一部のモダンブラウザでサポートされていますが、まだ広く対応しているわけではありませんので注意が必要です。対応情報はこちら。
DEMについて
個別要素法の詳しい内容については詳しく触れませんが、それほど難しい計算は行っていません。
GPUを使う上で工夫した点は、粒子の接触判定の際の隣接粒子探索のプロセスとメモリの削減です。
どちらも既往の研究でいろいろな方法がありましたが、今回のコードは私が理解できる範囲で、比較的簡単な方法を用いているので、最適ではないと思います。
今後より高速に計算するためには、このあたりの処理をうまく効率化することが重要になると考えています。
参考にしたサイトをページ下部に載せておきますので、そちらを参考にしていただければと思います。
実装コード
今回実装したコードは main.jsに集約されています。index.htmlは計算結果を表示するためのUIを提供します。ここでは、main.jsの主要な部分について解説します。
全体の流れ
main.jsの処理は大きく分けて以下の5つのステップで構成されています
- 初期化：WebGPUの利用準備として、デバイスやアダプタの取得と設定
- キャンバス設定：HTML上の描画領域とWebGPUの連携
- データバッファ準備：粒子情報（位置、速度、半径など）のGPUメモリ確保
- 計算処理：粒子の物理演算をコンピュートシェーダーで実行
- 描画処理：計算結果を視覚化するレンダリング
以下、コードをかいつまんで解説します。
初期化
グローバル変数・パラメータの読み取り
|
1 2 3 4 5 |
const params = new URLSearchParams(location.search); function parameter(name, def) { if (!params.has(name)) return def; return parseFloat(params.get(name)); } |
URLパラメータから設定値を読み取る仕組みです。これにより、URLを変えるだけでシミュレーションの条件を簡単に調整できるようにしました。
主な設定パラメータ:
min_radius/max_radius: 粒子の最小/最大サイズ(m)balls: シミュレーションする粒子の総数width/height: 描画領域のサイズ(px)width_length: シミュレーション空間の実際の幅(m)x1,y1,x2,y2: 固定線分の座標(m)
GPUリソースの初期化
|
1 2 3 4 5 6 7 |
if (!("gpu" in navigator)) fatal("WebGPU not supported."); const adapter = await navigator.gpu.requestAdapter(); const device = await adapter.requestDevice({ requiredLimits: { maxStorageBuffersPerShaderStage: 10, // 上限値の変更 }, }); |
このコードは、WebGPUの基本的なセットアップを行います：
- まずブラウザがWebGPUに対応しているか確認
adapterの取得（GPUとやり取りするためのインターフェース）deviceの取得（実際の計算リクエストを送るオブジェクト）
requiredLimitsでストレージバッファの上限値を引き上げて、より多くのストレージバッファを使えるようにしています。
キャンバス設定
|
1 2 3 4 5 6 7 8 9 10 |
const canvas = document.getElementById("webgpuCanvas"); const width = parameter("width", 1000); const height = parameter("height", 500); canvas.width = width; canvas.height = height; const ctx = canvas.getContext("webgpu"); ctx.configure({ device: device, format: navigator.gpu.getPreferredCanvasFormat(), }); |
index.htmlで準備された <canvas id="webgpuCanvas">に幅と高さを設定して webgpuコンテキスト（getContext("2d")のようなもの）を取得します。そして、canvas要素とWebGPUを連携させます
configureで実際に描画に用いるピクセルフォーマットを指定します（formatで現在のシステムで最適なものが適用される）。
この設定により、後ほどGPUで計算した結果をブラウザ上に表示できるようになります(私もこのあたりはよくわかっていない部分があります)。
データバッファの準備
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
const NUM_BALLS = parameter("balls", 12000); const BUFFER_COUNT = 8; const BUFFER_SIZE = NUM_BALLS * BUFFER_COUNT * Float32Array.BYTES_PER_ELEMENT; // 各粒子の初期値を定義 let inputBalls = new Float32Array(new ArrayBuffer(BUFFER_SIZE)); for (let i = 0; i < NUM_BALLS; i++) { inputBalls[i * BUFFER_COUNT + 0] = random(minRadius, maxRadius); // 半径 // ...省略... } // 粒子要素の初期化処理 function initializeInputBallsBuffer(){ return inputBalls; } // ...省略... // バッファの作成 const input = device.createBuffer({ size: BUFFER_SIZE, usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST, }); // ...省略... // inputバッファへ書き込み device.queue.writeBuffer(input, 0, initializeInputBallsBuffer()); |
GPUでの計算で用いる粒子の初期値の準備とバッファメモリの作成を行います。
最後に device.queue.writeBufferでバッファメモリに書き込みを行います。
このあたりがWebGPUを使う上でややこしくなってくるポイントだと思います。
基本的な流れは、
- バッファをあらかじめ準備（バッファの種類やサイズなど）
- 書き込みを行ってGPU内の計算で使用
です。
Float32Array(new ArrayBuffer(BUFFER_SIZE))について、
|
1 |
const BUFFER_SIZE = NUM_BALLS * BUFFER_COUNT * Float32Array.BYTES_PER_ELEMENT; // 粒子の配列に必要なバッファサイズ |
ですので、粒子数（12000個） ✕ 各粒子の情報（位置、速度、…）がここでは8個（BUFFER_COUNT） ✕ 情報のそれぞれが単精度の浮動小数点数（=4byte）=384000のバッファサイズの配列ということです。
具体的な配列の中身は[粒子番号1の位置x, 粒子番号1の位置y座標, …, 粒子番号nの位置x座標, 粒子番号nの位置y座標, …]と1次元の長い配列になっています。
計算処理
コンピュートシェーダー
|
1 2 3 4 5 6 7 8 9 |
const module = device.createShaderModule({ code: ` // WGSL での衝突判定と力学計算 @compute @workgroup_size(64) fn main(@builtin(global_invocation_id) global_id : vec3<u32>) { // 衝突判定 → 力の計算 → 位置更新 など } `, }); |
ここでWGSL（WebGPU Shading Language）によるコンピュートシェーダーが記述されています。
いわゆるGPU計算の部分です。
今回はDEMの計算を行うので、粒子同士の衝突の計算や、粒子と壁や線分の衝突の計算を行っています。ワーキンググループが粒子ごとに分かれており、1粒子ごとにmain関数が実行されているイメージです。
ワーキンググループについて
今回は粒子の数だけGPUでスレッド数を割り当てて計算しています。
以下のコードを見てください。
|
1 2 3 4 5 6 7 8 9 10 11 |
// 計算を更新する関数 function updateCompute(commandEncoder) { const dispatchSize = Math.ceil(NUM_BALLS / 64); const computePassEncoder = commandEncoder.beginComputePass(); // 「これから計算処理の指示を記録します」という宣言 computePassEncoder.setPipeline(pipeline); // どのシェーダープログラムを実行するかを指定 computePassEncoder.setBindGroup(0, bindGroup1); // シェーダーが使用するデータリソースを指定（粒子） computePassEncoder.setBindGroup(1, bindGroup2); // シェーダーが使用するデータリソースを指定（線分） computePassEncoder.dispatchWorkgroups(dispatchSize); // 「ワークグループ」と呼ばれる計算ユニットをいくつ起動するかを指定 computePassEncoder.end(); // 「計算処理の指示はここまでです」という宣言 return commandEncoder; } |
ワーキンググループは、dispatchWorkGroupsでその数が割り当てられます。この上限値は device.limitに記載されている maxComputeWorkgroupsPerDimension: 65535です。
さらに、先程のコンピュートシェーダで @workgroup_size(64)と指定しているので、各ワーキンググループが64個のスレッドに分けられます。
つまり、この例では、ワーキンググループを粒子数（12000個）/64個=188（小数点切り上げ）個作成し、1つのワーキンググループは64個に分けられるので、結果188✕64個=12032個のスレッドができるため、12000個の粒子を別々のスレッドで扱えることになります。以下イメージ図です。
ちなみに device.limitで maxComputeWorkgroupSizeX : 256 maxComputeWorkgroupSizeY : 256 maxComputeWorkgroupSizeZ : 64とあります。これは @workgroup_size()の上限を表しており、理論的には「(65535 × 256) × (65535 × 256) × (65535 × 64)」がスレッド総数の上限となる計算です。ただし、実際には使うハードウェアにより、制限があるようです。
diapatchWorkGroups()と @workgroup_size()は3次元で表すことができ、diapatchWorkGroups(188)は diapatchWorkGroups(188, 1, 1)、 @workgroup_size(64)は @workgroup_size(64, 1, 1)と等しいです。
今回は1次元（xのみ）を使っているので、main関数の global_id.xが何番目のスレッドかを表します。なので、その値をそのまま粒子番号として計算しています。
ワーキンググループの説明はこちらがわかりやすかったです。
コンピュートシェーダーへのバッファ受け渡しの流れ
@computeの main関数に device.queue.writeBufferで書き込んだ変数がわたるまでの流れは以下のとおりです。
createBindGroupLayoutでどのバインディング番号にどのバッファを割り当てるかを以下で定義します（binding=0→ storageバッファ、binding=3→ uniformバッファ など）。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
const bindGroupLayout1 = device.createBindGroupLayout({ entries: [ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage", }, }, // ...省略... ], }); const bindGroupLayout2 = device.createBindGroupLayout({ entries: [ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage", }, }, // ...省略... ], }); |
作成した bingGroupLayoutを createPipelineLayoutでパイプラインレイアウトに適用して、createComputePipelineで コンピュートシェーダーの main関数をエントリーポイントとして紐づけて、パイプラインを作成します。
|
1 2 3 4 5 6 7 8 9 |
const pipeline = device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [bindGroupLayout1, bindGroupLayout2], }), compute: { module, entryPoint: "main", }, }); |
まだ作成したバッファメモリの変数は対応付けられていません。
以下で、 createBundGroupで bindGroupLayoutと作成した変数を対応付けた（binding=0はinputバッファ、binding=1はoutputバッファ のように）bindGroupを作成します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
const bindGroup1 = device.createBindGroup({ layout: bindGroupLayout1, entries: [ { binding: 0, resource: { buffer: input, }, }, // ...省略... ], }); const bindGroup2 = device.createBindGroup({ layout: bindGroupLayout2, entries: [ { binding: 0, resource: { buffer: line_input, }, }, // ...省略... ], }); |
最後に、device.createCommandEncoderで作成したコマンドエンコーダーにパイプラインをセットし、setBindGroupで @group(0)と 作成した bindGroup1と @group(1)と bindGroup2を対応させるように指示します。
これで、main()関数に bindGroupで指定したバッファ変数が渡されます。
|
1 2 3 4 5 6 7 8 9 10 |
// GPUコマンドの作成 let commandEncoder = device.createCommandEncoder(); // 一度計算する commandEncoder = updateCompute(commandEncoder); // updateCompute関数の中身 // computePassEncoder.setPipeline(pipeline); // どのシェーダープログラムを実行するかを指定 // computePassEncoder.setBindGroup(0, bindGroup1); // シェーダーが使用するデータリソースを指定（粒子） // computePassEncoder.setBindGroup(1, bindGroup2); // シェーダーが使用するデータリソースを指定（線分） // ...省略... device.queue.submit([commandEncoder.finish()]); // commandEncoder.finish()でコマンドリストを完成させ、device.queue.submit()でGPUに送信することで、実際の処理が行われる |
描画処理
WGSLを用いて計算結果を canvas要素に描画します。
シェーダーモジュール
レンダリングのシェーダーモジュールは以下のように頂点（vertex）とフラグメント（fragment）に分かれています。
頂点シェーダーは描画する点の位置を記述します。フラグメントシェーダーは色を記述します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
// シェーダーモジュール const cellShaderModule = device.createShaderModule({ label: "Cell ader", code: ` // 構造体の定義 // 頂点シェーダーを定義する vertex fn vertexMain (@builtin(vertex_index vertex(Red, Green, Blue, Alpha)dex : u32, @built回転をズーム倍率を考慮してinstance_index省略...) instance_index: u32, @ // (X, Y, Z, W) w の値は、3 次元同 次座標系における頂点の 4 つ目の要素location(0) (Red, Green, Blue, Alpha): vec2<f32>) -> VertexOutput { // 粒子の位置と回転をズーム倍率を考慮して計算 // ...省略... return vertexOutp ut; // (X, Y, Z, W) w の値は、3 次元同(Red, Green, Blue, Alpha)系における頂点の 4 つ目の要素 } // フラグメント シェーダーを定義する（戻り値は0～1） @fragment fn fragmentMain(input: VertexOutput) -> @location(0) vec4f { return input.cellColor; // (Red, Green, Blue, Alpha) } ` }); |
頂点シェーダー @vertex
vertexMain関数について、
@builtin(vertex_index)は“今処理している頂点が何番目か”を示す変数
@builtin(instance_index)は“今処理しているインスタンスが何番目か（何番目の粒子か）”を示す変数
@location(0) pos: vec2<f32>は頂点バッファのレイアウトで location=0に対応する変数
となります。
今回の例では、locationでわたってくる円を表す頂点のデータを、コンピュートシェーダーで計算した粒子の位置を使って、適切な位置に表示させるように移動や回転をさせています。
頂点シェーダーへのlocationの変数の受け渡しの流れ
@locationに変数がわたる流れは以下のとおりです。
粒子に対応する頂点の位置を求めて、配列に格納します。粒子の形状を以下のように定義します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
// 頂点を定義する function createCircleVertices(centerX, centerY, radius, segments) { const vertices = []; for (let i = 1; i < segments; i++) { // 一つ欠けさせるために0番目はスキップ const theta1 = (i / segments) * 2.0 * Math.PI; const theta2 = ((i + 1) / segments) * 2.0 * Math.PI; const x1 = centerX + radius * Math.cos(theta1); const y1 = centerY + radius * Math.sin(theta1); const x2 = centerX + radius * Math.cos(theta2); const y2 = centerY + radius * Math.sin(theta2); // 頂点データに中心点と円周上の点を追加 vertices.push(centerX, centerY, x1, y1, x2, y2); } return new Float32Array(vertices); } // 頂点バッファを作成する const vertices = createCircleVertices(0, 0, 1, 10); |
ここでは、半径1の粒子を10の三角形との集合で表すようにしています。以下のようなイメージです。
1つ欠けさせているのは、粒子自体の回転を表現したかったからです（ただの円だと回転しているか見えづらいから）。
ここで作成した頂点を以下で、vertexBufferというGPUのバッファに書き込みます。
|
1 2 3 4 5 |
const vertexBuffer = device.createBuffer({ size: vertices.byteLength, usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.VERTEX, }); device.queue.writeBuffer(vertexBuffer, 0, vertices); |
頂点バッファのレイアウトを定義します。@location(0)に対応させるために、shaderLocation: 0としています。
|
1 2 3 4 5 6 7 8 9 10 |
const vertexBufferLayout = { arrayStride: 8, // vec2<f32>なので8バイト attributes: [ { shaderLocation: 0, // シェーダーの @location(0) に対応 offset: 0, format: "float32x2", }, ], }; |
レンダリングのパイプラインを作成するときに、vertexの buffersに上記の vertexBufferLayoutを設定します。まだ、vertexBufferLayoutに vertexBufferは対応付けられていません。
|
1 2 3 4 5 6 7 8 |
const cellPipeline = device.createRenderPipeline({ vertex: { module: cellShaderModule, entryPoint: "vertexMain", buffers: [vertexBufferLayout], // ここでレイアウトを関連付ける }, // ... }); |
最後に描画コマンドを作成するときに、以下のように setVertexBufferで @location(0)と vertexBufferを対応させるように指示します。
|
1 2 3 |
renderPassEncoder.setPipeline(cellPipeline); renderPassEncoder.setVertexBuffer(0, vertexBuffer); renderPassEncoder.draw( /* vertexCount */ , /* instanceCount */ ); |
これで頂点データがシェーダーの @location(0)に渡されます。
フラグメントシェーダー @fragment
続いて、fragmentMain関数について、fragmentMainの引数 inputは、vertexMainの戻り値である VertexOutputを受け取ります。
そして、fragmentMainの戻り値は色を表していて、頂点で形成された三角形のピクセル毎にこの @fragmentが呼び出されるようです。
今回は粒子番号や速度で切り替えられるようにしています。vertexMain関数のvertexOutput.cellColorで計算した色(r, g, b, a)をそのまま表示しています。
メイン関数
上記のコードを準備した後、以下のメイン関数をで計算をrequestAnimationFrameで回します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
// メインループ async function updateGrid(actualFps) { // 実際のFPSから反復回数を計算 const iteration = Math.ceil(1 / actualFps / dt); // コマンドを記録するエンコーダーを作成する let commandEncoder = device.createCommandEncoder(); // cpu側でループする for (let i = 0; i < iteration; i++) { commandEncoder.clearBuffer(gl_atomic); // 初期化 commandEncoder.clearBuffer(gl_output); // 初期化 // 計算更新 commandEncoder = updateCompute(commandEncoder); commandEncoder.copyBufferToBuffer(output, 0, input, 0, BUFFER_SIZE); // 粒子要素コピー commandEncoder.copyBufferToBuffer(gl_output, 0, gl_input, 0, GL_BUFFER_SIZE); // 粒子のインデックス配列をコピー } // END cpu側でループする場合 // レンダリング更新 commandEncoder = updateRender(commandEncoder); // 一連のコマンドをパッケージ化する const commands = commandEncoder.finish(); // コマンド送信 device.queue.submit([commands]); counter++; } |
requestAnimationFrameは、実行間隔は一定ではないため、前回からの時間間隔を記録しておいて、メイン関数に渡しています。
そこからその時間で計算する回数を計算し、その回数分GPUでの計算をコマンドに記載していきます。
計算が終わるとレンダリングを実行します。
DEMは設定パラメータによって、時間ステップが結構短くなる場合があるので、試行錯誤的に決定しています。
描画範囲の拡大・縮小・移動
一応実行結果を見やすくするために、描画範囲の拡大・縮小・移動をできるようにしています。
これらの処理はvertexMainの計算にも関わってくるので、少しめんどうでした。
動作確認
一応GPUで動いているか確認しました。
GPUの3D使用率が上昇していることが確認できました！
まとめ
今回のコードでは、WebGPUを用いた2次元のDEMを実装しました。コンピュートシェーダーを使って、並列計算をすることで、数万の粒子をリアルタイムで描画できるのは、これまでの経験上すごく魅力的な体験でした。
今後3次元化や描画方法の見直しを行うことで、より面白い計算が行えると思います。
あと最近はコードを調べるにしてもChatGPT等の生成AIに聞くことがほとんどで、今回のコードもだいぶ助けられました。なので、今回の記事のコードの解説などはそのうちなくなって行くんですかね。
参考サイト
DEM関係
- 内田吉彦, 伯野元彦: 粒状体シミュレーションによる岩屑流・土石流の解析
- 山本修一: 個別要素法による粒状体の力学的挙動に関する解析的研究（その1）
- 酒井幹夫: 粉体の数値シミュレーション, 丸善出版(2012)
WebGPU関係
コメント