学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

jh170020-NAJ

Optimisation of Fusion Plasma Turbulence
Code toward Post-Petascale Era II
前山伸也（名古屋大学）

概要 本研究課題は、プラズマ乱流シミュレーションコードの高度化とそれを利用した物理解析を行い、核融合プラズマ研究に資することを目的とする。コードの高度化に関して、前年度開発した GPU 版コードを CUDA/Thrust 系へ移植するとともに、通信演算オーバーラップを実装することで、~2 倍の性能向上を得た。また、アルゴリズム拡張として衝突項の陰解法化を行い、演算密度増大による計算効率の向上と、大きな時間刻み幅による求解時間（Time to solution）の短縮の両立を実現した。物理解析に関しては、
近年核融合分野で議論されている局所モデルによる輸送予測が装置周辺領域では過小評価になるという、Shortfall 問題に対して２つの異なる観点から解析した。第一に、磁力線方向への乱流自己相関を緩和した改良局所モデルを適用し、輸送増大が起こることを示した。第二に、大域的モデルを適用し、微視的乱流揺動と大域的新古典輸送の相互作用により輸送増大が起こることを示した。従来の局所モデル解析に含まれていなかったこれらの効果は、Shortfall 問題を解消する糸口になり得ると期待される。

1. 共同研究に関する情報担当する。青木尊之（東京工業大学）と荻野(1) 共同研究を実施した拠点名

正雄（名古屋大学）は、それぞれ TSUBAME

東京工業大学上での GPU への最適化および FX100 上での名古屋大学CPU への最適化支援を担当する。

2. 研究の目的と意義(2) 共同研究分野☑

超大規模数値計算系応用分野

本研究課題は、将来の基幹エネルギー源の

超大規模データ処理系応用分野候補の一つと期待される核融合プラズマの研

超大容量ネットワーク技術分野

究開発に資することを目的とし、シミュレー



超大規模情報システム関連研究分野

ションコードの高度化とそれを利用した物理解析を行う。
計算科学の観点からは、核融合プラズマの乱流解析に用いられる 5 次元ステンシル計算(3) 参加研究者の役割分担
代表者の前山伸也（名古屋大学）は核融合を如何にして最先端のメニーコアアーキテク

プラズマ中の乱流輸送現象に関する物理解析

チャ上で高速実行するかという問題があり、を担当する。副代表者の朝比祐一（フランス異なるアーキテクチャの特性を考慮したマル原子力庁）はシミュレーションコードのマルチプラットフォーム開発を行うことで、問題チプラットフォーム開発およびコードの最適の特性と対策を明らかにする。本年度は特に、化を担当する。課題参加者の Xavier Garbet

次世代のエクサスケール流体コード開発に資と Guillaume Latu（フランス原子力庁）は大することを目的とし、メニーコアアーキテク

域的プラズマ乱流シミュレーションに関した、チャ環境における転置通信の最適化手法につ渡邉智彦（名古屋大学）は局所的プラズマ乱いて調べた。

流シミュレーションに関した物理解析支援を

核融合学の観点からは、異なる特徴を持つ
1

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

コード（大域的コード GYSELA 及び局所的上が得られた。GPU 版 GKV コード向けに開コード GKV）を用いた相互補完的な解析によ発した技術の大部分は FX100 版にも適用可り、閉じ込め性能の予測モデルの高度化を目能であり、FX100 版の高速化も見込まれる。

指す。具体的には、今年度は従来局所的コー

また、本年度新たに FX100 において開発されドで予測が困難とされてきた、熱輸送のたノード間通信フリーの衝突項ソルバーは、

Shortfall 問題（局所的コードの予測は、炉心演算密度が高く GPU に極めて適した実装と部分について輸送レベルの実験値を再現する言える。

ものの、プラズマ周辺領域では実験値を大き

これらは、JHPCN の枠組みで提供されるく下回ることが知られ、問題となっている）

マルチプラットフォームの利点を生かした研に取り組んだ。

究と言える。

これらは計算科学および核融合学が直面している挑戦的課題であり、将来の HPC に対

4. 前年度までに得られた研究成果の概要

する核融合プラズマシミュレーションコードの適応性を明らかにするとともに、核融合エ

前年度は“Multi-platform Development of

ネルギー開発に資する物理的知見を得るといFusion Plasma Turbulence ”（ 課 題 番 号 ：う点から、学術的・社会的意義が高い。

jh160028-NAJ）の研究課題において、まず、プラズマ乱流シミュレーションコード GKVの TSUBAME2.5 への移植および性能評価を3. 当拠点公募型共同研究として実施した意義

行い、通信律速の性能劣化が課題として顕在

JHPCN では、FX100、TSUBAME3.0 など化した。次に、FX100 におけるアシスタント最先端のメニーコアアーキテクチャ利用環境コアを利用した通信・演算オーバーラップ手が複数提供されており、本研究で行うマルチ法の性能評価を行い、チャンクサイズ依存性プラットフォーム開発に最適の環境である。

を明らかにした。また、物理応用として、経マルチプラットフォームでアプリ開発を行い、済性に優れる高圧力プラズマで顕著になるとアーキテクチャの特性を活かした最適化手法考えられるマイクロティアリングモードの線を確立することは、ポストペタスケールでの形・非線形解析を行った。

研究展開の礎となる。また、国際協力に基づ以上の成果を受けて、今年度は GPU 環境

き複数コードで最適化を行うことで、最適化での通信最適化、アルゴリズムの拡張、およ手法ごとの利点・欠点などの比較を行うことび開発したコードを利用した物理解析を展開ができる。

した。

まず、GKV コードと GYSELA コードで実装されている転置通信のパターンを比較し、GYSELA コードで実装されている転置通信

5. 今年度の研究成果の詳細の方が通信コストの観点で効率的なことを明

今年度は、GKV・GYSELA コードの通信最らかにした。新たに開発された GPU 版 GKV適化（5.1, 5.2 節）、衝突項陰解法の高効率実コードでは、転置通信を GYSELA コードの装（5.3, 5.4 節）、物理解析への応用（5.5, 5.6

手法で実装し、通信コストを従来版の 1/3 程節）の 3 つに取り組んだ。

度に削減した。同様に、転置通信部分に関するオーバーラップを K N L 版の GYSELA コ5.1. GKV・GYSELA コードの通信最適化（前期）ードに適用したところ、最大 5%程度の性能向GPU や Xeon Phi KNL を始めとする次世2

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

グスケーリングの性能測定結果である。畳み込み積分において必要となる転置通信に関するスケーラビリティを調べるため、波数空間における並列度を変化させた。TSUBAME3.0におけるオーバーラップの効果は主に、ホストデバイス通信とノード間 MPI 通信のオーバーラップに起因する。TSUBAME3.0 においてはオーバーラップ版、非オーバーラップ版ともに高いスケーラビリティが得られている。FX100 における性能劣化は、転置通信のコストが最適化済みの GPU 版と比べ 3 倍程度であることに加え、ノード数が
図１．FX100 及び TSUBAME3.0 におけるストTSUBAME3.0 と比べ 4 倍であることによるロングスケーリング。性能は MLUPS (Megaと考えられる。

Lattice Update Per Second)で評価。
5.2. GKV・GYSELA コードの通信最適化（後期）
代のアーキテクチャにおける転置通信の最適

前期は、GKV の Cuda/Thrust への移植と化を行った。前年度開発した G P U 版 GKV

パイプライン化による通信演算オーバーラッ

コードを CUDA/Thrust 系へと移植するととプの並列 FFT への適用を行った。後期は、コ

もに、畳み込み積分、移流計算実行時の通信ア数に比べて小メモリであるメニーコア環境と演算のオーバーラップを実装した。高次元に対応するため、並列 FFT の最適化を進めた。

ステンシル計算であることに着目し、演算と具体的には、パイプライン処理において、メ無関係な次元についてパイプライン化を行う

モリ削減と高速処理を両立するために、ポイことで、従来よりも効率的にオーバーラップンタのスワップを利用して、メモリコピーをできることが明らかとなった。これに加えて行わずに作業配列を再利用するように改良を転置通信を最適化したことで、TSUBAME3.0

施した。これにより、作業配列のサイズをバ

では、FX100 と比べ 2 倍程度の性能向上に成ッチ数分の 1 へと削減でき、処理速度を損な功した。約 5 倍のプロセッサの演算性能差うことなくメモリの大幅な削減が実現された。

（ TSUBAME3.0 ～ 5TFlops, FX100 ～さらに、GYSELA コードにも同種の通信演

1TFlops）に対し低い性能向上率となってい算オーバーラップを実装した。この際、独立るのは、ノード間通信の性能差が 2 倍である次元方向へのパイプライン化処理だけでなく、ことによると考えられる（TSUBAME3.0～

タスクレベルでの通信演算オーバーラップも4x12.5GB/s,FX100～2x12.5GB/s）。特に試した。後者では、通信用コアが MPI 通信を

GPU 版においては、全体のコストのおよそ行っている間に、その他の演算用コアが独立50%がノード間通信コストとなっているため、な演算タスクを実行する。今回は簡易に手法

ノード間通信の性能がコード全体の性能を大の適用性を確認するために、通信・演算タスきく左右する。

クの組み合わせはコード中で静的に記述した。

図 1 は、FX100 及び TSUBAME3.0 におけ

通信演算オーバーラップの実装により、る、格子解像度 109 点(= 128 x 128 x 64 x 32

GYSELA コードにおいても演算性能の向上

x 16 x 2)、32～128 CPU/GPU でのストロンが実現された（表１）。しかしながら、並列度
3

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

表１．GYSELA コードのタスクレベル通信演算オーバーラップによるスピードアップ。
Nodes

Original

Overlap

Speedup6447.1 s

31.9 s1.48

12843.4 s30.4 s1.43

25628.7 s23.4 s1.22

51224.6 s22.1 s1.11に応じて MPI 通信や各演算タスクの比率が変化していくため、単一の静的コーディング
では必ずしも最良のタスクレベル通信演算オーバーラップが実現されないという課題も明

図２．時間刻み幅に対する誤差の収束性。赤が従らかとなった。このため、有向非巡回グラフ

来の陽解法、青が今回実装した陰解法によるシに基づくアルゴリズムや StarPU のようなタミュレーション結果。

スク処理フレームワークを利用してタスクの依存関係の解決を自動化することで、並列度へのスケーラビリティは確保しつつ、衝突に応じてタスク処理比率が変化しても、その項の陰解法による反復法ソルバーは MPI 通

都度ロードインバランスを最小化することが信フリーで計算できるというメリットがある。

望まれる。

図 2 は従来の陽解法（4th RKG）と新たに実装した陰解法（2nd Implicit）による、時間刻み幅に対する誤差を図示したものである。陽5.3. 衝突項陰解法の高効率実装（前期）

解法では衝突項の CFL 条件により、dt=8x10-次に、GKV コードのアルゴリズム拡張とし5R/vti 程度の小さな時間刻み幅までしか安定て、衝突項の陰解法化を行った。この拡張の動機は、特に Shortfall 問題で対象とするように解けないのに対し、陰解法では衝突項によな装置周辺領域では、プラズマの温度が低くる律速を超え、移流項による CFL 条件の限界なり、衝突周波数が高くなることに起因する。

dt=2x10-3R/vti 程度まで安定に解けている。ま

速度分布関数を扱う運動論的シミュレーショ

た、誤差のスケーリングから、今回設計した

ンでは、衝突周波数が粒子速度の 3 乗の逆数通りの 2 次精度の収束性が確認された。に比例するような依存性を持つため、低速のさらに、開発した手法を非線形乱流解析に電子の扱いにおいて衝突項の CFL 条件が時適用した際の演算性能や求解時間の評価を間刻み幅を律速する。この制約を解消するた

FX100 上で行った。その結果を表２にまとめめ、作用素分割法を利用して、無衝突項は従る。陰解法では、その反復法ソルバーのため来通りの陽解法で、衝突項は陰解法で解くアに多数の演算を必要とするが、前述したようルゴリズムの実装を行った。ここで、衝突項に、この反復演算は MPI 通信フリーに計算では速度空間に対する微積分オペレータであるきるため、演算密度が大幅に向上しており、ことに着目し、速度空間並列プロセスに対すFLOPS 値は約 3 倍に向上している。全体のるデータの転置通信を行い、実空間を細分化演算数が増えるために時間ステップあたりのし速度空間を大域的に持つようなデータの持処理時間は 8 割ほど増加した。しかしながら、ち替えを行う。これは、MPI 分割の部分集合

この例では時間刻み幅を 20 倍ほど大きくとに対する転置通信であるので、他の次元方向れるため、最終的な求解時間としては約 10 倍4

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

表２．FX100 上での非線形計算の性能評価。格子解像度は 1.8x109 点で、FX100 の 192 ノードを利用して計測。プラズマパラメータはShortfall 問題の解析に用いられる実験的パラメータとした。
4th Explicit

2ndImplicit

FLOPS

8.67

25.78(/PEAK)

TFLOPS

TFLOPS(4.01%)(12.5%)

図３．非線形乱流シミュレーションにおける熱1.76

3.13

輸送フラックスの時間発展。電子およびイオン

sec/step

sec/step

熱輸送について、緑破線、赤点線が従来の陽解法

Elapsedtimeper

step
Time

による結果、桃一点鎖線、青実線が今回実装した
step

5x10-5 R/vti

1x10-3 R/vti

陰解法による結果を示す。

size
Speed-up to 1

10solution

fastertimes

熱輸送フラックスは小さなずれ（時系列データの標準偏差内）にとどまっている。このことから、今回開発した衝突項陰解法が、計算の高速化が実現された。以上のように、衝突を～10 倍程度高速化しつつ、プラズマ乱流の項の陰解法実装は物理的にも計算効率的にも振る舞いを解析するうえで十分な精度を持つ

メリットがあることが示された。このようなことが確認された。

省 MPI 通信および高演算密度の実装はメニ
5.5. 物理解析への応用成果（前期）ーコアアーキテクチャにおいて特に有用と期

物理応用課題として、Shortfall 問題の原因待される。

究明に向けて、Flux tube train モデルと呼ばれる改良局所モデルを解析に適用した。Flux

5.4. 衝突項陰解法の高効率実装（後期）
衝突項陰解法の開発に関しては前期で一通tube train モデルでは、従来の局所モデルでり完了し、後期は開発した手法が非線形乱流想定される磁力線方向の乱流相関について緩の物理解析に適用可能であることを実証する

和して、磁力線方向に長いシミュレーションため、長時間の非線形シミュレーションを行領域をとることで、長い自己相関を持つ乱流い、乱流揺動により引き起こされる熱輸送フ揺動を扱うことができるモデルである。特に、ラックスの評価を行った。図３に示す通り、Shortfall 問題のパラメータ領域で不安定化

従来の陽解法と新たに実装した陰解法は、線される補足電子モードは磁力線方向に長い相形の範囲ではほぼ完全に一致していることが関を持ちうるため、Flux tube train モデルが示された。長時間の振る舞い（t>55R/vti）で有効となると期待される。表３は従来的な局はごく僅かにずれが生じるが、これは非線形所モデルによる解析と Flux tube train モデ

現象に見られる初期値鋭敏性のためであり、ルによる解析の比較を行った初期的結果であ時間積分法の違いによる数値誤差が長時間のる。Flux tube train モデルを用いることで、振る舞いのずれを生み出す。しかしながら、乱流輸送フラックスが 2 割程度増大している

期待される通り、乱流揺動の時間スケールやことが確認できる。この結果は、従来的な局5

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

表３．Flux tube train モデルによる Shortfall 問
題のイオン熱輸送フラックスの解析。トレインの長さ Nt=1 は従来的局所モデルに相当し、トレインを長くするほどイオン熱輸送が増大する傾向。
トレインの長さ

イオン熱輸送フラックス

図４．数値フィルタにより Convective cell を排Nt

Qi除しなかった場合（左図）、 排除した場合（右図）

1156.2の静電ポテンシャルのポロイダル断面図。

3155.65179.0

7190.2所モデルで仮定される磁力線方向の乱流相関が Shortfall 問題の要因の一つとなり得ることを示唆している。
図５．数値フィルタにより Convective cell を排5.6. 物理解析への応用成果（後期）
前期は Shortfall 問題に対して、GKV コー除しなかった場合（左図）、 排除した場合（右図）の新古典輸送の時空間プロット。

ドを利用し、局所モデルの改良という観点から議論した。後期は、GYSELA コードを利用し、大域モデルの観点から解析を行った。局駆動の Convective cell により新古典輸送が所モデルでは、乱流揺動のスケールが大域的増大する。このプロセスは、Convective cell

平衡のスケールに小さいという仮定の下でスの有無によってプラズマの分布関数が異なケール分離を行い、乱流輸送と大域的新古典る平衡状態へと緩和することに起因すると輸送を別々に取り扱うが、大域モデルでは両考えられる。本課題により、局所的コードで者の発展・相互作用を同時に取り扱う。

は表現できない物理過程による輸送増大

本課題では乱流が駆動する対流胞（Shortfall の解消）の可能性が示された。今（Convective cell）と新古典輸送との結合と後さらに、大域―局所モデル間の比較を進めいう、局所モデルでは表現できない物理過程ることで、Shortfall 問題に対する新しい視点に着目して解析した。図４は、数値フィルタを与えると期待される。

により Convective cell を排除した場合としなかった場合のポテンシャル構造の違いを示

6. 今年度の進捗状況と今後の展望す。基本的なモード構造に違いは無いが、フィルタによりポロイダル方向に長波長を持つ

当初の計画通り、前年度開発した GPU 版モードが排除出来ていることがわかる。続いGKV コードを CUDA/Thrust 系へと移植すて、図５に示すように新古典輸送レベルの比るとともに、畳み込み積分、移流計算実行時較をフィルタ有無の場合で行なった。図５右の通信と演算のオーバーラップを実装した。

図に見られるようにフィルタを適用した時点これによって、TSUBAME3.0 では、FX100(t = 350)で、新古典輸送が大きく低減されると比べ 2 倍程度の性能向上に成功した。
次に、Shortfall 問題で対象となる装置周辺ことが明らかとなった。言い換えると、乱流6

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

Lett. 119, 195002 (2017).

領域の解析に有用な衝突項の陰解法アルゴリズムの実装を行った。誤差評価と性能評価を行い、その物理的・計算効率的有用性を実証(2) 国際会議プロシーディングスした。
以上より、コード開発に関しては今年度の(3) 国際会議発表

当初計画以上の進展を達成した。さらに、今[3] （口頭，招待講演）S. Maeyama, “Cross-scale後の展望として、任意の問題規模・並列度にinteractions between electron and ion scale

おいて通信演算オーバーラップの効果を引き

plasma出すために、有向非巡回グラフに基づくアルstructures”, International Conference on

ゴリズムや StarPU のようなタスク処理フレ

Numerical Simulation of Plasmas, Leuven,ームワークを利用して、タスクの依存関係のBelgium, September 19, 2017.

turbulence[4] （ 口 頭 ） S.

解決を自動化し、ロードインバランスを最小

via

sub-ion-scale

Maeyama,“Multi-scale化することが有効となるだろうという開発指turbulence simulations of micro-tearing

針を得た。

modes and electron-temperature-gradientmode”,

物理応用に関しては、従来の局所モデルシ3rdUS-JapanJIFT

Exascale

ミュレーションで報告された Shortfall 問題

Computing Workshop, Kashiwa, Japan,（熱輸送の過小評価問題）に関して、拡張局

August 8, 2017.

所モデルにおける磁力線方向乱流相関の緩和[5] （ 口 頭 ） Y. Asahi, V. Grandgirard, Y.

や、大域モデルにおける乱流揺動・大域的新Idomura, X. Garbet, G. Latu, Y. Sarazin, G.

古典輸送相互作用といった、従来的局所モデ

Dif-Pradalier, P. Donnel, and C. Ehrlacher,

ルに含まれない効果が輸送増大を導くことを“Benchmarking of GT5D and GYSELA”,示した。どの物理効果がどの程度支配的になtheるか、Shortfall 問題を定量的に解消しうるかComputing Workshop, Kashiwa, Japan,

は今後の研究課題として残された。

August 8, 2017.

3rdUS-JapanJIFT

Exascale[6] （ポスター）S. Maeyama, T.-H. Watanabe,and A. Ishizawa, “Effects of electron-scale

7.

turbulence on micro-tearing modes”, 7th

研究成果リスト(1) 学術論文

Asia Pacific Transport Working Group[1] （投稿中）S. Maeyama, T.-H. Watanabe, Y.

International Conference, Nagoya, Japan,

Idomura, M. Nakata, and M. Nunami,June 6, 2017.

“Implementation of a gyrokinetic collision[7] （ ポ ス タ ー ） Y. Asahi, V.Grandgirard, Y.

operator with an implicit time integration

Idomura, X. Garbet, G. Latu, Y. Sarazin, G. Dif-scheme and its computational performance”,

Pradalier, P. Donnel, C. Ehrlacher, “Experience

submitted to Comput. Phys. Commun.

with GYSELA & Accelerators”, EoCoE workshop(2018).

F2F meeting, 29-30 November, Toulouse, France.

[2] S. Maeyama, T.-H. Watanabe, and A.
Ishizawa,
microtearing“Suppressionmodes

byofion-scale(4) 国内会議発表“マイクロティ[8] （口頭、招待講演）前山伸也，electron-scalenonlinear

アリングモードが乱流輸送に与える影響”，第interactions in tokamak plasma”, Phys. Rev.

23 回 NEXT 研究会，京都，
2018 年 3 月 1 日．

turbulenceviacross-scale7学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月[9] （口頭）前山伸也，國成希浩，渡邉智彦，“Fourier-Laguerre 展開による磁化プラズ
マ乱流中の位相空間カスケードの解析”，第73 回日本物理学会年次大会，東京理科大学，
2018 年 3 月 23 日．
[10]（ 口 頭 ） 前 山 伸 也 ，“ Cross-scaleinteractions between electron and ion scale
plasmaturbulencevia

sub-ion-scale，平成 29 年度閉じ込め・輸送研
structures”
究会，核融合研，2017 年 12 月 14 日．
[11]（口頭）前山伸也，“Effects of electronscale

turbulence

inmagnetic

fusion，実験室・宇宙プラズマ研究会「乱流・plasma”
輸送・粒子加速」，名古屋大学，2017 年 12 月
7 日．
[12]

（ポスター）前山伸也，渡邉智彦，石澤明 宏 ，“ Multi-scale interactions between
electron-temperature-gradient-modeandmicro-tearing-mode turbulence ”， Plasma
Conference 2017，姫路，2017 年 11 月 23 日．

(5) その他（特許，プレス発表，著書等）[13]

前山伸也, 第 12 回日本物理学会若手奨励賞（2018 年 3 月）．

8