先端的計算科学研究プロジェクト（ベンチマーク課題）成果報告書

非構造高次精度解析法を用いた圧縮性数値流体ソルバの大規模並列 GPU による高速化
多湖和馬 **

芳賀臣紀 *

堤誠司 *高木亮治 *

宇宙航空研究開発機構 第三研究ユニット

1 研究目的と研究成果
数値流体力学 (CFD) において，近年，不連続ガレルキン法（DG 法）に代表される非構造格子上で高い空間解像度を有する手法の研究が進展し，複雑形状周りで効率の良い Large Eddy Simulation (LES) が実現しつつある．これらの非構造高次精度スキームは，空間を分割する計算要素（セル）内に高次精度化の為の自由度を持つ為，同一のセル数であっても，次数をあげることで総自由度を増やすこと（高解像度化）が可能であることから，大規模解析を行う際に格子生成がボトルネックとならない．広く活用されている有限差分法や有限体積法を用いた流体計算は，メモリアクセス量に対して演算数が少なく，メモリ転送ネックになりやすい．
また，領域分割型並列化を行う際，高次精度化することによって MPI 通信するセル数 (ステンシル数) が増加する．一方，DG 法はセル内に自由度を持っていることからデータの局所性が高く，キャッシュの高速な転送速度の恩恵を得られる．また，領域分割間の通信はセル境界上の値だけであり，データ通信量が高次精度化の次数によらないことから，差分法と比較して大規模並列性能が期待できる．これらの利点に着目し，JAXA では DG 法と同等な精度を持ち，かつ，定式化が簡潔な Flux Reconstruction(FR) 法 [1] を使った CFD ソル
バ LS-FLOW-HO を開発してきた．

FR 法を利用した CFD ソルバとして有名な PyFR は，195 B 自由度を用いた低圧タービンの圧縮性 DNS 解析にて，18,000 台の NVIDIA Tesla K20X を利用し 13.7 PFlops の倍精度演算性能（理論性能に対して 58%の実行効率）を記録している．[2] これほど高い実行効率が得られる理由として，上述した FR 法のメリットに加え，最適化された BLAS ライブラリの DGEMM をうまく利用した実装を行っていることが挙げられる．
しかし，六面体要素に限った場合，行列の疎性から高次精度ではゼロ要素のみ計算する実装の方が高速なことを著者らは富士通 FX100 を用いて確認してきた [3]．そして，FX100 の 10 万コアまでを使用したストロングスケーリングなどを行い，本ソルバの高速化研究をこれまで行ってきた．
本研究では, GPGPU を用いて JAXA で開発した FR 法ソルバの高速化と大規模並列性能の評価を行い，LESにおける GPU の有用性を検証する. ソルバの GPU 化にはディレクティブベースでオリジナルのコードからの移植が容易な OpenACC を用いて実装した. これによる GPU と CPU との実行時間の比較を行った．さらに，ITO-B の 508 GPU を用いて並列スケーリング実行を行った．

1

2 Flux Reconstruction 法
Flux Reconstruction 法の概要を説明する．詳細は [1] を参照されたい．FR 法はヘキサ，テトラなどの一般的な非構造セルで使用可能であるが，本研究では六面体セルを対象として離散化を行う．空間の高次精度化のために，各要素には解の定義点として SolutionPoint(SP) を導入する．SP 値は支配方程式を解くことで保存，更新される自由度 (DOFs) となる．i 番目のセル内の保存量 Q は，各座標方向の離散的な K + 1 個の SP
K
値を用いて，K 次のラグランジュ多項式補間のテンソル積で近似される．Qi の多項式近似は QKi := I [Qi ]として与えられ，次のように与えられる．
K
QK
i := I [Qi ] =K
∑

K
K
K
QSP
i;j,k,l ϕj (ξ)ϕk (η)ϕl (ζ),(1)

j,k,l=0

ここで，i, j, k は各座標における SP のインデックスであり，ϕKj は K 次のラグランジュ多項式である．FR法は K の多項式を用いて，K + 1 次の空間精度を達成する．ここでは，次数 K の解多項式を用いた FR 法をpK スキームと呼ぶことにする．セル境界の数値 Flux を風上法などで計算することで，セル内の Flux 分布を修正する．ξ 方向を考えると，ξ 方向の非粘性 Flux Ê の微分は以下で計算される．

∂ Êi∂ Eˆiccom′com′=+ [Êi−1 − Êi (−1)]gL (ξ) + [Êi+ 12 − Êi (1)]gR (ξ)
2
∂ξ
∂ξ

(2)

時間積分には，3 段 3 次の TVD Runge-Kutta 法を用いる．

3 成果の詳細
3.1 OpenACC を用いた FR 法ソルバの GPU 化
GPU 化には今まで JAXA で開発を行ってきた非構造高次精度ソルバ LS-FLOW-HO を用いる．この ソ ル バ ー は Fortran で 実 装 さ れ て お り ，Fujitsu FX100 用 に チ ュ ー ニ ン グ さ れ て い る ．GPU 化 にはオリジナルコードの改変を最小現にするために OpenACC を用いた．コンパイラには PGI FortranCompiler(ver.19.4) を利用した．本研究ではプログラムの実行時にほとんどの時間を占める時間積分の処理を OpenACC を用いて GPU 上で実行可能にした．この時，ホストからデバイスへのメモリ転送は時間積分に入る前に行い，時間積分の計算中は結果の書き出し以外ではホストとデバイス間のデータ通信は極力行わないように実装を行う．プロセス間通信には，PGI コンパイラに同伴の OpenMPI 3.1.3 を使用
し，GPU 間データ転送には CUDA-Aware MPI を使用する．PGI コンパイラのコンパイルオプション
は -Mpreprocess -Minfo=all,ccff -Mneginfo -mcmodel=medium-lblas -acc -ta=tesla:lineinfoを使用する．

3.2 CPU と GPU の性能比較
CPU と GPU で実行時間を比較した結果を図 1 に，比較に使用した各プロセッサの緒言を表 1 に示す．
計算には p7 スキームを使用し，総自由度 5,451,776 を用いて Taler Green Voltex モデルを計算する．GPUの実行時間の比較として Intel Xeon E5-2698v3 (Haswell) と Intel Xeon Gold 6154 (Skylake) を使用する．

Haswell と Skylake ではそれぞれ 16，18 スレッド並列による実行を行った．NVIDIA Tesla P100 による実2

行時間のスピードアップは Haswell と比較して 7.3 倍，Skylake と比較して 4.3 倍となった．各プロセッサの理論 DP FLOPS 値を考えると，P100 は Haswell に対しては 9 倍, Skylake に対しては 4 倍となっている．

Skylake と P100 の理論 FLOPS は約 4 倍ほどなので, Skylake に対しては FLOPS 値の増加分速くなっており理想的である．

図 1 Haswell を基準とした実行時間の比較

Intel Xeon E5-2698v3(Haswell)

Skylake(Intel Xeon Gold 6140)

NVIDIA Tesla P100(Pascal)

コア数

1618

-

周波数 GHz (turbo)

2.3(3.6)

2.3(3.7)-理論 DP GFLOPS

588.8

1324.85304.3

メモリ帯域 GB/s

68.0

127.9

732.0

B/F

0.11

0.090.13

表1

各プロセッサの性能緒言

3.3 小規模ストロングスケーリング性能
GPU の大規模並列化を行うにあたり，小規模で GPU の並列性能の確認を行った．小規模では，同一ノード内の 1 つの GPU から 4 GPU までのストロングスケーリング, 及びノード 1 GPU から 4 ノード 4 GPUまでのストロングスケーリング (ノード間並列) の 2 パターンを計測した．ITO-B は 1 ノードに 2 ソケットあり，18 コア CPU を搭載しているので，実行条件として両ケースともホストでは 1 プロセス 9 スレッドを割り当てる．実行時オプションとして --mca btl_openib_want_cuda_gdr 1 を指定して GPU Directを有効にし，ホストへのメモリコピーを抑制した．計算にはセル数 12,950 の ANSS 5th Aerodynamics

Prediction Challenge(APC-5) の高揚力装置翼型の LES ケースを使用し [4]，p7 スキームで計算を行った．
ノード間比較のプロセス割付は mpirun オプションの--map-by ppr:1:node, ノード内でのプロセス割付は--mca-by ppr:2:socket(4GPU プロセス並列の場合) などで設定している．計測結果を表 2 に示す．計
測値は 10 回実行した値の平均値を使用してる．

3

ノード内 GPU 並列

ノード間 GPU 並列# of GPU

12

3

4

1

2

3

4

elapse time

3.94E+12.08E+01

1.48E+01

1.20E+01

3.94E+01

2.03E+01

1.43E+01

1.11E+01

speedup rate

1.0

1.9

2.7

3.3

1.0

1.9

2.83.6

表2

ノード内の GPU を用いた並列性能とノード間の GPU を用いた並列性能

ITO-B の GPU-GPU 間接続は NVLink により 20 GB/sec × 1 or 2，ノード間接続は InﬁniBand EDR 4xにより 12.5 GB/sec × 2 となっており，理論性能から見ればノード内 GPU 並列のほうが速いと予想されたが，計測結果ではノード間 GPU 並列のほうが速い結果となった．4 GPU 並列のとき，mpi wait を呼び出してデータパッキングを行うルーチンの実行時間の合計は，ノード内並列実行では 2.87 s(実行時間に占める割合は 23.9 %)，ノード間並列実行では 2.94 s(実行時間に占める割合は 26.5 %) となっており，通信の待ち時間の絶対値はノード間並列のほうが若干だが大きい．次にノード間並列とノード内並列で GPU 処理の実行時間が異なる可能性が考えられるので，4 GPU 並列実行時の GPU 実行時間の内訳を NVIDIA のプロファイラツールを用いて取得した．結果を図 2 に示す．within node はノード内並列，inter-node はノード間並列を表し，各プロセスごとに GPU 処理を行った箇所の実行時間の積である．各処理でノード間，ノード内でのGPU の実行時間に大きな違いは見られなかった．

図 2 GPU 区間の実行時間内訳

このことから，ホスト側の処理で実行時間の差が発生していると予測される．表 2 の実行時間は時間積分区間を計測した値になっており，ホスト側の処理時間も含まれている．OpenACC の特徴としてループ区間を対象に GPU 化を行うことから，ループ処理以外の制御はホスト側で行うことになる．ホスト側の処理として差が現れそうな処理としてはメモリの確保などが考えられるが，本調査では原因箇所の特定にいたれなかった．

3.4 大規模ストロングスケーリングセ ル 数 を 一 定 に し た 計 算 格 子 を 用 い て ，使 用 す る GPU 数 を 増 加 さ せ た と き の ス ト ロ ン グ ス ケ ーリ ン グ の 計 測 を 行 う ．計 算 ケ ー ス は ANSS 5th Aerodynamics Prediction Challenge(APC-5) の 高
揚 力 装 置 翼 型 の LES ケ ー ス を 使 用 し た [4]．実 行 し た 4 つ の ケ ー ス に 関 し て 表 3 に 示 す ．GPU 並

4

列 数 は 16, 64, 256, 508 と 増 加 さ せ た ．こ の と き ，領 域 分 割 に は グ ラ フ 分 割 ツ ー ル の metis を 使 用し ，領 域 が 均 等 に な る よ う に 分 割 を 行 っ た ．case1 と case3 は ，最 小 の 16 GPU 並 列 の と き に デ バ
イ ス メ モ リ の 容 量 の 限 界 を 想 定 し て セ ル 数 を 設 定 し て い る ．並 列 実 行 コ マ ン ド の オ プ シ ョ ン と し て--mca btl_openib_want_cuda_gdr 1 を指定し，1 ソケットに 2 プロセス--map-by ppr:2:socketを割り当て，GPU Direct RDMA を有効にした．各計測は case1, 4 は 3 回, case2, 3 は 5 回行い，その平均値を使用した計測結果を図 3 に示す．16 GPU での実行時間に対して 508 並列の case1 は 14.3 倍，case2 は5.5 倍，case3 は 10.1 倍，case4 は 3.9 倍であった．p7 を使用した case1 と case2, p4 を使用した case3 とcase4 を比較すると，16 GPU 並列時点でセル数が大きいケースのほうがストロングスケーリングが良い結果を示している．
並列実行での通信コストを評価するため，セル数が同じである case1(p7) と case4(p4) の mpi wait の時間を計測し，実行時間に占める割合 (各プロセスの平均) を表 4 に示す．16 並列では case4(p4) と比べて case1(p7)のほうが mpi wait の占める割合が大きいが，64, 256, 508 と並列数が増加すると case4(p4) のほうが通信コストの割合が大きくなっている．FR 法では，p スキームの次数が大きいほどセル内の SP 数が多くなり，演算数は増加する．さらに 1 セル当たりのデータ量も多くなることでプロセス間で通信されるセル境界のデータ量も多くなる．演算数が多いとき，演算による通信時間の隠蔽が行われやすいと予想される．case1(p7) はcase4(p4) と比較して演算数が多く，mpi wait の待ち時間の隠蔽に対して有効だと考えられる．case4(p4) は
case1(p7) と比較して，プロセスごとのデータ通信量が少ない利点はあるが，並列数が増えると演算リッチなcase1(p7) のほうが実行時間に占める通信待ち時間の割合は下がることがわかった．
case1

case2

case3

case4

スキーム

p7

p7

p4

p4

セル数

194,250

97,125

796,425

194,250

DOFs

99,456,000

49,728,000

99,553,125

24,281,250

表3

ストロングスケーリングで使用した格子数と自由度

図3

ストロングスケーリングの計測結果

case1 の 16 GPU 並列と 508 GPU 並列の 1 プロセスにおける実行時間の内訳を図 4, 5 に示す．値は nvprof を用いて取得した．2 つの図から 16 並列と 508 並列ではコストが大きい処理の順位に変化が5# of GPU

1664256508

p7 case1(%)

25.3

28.2

33.0

33.2

p4 case4(%)

21.2

29.2

37.2

38.4

表 4 実行時間に占める mpi wait の実行時間の割合あることがわかる．16 GPU ではコスト上位に来ていなかった isrncomm ho face grad2 openacc, isrn-comm ho face2 openacc, computeprimitivevalues openacc の 3 つが 508 並列では上位に来ている．isrncomm のルーチンは MPI 送受信のパッキング処理を行っており，ループに reduction を用いる処理となっている．この 3 つのルーチンが 508 GPU のときにボトルネックとして現れているので，チューニングを行うことでスケーリングが改善する可能性が考えられる．

図 4 16 GPU 並列実行時の 1 プロセスにおける GPU 実行時間の内訳

図 5 508 GPU 並列実行時の 1 プロセスにおける GPU 実行時間の内訳

6

4 まとめ本研究では，FR 法ソルバーの GPU における性能評価を行った．まず，CPU と GPU での実行時間の比較を行った．この結果から，GPU 単体で，FLOPS 分の性能向上を確認でき, Haswell に対して 7.3 倍，Skylakeに対して 4.3 倍ほどの高速化であった．次に大規模 GPU を用いたストロングスケーリングの計測を行い，GPU 数 31.75 倍に対して最大で 14.3 倍のスケーリングを確認できた．また，ストロングスケーリングの疎外要因として，特定の GPU ルーチンの実行時間のスケーリングが悪い点が考えられる．今後はこのルーチンのチューニングと通信によるスケーリングの影響を他のスパコンシステムと比較・評価を行いたい．

参考文献[1] Hung T Huynh. A ﬂux reconstruction approach to high-order schemes including discontinuousgalerkin methods. In 18th AIAA Computational Fluid Dynamics Conference, p. 4079, 2007.
[2] Peter Vincent, Freddie Witherden, Brian Vermeire, Jin Seok Park, and Arvind Iyer. Towards green
aviation with python at petascale. In SC’16: Proceedings of the International Conference for High
Performance Computing, Networking, Storage and Analysis, pp. 1–11. IEEE, 2016.
[3] Kazuma Tago, Takanori Haga, Seiji Tsutsumi, and Ryoji Takaki. Performance improvement of ascalable high-order compressible ﬂow solver on unstructured hexahedral grids. In Proceedings of theInternational Conference on High Performance Computing in Asia-Paciﬁc Region, HPCAsia2020, p.
227–236, New York, NY, USA, 2020. Association for Computing Machinery.
[4] The Japan Society for Aeronautical and Space Sciences Aerodynamics Division Committee. Fifth
aerodynamics prediction challenge(apc-v). https://cfdws.chofu.jaxa.jp/apc/apc5/.

7