はじめに
格子ボルツマン法による 2 次元流体計算をブラウザにさせるウェブアプリを作った。ユーザーが自由に描いた物体を一様流中に置き、流れの様子を見守ることができる。流線は LIC 法によって可視化した。アプリの UI とお絵かきの部分は JavaScript、流体計算と可視化の部分は WebAssembly (Rust の wasm-bindgen) によって実装されている。この記事では、実装のために必要だった知見 (を得るために必要な情報) を忘備録としてまとめる。
今回は、アプリで用いた可視化手法についてまとめる。具体的には、HTML の canvas の使い方と流線の計算方法 (LIC) である。流体計算と WebAssembly については、その1の記事にまとめたので、併せて読んで欲しい。
実装されたアプリは以下のページに置かれている。楽しくて癒されるアプリなので、ぜひ遊んで欲しい。高スペックのパソコンやタブレット推奨。
気に入っていただけたら、ぜひ以下のツイートも拡散して欲しい。
2次元流体計算をブラウザに計算させるウェブアプリを作りました。格子ボルツマン法を用いています。ユーザーが自由に描いた物体を一様流中に置いて、流れの様子を見守ることができます。 いろいろイジって遊んでください。 #draw_in_flowhttps://t.co/5mac8evgHF pic.twitter.com/IXq6slXEnh
— みね@太陽の科学館 (@solarphys_info) February 16, 2023
ソースコードは次のレポジトリに置かれている。
canvas の使い方
アプリの図の部分には、HTML の canvas
タグを用いている。これは、JavaScript によってピクセルレベルで色を指定して、その通りに画面に描画させる機能である。ブラウザによっては png 画像としてダウンロードできる。今回は用いていないが、GPU も用いることができる。基本的な使い方はこのページを読めば分かる。
JS では、id か何かで canvas 要素を取得し、そこから "2d"
コンテキストを取り出す。描画は基本的に、このコンテキストが持つ要素や組み込みメソッドを通じて行う。2D コンテキストの持つ機能は上と同じサイトの次のページに分かりやすくまとめられている。
<canvas id="canvas-id"></canvas>
const canvas = document.getElementById("canvas-id");
const ctx = canvas.getContext("2d");
canvas 周りでメモっておきたいことは次の通り。
- canvas の大きさ (高さ & 幅) には次の 2 種類がある。1 つ目は、描画スペースとしての大きさであり、言い換えると、JS で色を指定するときのピクセル数のこと。これは、
canvas.height
,canvas.width
に値を代入することで設定できる。2 つ目は、実際にブラウザに表示される大きさ (CSS と同じ意味での論理ピクセル数) である。これは、canvas.clientHeight
,canvas.clientWidth
を参照することで、取得できる。 - canvas はデフォルトでは左上を原点にして座標を指定することになる。左下を原点にしたいときは、原点を移動させると同時に、Y スケールを反転させる。
ctx.translate( 0, canvas.height );
ctx.scale(1, -1);
- RGB 情報を 256 の整数値として格納した変数によって色を指定したい場合は、次のように、文字列の連結として描画スタイルを指定してから、
ctx.fillRect()
で描画する。
// 0 <= r, g, b <= 255
ctx.fillStyle = 'rgb(' + r + ',' + g + ',' + b + ')';
カラーマップ
計算結果はカラーマップとして図示される。カラーマップには K. Moreland 氏によって開発されたものを用いた。これは Matplotlib の "coolwarm" として用いられているものである。
カラーマップの計算は WebAssembly で行っている。wasm についてはその1の記事で説明しているが、今回は Rust/wasm-bindgen を用いて wasm を生成している。
上のカラーマップを 256 段階に分割し、それぞれの段階における RGB 情報を Rust の u8
整数型 (0 から 255 の値を取る) で格納した。つまり、[u8; 256]
配列 3 個 r
, g
, b
を用意した。あるピクセルに対して計算結果 result
とプロットの最大値 vmax
、最小値 vmin
が与えられれば、次のようにして、プロットのための RGB 情報を得る。
let idx: u8 = ((result - vmin) / (vmax - vmin) * 256.0) as u8;
("RGB 情報") = (r[idx as usize], g[idx as usize], b[idx as usize]);
LIC 法の概要
流線は LIC (line integral convolution) 法を用いた。これは、種となるノイズ的な画像に対して、流線に沿って畳み込み線積分を行うものである。流線が一本の線として可視化されるわけではないが、全体としてなんとなく流線が浮かび上がる。ベクトル場の可視化以外にも、画像を鉛筆画やゴッホ風タッチに変換する芸術的用途でも用いられるようだ。例えば、次のサイトが詳しい。
因みに、LIC 以外の流れ場の可視化手法については、次の文献が詳しい。
まず、種となる画像 $F(\vec{x})$ を用意する。$\vec{x}$ はプロット座標である。今回は、生成されたテクスチャを上述したカラーマップに (RGB に対して一様に) 乗算することで、最終的なプロット画像を導出する。よって、$F(\vec{x})$ はスカラー関数 (モノクロ画像) である。$F(\vec{x})$ にはノイズのような、不規則な画像を用いることになる。これの与え方によって、出来上がる LIC テクスチャの雰囲気がだいぶ変わる。
流れ場 $\vec{v}(\vec{x})$ と初期位置 $\vec{x}_0$ が与えられれば、次の微分方程式を解くことで、$\vec{x}_0$から伸びる流線が媒介変数表示 $\vec{x}(s; \vec{x}_0)$ として求まる。
\frac{d\vec{x}}{ds} = \frac{\vec{v}(\vec{x})}{|\vec{v}(\vec{x})|}
各点 $\vec{x}_0$ から伸びる流線に沿って、次のように線積分することで、LIC テクスチャ $G(\vec{x})$ が生成される。
G(\vec{x}_0) = \left. \int_{s_0}^{s_1} k(s)F[\vec{x}(s; \vec{x}_0)] ds \ \right/ \ \int_{s_0}^{s_1} k(s) ds
線積分に用いるカーネル $k(s)$ と、積分範囲 $s_0$, $s_1$ にも選択の自由度があり、テクスチャの雰囲気を変える。具体的なアルゴリズムは次の文献にまとまっている。
今回用いたパラメータは次の通り。種画像とカーネルの選択には次の文献を参考にした。
- 計算結果は 179 * 89 の等分割格子上で求まるが、流線の可視化は、その 2 倍 (面密度にして 4 倍) の解像度で行った。
- 種画像には、次のように、点を適度に間隔を保ちつつ、不規則に並べたものを用いた。
- 積分範囲について、初期位置から流れ場に沿った方向に積分した。積分する長さは初期位置の流速に比例させた。ただし、無限遠での流速を 1 として、 流速が 2 を超える場合は、 2 で頭打ちにさせた。
- カーネルには次の周期関数を用いた。
- 描画のたびにカーネルの位相を変化させることにより、まるで流線が流れているかのような印象を与えるアニメーションを作り出した。
アニメーション
アプリのアニメーションは再帰関数と requestAnimationFrame()
を用いて実現している。再生ボタンが押されると、次の renderLoop()
が呼び出される。
let animationId = null;
let startTime = new Date().getTime();
const maxFps = 8.0;
const minInterval = 1000.0 / maxFps;
const renderLoop = () => {
fluid.iterate("some_arg"); // wasm に時間積分と可視化を実行させる
drawCells(); // canvas を描画
// 描画の速度が maxFps を超えてしまいそうな場合は無限ループ
let interval;
let endTime;
do {
endTime = new Date().getTime();
interval = endTime - startTime;
} while (interval < minInterval);
startTime = endTime;
animationId = requestAnimationFrame(renderLoop); // 自分自身を呼び出す
};
一時停止や初期化ボタンが押されると、次の文が実行され、描画の連続が止まる。
cancelAnimationFrame(animationId);
animationId = null;
animationId
が null
かどうかで、現在アニメーションが止まっているのかどうかを判別できる。