Mercury-Opal: the GPU-accelerated version of the n-body code for planet formation Mercury-Arxes

技術要約：Mercury-Opal

問題提起
完全に相互作用する天体を含む大規模な$n$体シミュレーションは、計算惑星形成研究における重要なフロンティアである。汎用インテグレータであるMercury (Chambers, 1999) およびその特化した並列拡張版であるMercury-Arχes (Turrini et al., 2019, 2021) は、惑星形成の詳細なモデリングを可能にしてきたが、これらは計算力学アルゴリズムと利用可能なハードウェアリソースによる制約を受けている。現在のCPUベースの並列化（例：OpenMP経由）では、$\sim 10^5$個の天体を$\sim 10^7$年間にわたってシミュレーションすることが可能であるが、$\sim 10^5$個の天体による完全な相互作用セットを管理することは、$O(N^2)$のスケーリングのため依然として困難である。GPU加速コンピューティングへの移行は、計算密度の向上とエネルギー効率の向上をもたらすが、GPU内部RAMの制限やホスト・デバイス間の帯域幅を克服するための、注意深いデータ管理とコードの再構築を必要とする。

手法
著者らは、OpenACC 3.3仕様を用いて開発された、Mercury-ArχesのGPU加速版であるMercury-Opalを提示する。OpenACCの選択は、パフォーマンスの向上を維持しつつ、大幅なコードの再構築を最小限に抑え、クロスプラットフォームのサポートを確保する必要性に基づいている。

主な手法ステップは以下の通りである：

1. プロファイリング： 大質量天体の数（$N_{big}$）を変化させながら、Mercury-Arχesのシリアル実行をプロファイリングした。分析の結果、mdt_hyサブルーチン（重力および非重力力の処理）が主要な計算ボトルネックであり、大きな$N$に対して実行時間の最大68%を消費していることが特定された。
2. GPU移植戦略：
   • 並列化： mdt_hyルーチンをGPUに移植し、力計算（重力およびガス抗力）、ケプラー・ドリフト、および衝突チェックを包含させた。
   • Bulirsch-Stoerの除外： Bulirsch-Stoer積分による近接遭遇を処理するmdt_hkceルーチンは、CPU（ホスト）側に保持した。この決定は、近接遭遇が少数の不規則な天体サブセット（<1%）に関与するため、ロードインバランスと並列化の利得の制限が生じること、およびデータ転送量を低く抑えることを考慮して行われた。
   • コードの再構築： OpenACCのスタックサイズ制限（512 kB）に対応するため、大きなルーチン（mfo_drct, arxes_gas, mce_snif）を解体し、mdt_hy内に直接移動させた。
   • 言語変換： 統合されたWHFASTライブラリは、Fortran-Cインターフェースにおけるコンパイルエラーを解決し、シリアル実行効率をわずかに向上させるために、C99からFortranへと変換された。
3. データ管理： mdt_hyを包含するデータ領域を定義した。スカラーの円盤パラメータは一度だけコピーされる。位置、速度、および質量の配列は、領域の開始時と終了時に転送される。GPU上で近接遭遇が検出された場合、相互作用する天体のインデックスを含む特定の配列がCPUベースの積分のためにホストへ転送され、更新された位置が送り返される。
4. テスト環境： 消費者向けハードウェア（Intel Core i5-10200H CPUおよびNvidia GeForce RTX 3060 mobile GPU）を用いて、CPUに有利な条件（具体的には、GPUの低い倍精度性能）の下でGPUコードのストレス・テストを実施した。

主な貢献

• Mercury-Opal： 完全に相互作用する天体を用いた惑星形成のシミュレーションを可能にする、Mercury-Arχesコードの初のGPU加速実装である。
• OpenACCの実装： 元のコードの物理的レシピ（ディスクと惑星の相互作用や空気力学的抗力を含む）を維持しながら、最小限の再構築で複雑なマルチフィジックス$n$体コードを効果的に移植できることを実証した。
• ハイブリッド実行モデル： $O(N^2)$の力計算の大部分はGPUで実行し、近接遭遇積分の複雑で不規則なロジックは、データ転送のボトルネックとロードインバランスを回避するためにCPUに残すというハイブリッドアプローチを確立した。

結果

• パフォーマンス：
  • 天体数が少ない場合（$N \lesssim 100$）、カーネル起動のオーバーヘッドとデータ転送コストにより、GPU実行はシリアルCPU実行よりも遅くなる。
  • $N=1000$において、GPU並列実行はシリアルCPU実行に対して**2.2倍の加速（$S_w$）**を達成した。
  • 1体の場合に正規化すると（相対的な加速 $S_1$）、$N=1000$においてGPUは16倍高速である。
  • Nvidia Nsight Systemsによるプロファイリングでは、大きな$N$においてデータ転送操作が総ウォールタイムの1%未満であり、主な効率低下はCPUバウンドのmdt_hkceルーチンに起因することが示されている。
• 精度と安定性：
  • $10^5$年にわたるエネルギー保存テストでは、シリアル実行とGPU実行の両方において相対的なエネルギー振動が$\sim 10^{-13}$以内であり、オリジナルのMercury-Arχesの性能と一致している。
  • $N=10$および$N=100$の系における軌道要素（離心率および軌道長半径）の進化は、シリアル実行とGPU実行の間で無視できる程度の偏差しか示さず、GPUポートが重大な数値誤差を導入しないことを確認した。

意義
本論文は、Mercury-Opalが、特にESAのArielミッションの準備を支援する、物理的に現実的かつ効率的なインテグレータを惑星科学コミュニティに提供すると主張している。GPU加速を活用することで、本コードは、CPUのみの並列化では以前は不可能であった、より多くの数の完全に相互作用する天体のモデリングを可能にし、$O(N^2)$の相互作用という計算上のボトルネックに対処する。著者らは、OpenACCの使用により、低レベルのGPUプログラミングに関する深い専門知識を必要とせずに、コードのアクセシビリティとポータビリティが維持され、より広いコミュニティへの採用が容易になることを強調している。本研究は、倍精度性能が最適でない消費者向けハードウェアであっても、GPU加速が大規模な$n$体シミュレーションにおいて具体的な利点を提供することを検証している。