技術概要：非保存項を有するエントロピー安定性 discontinuous Galerkin 法によるオイラー方程式ソルバの GPU パフォーマンス

問題提起

大気モデルのための正確で安定かつ効率的なダイナミカルコアの構築は、依然として重要な課題である。discontinuous Galerkin（DG）法は高次精度と優れた並列化の可能性を提供するが、地球流体力学への応用には、高次における安定性と計算効率という 2 つの主要な障壁が存在する。エントロピー安定性 DG（ESDG）法は、エントロピーを保存または散逸させる数値フラックスを構築することで安定性問題を解決するが、2 点フラックス評価の必要性により計算コストが高くなる。さらに、浮力を圧縮性オイラー方程式に組み込むと、保存則の集合から非保存項を有する平衡則の系へと変換され、フラックスの定式化が複雑化する。

オイラー方程式および浅水方程式に対する GPU 実装に関する文献は豊富であるが、特に浮力を伴うオイラー方程式に対する ESDG 法の GPU 最適化に関する公表された研究には大きなギャップがある。高い計算コストに対処しなければ、ESDG 法は大規模な大気モデルの実用化において非現実的となるリスクがある。本論文は、現代のグラフィックス処理ユニット（GPU）アーキテクチャ上で効率的な ESDG 法を構築する必要性に応えるものである。

手法

著者らは、NVIDIA A100 GPU 上で浮力を伴う圧縮性オイラー方程式に対するエントロピー安定性 discontinuous Galerkin ソルバを実装した。数値的手法は、以前の研究 [55] で詳述された ESDG 枠組みに従い、GPU 実行向けに適応されている。

支配方程式と離散化

ソルバは、地ポテンシャル（浮力）およびコリオリ項を含む時間依存の圧縮性オイラー方程式を処理する。計算領域は、重なりを持たない六面体要素に離散化され、これらは $(N+1)^3$ 個のルジャンドル・ガウス・ロバット（LGL）ノードを含む参照要素に写像される。半離散形式は、メモリ読み込みを最小化するためにテンソル積構造を利用するフリーマトリックスアプローチを採用している。エントロピー安定性は以下の手段で強制される：

1. 体積項： エントロピー保存型の 2 点フラックスを使用。
2. 表面項： 行列散逸（ローのリーマンソルバに類似）を付加したエントロピー保存型フラックスを使用。
3. 時間積分： 低記憶型の 4 次ルンゲ・クッタ法を使用。

GPU 最適化戦略

本論文の核心は、2 点フラックスおよび非保存項の計算オーバーヘッドを軽減するための一連の最適化にある：

1. スレッドマッピング： 各要素がスレッドブロックに割り当てられ、スレッドは自由度（DOF）を表す。
2. 演算削減：
   • 事前計算： 原始変数（例えば運動量からの速度）および対数項は、各スレッドごとに 1 回のみ事前計算され、各フラックス評価ごとに再計算されない。
   • 対数平均： 著者らは、除法を削減し、$\log(a)$ と $\log(b)$ を個別に計算して評価間で再利用可能にするために、ホーナー法を用いた効率的な定式化 [42] を採用する。
3. 対称性の利用：
   • 標準的なオイラー方程式では、フラックス対の対称性により計算量が半分になる。
   • 非対称な重力処理： 浮力項は非対称な成分を導入する。著者らは運動量フラックスを対称部分と重力寄与に分解する。重力項は対ごとに 1 回計算され、パートナーノードに対してスケーリングされることで、ほぼ 2 倍の高速化を維持する。
4. 負荷分散： スレッドがフラックス対（$j > i$）を計算する際の負荷不均衡を防ぐため、特定のインデックス付け戦略を採用する。偶数次多項式の場合、スレッドは均等な作業を行う。奇数次の場合、重み付けアプローチによりループ反復を均一化し、同期の問題を回避する。
5. メモリとカーネル融合： 状態変数および事前計算データは、高速アクセスのために共有メモリに格納される。3 つの空間方向は、冗長な共有メモリ読み込みを排除するために単一のカーネル起動に融合される。共変フラックスは、中間行列を形成・縮約することを避けるために直接計算される。

主要な結果

単一 GPU 上のパフォーマンス

• 計算効率： NVIDIA A100 ハードウェアにおいて、最適化された体積カーネルは、多項式次数 $N \ge 3$ に対して 64 ビット浮動小数点ピークパフォーマンス（ルーフライン天井）のほぼ**70%**を達成する。カーネルはメモリーバウンドではなく計算バウンドであるため、GPU 加速に非常に適している。
• 最適化の影響： 事前計算、効率的な対数平均、対称性の利用といった累積的な最適化により、ベースライン GPU カーネルに対して7.3 倍の総高速化が達成された。
• カーネルのボトルネック： 最適化前には体積カーネルが実行時間を支配していた。最適化後、表面カーネルが主要なコストとなり、今後の最適化努力はメモリーバウンドの表面積分を対象とすべきであることを示している。
• スケーリング： ソルバは強力な飽和能力を示す。高次多項式（$N=6$）の場合、わずか 4,000 個の要素で GPU は完全に飽和する。

複数 GPU へのスケーリング

• 強スケーリング： 上昇熱気泡テストケース（$N=4, L=6$）を使用すると、ソルバは 32 GPU までスケーリング可能である。通信オーバーヘッドにより、32 GPU での並列効率は 60% に低下する。これは、問題サイズがプロセッサ数と一致するときに最適なパフォーマンスが達成されることを示唆している。
• 弱スケーリング： 問題サイズとプロセッサ数が比例して増加する場合、ソルバはステップあたりの時間をほぼ一定に維持し、理想的な弱スケーリング動作を示す。

GPU と CPU の比較

ソルバは、AMD EPYC 7763 Milan ノード上で実行される高度に最適化された CPU コード（NUMA モデル）と比較された。

• 速度： GPU 実装は、同じ問題サイズに対して CPU コードより約10 倍高速である。
• エネルギー効率： GPU 実装は CPU コードよりも13 倍以上エネルギー効率が優れている。
• 精度のトレードオフ：
  • 上昇気泡： 64 ビットから 32 ビット精度に切り替えると、大規模な対称性の最小限の損失で2.16 倍の高速化が得られる。
  • バロクリニック不安定： 問題サイズが小さい場合、計算に対する通信オーバーヘッドの相対的な大きさにより、32 ビット精度で1.5 倍の高速化となった。しかし、問題サイズを増やすことで2.16 倍の高速化が回復した。32 ビット精度の視覚的結果は、64 ビットと同様のダイナミクスを捉えており、ノイズは手法の安定性ではなく、総エネルギー変数の大きさによるものであると帰属された。

意義と主張

本論文は、GPU 最適化に特化したオイラー方程式に対する ESDG に関する最初の研究であると主張している。その主な意義は、採用の障壁となることが多いエントロピー安定性手法の計算コストが、現代の GPU ハードウェア上で効果的に軽減できることを実証した点にある。

主な主張は以下の通り：

1. 実現可能性： 浮力駆動流に対する ESDG 法は GPU 上で実現可能であり、高い演算強度とピークパフォーマンスに近い性能を達成する。
2. 最適化ロードマップ： 本論文は、非保存項の対称性処理、負荷分散、演算削減といった具体的な戦略の詳細なロードマップを提供し、効率を大幅に向上させる。
3. パフォーマンスの向上： 本実装は、最先端の CPU コードに対して 10 倍の高速化と 13 倍以上のエネルギー効率の向上を達成し、高次・エントロピー安定性の大気モデルを実用的にする。
4. 精度の柔軟性： ソルバは、地球物理学的問題の解の質を維持しながら大幅な高速化（最大 2 倍）を実現する 32 ビット精度をサポートしており、大規模シミュレーションにおける実用的なトレードオフを提供する。

著者らは、現在の作業は陽的時間積分に焦点を当てているが、実証された空間離散化のパフォーマンスは、陰的時間ステップ法を取り入れる可能性のある将来の実用的な大気モデルの基盤を築くと結論付けている。