J. Plasma Fusion Res. Vol.89, No.3 (2013)１
７１‐１７９講座 核融合プラズマシミュレーションの技法−大規模並列計算環境の活用−

４．Vlasov シミュレーションのコーディング技法
渡 邉 智 彦，井 戸 村 泰 宏１）
核融合科学研究所・総合研究大学院大学，１）日本原子力研究開発機構（原稿受付：２０
１
３年２月８日）本章では，磁場閉じ込め核融合プラズマの乱流輸送解析で用いられるジャイロ運動論的シミュレーションにおけるコーディング技法について解説する．前半は，局所フラックス・チューブ・モデルを用いた GKV コードにおける並列化方針と集団通信について述べる．後半は，グローバル配位において全分布関数の時間発展を求める
GT5Dコードを例に，核融合プラズマの大規模シミュレーション研究で用いられている最先端の技法を概説する．
Keywords:
gyrokinetic simulation, turbulence, parallel computing

４．
１ ジャイロ運動論的シミュレーションについてよって扱う GT5D コード［５］での技法について概説する．

数 keV にも及ぶ高温の磁場閉じ込め核融合プラズマで第１章で簡単に紹介されたように，５次元位相空間上のは，粒子の平均自由行程は装置サイズの数百倍にも及び，分布関数を取り扱うジャイロ運動論的シミュレーションは

通常の流体近似の限界を超えて，運動論的な取り扱いが不多大な計算コストを必要とする．このため，京コンピュー可欠である．特に，異常輸送の原因となる微視的乱流につタや国際核融合エネルギー研究センターの Helios，核融合

いては，ジャイロ運動論［１］に基づいたシミュレーション

科学研究所（核融合研）のプラズマ・シミュレータなどにが輸送フラックスの定量的な評価に必須の手段となっていおいて大規模なシミュレーションが進められている．ここる．このため，大規模なジャイロ運動論的シミュレーショ

では，そうした大規模シミュレーションで使われている手ン研究が精力的に進められ，輸送機構の解明と理解，さら法の一端を紹介したい．読者諸賢の今後の研究の一助となに閉じ込め性能の予測に関する研究が世界各国で展開されれば幸いである．

ている［２］．磁場閉じ込め核融合プラズマにおける微視的

４．
２ 局所ジャイロ 運 動 論 的 Vlaosv シ ミ ュ レ ーションの技法乱流シミュレーションについては，学会誌の講座としてもまとめられているので，ご参照いただきたい［３］．

４．
２．
１ フラックス・チューブ・モデルとその数値解法

さて，ジャイロ運動論的シミュレーションでは，実空間３次元と速度空間２次元の計５次元からなる位相空間上のトカマクやヘリカル型の磁場閉じ込め装置において発生分布関数の時間発展方程式を数値的に解くことになる．そする微視的乱流は，熱や粒子の異常輸送を引き起こす．

の解法には，粒子（マーカー）を用いた PIC（Particle-in-

ジャイロ半径 !
%"程度の大きさをもつ乱流渦は，装置サイ

Cell）法に基づくものと，分布関数の発展方程式を数値格ズ（!）に比べ十分小さいと考え，局所的な乱流輸送シミュ

子や関数展開などを使って離散化して近似的に解く，いわレーションが行われている．そのスケール比は，乱流揺ら

ゆる Vlasov シミュレーション法とがある．本章では，後者

ぎの磁力線平行方向 !
&"と垂直方向の特性長の比と同じの手法，すなわちジャイロ運動論的 Vlasov シミュレーショ

オーダーとなり，実空間において強い非等方性をもつ揺動

ンのコーディング技法について解説する．

成分を取り扱う必要がある．そのために，磁力線に沿って分布関数を連続体として取り扱う Vlasov シミュレ ー

シミュレーション領域を設定した，フラックス・チュー

ションは，数値手法の観点から３章で述べられた流体・

ブ・モデル［６］が有用である．文献［４］に挙げた GKV コー

MHD シミュレーションと多くの共通点を持つ．実際，差ドでは，実空間座標を !
$
!
%
!
&
"とし，さらに速度空間につい分に伴う袖通信や FFT など共通する課題が多い．一方，５て!
##!
$"の座標を選ぶ．ここで，$は動径方向，%は磁力線次元という高次元性や，速度空間積分を含む微積分方程式

ラベル方向，##は磁力線平行方向速度，$は磁気モーメン

であることなどから，コーディング技法としてユニークなトを表す．

トピックも含んでいる．これらを踏まえ，次節では局所的

ここでは５次元位相空間 !
$
!
%
!
&
!
##!
$"上で定義されたな乱流輸送を扱う GKV コード［４］を用いたフラックス・

ジャイロ中心の揺動分布関数 #"の時間発展を，静電近似

チューブ・シミュレーションで使われる手法について述（低 "近似）のもと，ジャイロ運動論的方程式

べ，４．
３節ではグローバルな乱流輸送問題を全分布関数に
4. Coding Techniques for Vlasov Simulations
WATANABE Tomo-Hiko and IDOMURA Yasuhiro

corresponding author’s e-mail: watanabe.tomohiko@nifs.ac.jp

１７１!2013 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.89, No.3 March 2013

& '!
(
( #
%
%
#
#
(
(
'!%(
+
" *
%
'"%#%
%
'") "
! "-.#!!
(,(-.

割と並列化の指針をまとめてみよう．
低並列度における実装'''%%$!%#!-.#
# ""(
#
(
)
%
) （１）
&(
%$ "

SX-7 は，１ノードが３２台のベクトル・プロセッサで構成
さ れ，そ れ ら を 一 つ の 筐 体 に 搭 載 し た 計 算 機 で あ り，にした が っ て 求 め る 場 合 を 考 え よ う．こ こ で，左 辺 の１ノードあたり約２８２GFlops という当時としては強力な演/
20内の項は線形微分演算からなり，*!+は "#!ドリフト

算性能を有していた．このSX-7の５ノード構成をフルに利による非線形対流項をポアッソン括弧式の形式で表したも用して，トカマクにおけるイオン温度勾配（ITG）乱流ののである．右辺の "(#()は速度空間座標についての２階微シミュレーションが行われた［４］．SX‐７では，ノード内の分を含むモデル衝突項を意味する．また，マクスウェル分

３２台の CPU でスレッド（共有）並列を行うことができ，ま布 #" に比例する項は，#(についての式としては非斉次項

た，コンパイラによる自動並列化機構が有効に働いていにあたり，ジャイロ半径効果を含むポテンシャル揺動 'た．この環境を利用するため，GKV コードでは，分散並列に依存する．このポテンシャル揺動は，分布関数揺動 #(化は５つのノードそれぞれに１つの MPI プロセスを割りを速度空間上で積分したジャイロ中心密度の揺動成分

当てるのみにとどめ，ノード内は自動並列を活用するという方針をとった．わずか５つの領域に分割するだけなの

-.#%
&*)
#
(
#
+$ &#$!($!#（２）

で，１次元領域分割を磁力線平行方向（0）座標に適用してコーディングを行った．同様の手法は，やはり比較的少数から準中性条件を使って求められる．ここで下付き添字の強力な計算ノードからなる SX‐９においても有効である．

$は，(
.
!
/)座標に関してフーリエ変換を行った揺動の波数

この場合に必要となる MPI 通信は，磁力線平行方向の移ベクトル $&(*/)を意味する．ゼロ次ベッセル関数 $*.!
!

流項を扱う差分演算にともなう袖通信と，式（３）の磁気面

は，揺動波数に対する有限ジャイロ半径効果を表す．また，平均における 0方向の総和演算に伴う集団通信である．前

磁気面上で一定のポテンシャルをもつ，いわゆるゾーナル

者については，第３章の MHD シミュレーションで解説さフロー成分に対して電子の密度揺動がゼロになることを取れたものと同様に行えばよい．後者については，数値積分り扱うために，磁気面平均を求めることが必要な場合があにおいて現れる総和演算（リダクション演算の一種）を領る．これは，ポテンシャル揺動の */&!成分に関して磁力域分割された変数に対して行うため，MPI_Allreduce

線平行方向 (
0
)座標について次式の積分操作を施すことでを利用する（具体的なコーディング例については４．
２．
３小

得られる．

節で解説する）．
ノード内の自動並列化については，ほとんどのdoループ$,$.&#1)$$.!
0
$/&!#
#1)#0

に対してソースコードの変更なしに動作した．スレッド並（３）

列は，主に最外側ループとなる%（磁気モーメント）方向に対して行った．ただし，２次元 FFT を構成する個々の１次

次に GKV コードで用いる数値手法を簡単にまとめよう．

元 FFT の手続きをまるごと自動並列で呼び出すための指上記の基礎方程式の性質から，(.
!
/)座標について波数空示 行 を 利 用 す る．も ち ろ ん，同 様 な ス レ ッ ド 並 列 は，

間(*/)へとフーリエ変換し，スペクトル法を用いるの*.!

OpenMP を用いても可能である．

が自然である．そのため，非線形項については，実空間へ

高並列度における実装と２次元 FFT を用いて逆変換した上で２つの項の積を取

地球シミュレータやプラズマ・シミュレータ（日立 SR

り，それを再度波数空間へと FFT により変換する．一方，

１６０００）においては，SX‐７／５M に比べて演算性能は数十

線形項については，(.
!
/)座標に関する空間微分は波数ベ倍程度増大し，これを活用して大型ヘリカル装置におけるクトルを乗ずるだけで簡単に計算できる．磁力線平行方向

ITG 乱流輸送のシミュレーション［７］を実行することがでの 移 流 項 と ミ ラ ー 項，さ ら に 衝 突 項 に は，座 標 成 分

きた．この時，利用できるノード数は２５６から５１２程度へと0!
-.!
%)に関する微分が現れるが，これらは差分近似に(増大し，個々のノードには８から３２のスレッド並列を行う

よって評価する．時間積分には，４次精度をもつルンゲ･

CPU または演算コアが実装されていた．領域分割数が増大

クッタ･ギル法を適用し，ジャイロ運動論的方程式を数値するのに伴い，一般に通信コストも上昇する．ただし，袖的に解く．こうして次の時間ステップで得られた #($ を数通信に関しては，分割領域の表面積を小さくすることで，値積分してジャイロ中心密度を求め，ポテンシャル揺動を

MPIプロセスあたりの計算量に対する通信負荷を減らすこ

計算するという手順を取っている．以下では，大規模並列とができるので（通信回数の増加による起動時のオーバー

計算機における実装技法についてみていこう．

ヘッドが問題にならない範囲では），多次元の領域分割が

４．
２．
２ GKV コードでの並列化方針有利となる．そこでこれらの計算機に対しては，(
0
!
-.!
%)

オリジナルの GKV コードは，核融合研で稼働していたの３次元について領域分割を行うことにした．一方，２次

SX-7/5M（NEC）において開発され，その後，地球シミュ

元の波数空間 (*.!
*/)については分割せず，同一の MPI

レータや核融合研のプラズマ・シミュレー タ（日立 SR

プロセスに納めるようにした．これは第３章で議論された

16000）をはじめ京コンピュータにいたるまで，様々な超並

２次元 FFT に伴う転置処理とその通信コストの負荷を避列計算機に移植されてきた．まず，それらにおける領域分けるためである．
１７２

Lecture Note

4. Coding Techniques for Vlasov Simulations

T.-H. Watanabe and Y. Idomura

速度空間座標 !
!#!
""についても分割されたため，それらに関する微分項（ミラー力と衝突演算子）の部分には袖通信が発生する．加えて，式（２）の速度空間積分に伴う総和演算を MPI プロセス毎に並列で行った後，リダクション演
算と集団通信を行う必要がある．すなわち，磁気面平均と速度空間積分という２つの異なるリダクション演算を実行することになる．また，"方向も領域分割するためにループ長が短くなることから，スレッド並列の並列度は４ない
し８程度が適度な設定である．
４．
２．
３ 積分演算の実装法
この小節では，上述した数値積分に伴うリダクション演
算の実装についてみていこう．並列化していない場合，台形公式に基づいた数値積分の核となる部分は
図１

int_phi(:) = (0.d0, 0.d0)
do i = 1, global_nz__do m =1, nm____int_phi(m) = int_phi(m) + phi(i,m) * dz__end doend do領域分割とコミュニケータの概念図．z 方向への分割された領域は異なるパターンで示され，それぞれが２次元速度
空 間 積 分 を 計 算 す る た め の コ ミ ュ ニ ケ ー タ（v_comm_world）を作る．一方，磁気面平均は，z 方向のパターンの異なる立方体一つずつからなるカラムでできたコミュニ
ケータ（z_comm_world）内のリダクション演算で計算される．また，それぞれの粒子種が一つのコミュニケータを使って定義されている．

の do ループからなる．ここで global_nz は，"方向の全のであり，特にその中での集団通信を行うのに有効であ

格子点数を，nm は他の次元についての配列サイズを表す．

る．今回の例では，磁気面平均と速度空間積分について異"方向を nprocz 個に領域分割したとすると，MPI プロセ

なるコミュニケータを利用することになる．４．
２．
２小節前

ス内の和を

半で述べた "方向への１次元領域分割の場合，上の例に挙げた z_comm_world は，全 MPI プロセスからなる既定のint_phi(:) = (0.d0, 0.d0)

コ ミ ュ ニ ケ ー タ MPI_COMM_WORLD に 等 し い．一 方，

do i = 1, local_nz""での３次元領域分割を行った場合，速度空間積分
"
!
!#!
!

__do m =1, nm

用のコミュニケータも別に用意する必要がある（図１参____int_phi(m) = int_phi(m) + phi(i,m) * dz照）．それを v_comm_world としよう．それぞれのコミュ__end doend doニケータは，全ランク番号（rank），"方向へのランク番号（zrank）と "方向の分割数（nprocz）から，次のように定義できる．

と し て 求 め た 後（こ こ で，global_nz=local_nz*
nprocz とする），nprocz 個のプロセスについて次のよう

rankz = mod( rank, nprocz )に MPI ライブラリを呼び出して int_phi について総和を

zcolor = rank / nprocz

とる．

vcolor = rankz
call MPI_Comm_split( MPI_COMM_WORLD, &call MPI_Allreduce( int_phi, int_phi_all, &zcolor, rank,z_comm_world, ierr )nm,MPI_DOUBLE_COMPLEX, MPI_SUM, &call MPI_Comm_split( MPI_COMM_WORLD, &z_comm_world, ierr )ここで，総和の結果は int_phi_all に格納され，"方向分

vcolor, rank,v_comm_world, ierr )割したすべての MPI プロセスにおいて同じ値となる．MPI

ここで，それぞれのコミュニケータは，zcolor または_DOUBLE_COMPLEX は，int_phi，int_phi_all が倍精

vcolor が同じ値をもつ MPI プロセスのグループから構成

度複素数型の変数であることを意味する（もちろん，int

される．このようにして，３次元領域分割された変数の数_phi の型に応じてここは変更する）．ここでは，MPI_SUM

値積分を簡単に計算することができる．

というパラメータ変数を指定することで，種々のリダク

４．
２．
４ より進んだ計算モデルとその並列化

ション演算の中から総和を選んでいる．z_comm_world

４．
２．
１小節では，静電的なジャイロ運動論的方程式にもは，"方向の分割に伴って定義したコミュニケータである．

とづくシミュレーション技法について説明した．現在さらこうして領域分割された座標方向への積分を求めることがに，電磁揺動，多粒子種，マルチ・フラックス・チューブできる．

などを取り入れた拡張が進められている．これらの拡張で必要となる技法を以下で簡単にみていこう．

ここで，領域分割の方法に応じて，どのようにコミュニケータを使い分けるかをもう少しみていこう．コミュニ

ジャイロ運動論ではゆっくりした電磁場変動のみを扱う

ケータは，全 MPI プロセスのなかの部分集合を定義するもため，変位電流を無視したアンペールの式を解く．そこで
１７３

Journal of Plasma and Fusion Research Vol.89, No.3 March 2013

は，揺動分布関数の１次モーメントで与えられる電流密度#$$ "+.$,"&!"(+"#

をソースとして，磁場揺動のベクトルポテンシャルを求めることになる．コーディング技法としては，静電場を解く（６）を用いて与えられる．ここでは前節と同様に静電近似，断場合と本質的な違いはない．
一方，ジャイロ（またはドリフト）運動論的電子や，多熱的電子応答のイオン系乱流を考えており，!")は粒子

成分イオンなどを扱う場合には，２つの方法が考えられ

位置，(と + はイオンの電荷と質量，,!)は平衡電子密度，

る．一つは，分布関数の配列の定義において，粒子種を新).+は熱速度で評価した軌道半径，%&)，%&+は電子とイオン

たな次元として追加する方法である．この場合，粒子種にのデバイ長，+（$*）は軌道平均した静電ポテンシャル，$,"

ついては並列化が難しくなるとともに，コード全体にわた（$+は磁気面平均した静電ポテンシャルをそれぞ
$,
$,
+
*
/）
*

る改訂が必要となる．もう一つは，粒子種についても分散れ示す．式（４）は衝突項"および熱源%-,(がない極限では，並列化する方法である．この場合，MPI のランク番号にリウヴィル方程式となり， )はジャイロ中心の軌道に沿っよってどの粒子種を担当するかを決めておき，さらに，そて保存する．式（５）は準中性条件から導かれるポアッソン

れぞれの粒子種をコミュニケータによってグループ分けす

方程式であり，左辺第１項は軌道効果によるイオン分極密ればよい．こうすることで，コードの大幅な変更は不要と度，第２項は電子断熱応答を示す．

なり，かつ，並列性能を損なうことなく，分散可能な並列

大域的モデルではフラックス・チューブ・モデルで仮定数を大きくすることができる．この手法は，また，径電場

している磁力線に垂直な平面での背景プラズマ分布，背景効果を取り入れるためにヘリカル系プラズマのシミュレー

磁場パラメータの一様性と周期境界条件を利用できないた

ションにおいて進められている，マルチ・フラックス・め，様々な問題を検討する必要がある．まず，準定常状態

チューブ・シミュレーションに対しても有効である．

のプラズマ分布における揺動分布関数 #)と背景分布関数

これまでの例では，２次元 FFT の効率的な計算のため，（〜#の時空間スケールの分離を仮定し，前者のみを発)!
'）波数空間は分割せず一つの MPI プロセスに納めていた．し展させる #)モデルを計算する場合，フラックス・チュー

かし，イオンおよび電子温度勾配乱流が共存する場合な

ブ・モデルでは，周期境界条件の下で式（１）のように解析

ど，より幅広い波数空間を扱う際には，波数空間について的に #)を分離できるのに対し，大域的モデルでは乱流輸送

も領域分割が必要となる．この際に問題となる波数空間内

による背景プラズマ分布の緩和や時間発展を打ち消して準でのデータ転置にともなう通信コストを隠蔽するため，定常状態を維持する熱源モデルが必要となる．このようなOpenMP を利用した通信部分と演算部分の並列実行手法

分布固定熱源は，分布補正の時定数の任意性や実験とは異が新たに開発され，京コンピュータのように１０万コアを超なるソース・シンク分布（例えば，プラズマ中に現れるシえる最先端の超並列環境においても，高い実行性能を達成

ンク）など，その取り扱いに注意が必要である．一方，実することができるようになった［８］．この新技法を用いた

験条件と同様に一定の熱源とプラズマ境界におけるシンク

大規模なジャイロ運動論的シミュレーションにより，イオを与えてプラズマ分布を自由に発展させる full- )モデルのンスケールから電子スケールまでを同時に含んだ乱流揺動

場合にはプラズマ分布が準定常状態に達するまで長時間ののスペクトル構造が明らかにされ，異常輸送現象におけるシミュレーションが必要になる．
full- )モデルの長時間乱流シミュレーションでは粒子数

複数の物理機構の協調および競合過程を詳らかにすること

保存の誤差がシミュレーションに大きな影響を与えるのができるであろう．

で，全分布関数 )の高い保存精度が要求される．一方，高温

４．
３ 大域的ジャイロ運動論的 Vlaosv シミュレーションの技法の核融合プラズマは粒子衝突による散逸効果が小さく非線

４．
３．
１ 大域的モデルとその数値解法ことが要求される．フラックス・チューブ・モデルの場合形性が強いため，非線形移流項の数値的安定性を維持する

前節で紹介したフラックス・チューブ・モデルと異なには擬スペクトル法の採用によってこの問題を解決していり，大域的モデルはトーラスプラズマ全体を計算領域とるが，大域的モデルでは背景プラズマ分布と背景磁場の非し，５次元位相空間'!!
.-!
&(上で定義されたジャイロ中心の全分布関数 )をジャイロ運動論方程式系
,
)
,
) /
,
)
/
)
!
$*
)
")
$ "!
#
&)
".
$"'
(
"%-,(
,
,
,
.-

（４）

#

).+
" $ $
&%$" # .
!+ ,/$!&%# #
*
%&+
%&+
!")/
'(#)
'%&!,!)"
$
.
!"(
! #'
（５）
に従って発展させる．式（４）の左辺は５次元位相空間のポ
図２

アッソン括弧式 )!*とジャイロ中心のハミルトニアン

１７４

大域的モデルで使用する円柱座標系(R , +, Z )および磁束座
標系( , (, -).

Lecture Note

4. Coding Techniques for Vlasov Simulations

T.-H. Watanabe and Y. Idomura

一様性（軸対象トカマクの場合にはトロイダル角について

式（５）は磁束座標系において取り扱い，軸対称トカマクにの対称性は存在）および非周期的な境界条件のために，スおける演算子のトロイダル角 #についての対称性を利用しペクトル法を用いることができない．このため，大域的モて，トロイダル方向のフーリエモード展開，および，ポロデルの場合には精度と安定性を両立する数値手法の選択がイダル断面上で２次元有限要素近似を組み合わせて静電ポ

きわめて重要になる．

テンシャルを解く．

フラックス・チューブ・モデルとのもう一つの大きな違

時間積分には半陰的ルンゲ・クッタ法を適用し，式（４）いは，磁気軸の取り扱いが必要になることである．これににおいて磁力線方向の熱運動を円柱座標系で表現すること

ついては磁束座標系（図２参照）の構造格子を用いる多くに起因するスティフな線形移流項を陰的ステップとして取のコードにおいて，磁気軸近傍のクーラン条件の発散を回り扱う．この陰的ステップの連立一次方程式は〜１０１０次元避するために磁気軸近傍を解かずに境界条件を設定する中の巨大な非対称疎行列となるため，非対称行列に対する基

空トーラス配位を採用している．一方，円柱座標系を採用

礎的なクリロフ部分空間解法として知られている共役残差して軸の問題を回避するコードも存在する．前者において

法［１１］に基づく反復法ソルバーを適用する．ここで，フは，プラズマ中心領域に設定された人工的な境界条件の影

ラックス・チューブ・モデルでは磁力線方向の座標 "で熱響に注意が必要であるのに対し，後者の場合には計算格子運動を表現することにより陽解法で時間積分を行えたのにが磁気面に沿っていないので，磁気面量や磁気面平均の取対し，大域的モデルではトロイダル角 "で熱運動を表現しり扱いに注意が必要である．

たために陰的ステップで大規模行列を取り扱うというデメ

４．
３．
２ GT５D コードの数値解法

リットが発生する．一方，円柱座標系の採用には，磁気軸次に GT5D コード［５］で用いる数値手法を簡単にまとめの問題を回避できる，トロイダル角 "（"!#）に中心差分

よう．式（４）は位置 !に関しては円柱座標系において取りを適用することによって離散的なトロイダル対称性が満た扱い，左辺の移流項に４次精度の森西差分スキーム［９］をされてトロイダル角運動量が厳密に保存される，あるい５次元一般曲線座標として与えられる案内中心座標系に拡は，トロイダル対称性を生かしたフーリエモード展開およ

張したスキーム適用する．森西差分スキームはポアッソン

びその並列化が可能になるといったメリットもある．この!
という２つの保存則を同時
括弧式が持っている !
および !

ように，大域的モデルの設計においては解きたい問題の特に満たすように中心差分を修正することによって，数値的性を考慮した上で各手法のメリット，デメリットを総合的な散逸無しに数値安定性を保証し，精度と安定性を両立すに勘案し，現実的な計算コストの範囲内で最適な座標系やる．一方，右辺の衝突項は線形フォッカー・プランク演算

数値手法を選択することが重要である．

子に６次精度の中心差分を適用し，かつ，衝突演算子の基

４．
３．
３ GT5D コードの並列化方針

本的特性である粒子数，運動量，エネルギーの保存則を厳GT5D は 原 子 力 機 構 に 導 入 さ れ た Altix3700Bx2（Ita-

密に満たす手法［１０］を採用する．このように，GT5D では

nium2，2048 コア，１３TFlops），あるいは，BX900（Nehalem

full- !モデルの長時間シミュレーションを実現するために-EP，１
７１９２コア，２００TFlops）といった超並列スカラー計

式（４）全 体 を 保 存 型 ス キ ー ム で 実 装 し て い る．一 方，算 機 で 開 発 さ れ，最 近 で は，Helios（SandyBridge-

図３

GT5D における Gyrokinetic solver（式（４）左辺），Collision operator（式（４）右辺），および，Field solver（式（５））の並列化モデル．(R , Z , !)方向の３つの MPI コミュニケータによって構成される２階層の MPI ネットワークとマルチコアのスレッド並列によるSMP層から構成される３階層の通信ネットワーク上で各演算子の対称性を利用した３次元ないし２次元の領域分割を実装している．
図は !方向４並列，R 方向２並列，Z 方向１並列，SMP2 並列の例を示し，通信処理を行うコミュニケータは色のついたボックスで示している．ここで，色の異なる通信処理は独立に処理される．

１７５

Journal of Plasma and Fusion Research Vol.89, No.3 March 2013

EP，７０５
６０コア，1.5 PFlops）
や京（SPARC６４VIIIfx，７０５０２４

ソルバー間のデータ転置処理を１階層の部分空間に閉じ込

コア，11.3 PFlops）といったペタスケール計算機にも移植めている点である．これによって，通信処理で使用するコされている．このような超並列環境における並列化で最もミュニケータのサイズ（通信に参加するノード数）が一桁重要になるのが，問題の対称性を利用して階層的な並列化以上小さくなり，通信コストが大幅に削減されている．

モデルを構築することである（図３）．例えば，式（４）左辺

４．
３．
４ 通信マスク手法の移流演算子（Gyrokinetic solver）は断熱不変量 $をパラ

上記の階層的並列化モデルによって BX900 上で数千〜メータとする４次元$!!
%&%の移流項となり，式（４）右辺の１万コアを用いた並列処理が実現し，JT-60U 規模の問題

衝突演算子（Collision operator）は位置 !をパラメータとサイズ（〜１０１０格子）に対して 99.9969％ という並列化率がする２次元 $%&!
$%の移流‐拡散項となる．一方，トロイダ達成された．ここで，並列化率 #はアムダールの法則におル方向にフーリエモード展開を適用すると，式（５）（Fieldいて全体処理に占める並列化されている部分の割合と定義

solver）はトロイダルモード数 $をパラメータとする２次され，１コアによる逐次処理に対する $コアの並列加速率

元ポアソン方程式となる．これらの演算子を対称性パラ!$ はメータについて分割すると，分割後の部分空間の演算は完!$ #全に独立に実行できる．このため，まずはこれらの対称性
パラメータを用いて粒度の粗い並列化を行うことで，並列!
$!!#"#"（７）と与えられる．逆に，並列化率 #は１コアと $コアでの処処理における通信コストを大幅に削減する．ただし，上記のように演算子毎に異なる対称性パラメータを用いて並列

理時間 "! および "$ で与えられる並列加速率 !$ #"!""$

化を行う場合には，ソルバー間でデータ転置処理が発生すから算出できる．ただし，メモリサイズや処理時間の制約るため，各ソルバーの演算コストと転置処理コストのバラにより１コアと $コアで同じ問題を処理して並列加速率!$を求めることが難しい場合には，$コアと#コアでの処

ンスに注意が必要である．

理 時 間 "$ お よ び "# か ら 得 ら れ る 並 列 加 速 率 の 比次に検討する点が対称性パラメータで分割した部分空間の問題に対する細かい粒度の並列化である．この並列化は!$"!# #"#""$ に式（７）を代入して並列化率 #を算出す

通信処理を伴うため，演算量と通信量のバランスを考慮し$る．式（７）はコア数 $を増やしていくと，並列化部分 #"た分割軸と並列度の選択が必要となる．差分法のようなスとそれ以外 !!#が釣り合う点 $"##"!!#%が現れ，それ$

テンシル演算の場合，演算量と通信量は領域分割した部分を超えると並列加速率はコア数無限大で得られる上限値

空間の体積と表面積にそれぞれ比例するため，できるだけ"!
!!#%に向かってゆっくりと漸近していくこ とを示$

多次元の並列分割軸を利用し，かつ，分割された領域が球す．このため，あるプログラムおよび問題サイズで実質的（あるいは，直方体）になるべく近くなる並列度の選択を行に利用可能なコア数の目安は $"となる．このことは，例えうことで与えられた並列度に対する通信量を最小化でき

ば，１００万コアの並列処理を実現するには，プログラムの並る．ただし，多次元並列化を行う場合には通信相手の数が列処理性能として 99.9999％（six nine）の並列化率が要求増大するため，通信ネットワークのハードウェア機構にさ れ る こ と を 示 し て い る．上 の GT5D の 例 で は，##

よっては処理コストが増大する場合もある．

99.9969％ に対して $"#32,257 となるが，実際に 16,384 コア

上記の考察に基づいて GT5D では各ソルバーに３次元なでの処理コスト分布をみると，通信処理のオーバーヘッド

いし２次元の領域分割モデルを適用し，第２章で説明したが３〜４割を占めて並列加速率が頭打ちになりつつあっ

MPI 並 列 化 と OpenMP に よ る SMP（Symmetric Multi-

た．このボトルネックを解決し，１０万コア以上のペタス

Processing）並列化を組み合わせて構成する３階層の通信

ケール計算に十分な並列化率を実現するために，GT5D でネットワーク上に実装している．ここで重要な点は，２階は通信マスク手法の開発を行った．

層の MPI ネットワークの一方の分割軸に対称性パラメー

図４に領域分割した差分演算の例を示す．通常の実装方

タを適用することにより，ソルバー内の通信処理，および，

法は袖領域のデータを隣接ノードから転送し，通信処理の

図４

領域分割した差分演算における通信マスク手法の概念図．通常は袖領域通信を実行してから差分演算を実行するが，通信マスク手法
では袖領域通信と内部領域の差分演算を同時処理することによって通信処理のオーバーヘッドを隠蔽する．

１７６

Lecture Note

4. Coding Techniques for Vlasov Simulations

T.-H. Watanabe and Y. Idomura

完了後に差分演算を実行する．一方，通信マスク手法の場

手法 C では，OpenMP のマスター構文を使用して，通信処

合には，袖領域通信と袖領域を参照しない領域の差分演算

理の完了までマスタースレッドを通信専用に割り当てるこを同時実行することによって通信コストのオーバーヘッドとによって通信処理が中断する問題を回避する．手法 Bを隠蔽する．このような演算と通信の同時処理は，従来，では演算と通信の同時処理が実現するのに対し，手法 C

MPI_Isend/MPI_Irecv のような非同期通信関数を用いではマスタースレッドが通信に占有されるため，通信処理ることで実装可能と考えられてきたが，MPI ライブラリやのオーバーヘッドがスレッド間で共有されて，例えば，８計算環境によってはこのような標準的な実装が機能しない

スレッド中１スレッドを通信に割り当てる場合には通信の場合があるので注意が必要である．多くの MPI ライブラリ

オーバーヘッドが 1/8 に減少する．現在のプロセッサは８ではある程度のデータサイズ以上の転送に RendezVouz

コア，１６コアとコア数が増加する傾向にあるため，この方プロトコルと呼ばれる通信方式を採用しており，ライブラ法でも十分にオーバーヘッドの削減が見込める．また，手リの内部ではデータ転送の開始前に通信相手ノードと制御

法 C は 非 同 期１対１通 信 だ け で な く 同 期 通 信，さ ら に通信を確立するような処理を実行している．しかしながは，集団通信にも適用できるため汎用性が高い．

ら，通信マスク手法では MPI 関数のコール直後に演算処理

実際，GT5D では図３の Collision operator を計算する前を開始するため，通信相手のノードが演算処理に占有され後に発生する MPI_Alltoall によるデータ転置にもこのて制御通信に応答できない状況が発生しうる［１２］．このよ手法を適用し，通信処理のオーバーヘッドを削減した．

うな場合，演算終了まで通信処理が中断するため，演算と図６に示す実装例では，データ転置前後の分割軸に関係し通信の同時処理が行われない．GT5D の開発では，BX900，

ない "方向に処理全体を分割し，さらに，データの転置とFX1，Helios，京等の多くの計算環境でこの問題に直面し逆転置の処理順序を入れ替えて，データ依存性のない "方

たため，これを回避する２つの方法を開発した［１３］．

向インデックスの転置，差分演算，逆転置を組み合わせて

図５は図３における Gyrokinetic solver の差分演算と袖同時処理することでデータ転置のオーバーヘッドを削減す

領域 通 信 の 実 装 例 を 示 す．手 法 B で は，演 算 開 始 後 にる．この手法は，前述の GKV の並列化 FFT のデータ転置

MPI_Test（非同期通信の状態をチェックする関数）のよ処理にも応用されている［８］．これらの通信マスク手法のうなダミー MPI 関数を定期的にコールしてノードを一時開発によって京において１０万コア以上を用いるペタスケー

的に MPI 処理に復帰させることによってRendezVouzプロ

ル計算が実現し，ITER 規模の問題サイズ（〜１０１１格子）のトコルの制御通信を確立し，通信処理を開始させる．一方，計算で !!99.9998％（!!!500000）を達成した．

図５８コアの環境において Gyrokinetic solver の差分計算と袖領域通信を実装した例．A は同期通信を用いるオリジナルの実装．B は演
算中に MPI_test をコールすることによって非同期通信を促進する通信マスク手法．C はマスタースレッドで通信処理を実行する通信マスク手法．

図６８コアの環境において Collision operator の差分演算とデータ転置を実装した例．データ転置前後の分割軸に関係しない "方向に処理全体を分割し（中図），データの転置と逆転置の処理順序を入れ替えて，データの依存性のないインデックスの転置，差分演算，逆転置を組み合わせて通信マスク手法を構築する（右図）．

１７７

Journal of Plasma and Fusion Research Vol.89, No.3 March 2013４．
３．
５ マルチコアプロセッサにおける最適化
第２章で現在のスカラー計算機におけるメモリの壁の問題が指摘されたが，この問題はしばしばルーフラインモデ
ル［１４］に基づいて議論される．このモデルでは，理論演算性能 "（Flops），理論メモリ帯域幅 !
（Byte/s）の計算機上で演算数 （Flop），メモリ参照数 $（Byte）のプログラムを%実行する際の実効演算性能 #（Flops）
が
#"

%
%
"
"!$"
!

（８）

と与えられる．この関係式を変形するとプログラムの対ピーク性能比 #""がプログラムの演算密度 %"$および計算機の !""比によって
%
"
$
#
"
" %
"
$!""
!

図７

（９）

ルーフラインモデルにおける対ピーク性能比 S/F のプログ
ラムの演算密度Flop/Byteおよび計算機のB/F比に対する依存性．図中に示すアプリケーション特性データは文献［１４］に基づいている．

と与えられる．図７に式（９）を示すが，この関係式は %"$が低い領域ではメモリ帯域幅 ! が性能を律速し，%
"
$が高

enddo

い領域では演算性能 "によって性能が頭打ちになる状況を

enddo

示 し て い る．図 中 で !"""#は 地 球 シ ミ ュ レ ー タ，

enddo

""!
!
!"
"""は Altix3700Bx2，!"$は BX900 や京に相当する性能であるが，ステンシル計算の演算密度の領域に着ここで，a，b，c はスキーム，格子幅等によって決まる定

目すると，ベクトル・プロセッサでは５０％以上を誇ってい数を示す．この差分カーネルの演算数はソースコードよりた CFD（数値流体力学）計算の対ピーク性能比がシングル

25 Flop となる．一方，メモリ参照数の見積りは少し複雑な

コアスカラープロセッサでは２０％程度まで低下し，現在の計算を必要とする．

マルチコアスカラープロセッサでは１０％程度といわれる実第２章でも紹介されたように，ほとんどのスカラー機は効処理性能の変遷を理解することができる．

キャッシュを搭載した階層的なメモリ機構をもっている．

このモデルから，現在のマルチコアスカラー機の性能を

例えば，京の場合には８コアで共有する６MByte の L２

引き出すには演算密度を向上する最適化がきわめて重要で

キャッシュを搭載しており，メインメモリのデータ転送速あることがわかる．本節では，そのようなマルチコアプロ

度が６４GByte/s となるのに対し，L２キャッシュのデータ

セッサにおける演算密度を向上する最適化技術を１つ紹介

転送速度は２５６GByte/s となる．このため，メインメモリへしたい．ここで，以下のような４次元の４次精度中心差分のデータ参照コストに比べてキャッシュへのデータ参照はを考える．

ほぼ無視できると仮定して上記差分カーネルのデータ参照
数を見積もる．ここでは，問題規模として nx=ny=nz=nv

real*8 f(-1:nx+2,-1:ny+2,-1:nz+2,-1:nv+2)＝４０のケースを考える．Fortran における配列のメモリア

real*8 df(-1:nx+2,-1:ny+2,-1:nz+2,-1:nv+2)

ド レ ス は N(i,j,k,l)=i+1+(j+1)*(nx+4)+(k+1)*!$OMP DO ←デフォルトではブロック分割となる(nx+4)*(ny+4)+(l+1)*(nx+4)*(ny+4)*(nz+4)のよ

do l=1,nv

うに変化する．このため，i 方向差分では N(i,j,k,l)-4

do k=1,nz〜N(i,j,k,l)+4 を参照するのに対し，j 方向差分では

do j=1,ny

N(i,j,k,l)-2*(nx+4)〜N(i,j,k,l)+2*(nx+4)を参

do i=1,nx

照し，k 方向，l 方向ではさらに参照するアドレスの範囲が

fx=a*(f(i+1,j,k,l)-f(i-1,j,k,l))+ &拡大する．このデータサイズを計算すると，i 方向：40b*(f(i+2,j,k,l)-f(i-2,j,k,l))Byte，j 方 向：1.68 kByte，k 方 向：70 kByte，l 方 向：

fy=a*(f(i,j+1,k,l)-f(i,j-1,k,l))+&2.96 MByte となる．一方，コアあたりの L2 キャッシュサイ

b*(f(i,j+2,k,l)-f(i,j-2,k,l))ズは 6 MByte/8 コア＝７
５０kByte となるため，k 方向までは

fz=a*(f(i,j,k+1,l)-f(i,j,k-1,l))+&オンキャッシュとなるのに対し，l 方向はメモリアクセス

b*(f(i,j,k+2,l)-f(i,j,k-2,l))となる．この議論に基づいて上のソースコードではキャッfv=a*(f(i,j,k,l+1)-f(i,j,k,l-1))+&

シュからロードするデータを実線，メモリからロードする

b*(f(i,j,k,l+2)-f(i,j,k,l-2))

データを破線，メモリにストアする（正確にはロードしてdf(i,j,k,l)=f(i,j,k,l)+ &からストアする）データを点線で表している．この判定条

c*(fx+fy+fz+fv)

件でメモリ参照数を見積もるとロード５，ロード１/スト

enddo

ア１となり，倍精度配列で 8 Byte×7＝56 Byte のメモリ参
１７８

Lecture Note

4. Coding Techniques for Vlasov Simulations

T.-H. Watanabe and Y. Idomura

照となる．したがって，このカーネルの演算密度は "!
!$のに対し，このケースではサイクリック分割（$ 番目のス

25 Flop / 56 Byte=0.45 となる．

レッドに #$%%!!&#""$ を割り当て）となる．このため，
隣接スレッドがロードしたキャッシュ上のデータを再利用

これに対して OpenMP のループ分割方法をサイクリッ

することによって，キャッシュ上のデータだけで l 方向のク分割に変更したのが，以下の例である．

差分を実行できる．この結果，差分カーネル全体のメモリ
real*8 f(-1:nx+2,-1:ny+2,-1:nz+1,-1:nv+2)

参照数は 24 Byte にまで削減され，演算密度も "!
!$25 Flop

real*8 df(-1:nx+2,-1:ny+2,-1:nz+1,-1:nv+2)/ 24 Byte=1.04 に向上する．この手法は共有キャッシュを有するマルチコアプロセッサに特有の最適化技術である．

!$OMP DO schedule(static,1)←サイクリック分割を指定
do l=1,nv

GT5D でもこの手法によって Helios および京における差分

do k=1,nz

カーネルの演算密度および実効演算性能を大幅に向上し

do j=1,ny

た．

do i=1,nx

参考文献

fx=a*(f(i+1,j,k,l)-f(i-1,j,k,l))+ &［１］洲鎌英雄：プラズマ・核融合学会誌 79, 107 (2003).
［２］X. Garbet et al., Nucl. Fusion 50, 043002 (2010).
［３］渡邉智彦 他：プラズマ・核融合学会誌 81, 534 (2005);
ibid 581 (2005); ibid 686 (2005).
［４］T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
［５］Y. Idomura et al., Comput. Phys. Commun. 179, 391 (2008).
［６］M.A. Beer et al., Phys. Plasmas 2, 2687 (1995).
［７］T.-H. Watanabe et al., Phys. Rev. Lett. 100, 195002 (2008).
［８］S. Maeyama et al., submitted to Int. J. High Perform Comput. Appl. (2012).
［９］Y. Morinishi et al., J. Comput. Phys. 143, 90 (1998); J. Comput. Phys. 197, 686 (2004).
［１０］S. Satake et al., Plasma Fusion Res. 3, S1062 (2008).
［１１］Y. Saad, Iterative Methods for Sparse Linear Systems, Second
Edition (SIAM, 2003).
［１２］S. Sayantan et al., in Proceedings of PPoPP06 (2006).
［１３］Y. Idomura et al., submitted to Int. J. High Performanceb*(f(i+2,j,k,l)-f(i-2,j,k,l))fy=a*(f(i,j+1,k,l)-f(i,j-1,k,l))+ &b*(f(i,j+2,k,l)-f(i,j-2,k,l))
fz=a*(f(i,j,k+1,l)-f(i,j,k-1,l))+ &b*(f(i,j,k+2,l)-f(i,j,k-2,l))fv=a*(f(i,j,k,l+1)-f(i,j,k,l-1))+ &b*(f(i,j,k,l+2)-f(i,j,k,l-2))df(i,j,k,l)=f(i,j,k,l)+ &c*(fx+fy+fz+fv)
enddo
enddo
enddo
enddo
前 の ケ ー ス で は OpenMP（あ る い は，多 く の 自 動 並 列Comput Appl. (2012).
［１４］S. Williams et al., Commun. ACM 52, 65 (2009).

化）でデフォルトのブロック分割（%$１，２，…，５として$ 番目のスレッドに #$%$!!&#""%を割り当て）となっていたため，各スレッドで全く独立な領域を計算していたひこ

渡 邉 智 彦
気がつけば，シミュレーションをやるようになってからの人生の方が，その前よりも長くなりました．次はどんなスーパーコン
ピュータを使えるようになるのか楽しみに，そろそろ，また新しい研究を始めたいと思っています．

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!

!!!!!!!!!!!!

とも

い

ど

むら

やす

ひろ

井 戸 村 泰 宏

日本原子力研究開発機構システム計算科学センター所属．ジャイロ運動論シミュレー
ションによる乱流輸送研究，および，関連する HPC 技術の開発に取り組んでいます．
この分野に従事して１０年以上になりますが，計算機の目覚ま
しい進歩にいつも驚かされます．シミュレーション手法やモデルも計算機パワーに負けないように高度化できればと考えています．

!!!!!!!!!!!!!!!!!!!!!!!!!!

なべ!!!!!!!!!!!!!!!!

わた!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

１７９!!!!!!!!!!!!!!!!!!!!!!!!!!