17:15 〜 19:15
[SSS09-P02] CPUとGPUの双方で動作する高速並列地震波動伝播計算コードの開発
キーワード：地震波伝播、数値シミュレーション、GPGPU、並列計算、差分法
科学技術計算の多くの分野で，Graphics Processing Unit (GPU) による計算（General-Purpose computing on Graphics Processing Unit; GPGPU）が活用されつつある．しかし従来のGPGPUでは，C/C++やFortranの代わりに専用言語である Compute Unified Device Architecture (CUDA)を用いる必要があり，その実装はやや煩雑であった．しかも，CUDAを用いたプログラムはGPU環境でしか動かすことができなかった．現代でもCPUを中心とした大規模スーパーコンピュータや計算機クラスタも数多く活用されているなか，CPUおよびGPU双方で動作する大規模数値シミュレーションコードを維持することは難しい問題であった．しかし近年，CUDAを使わずとも従来のコードに数行のOpenACC指示子を付与，あるいは最新の並列性を担保する言語規格に基づいた命令を使うことで，CPU環境での動作はそのままにGPU環境ではGPGPUの恩恵を受けることが可能となりつつある．本研究は，これまで我々がCPU環境において開発してきた地震波動伝播コードのアルゴリズムをGPUにも対応させ，GPUとCPUの両方の環境で動作する地震波動伝播コードを開発する．
本研究で検討対象とするのは，地震波動伝播のコミュニティ・コードOpenSWPC (Maeda et al., 2017)から抽出したアルゴリズムである．このコードは粘弾性体の運動方程式と構成関係式についてメモリ変数法を用いて多変数連立偏微分方程式として扱い，それを時間空間それぞれ等間隔グリッドの並列差分法で解くものである．本研究では，OpenSWPC全体ではなく，その主要な演算カーネル部分を抽出したものを利用し，書き換えやGPU化による効率化を評価した．
本研究では，CPUで動作していたコードをGPUに対応させる手法として，以下の3種類を検討した．1つめはFortran言語標準の並列化機能を用いるものである（方法A）．この方法では，Fortranにおけるdo文による多重ループをFortran2008の新しい命令である do concurrent 命令に置き換える．2つめは最低限のOpenACC指示子を用い，GPUで計算する領域のみを指定する方法である（方法B）． 3つめは詳細なOpenACC指示子の指定により，CPUとGPUの間のデータ伝送まで含めてすべて明示的に記載する方法である（方法C）．
3種類の方法に基づく計算コードをそれぞれ作成し，CPUおよびGPU環境で動作確認と速度に関する検討を行った．まず単一ノードにおける実行において，方法Aは最も単純であるがCPUで顕著な速度低下が見られた．一方，単一ノードでMPI通信なしで実施する場合には，方法BがCPUコードを変更しないままGPUでの十分な加速の恩恵が受けられることがわかった．方法Cは最善のチューニングにより方法Bをさらに上回る速度が出るが，データ伝送をもれなくOpenACC指示子として記載しなければならず，不十分な記載では方法Bよりも低速となる場合もあった．次に並列計算の場合には，方法AおよびBでは一旦CPU側にデータを書き戻してから通信を行う必要がある．これは，方法AおよびBのメモリ管理法ではGPU間の直接伝送機能（GPU-aware MPI）を使えないという制約によるもので，ノード間並列計算のボトルネックとなってしまう．そのため，本研究で対象としたコードでは，最終的に方法Cが最適であることがわかった．ただしその場合でも，まず方法BにおいてGPUの単一ノードでどの程度加速するかを見極めたうえで，その性能を下回らないようOpenACC指示子を追加することが開発プロセスとして肝要であった．
開発した計算コードの評価には東京大学情報基盤センターのMiyabi スーパーコンピュータシステムを用いた．Miyabiは2025年1月に稼働を開始し，汎用CPUノード（Intel Xeon Max 9480）とGPU（NVIDIA GH200 Superchip）の演算加速ノードとを持つ．本研究で開発した方法Cによるコードは，同数のCPUとGPUを利用した場合で，GPUのほうが8〜10倍程度高速に計算を実行できた．高並列化における性能向上はCPUでも十分に見られたが，GPUはそれを更に上回り，256GPUまで十分な高速化を実現した．この知見は近日中にコミュニティ・コードOpenSWPCに反映・公開することを計画しており，PCによるプロトタイピングからCPUクラスタ，そしてGPUスパコンまで一貫して利用できる，使いやすい数値シミュレーションコードが実現すると期待される．また，本研究で開発したCPU-GPU双方で動作するカーネルコード例も，講演当日にGitHubを通じて公開予定である．本研究で得られたGPU化の知見は地震学におけるさまざまな解析・計算コードの高速化に対しても有効であると期待される．
本研究で検討対象とするのは，地震波動伝播のコミュニティ・コードOpenSWPC (Maeda et al., 2017)から抽出したアルゴリズムである．このコードは粘弾性体の運動方程式と構成関係式についてメモリ変数法を用いて多変数連立偏微分方程式として扱い，それを時間空間それぞれ等間隔グリッドの並列差分法で解くものである．本研究では，OpenSWPC全体ではなく，その主要な演算カーネル部分を抽出したものを利用し，書き換えやGPU化による効率化を評価した．
本研究では，CPUで動作していたコードをGPUに対応させる手法として，以下の3種類を検討した．1つめはFortran言語標準の並列化機能を用いるものである（方法A）．この方法では，Fortranにおけるdo文による多重ループをFortran2008の新しい命令である do concurrent 命令に置き換える．2つめは最低限のOpenACC指示子を用い，GPUで計算する領域のみを指定する方法である（方法B）． 3つめは詳細なOpenACC指示子の指定により，CPUとGPUの間のデータ伝送まで含めてすべて明示的に記載する方法である（方法C）．
3種類の方法に基づく計算コードをそれぞれ作成し，CPUおよびGPU環境で動作確認と速度に関する検討を行った．まず単一ノードにおける実行において，方法Aは最も単純であるがCPUで顕著な速度低下が見られた．一方，単一ノードでMPI通信なしで実施する場合には，方法BがCPUコードを変更しないままGPUでの十分な加速の恩恵が受けられることがわかった．方法Cは最善のチューニングにより方法Bをさらに上回る速度が出るが，データ伝送をもれなくOpenACC指示子として記載しなければならず，不十分な記載では方法Bよりも低速となる場合もあった．次に並列計算の場合には，方法AおよびBでは一旦CPU側にデータを書き戻してから通信を行う必要がある．これは，方法AおよびBのメモリ管理法ではGPU間の直接伝送機能（GPU-aware MPI）を使えないという制約によるもので，ノード間並列計算のボトルネックとなってしまう．そのため，本研究で対象としたコードでは，最終的に方法Cが最適であることがわかった．ただしその場合でも，まず方法BにおいてGPUの単一ノードでどの程度加速するかを見極めたうえで，その性能を下回らないようOpenACC指示子を追加することが開発プロセスとして肝要であった．
開発した計算コードの評価には東京大学情報基盤センターのMiyabi スーパーコンピュータシステムを用いた．Miyabiは2025年1月に稼働を開始し，汎用CPUノード（Intel Xeon Max 9480）とGPU（NVIDIA GH200 Superchip）の演算加速ノードとを持つ．本研究で開発した方法Cによるコードは，同数のCPUとGPUを利用した場合で，GPUのほうが8〜10倍程度高速に計算を実行できた．高並列化における性能向上はCPUでも十分に見られたが，GPUはそれを更に上回り，256GPUまで十分な高速化を実現した．この知見は近日中にコミュニティ・コードOpenSWPCに反映・公開することを計画しており，PCによるプロトタイピングからCPUクラスタ，そしてGPUスパコンまで一貫して利用できる，使いやすい数値シミュレーションコードが実現すると期待される．また，本研究で開発したCPU-GPU双方で動作するカーネルコード例も，講演当日にGitHubを通じて公開予定である．本研究で得られたGPU化の知見は地震学におけるさまざまな解析・計算コードの高速化に対しても有効であると期待される．