# Beyond Exascale: Dataflow Domain Translation on a Cerebras Cluster > [!abstract] 概要 > 物理システムのシミュレーションは科学および工学の諸分野において不可欠である。一般的に用いられている領域分割法（Domain Decomposition Methods）は、ネットワーク接続された計算環境において高いシミュレーション速度と高い計算資源利用率を同時に達成することができない。とりわけエクサスケールシステムであっても、これらのワークロードに対してピーク性能のわずかな割合しか発揮できていない。本論文では、これらの限界を克服するために設計された新しい「領域平行移動（Domain Translation）」アルゴリズムを導入する。64台のCerebras CS-3システムからなるクラスタ上で本手法を用い、多岐にわたる指標において前例のないクラスタ性能を実証する。毎秒160万タイムステップを超えるシミュレーション実行を示すとともに、ピーク性能の88%において完全な弱スケーリング（weak scaling）を実証する。このクラスタ規模において、本実装は電力無制約環境で112 PFLOP/s、電力制約環境で57 GFLOP/Jを提供する。本手法の有用性を示すため、惑星規模での小惑星衝突に伴う津波を460メートル解像度の浅水方程式（Shallow Water Equations）を用いてモデル化・シミュレーションを行う。 ## 問題設定 偏微分方程式（PDE）の離散化に基づく数値シミュレーション（有限差分法、有限要素法、有限体積法など）は、流体、熱伝導、波動、電磁気学など物理現象の解析に不可欠である。しかし、従来のフォン・ノイマン型スーパーコンピュータによる分散クラスタ環境では、次の根本的な課題が存在する。 - **メモリウォールと強スケーリングの停滞**: 弱スケーリングによって空間解像度や格子規模を拡大することは可能であるが、フォン・ノイマン型アーキテクチャのメモリ帯域制約（メモリウォール）により、時間発展の進行速度（強スケーリング）は頭打ちとなっている。典型的な地球システムモデルの実効性能はピーク比5%未満に留まり、ペタ・エクサスケール級の大規模マシンでも実効1.2〜8 PFLOP/s（最高でも25.96 PFLOP/s）しか達成できていない。 - **従来の領域分割法（Domain Decomposition）の限界**: - **静的分割（Static Partitioning）**: 格子を各ノードに固定配置する方式では、境界に隣接する格子点が毎タイムステップごとにノード間ネットワーク遅延 $\tau$ の影響を受ける。この結果、シミュレーション全体の反復速度は最も遅いリンクの遅延 $1/\tau$ に律速される（**Figure 1b**）。 - **ゴーストセル重複法（Overlapping / Ghost Cell Method）**: ノード境界に $g$ 層のゴーストセルを配置して遅延を償却する手法では、目標反復周波数 $f$ に対し $g \propto \tau f$ が要求される。$d$ 次元領域においてゴースト領域が占める体積割合は $1 - ((n-2g)/n)^d$ となり、計算資源利用率が $(\tau f)^d$ の多項式で急速に悪化する。時間ステップ速度を高めようとすると冗長計算と電力消費が爆発するトレードオフを抱える。 **Table 1: 記号と変数** | 記号 | 単位 | 説明 | |---|---|---| | $\tau$ | s | ネットワークリンクの遅延時間（Time of flight） | | $n$ | pts | 各ノードが担当する格子点の線形スパン | | $g$ | pts | ゴースト領域の厚さ（重複幅） | | $p$ | pts | ステンシル半径（マンハッタン距離） | | $t$ | s | 1タイムステップの反復計算時間 | | $f$ | Hz | タイムステップ反復周波数 | | $c$ | s | 1格子点あたりのシリアル更新時間 | | $d$ | - | 空間の次元数 | | $w$ | cores | コア配列の線形スパン（ファブリックサイズ） | | $\phi$ | - | 緯度（Latitude） | | $\lambda$ | - | 経度（Longitude） | | $\ell$ | - | コリオリ力係数（Coriolis force） | **Figure 1: 1次元3点ステンシルによる領域分割の比較** ![[_attachments/arxiv-2511.11542/fig01-stencil-domain-translation-concept.png]] - **Figure 1a**: 単一ノードでは外部依存がなく各時間ステップが1単位時間で進む。 - **Figure 1b**: 固定領域分割ではノード境界をまたぐ通信遅延（例：10単位時間）が毎ステップ加算される。 - **Figure 1c**: 分割境界を毎ステップ1格子点ずつシフトさせることで、境界をまたぐ依存関係が単方向化され、ノード間遅延は1回のみ印加される。 ## 提案手法 本論文は、ノード間ネットワーク遅延を完全に隠蔽する並列アルゴリズム「**Domain Translation（領域平行移動）**」を提案し、ウェーハスケールエンジン（WSE）クラスタ上に実装する。 ### 1. 領域平行移動アルゴリズム（Domain Translation） - **格子・プロセッサ対応関係の動的シフト**: ノード間をリング（またはトーラス）状に接続し、反復ごとに格子点からプロセッサへのマッピングをステンシル幅 $p$ だけ平行移動させる。 - **単方向通信トラフィックへの変換**: ステンシル計算本来の双方向通信（$\pm p$）と移動処理を合成することで、ネットワークリンク上のデータ転送が移動方向への $2p$ の**単方向トラフィック**へと変換され、逆方向への転送が完全に消滅する（**Figure 1c**）。 - **遅延の線形加算の排除**: 各格子点はノードの担当領域全体を横断した後にのみネットワーク遅延を経験するため、遅延の影響がサブドメイン幅 $n$ 全体へと償却される。サブドメインが一定の臨界サイズを超えると、ネットワーク遅延による利用率損失はゼロとなる。 **Figure 2: 領域平行移動法における時空状態遷移の定常サイクル** ![[_attachments/arxiv-2511.11542/fig02-spacetime-domain-translation-duty-cycle.png]] 格子点（x軸）とタイムステップ（y軸）の関係を示す。ノードがパッケージ（青）を受信すると、保持データ（黄）を用いて計算スイープを実行し新たな格子点状態（赤）を生成する。三角形の時空図が平行四辺形へと伸長し、受信データが到着するまでの間に $n/(2p)$ ステップ分の計算を先行実行して遅延を完全に隠蔽する。 **Table 2: 領域分割手法の比較** | 手法 | 利用率乗数 | 周波数限界 | 最大速度時の利用率 | |---|---|---|---| | 静的領域分割（Static Domain Decomposition） | 1 | 遅延（Latency） | $c/\tau \ll 1$ | | 重複領域分割（Overlapping Domain Decomposition） | $(1 - f\tau/n)^d$ | 帯域幅（Bandwidth） | $1/(f\tau)^d \ll 1$ | | **領域平行移動（Domain Translation, 本手法）** | **1** | **帯域幅（Bandwidth）** | **1** | ### 2. ハードウェア制約と性能限界 ノードの性能限界は計算限界、遅延限界、帯域限界の3要素によってモデル化される（**Figure 3**）。 - **計算限界（Compute Limit）**: 1点あたりの計算時間を $c$ とすると、ノードは $1/(cn)$ のレートでデータを生成・消費する。 - **遅延限界（Latency Limit）**: パイプライン内のパッケージ間隔から、遅延限界リンクの供給レートは $n/(2p\tau)$ となる。完全な計算バウンド利用率を達成するための条件は次式で与えられる。 $n^2 > \frac{2p\tau}{c}$ すなわち、隠蔽可能な遅延時間は**ノードが保持するメモリ量（格子点数 $n$）の2乗に比例**して拡大する。 - **帯域限界（Bandwidth Limit）**: 境界格子点数 $n^{d-1}$ に比例して伝送レートが制限される。格子点数 $n$ を十分に大きく確保することで、遅延および帯域幅の影響を完全に覆い隠し、ノードを100%の計算カーネル効率で動作させることができる。 **Figure 3: 格子点数とタイムステップ反復周波数の関係** ![[_attachments/arxiv-2511.11542/fig03-performance-limits-timestep-rate.png]] 格子点数が小さい領域ではネットワーク遅延が支配的であるが、格子点数の増加に伴い帯域限界を経て計算限界へと移行し、フルカーネル利用率に到達する。 ### 3. WSE クラスタへのマッピングとDSLフレームワーク - **空間データフロー実行モデル**: WSEはフラットなPEグリッド（各PEに局所SRAMとルータ）で構成され、グローバル同期なしに自律的データフローで駆動される。計算平面を時空上で45度傾斜させることで（相対論のペンローズ図に類似）、すべてのデータ移動が時間方向2ホップ、空間方向1ホップ以内に収まり、常に局所SRAM内で演算が完結する。 - **Tungsten DSL フレームワーク**: 汎用ステンシルフレームワークを約1,000行のTungsten言語コードで実装。コア間・ノード間通信と全カーネルを統一的に記述する。 **Table 3: 5点ステンシル熱伝導方程式のメインタイムステップループ（Tungsten言語実装）** ``` function sendrecv() { parallel { // Send right and up ∀ i ∈ [w, n + w) ∀ j ∈ [n, n + w) right[] ← x[i][j]; ∀ i ∈ [n, n + w) ∀ j ∈ [w, n + w) up[] ← x[i][j]; // Receive from left and from below ∀ i ∈ [w, n + w) ∀ j ∈ [0, w) x[i][j] ← left[]; ∀ i ∈ [0, w) ∀ j ∈ [w, n + w) x[i][j] ← down[]; } parallel { // Send corner data received from left to above ∀ i ∈ [n, n + w) ∀ j ∈ [0, w) up[] ← x[i][j]; // Receive corner data from below ∀ i ∈ [0, w) ∀ j ∈ [0, w) x[i][j] ← down[]; } } function compute() { let i ∈ [0, n); let j ∈ [0, n); ∀ i ∀ j y[i + w][j + w] ← avec[0]*x[i + p][j + p - 1]; ∀ i ∀ j y[i + w][j + w] += avec[1]*x[i + p - 1][j + p]; ∀ i ∀ j y[i + w][j + w] += avec[2]*x[i + p + 1][j + p]; ∀ i ∀ j y[i + w][j + w] += avec[3]*x[i + p][j + p + 1]; ∀ i ∀ j y[i + w][j + w] += avec[4]*x[i + p][j + p]; ∀ i ∀ j x[i + w][j + w] ← y[i + w][j + w]; } function innerloop(sp niter) { sp iter = 0.0; while (iter < niter) { sendrecv(); compute(); iter += 1.0; } } ``` **Table 4: 対象ワークロードのハイパーパラメータ** | システム | 場変数数（Field Vars） | FLOPs/点/ステップ | |---|---|---| | 熱伝導方程式（HE, 5点ステンシル） | 1 | 9 | | 熱伝導方程式（HE, 9点ステンシル） | 1 | 17 | | 浅水方程式（SWE, 偶数ハーフステップ） | 3 | 94 | | 浅水方程式（SWE, 奇数ハーフステップ） | 4 | 61 | | 浅水方程式（SWE, 完全タイムステップ） | 7 | 155 | ### 4. クラスタスイッチトポロジとルーティング - **スイッチトポロジ**: 各CS-3ノードは12本の100Gbpsポートを備え、12台の64ポートEthernetスイッチ群を介して64ノードが全結合される（**Figure 4**）。 - **鏡像チェッカーボード配置**: ハードウェアの物理ポート配置とアプリの送受信方向の不一致を解決するため、Y軸方向の鏡像コードを交互に配置するチェッカーボードパターンを採用（**Figure 5**）。 - **オンウェーハルーティング**: 水平トラフィック（紫）と垂直トラフィック（濃緑・淡緑・青）を分離し、8個の仮想チャネル（VC）を用いてチップ内デイジーチェーンを構築（**Figure 6**）。 **Figure 4: 任意数のWSEを接続するスイッチトポロジ** ![[_attachments/arxiv-2511.11542/fig04-switch-topology-cluster.png]] **Figure 5: クラスタトポロジとY軸鏡像チェッカーボード配置** ![[_attachments/arxiv-2511.11542/fig05-cluster-topology-mirroring.png]] **Figure 6: オンウェーハルーティングとデイジーチェーン機構** ![[_attachments/arxiv-2511.11542/fig06-on-wafer-routing-daisy-chain.png]] ## 新規性 1. **WSEクラスタによる世界初の分散偏微分方程式（PDE）ソルバー**: 単一ウェーハに閉じていたデータフロー空間計算を、64ノード・5,760万コア規模の分散クラスタ環境へ拡張。 2. **ノード間遅延を完全に無効化する領域平行移動アルゴリズム**: 従来のゴーストセル方式のように利用率を犠牲にすることなく、約10μsのEthernet遅延を100%隠蔽して計算バウンドなスケーリングを達成（**Table 2**）。 3. **ステンシル計算における前例のない利用率（88%）と電力効率（57 GFLOP/J）**: フォン・ノイマン型クラスタが5%未満で喘ぐ中、88%の実効利用率と毎秒160万ステップ超の時間発展速度を実証。 ## 実験 ### 実験設定と性能モデル 64台のCerebras CS-3クラスタ（クロック750MHzおよび1.2GHz）を用い、弱スケーリング、強スケーリング、問題サイズの網羅的測定を実施した（**Table 5**）。 **Table 5: 性能特性評価に用いたグリッドスイープパラメータ設定** | パラメータ | 次元 | 値 | パラメータ範囲 | |---|---|---|---| | 弱スケーリング（Weak Scale） | ノード数 | 1, $2^2$, $4^2$, 4×6, 2×30, 6×10, $8^2$ | 1〜64 ウェーハスケールノード | | 強スケーリング（Strong Scale） | 点/コア | $2^2$, $4^2$, $8^2$, $16^2$, $32^2$, $64^2$ | 4〜4,096 格子点/コア | | ノード規模（Node Size） | コア/ノード | $720^2$, 744×1116 | 50万〜90万コア/ノード | | ノードクロック（Node Clock） | GHz | 0.75, 1.2 | 0.8〜2 PFLOP/s ピーク性能 | **Table 6: 各種PDEシミュレーションの性能モデル** | コード | 性能モデル $f(n)$ [cycle/step] | Flops/step | 漸近利用率（$n \to \infty$） | 実測利用率 | I/O [words/step] | |---|---|---|---|---|---| | 5点熱伝導方程式 | $105 + 3.74n + 6.72n^2$ | $9n^2$ | 67% | 67% ($n=64$) | $4n + 4$ | | 9点熱伝導方程式 | $97 + 3.5n + 9.37n^2$ | $17n^2$ | 91% | 88% ($n=64$) | $4n + 4$ | | 浅水方程式（SWE） | $1026 + 183.2n + 137.6n^2$ | $155n^2$ | 56% | 53% ($n=24$) | $64n + 80$ | ### 1. 5点熱伝導方程式の弱スケーリングと強スケーリング - **弱スケーリング**: 4台から60台のノード数拡大において、98.8%（$n=2$）から99.9998%（$n=64$）というほぼ完璧な弱スケーリングを記録（**Figure 7**）。 - **性能モデルとの一致**: 1コアあたり256点未満ではI/Oバウンド、256点以上で計算バウンドへと移行し、ピーク比66%（1.32 FLOPs/cycle/コア）を達成（**Figure 8**）。 - **強スケーリング**: 各種問題規模（6.8B点〜54B点）において、64ノードまで完璧にスケーリングし、毎秒160万タイムステップ超を記録（**Figure 9 Top**）。 **Figure 7: 5点熱伝導方程式の弱スケーリング（4〜60 CS-3ノード）** ![[_attachments/arxiv-2511.11542/fig07-heat-eq-weak-scaling.png]] **Figure 8: 性能モデル予測と実測値の比較（I/Oバウンドと計算バウンドの遷移）** ![[_attachments/arxiv-2511.11542/fig08-performance-model-vs-measured.png]] **Figure 9: 熱伝導方程式の強スケーリング（上）および弱スケーリング（下）** ![[_attachments/arxiv-2511.11542/fig09-strong-and-weak-scaling-heat-eq.png]] 上段：750MHz動作における強スケーリング。下段：1.2GHz動作における9点熱伝導方程式の弱スケーリング。ソフトウェア電力最適化版で84.7 PFLOPS、強化電源版で112 PFLOPS（利用率88%）を達成。 ### 2. 9点熱伝導方程式とピーク性能・電力効率 - **電力最適化版（1.2GHz動作）**: 高クロック時の電力スロットリングを抑えるためコア内低電力メモリを活用したコードを開発。64ノード（1,920億格子点）で**84.7 PFLOPS**、電力効率**57 GFLOP/J**（57 GFLOP/W）を達成（**Figure 9 Bottom**）。Green500首位（JEDIの密行列72.7 GFLOP/W）に迫る効率を、疎なステンシル計算で記録。 - **強化電源ノードによるフル稼働**: 電力スロットリングのない特別仕様ノードにおいて、1コアあたり2.1 GFLOPS（ピーク比88%）を達成。64ノード換算で**112 PFLOPS**の性能に相当。 ### 3. 地球規模浅水方程式（SWE）による小惑星衝突津波シミュレーション - **実地球地形データセット**: GEBCO 2024 Grid（15秒角、解像度462m）を使用。 - **衝突津波の伝播**: 衝突運動エネルギーの90%（TNT換算240万トン相当）が海洋のポテンシャルエネルギーに変換された条件（水塊隆起200m、範囲3万km²）をモデル化。衝突後14時間にわたる惑星規模の津波伝播をシミュレート（**Figure 10**）。サンフランシスコ湾岸への局所波浪到達を高解像度に再現（**Figure 11**）。 - **性能**: 浅水方程式においてもほぼ完全な弱スケーリングを維持し、計算バウンド領域でピーク比53%の実効性能を達成。 **Figure 10: 小惑星衝突14時間後の惑星規模津波伝播シミュレーション** ![[_attachments/arxiv-2511.11542/fig10-planetary-tsunami-wave-propagation.png]] **Figure 11: サンフランシスコ湾岸への津波波浪到達の局所ズーム** ![[_attachments/arxiv-2511.11542/fig11-tsunami-sf-bay-impact.png]] ## 考察 - **時間領域における強スケーリングの獲得**: 従来の空間弱スケーリング一辺倒から脱却し、時間発展方向（強スケーリング）を毎秒160万ステップにまで加速できたことで、超長期物理現象の解析、不確実性定量化（UQ）、設計最適化、実時間デジタルツインが現実的な時間枠で可能になる。 - **都市間WANをまたぐ超分散エクサスケールシミュレーションの可能性**: 領域平行移動アルゴリズムの遅延隠蔽能力がメモリ量の2乗に比例する（式1: $n^2 > 2p\tau/c$）ため、クラスタ全体のメモリを集約すればミリ秒単位の広域ネットワーク遅延をも克服できる。地理的に離れた複数のエクサスケールスパコンを接続した協調シミュレーションの道を開く。 - **次世代数値気象予報・地球システムモデルへの波及**: 現代の大気・海洋大循環モデル（CESM, E3SM, FV3, MPASなど）は多層浅水方程式（Stacked SWE）を中核ボトルネックとしている。本成果は、地球シミュレーションのスループットを1桁、エネルギー効率を1.5桁向上させる基盤となる。 ## 出典 - [[.raw/papers/arxiv-2511.11542.pdf]] - [[.raw/papers/arxiv-2511.11542.txt]]