🧱 AMR — how numerical relativity actually computes (C++ / WASM)

これはアインシュタイン方程式のデモではありません。数値相対論のコードを組み立てている「機械」のデモです。 機械の方は解く方程式と無関係なので、いちばん単純な波動方程式に載せれば、本物と同じものが動かせます。
右辺を BSSN に差し替えれば、これはそのまま数値相対論のコードになります。

部品何をしているか
method of lines空間 4 次中心差分 → RK4。特性速度から CFL。BSSN でも同じ
Berger–Oliger の多層格子細分化率 2 の入れ子。時間も細分化:レベル l は粗い 1 歩あたり 2l
prolongation細格子のゴーストを親から補間。空間 4 次+時間 2 次。細格子は親が一度も訪れていない時刻の値を要求する ―― 本物のコードがいちばん間違える場所
restriction細格子が追いついたら、その内部を親へ書き戻す
入れ子条件とバッファ各箱は親の内部に厳密に収まる。境界の誤差が次の regrid までに機能へ届かない幅を確保
regridding箱が機能を追って動き、新しく現れたセルを prolongation で埋める。moving punctures と同じ仕組み
拘束条件の監視s − ∂xf = 0 は解析的に真で、進化方程式は課していない。H = 0 と完全に同じ立場なので同じように見張る
読み込み中…

画面は上から、①領域全体の箱(パルスを追って入れ子で動く)、 ②パルス拡大(白=厳密解、青の濃淡=各レベル、赤=粗い格子 1 枚だけの計算、軸下の点=最細レベルの格子点)、 ③subcycling の実績(1 列 1 ステップ。下のレベルほど密=2 倍ずつ多く踏んでいる)。

機械が正しいことの証明は収束試験です。レベルを 1 枚足すと最細の h が半分になるので、 4 次精度なら誤差は 16 分の 1 になるはず ―― subcycling も prolongation も restriction も、 毎ステップの regrid も全部オンのまま:

レベル数最細 h誤差下がった倍率次数
50.06253.575×10⁻²
60.031252.498×10⁻³14.313.84
70.0156251.559×10⁻⁴16.024.00
80.00781259.907×10⁻⁶15.733.98

拘束条件の出方も本物と同じです。レベル 0(h=1、パルス幅 0.4 なので幅あたり 0.4 点=表現不可能)では 1.77 という O(1) の違反、最細レベルでは 2.3×10⁻⁷。 実際の数値相対論でも H の違反は「解像度がぎりぎりの場所」に集中します。

コスト。既定の 8 レベルで格子点 3280、最細 h の一様格子なら 30721 点。 時間刻みも粗いレベルは 128 倍大きく取れるので、1 次元では 35.6 倍の節約です。
ここが肝心で、同じ階層を 3 次元でやると節約は 178196 倍になります(点数が 3 乗で効く)。 これが「AMR は便利」ではなく必須である理由で、一様格子なら 2×10¹⁵ 点という始めることすら不可能な数になります。

作る途中で踏んだバグを 2 つ、隠さずに書いておきます。どちらも「収束していないコード」に見せる種類のものでした。
Kreiss–Oliger 散逸の階数不足。5 点(O(h³))フィルタは 4 次スキームの打ち切り誤差より大きく、 モード k への 1 ステップ倍率が 1−16σsin⁴(kh/2)。最細レベルは 2⁷ 倍のステップを踏むのでこれが積もり、 パルスの振幅を 24% 削っていました。本物の 4 次コードは 7 点(O(h⁵))を使います。
時間補間の履歴管理。あるレベルが最初の 1 歩を踏んだ直後、履歴は t₂ = t₁ なのにフラグが立っていて、 2 次ラグランジュの分母がゼロ → ∞−∞ → NaN。通常の進化では子が親の履歴を読むのがフラグ設定より前なので隠れていて、 regrid だけがフラグ設定後に読むため、そこでだけ発火しました。

②の方が、数値相対論が「重い」ことの正体だと思います。計算量ではありません ―― これはブラウザで 60fps で回ります。 1990 年代の Grand Challenge がスパコンを持っていて合体を出せなかったのは、この階層の話です。
→ ゲージの側で同じことが起きるのは gaugeshock_os、 座標が観測量に効かないことは radar_os

ソース: github.com/yomei-o/universe_cpp