RISTニュース No．
５７（２０１４）

ハイパフォマンスコンピュティングの進展によるプラズマ乱流研究のブレクスル      

（独）日本原子力研究開発機構
システム計算科学センター高度計算機技術開発室
前山 伸也
磁場閉じ込め型核融合装置におけるプラズマの閉じ込め特性はプラズマ中で生じる乱流輸送により大きく左右される。近年の超並列計算機の性能向上と並列計算技術の開発により、プラズマ乱流シミュレーション研究は大きく進展してきている。本稿では、京コンピュータにおけるプラズマ乱流シミュレーションコードの最適化手法と並列演算性能向上に対する効果について解説する。さらに、コードの大幅な性能向上により実現された、電子／イオン系マルチスケール乱流シミュレーション解析について紹介する。
１．はじめに

ラズマ乱流について簡単に触れ（第２節）、プ

磁場閉じ込め型核融合炉開発においては、ラズマ乱流の数値シミュレーションについて

炉心プラズマ物理、超伝導技術、熱流体、炉説明する（第３節）。その後、京コンピュータ

壁材料、システム技術等の多様な物理的・工における並列計算の最適化手法（第４節）、パ学的研究課題が存在するが、中でも大規模なフォーマンス向上（第５節）について解説し数値計算を必要とするのがプラズマ乱流シた後、コードの大幅な性能向上により新たにミュレーションである。そのため、プラズマ実現された電子／イオン系マルチスケール乱乱流シミュレーション研究では、物理モデル

流シミュレーション解析について紹介するの拡張とともに並列計算手法の開発が精力的（第６節）。に進められている。近年では、スーパーコンピュータ「京」を中核としたハイパフォーマ

２．磁場閉じ込め型核融合装置とプラズマ乱流ンスコンピューティングインフラストラク

核融合エネルギーは、燃料資源が豊富であチャー（）の活用を目的とした
戦ることや高レベル放射性廃棄物が少なく環境略プログラム分野４「次世代ものづくり」

負荷が低いことから、長期的な代替エネルギ（２０１１年度〜）や、核融合シミュレーションコー源として期待されている。最も研究開発のードのエクサスケールに向けた大規模並列進んでいる核融合エネルギーへのアプローチ

化・高性能化を多国間協力関係の下で進めるは、強力な磁場によって高温・高圧のプラズ

８多国間国際研究協力事業（２０１１

マを環状装置に閉じ込めて、熱核融合反応を年度〜２０１３年度）といった、計算科学との連引き起こす方法である（図１）。このような磁携も活発に行われており、その研究成果はプ場閉じ込め核融合装置の閉じ込め性能は、プラズマ・核融合分野のみならず、計算科学分ラズマ中の電磁的な揺らぎを伴う乱流による野においても高く評価されている［１］。

輸送に大きく影響される。そのため、プラズ

以下では、磁場閉じ込め型核融合装置とプマ乱流は核融合研究の中心的課題の一つであ−−RISTニュース No．
５７（２０１４）

図１

磁場閉じ込め核融合装置概念図。環状に配置した超伝導コイルの作る強力な磁場でプラズマを閉じ込め、外部加熱（または核融合反応による内部加熱）により高温・高圧に保つことで熱核融合反応を引き起こす。プラズマ中で生じる電磁的揺らぎが作る乱流輸送によって閉じ込め特性が左右される。

る。高温・高圧のプラズマ内部で起こる微細えて、粒子速度２次元を考慮しなければならな乱流構造の直接計測は困難であるため、数ないため、通常の３次元流体シミュレーショ

値シミュレーションによる乱流輸送の評価・ンと比較すると大規模な数値計算が必要とな予測が求められる。

る。

核融合プラズマは装置サイズ（数）に比プラズマ乱流の空間スケールは荷電粒子のべて格段に長い平均自由行程（１０）を持旋回半径と同程度であり、燃料水素イオンでつため粒子軌道や波・粒子共鳴といった流体は数程 度、電 子 で は０
１程 度 で あ

モデルでは表せない粒子の運動論的効果が重１２る。両者は質量比の平方根（）４３（水要となる。そのため、プラズマ乱流を記述す素イオンの場合）倍だけスケールが離れている第一原理的方程式として、理論・数値解析

るため、従来はスケール分離を仮定した単一に広く用いられているのがジャイロ運動論でスケール乱流解析が主に行われてきたが、両ある。ジャイロ運動論は、プラズマ乱流の時

者は中間スケールの波との非線形結合を介し間スケールが強力な閉じ込め磁場のローレンて相互作用しうる。そのため、プラズマ乱流ツ力による荷電粒子の旋回運動（サイクロトは本質的に電子スケールとイオンスケールをロン運動）の周期よりも非常に長いことを利同時に含むマルチスケール物理である。しか用して、旋回位相について平均化した方程式しながら、電子とイオンの時間・空間スケーである。粒子速度は磁力線平行方向速度v

ルを同時に解像したマルチスケールシミュレと磁気モーメントμ（磁力線垂直方向速さにーションは、単一スケールシミュレーション関連した量）のみで記述されるように簡約化と比べて大幅に計算コストが増大する。

される。同時に、旋回位相やサイクロトロン

上述のように、プラズマ乱流を理論的・数周波数の依存性が方程式系から除去されるた値的に解析する際には、運動論効果とマルチめ、数値計算コストが削減されるという利点スケール性が課題となる。これらを解き明かがある。しかしながら、粒子位置３次元に加す一つの直接的手段であるジャイロ運動論に−−

RISTニュース No．
５７（２０１４）

図２ （）

的（粒子的）
描像：様々な位置（）にいる様々な速度（，μ）を持った電子やイオン（）の運動としてプラズマの挙動を扱う。
（）
的（連続体的）描像：ある位置、ある速度の粒子の存在確率を表す分布関数
（μ）の変化としてプラズマの挙動を扱う。

基づくマルチスケール乱流シミュレーション

その基礎方程式は５次元位相空間におけるは、ハイパフォーマンスコンピューティング

分布関数揺動（μ）の時間発展を表す

方程式である。

としても挑戦的な課題である。マルチスケー
ル性による空間高解像度に加え、速い電子の運動と遅いイオンの応答を同時に取り扱うた（１）

めに大きな時間ステップ数を必要とするため、現実的な時間で計算を行うためにはスト

ここで、
 はジャイロ中心速度を表し、プラズマは電磁揺動による力を受けて磁力線垂ロングスケーリングの向上が極めて重要となる。

直方向にドリフト運動する。
はそれぞれ磁力線平行方向加速度、粒子間のク３．プラズマ乱流の数値シミュレーション

ーロン衝突を表す衝突項、平衡分布による寄ジャイロ運動論に基づくプラズマ乱流シ与を表す。電磁揺動ポテンシャルは、分ミュレーションは、歴史的には的

布関数揺動の速度モーメントとして与えられ（粒 子 的）描 像 に 基 づ い た





る電荷密度・電流密度揺動（）法によるシミュレーション［図２（）］て、以下の





が先行して開発された。その背景として

無撞着に決定される。によっ方程式から自己法は比較的小さなメモリでも計算が可能であったことが挙げられる。一方、近年のメモ（２）

リ容量の増大により、
的（連続体的）描像に基づいた格子法によるシミュレーション［図２（）］が可能となった。粒子ノイズ（３）除去や計算精度の安定性向上などの理由か
ら、現在では後者が主流となっている。

（ここで、
−−

はジャイロ中心位置と粒子位

RISTニュース No．
５７（２０１４）

置のずれによる分極密度である。
）式（１３）プロセッサコアを内蔵したと１６の共によって閉じた方程式系が構成される。物理有メモリ、ノード間データ通信用インターコ

的には、平衡分布における密度勾配や温度勾ネクトコントローラで構成されている。この配といった熱力学的な力がプラズマ中の揺らような階層的メモリ構造を持った並列計算機ぎの駆動源として働く。そうして駆動されたにおいては、ノード間の分散メモリ並列には電磁揺動による分布関数揺動の非線形移流過





  （）を、ノ程によって、揺らぎの非線形混合が起こり乱ー ド 内 の 共 有 メ モ リ 並 列 に は 流状態へと発展する。



 （）を用いる方法が効果

数値計算の観点からは、式（１）は各粒子的である。コードでも ハ種（）についての５次元位相空間（

イブリッド並列を実装し、計算ノードあたりμ）における数値流体力学（）計算と１つのプロセスを割り当て、さらに計算みなすことができる。磁力線垂直方向への局ノード内で８つのスレッドを立ち所性や統計的一様性を仮定することで、中性上げることで８つのプロセッサコアで並列処流体と同様にフーリエスペクトル法による高理を行う。

精度計算が適用される。乱流のダイナミクス
４．
２ 多次元領域分割が重要となる磁力線垂直方向にはフーリエスペクトル法を、その他の次元方向には差分法

を用いた並列化には、各プロセスを用いてプラズマ乱流の数値解析を行うのが異なる計算領域を担当するように領域分割が、ジャイロ運動論的シミュレーションコーを用いる。多次元・多粒子種の問題であるこド［２，３］である。主要な計算は磁力線とを活用し、図３のように
（または），


垂直方向座標に関するフーリエスペクトμについての５次元領域分割を適用する。

ル法、磁力線平行方向座標
および速度空間複数の次元で領域分割することは、各ノード座標μに関する４次または５次精度差分間の通信量を低く抑えるとともに、通信に参法、電荷密度・電流密度を求めるための速度加する
プロセスの部分集合の要素数が空間・粒子種μに関する積分であり、こ小さくなるため
通信におけるレイテン

れ と ４ 次 精 度 陽 的

 

法を組

シ（データ転送を要求してから実際に転送さみ合わせて乱流揺動の時間発展を計算する。

れるまでの遅延時間のこと。通信に参加する
プロセス数に比例。
）を抑える観点から４．コードの並列計算モデルと京コンも効果的である。

ピュータにおける最適化

５次元領域分割を適用するにあたり、フーリエスペクトル法における
方向の並列高

４．
１ 分散・共有メモリ型並列計算機にお
けるハイブリッド並列モデル

速フーリエ変換（）のための転置通信、




μ方向の差分演算のための一対一通信、

京コンピュータは８２９４４の計算ノードから


μ


方向の積分のための総和通信といっなる大規模分散・共有メモリ型並列計算機でた複数の通信が必要になる。
通信

ある。計算ノード間はインターコネクは領域分割により生じた余剰の処理なので、トと呼ばれる６次元メッシュ トーラストポ演算に比べてできるだけ低く抑えたい。しかロジーを持ったネットワークで接続されておし、並列数の増加に伴い、通信コストと演算り、適応性の高い３次元トーラスネットワーコストの比は一般的に増大していく。そのたクをユーザに提供する。計算ノードは８つのめ、数万ノード規模の並列スケーラビリティ−−RISTニュース No．
５７（２０１４）

図３５次元位相空間における領域分割。それぞれの粒子種
について、
，，
μ方向は３次元ボックス状に分割し、それぞれのボックス内でさらに（または）方向への１次元領域分割を行う。全体として粒子種１次元＋空間４次元の５次元領域分割となる。必要となる通信として、並列に関する転置通信（赤矢印）、差分演算に関する１対１通信（緑矢印）、電化密度・電流密度を求める総和通信（橙枠）を代表的に示す。

を実現するには通信のコストを削減すれの転置通信を各セグメント内に閉じ込めるることが重要となる。

ことができ、通信帯域を最大限に活用して通信コストを最小化できる。さらに



μ方

４．
３ ３次元トーラスネットワークにおける向の一対一通信も隣接するセグメント間のみで行い、かつ、
μ方向の総和通信も３次
プロセス配置最適化
格子計算を行うアプリケーションを元ネットワークの各平面内に収まるように数万ノード上で並列実行する際に、メッシュプロセス配置を設計することで、これら／トーラストポロジーを持つネットワーク構の通信コストも極小化する。これにより通信造はよいスケーラビリティを得やすい。インコストを大幅に減少させることが可能となっターコネクトの性能を最大限に引き出すためた。にはネットワークトポロジーに留意した最適化が必要となるが、その効果は具体的なアプ

４．
４ 通信と演算のオーバーラップリケーションに依存する。ここでは、多次さらに、通信コストを実効的にマスクする

元・多粒子種問題における通信コストを３次ためには、通信と演算のオーバーラップが有元トーラスネットワーク上で最小化すべく、効である。ここでは、図５に示すように、セグメント化プロセス配置を考案した（図の１スレッドを通信用スレッドと４）。キーとなるアイデアは、最も負荷の大して利用する方法［４］を用い、並列やきい方向の転置通信に関与する
プロ

差分演算における通信と演算の同時処理のたセスを３次元ネットワーク上で個別に密集さめに、異なる次元方向に関する計算の独立性せることである（その集団の単位をセグメンを用いて、パイプライン化オーバーラップをトと呼ぶことにする）。これにより、それぞ適用した［５］。通信スレッドを利用して、通−−RISTニュース No．
５７（２０１４）

図４ （）３次元トーラスネットワーク構造と（）セグメント化プロセス配置最適化。並列に関する転置通信（赤矢印）、差分演算に関する１対１通信（緑矢印）、電化密度・電流密度を求める総和通信（橙枠）を代表的に示す。

信と演算のオーバーラップを行うというアイ加やのロードインバランス、スケデア自体は従来も研究されてきた。しかし、ジューリングオーバーヘッドを抑制しつつ、効果的に通信コストをマスクするためには注効率的な通信と演算の同時処理が可能となっ意深い実装が必要である。具体的には、通信た。

コストをマスクするのに十分な量の演算処理を確保するために演算カーネルの再配置を行うこと、パイプライン数を適切に選ぶこと（細かいパイプライン化はオーバーラップできる演算量を増やすが、同時にレイテンシの増加も招く）、スレッド並列におけるタスクの粒度を適切に調整すること（小さな粒度への細分化はロードバランスを向上させるが、同時にスケジューリングに関わるオーバーヘッドの増加も招く）などに注意が必要で
図５ 通信スレッドを利用した通信と演算の同時処理。

ある。コードにおいては、５次元空間に対応した多重ループに対して、パイプライン数・タスク並列化の粒度の調整を行った。
さ ら にのス ケ ジ ュ ー リングを用いることで、通信スレッドも通信完了後に演算処理に加わるようコーディングを行った。これにより、のレイテンシの増−−RISTニュース No．
５７（２０１４）

京コンピュータにおけるコードのストロングスケーリング。（問題規模：１０２４×１０２４×９６×９６×３２×２
６２００億格子点）

図６ 最適化なし、プロセス配置最適化、通信と演算の同時処理の場合の性能比較。
（問 題 規 模：２５６×２５６×３２×３２×３２× ２

４３億格子点、並列数：４０９６コア）

図７

５．パフォーマンスのオーバーラップにより、京コンピュータのマルチスケール乱流計算の性能評価を行う

フルシステム規模の約６０万コアまでの良好な前に、スケールダウンしたテストケースで、ストロングスケーリングが得られた。６０万並最適化手法による性能向上の効果を検証し列実行における並列化効率は７６８、演算性た。通信と演算の同時処理は分散・共有メモ

能 は７８６４、理 論 ピ ー ク 性 能 比 はリ型並列計算機における共通技術のため、京１０
７９８
２９で推移し、アムダール則より見積

コ ン ピ ュ ー タ の み な ら ず、原 子 力 機 構 のもられる実効並列化率は９９
９９９９４となり、９００（２１３４計算ノード、４プロセッサコア１００万コア級並列への適用可能性を実証した。

／計算ノード、ファットツリーインフィニバ

こうして、大幅な並列性能向上により、電子ンド）においても性能向上が得られることが／イオン系マルチスケール乱流シミュレー

確認された。しかしながら、並列を行うションを実行する見通しを得た。

，方向に高解像度を必要とするマルチスケール乱流計算を９００上で実行した場合、６．マルチスケール乱流シミュレーション

転置通信のコストがの演算コストに比べこれまでの多くのプラズマ乱流シミュレーて大きく、そのため通信コストを十分にマスションでは、スケール分離を仮定して電子まクできないという結果になった。一方、３次たはイオンの単一スケールのみの解析が行わ元トーラスネットワークを持つ京コンピューれてきた。これは、乱流を駆動する微視的不タ上では、図６に示すようにプロセス配置の安定性の線形解析において、イオンスケール最適化により通信コストが大幅に削減されで起こるイオン温度勾配不安定性と電子スケた。そのため、演算とのオーバーラップによールで起こる電子温度勾配不安定性は時間・り通信コストは実効的にマスクされ、演算性空間スケールとも質量比の平方根（）能の向上が実現された。

１２
４３倍だけ離れているためである。しかし

コードの京コンピュータにおけるスながら、両者の間には非線形性を介したマルトロングスケーリングを図７に示す。多次元チスケール相互作用が存在しうるため、スケ領域分割、プロセス配置最適化、通信と演算ール分離の仮定は自明ではない。マルチスケ−−RISTニュース No．
５７（２０１４）

図８

電子／イオン系マルチスケール乱流シミュレーションにおける静電ポテンシャル揺動のスナップ
ショット［（）初期飽和時と（）定常時］。図はトーラス垂直断面の局所領域（イオン旋回半径の１００倍程度）と、それをさらに１０倍に拡大した領域を表示している。
（）初期は成長率の大きい電子スケール乱流が卓越するが、（）定常状態ではイオンスケール乱流が支配的になる。

ール乱流の直接シミュレーションでは、大き子スケール乱流とイオンスケール乱流においくスケールの離れた２つの乱流駆動源を同時てもマルチスケール相互作用が重要となり、に解像しつつ時間発展を解く必要があるた

イオンスケール乱流による電子スケール乱流め、京コンピュータ上で最適化されたの抑制が起こることが示された。

コードにより初めて解析が可能となった。典
型的な計算量は２×１０１１格子点、３×１０５時間

７．おわりにス テ ッ プ 程 度 で あ り、１２，２８８計 算 ノ ー ド

超並列計算機と計算手法の進展により、従（９８，３０４プロセッサコア）を用いて１２０時間程来は実現不可能であった大規模シミュレー度の計算コストがかかる。図８にマルチスケション解析が核融合プラズマ研究にブレークール乱流シミュレーションにおける静電ポテスルーをもたらしている。しかし、演算性能ンシャル揺動のスナップショットを示す。初を引き出すためには高度な並列アプリケー

期は成長の速い電子温度勾配不安定性によりション開発が求められ、その実現のためには電子スケールの乱流が卓越する。このとき、物理モデル・計算機の両方に対する深い理解半径方向に長く伸びた渦構造を形成し、電子が必要となる。ここでは、プラズマ乱流の熱輸送を引き起こす。この構造はストリーマ

ジャイロ運動論的シミュレーションという多と呼ばれ、従来、電子熱輸送を担う物理過程次元問題に合わせた５次元領域分割モデル、と考えられてきた。しかし、成長の遅いイオ

京コンピュータの３次元トーラスネットワー

ン温度勾配不安定性によるイオンスケールのクに最適化したセグメント化プロセス配置、乱流が徐々に成長し、定常状態ではイオンス分散・共有メモリ型並列計算機に適合した通ケール乱流が卓越する。同時に、電子スケー

信スレッドを用いた通信と演算のパイプライ

ルのストリーマが抑制されることがマルチスン化オーバーラップといった並列計算手法をケール乱流シミュレーションにより明らかに実装し、約６０万コアまでの良好なストロングなった。そのため、最終的な定常状態ではイスケーリングと９９
９９９９４という高い実効並オンスケール乱流が電子熱輸送を担う。こう列化率を達成した。ストロングスケーリングして、大きくスケールの離れた実質量比の電の向上は、多数の時間ステップを計算する必−−

RISTニュース No．
５７（２０１４）

要があるマルチスケール乱流シミュレーショ

参考文献

ンの実現のために無くてはならない要素で［１］   


あった。こうして高速化されたコード





により、実質量比の電子・水素系におけるマルチスケールプラズマ乱流シミュレーション 


を実現し、開発した計算手法の有効性を実証  した。また、このシミュレーションの結果か     ［２］  


ら、大きくスケールの離れた実質量比の電 

子・水素系においても両者の相互作用が重要
 

［３］
   


であること、イオンスケール乱流によって電子スケールの乱流が抑制され、電子熱輸送に   

おいてもイオンスケール乱流が支配的になる   

ことを新たに明らかにした。こうした電子／      ［４］

イオン系マルチスケール乱流の直接シミュレ     ーションはマルチスケール相互作用のモデル

 


化や電子熱輸送評価に貢献する貴重なデータ  を提供すると期待される。

   謝辞
本研究は
戦略プログラム分野４「次  



［５］


   



世代ものづくり」および８多国間国際研究





協力事業の支援を受けました。ま






た、筆者は科研費２６８００
２８３の助成を受
  けました。本研究の主要な計算は理化学研究所の京コンピュータにおいて行われ、通信と演算の同時処理の汎用性検証には日本原子力機構の９００を、マルチスケール乱流の線形解析には国際核融合エネルギー研究センター
計算機シミュレーションセンター（）の
を利用しました。

−− 