マントル対流をWeb Workerでリアルタイムに計算する

ブラウザで流体シミュレーションを動かすと聞いて、最初に心配になるのは重さだと思います。毎フレーム格子を解いたらUIが固まるのではないか、と。実際に手を動かすと、同じくらい厄介な壁がもう一つあります。数値の発散です。「地球地下シミュレータ」のマントル対流は教育用の公開アプリの一機能で、訪れた人のブラウザで動き続けることが前提です。発散して画面が壊れることは許されません。この記事では、渦度と流れ関数によるマントル対流をWeb Workerで解き続けるための、解法の選択、計算とUIの分離、壊れたときの自動復帰という三段構えを紹介します。

慣性を持たないマントルの支配方程式

マントルは固体ですが、地質学的な時間スケールでは流体として対流します。粘性が極端に高くて慣性が無視できる流体であり、この極限を無限プラントル数の近似と呼びます。慣性がないと、流れは時間発展しません。各瞬間の温度分布から、その場でつり合う流れが一意に決まります(準静的なストークス流)。時間微分を持つのは温度だけです。

計算領域は断面表示に合わせた2D極座標の環状領域で、マントルの深さを1に正規化し、内半径1.219(核マントル境界相当)から外半径2.219(地殻直下相当)までとしました。半径比は実地球の3480/6336 kmを再現しています。支配方程式は次の4本です。

∇²ω = Ra · (1/r) ∂T/∂θ    (浮力による渦度の生成)
∇²ψ = −ω                   (流れ関数)
u_r = (1/r) ∂ψ/∂θ、 u_θ = −∂ψ/∂r
∂T/∂t + u·∇T = ∇²T         (温度の移流拡散)

ωが渦度、ψが流れ関数です。温度の水平勾配が浮力として渦度を生み、渦度から流れ関数を介して速度場が決まり、その速度が温度を運ぶ、という循環になっています。対流の活発さを決めるのがレイリー数Raで、本実装では1.5e5。格子は周方向192、半径方向48です。初期条件は線形の伝導プロファイルにごく小さな乱数と波数5の摂動を加えたもので、対称性が破れて約5本のプルームが立ち上がります。境界条件は、底面が温度1(核による加熱)、上面が温度0(地表への冷却)の固定です。

「時間発展するのは温度だけ」という構造は、安定性に効きます。速度場は毎ステップ温度から診断し直すので、速度の誤差が蓄積して暴走する経路が最初からありません。

発散しない数値解法の組み合わせ

残るのは温度の時間発展です。素朴に陽解法(前進オイラーと中心差分)で書くと、時間刻みが格子幅と流速で決まる条件を超えた瞬間に発散します。公開物として長時間動き続ける以上、条件を踏み外すと壊れる解法を選ぶわけにはいきません。そこで、時間刻みによらず爆発しない構成として知られるStable Fluids系の組み合わせを使いました。

• セミラグランジュ移流：各格子点から速度場に沿って時間を遡り(後退トレース)、遡った先の温度をバイリニア補間で拾って新しい値にする。補間は既存の値の重み付き平均なので、元の場になかった極値を新しく作らない
• 陰的拡散：(I − dt∇²) T_new = T_adv をJacobi反復16回で解く。陽的な拡散計算と違い、時間刻みの制約を受けない
• SORによるポアソン2連：渦度と流れ関数の2本のポアソン方程式を逐次過緩和法(緩和係数1.7)で解く

移流の後退トレースは、コードにするとこれだけです。

const rBack = r - this.uR[row + i] * dt;
const thetaBack = i * dTheta - (this.uTheta[row + i] / r) * dt;
this.tScratch[row + i] = this.sampleT(rBack, thetaBack);

引っかかるのはSORの反復回数だと思います。1ステップにつき各方程式4スイープしか回しません。ポアソン方程式が4回の反復で解けるはずがない、という感覚は正しいです。解き切ってはいません。代わりに、前ステップの解を初期値に使うwarm startで追随させています。温度場は1ステップでわずかしか変わらないため、前ステップの解と今欲しい解の差は小さく、少ないスイープでも解が準静的についていきます。

場を所有するWeb Worker

解法がどれだけ軽くても、メインスレッドで毎フレーム解けば描画の予算を食います。このアプリは3D地球の描画で60FPS維持を設計要件にしているため、対流計算はWeb Workerに完全に分離しました。分離で決めたことは二つです。場はWorkerが所有すること、通信は間引くことです。

Workerはタイマーで約30Hzに自走し、1回のtickで6サブステップ(無次元時間刻み2e-4)進めます。メインスレッドへの送信は3tickに1回、約10Hzです。

function tick(): void {
  if (!sim) return;
  sim.stepMany(SUBSTEPS_PER_TICK);
  totalSteps += SUBSTEPS_PER_TICK;
  tickCount++;
  if (tickCount % SEND_EVERY_TICKS === 0) {
    sendFrame();
  }
  timer = setTimeout(tick, TICK_MS);
}

送るのは全解像度(192×48)の温度場と、64×24に間引いた速度場です。速度場はシェーダでノイズテクスチャを流れ方向に歪ませる用途にしか使わないので、粗くて足ります。10Hzという頻度は描画の60FPSに比べるとずいぶん疎ですが、マントル対流は見た目のゆっくりした現象で、受信側の断面ヒートマップがノイズの2位相スクロールで更新の合間を埋めるため、離散的な更新は目立ちません。受信時はテクスチャを書き換えるだけで、Reactのstateには載せません。

場の所有権をWorkerに置いた効果は、一時停止に表れます。pauseメッセージはタイマーを止めるだけで場を捨てないので、resumeすると暖まった対流がそのまま続きます。

数値事故からの自動復帰

無条件安定な構成にしても、保証されるのは「爆発しにくい」ことまでです。シミュレーションのパラメータはWorkerへのメッセージで差し替えられる作りで、想定外の値が入る余地があります。公開物なので、防御は多層にしました。

• コンストラクタで時間刻みを5e-4、SORの緩和係数を1.8に上からクランプする
• 温度を毎ステップ[−0.05, 1.05]にクランプし、ブシネスク近似の有界性を強制する
• 100ステップごとにNaNを検査し、見つけたら場を初期状態にリセットする

最後の砦がNaN検査です。

step(): void {
  this.updateFlow();
  this.updateTemperature();
  this.stepCount++;
  if (this.stepCount % NAN_CHECK_INTERVAL === 0 && this.hasNaN()) {
    // 数値事故からの自己回復(公開教育物なので絶対に止めない)
    this.reset();
  }
}

リセットが起きれば対流が最初からやり直しになり、見ている人には不自然に映るかもしれません。それでも、固まった画面や真っ黒なヒートマップよりはましだ、というのが教育用アプリとしての割り切りです。

防御とセットで用意したのが再現性です。初期条件の乱数はmulberry32というシード付きの決定論RNGで生成しており、同じシードなら何度実行しても同じ場に発展します。ユニットテストは10ケースで、静止状態から対流が立ち上がること、1000ステップ回してNaNが出ず温度が有界に保たれること、境界条件と周方向の周期性が維持されること、同じシードで120ステップ後の温度場が全要素一致すること、本番解像度で8サブステップが150ms未満に収まること(テスト上の予算として設定)などを検証しています。決定論RNGを選んだのは、この種のテストを書くためでもあります。

対流から噴火強度への変換

このアプリでは、マントルの上昇流が地表の火山の噴火につながります。火山の位置で上昇流をサンプリングし、マントル温度、マグマ圧、ガス、地殻応力と合わせて重み付き線形結合(重みは順に0.32、0.26、0.22、0.12、0.08)を取り、しきい値を引いた残りを幅0.25のsmoothstepに通して0から1の噴火強度にします。実際の噴火予測とは別物の、教育用の簡略化です。

ただし、噴火に使う上昇流は、Workerで解いている物理シミュレーションからは取っていません。流れ関数を閉じた式で与えた解析的な対流場(対流セル5対)からサンプリングしています。噴火連動は決定論的なテストが要る経路であり、非同期に自走するWorkerの場は表示専用の系統だからです。同じ「マントル対流」でも、見せるための場とロジックを駆動するための場を、意図的に分けています。

ブラウザで流体を動かし続けるための優先順位

冒頭の二つの壁に戻ると、重さにはWorkerへの分離と約10Hzへの通信間引きで、発散には無条件安定な解法とクランプとNaN自動リセットの多層防御で答えたことになります。振り返って効いたのは、「速く正確に解く」より「何があっても止まらない」を優先順位の先頭に置いたことでした。SORを4スイープで打ち切る判断も、NaNで場を作り直す割り切りも、この優先順位から自然に決まっています。研究用の数値計算とは逆向きの妥協ですが、公開アプリの流体は、まず生き続けることに価値があります。

なお、表示される対流の速度は教育目的で大幅に誇張しており、その旨はUI上で注記しています。同じアプリの地震波シミュレーションについては、「地震波の波線追跡でシャドウゾーンを再現する」で扱っています。