J. Plasma Fusion Res. Vol.89, No.2 (2013) 119‐127

講座 核融合プラズマシミュレーションの技法−大規模並列計算環境の活用−

３．MHD シミュレーションのコーディング技法
三 浦 英 昭，藤 堂

泰，後 藤 俊 幸１）

核融合科学研究所，１）名古屋工業大学（原稿受付：２
０１
２年１月６日）

この章では，MHD 方程式のシミュレーションを行うにあたっての並列化のポイントについて簡単に紹介する．取り上げる話題は，差分法によるコードの並列化と，３次元高速フーリエ変換の高速化である．
Keywords:
MHD, finite difference, FFT, parallelization

３．
１ はじめに

値シミュレーションなどに用いられる．
（直接数値シミュ

MHD シミュレーションの一番簡単な形は，速度や電場

レーションとは，支配方程式を人工粘性などに頼らず極力などの物理量を空間上の離散点（格子点）の集合として表方程式を誠実に解く場合に用いられる用語であり，必然的

現し，空間微分をなんらかの形で近似，これを時間方向にに数値シミュレーション手法には高精度・高解像度である

積分することで数値的な近似解を得るものである．このたことが要請される．
）また，トーラスプラズマのトロイダ

め，数値シミュレーションコードの作成の際には，如何にル方向およびポロイダル方向もフーリエ変換可能であるた微分を近似するか，その近似を効率よく行うためにデータ

め，動半径方向は差分法で，トロイダル・ポロイダル方向をどのように主記憶に配置するかが重要となる．最も古典はスペクトル法を使うこともある．ここで重要な点は，的な（しかし今でも頻繁に使われる）数値解法として中心

フーリエ変換の演算には格子線に沿ったすべての格子点上差分法がある．最も粗い近似の中心差分法である２次精度のデータを使わざるを得ないということである．フーリエ中心差分法を例にとると，ある物理量 %(&の空間微分を空%変換を行う方向のデータが複数の計算ノードに分散して置間上の格子点 (!
!
!
"#'（!，" はそ!""&&##&，&#!
&#%

かれている場合にはノード間通信が発生する．３．
３節では，れぞれ空間の全長と，格子点の総数）でこのようなノード間通信の発生を念頭に，３次元高速フー
リエ変換（Fast Fourier Transform，FFT）における並列化

!
$%!
""!%
!"
&
&
!# %
##&
"#%
#
$(!
#
&

の方法について述べる．また，近年は優れた FFT ライブラ
リーのソースコードが公開されていることから，これにつと近似する．このような差分法の特徴は，考えている格子

いても同節で触れる．

点($"の情報があれば微分が近似できる&とその隣接格子 (&３．
２ MPI 並列化のための領域分割点である．（より高次の差分法では必要とする格子点数が増加する．）このため，格子点データが複数の計算ノード

３．
１で述べたように，計算法として差分法を用いて，MPIにまたがっていても，２．
２節で紹介した MPI ライブラリに通信を伴う並列化 MHD シミュレーションを行うことを考よるノード間通信を行って比較的簡単に計算が可能であえる．MPI プロセスと分割領域は一対一に対応し，各 MPI

る．このような場合についての並列化を３．
２節で説明する．

プロセスが一つの分割領域について同一のプログラムを実中心差分法が格子点 (&の周辺の情報だけを利用するの行するものとする．このとき境界値計算のためのデータ交に対して，空間全体の情報を利用する手法もある．その代換等のため他のプロセスとの通信が必要となるが，差分法

表的な例が，フーリエ変換によるスペクトル法である．ス

で参照されるデータは隣接する格子点に限定されているのペクトル法では，物理量 %(&の空間微分をそのフーリエ係%%数%' を用いてで，プロセス間の通信は基本的に１対１通信である．並列化によって短縮される計算時間と比べて通信時間が十分に

! ""#
$%!
&
'
(
!# " &
%
'
%
'$
$(!
#
& '#!""

小さければ，プログラム全体の実行時間を短縮することが&できる．（このような差分法のプログラムでは，プログラ
ム作成者やコンピュータのアーキテクチャーが異なってもと近似する．この手法は空間に対する周期性を利用する，

領域分割の様式は比較的似通った形になりやすい．これ

空間３次元が周期的な一様乱流，一様 MHD 乱流の直接数

は，プログラム作成者の考えやコンピュータのアーキテク3. Coding Techniques of MHD Simulations
MIURA Hideaki, TODO Yasushi and GOTOH Toshiyuki

corresponding author’s e-mail: miura.hideaki@nifs.ac.jp

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

Journal of Plasma and Fusion Research Vol.89, No.2 February 2013

チャー次第で デ ータ の 配 置 が 大 き く 変 わ り 得 る３次元

セスが担当する計算量は，正味の格子点数 lx*ly=(lx_FFT（次節）と対照的である．
）

global*ly_global)/nprocess に比例するので，全体

ここではプログラム言語として Fortran を使用する．モのプロセス数 nprocess に反比例して減少する．一方で通ジュール mpi を使用し（use_mpi），その中で定義されて信量は，上下左右のプロセスとデータを交換するので，

いな い 変 数（以 下 の 例 で は mpi_err，nprocess，my_

2*m*(lx+ly)=2*m*(lx_global/nprocess_x + ly_

rank など）は特に断らない限り整数型（integer）で宣言す

global/nprocess_y)に 比 例 す る．通 信 量 はるものとする．MPI 並列計算では，まず最初にnprocess_x と nprocess_y の積である nprocess ではなく，nprocess_xまたはnprocess_yに反比例して減少

call mpi_init (mpi_err)

するのである．通信量と計算量の比，（通信量）／（計算量）
は 2*m*(1/lx + 1/ly)に比例することになる．実際の計を実行し，計算の最後には算における通信時間と計算時間の比は，通信量と計算量のcall mpi_finalize（mpi_err）

比に加えて，計算機の演算性能および通信性能にも依存するが，lx と ly が十分に大きければ通信時間は無視できて，を実行しなければならない．
次に，MPI 並列計算では全プロセスが同一のプログラム

全体の計算時間（＝計算時間＋通信時間）を nprocessを実行するので，各プロセスは自分が何番目のプロセスに反比例して短縮することが可能となる．このため，差分（ランク，以下の例では my_rank）であるかをプログラム

法は大規模並列計算に適しているといえる．領域分割数が側で認識する必要がある．ランクと計算全体のプロセス数増大して lx または ly が小さくなると通信時間が計算時間（以下の例では nprocess）は，以下で与えられる．

と比較して無視できなくなり，プロセス数を増やしても期
待されるほど全体の計算時間は短縮されなくなる．また，

call mpi_comm_rank(mpi_comm_world,

nprocess 一定の条件下で（通信量）／（計算量）を極小化す

my_rank,mpi_err)

るには（同時にlx*ly=lx_global*ly_global/nproc-

call mpi_comm_size(mpi_comm_world,

ess も一定であることに注意），1/lx + 1/ly=(lx + ly)

nprocess,mpi_err)/(lx*ly)であるから，lx+ly を極小化すればよいことが

ここで，my_rank がとりうる値の範囲は０から nprocess−1

わかる．２次元の問題では lx=ly となる正方形に近い分割である．プログラム側ではコミュニケータ mpi_comm_world

領域が高速計算には有利であるといえる．

を定義する必要はない．コミュニケータは各通信に関係す

次に考慮しなければならないのは，分割領域格子と全体るプロセスを指定しており，mpi_comm_world はこの通信にランク０から nprocess−1 のすべてのプロセスが関格子の対応である．まずは各分割領域が x，y 方向に何番目_ my_rank
かを my_rank_x，my_rank_yで定義しよう(0 <係することを意味している．コミュニケータとして mpi_

_ nprocess_x−1,0 <_ my_rank_y <_ nprocess_y−1)．
_x <comm_world を使用すればあらゆる通信を行うことがで

ここでは分割領域を x 方向に左から順番に並べ，右端に達きるが，通信の効率を考えるとプログラム側でコミュニ

すると次は y 方向に一つ移動して左から順番に並べるもの

ケータを定義して対象範囲を必要なプロセスに限定する方とする．こ の と き my_rank = my_rank_x + my_rank_yがよい．その例については後述する．

*nprocess_xの関係が成り立ち，Fortranプログラムでは

さて，ここでは２次元座標系!
!
!
""の問題を考えよう．全my_rank_x，my_rank_y は以下のように計算できる．

体でそれぞれ lx_global，ly_global 個の格子点を使用

my_rank_y = my_rank/nprocess_x

し，それぞれを nprocess_x，nprocess_y 個の領域に均

my_rank_x = my_rank − my_rank_y*nprocess_x

等に分割するものとする．一つの領域が一つのプロセスに対 応 す る の で，nprocess=nprocess_x*nprocess_y

x 方向の境界条件は周期境界としてみよう．このとき物の関係がある．各プロセスが分担する分割領域の格子点数

理的な境界条件として，全体格子上の配列 G が G(1,:)=( lx, ly ) は lx = lx _ global / nprocess _ x ， ly = ly _G(lx_global+1,:)を満たし（全体格子上に実際に配列

global/nprocess_yである．差分法では近隣の格子点の情報を参照して差分計算を行う．そのため分割領域の端点を定義する必要はない），各分割領域のx方向にi番目の格_i<_ lx+2*m)全体では iglobal 番目の格子点
子点が(1 <

では差分計算に参照する格子点が近隣プロセスに存在するであるとする．my_rank_x=0 の分割領域の x 方向に m+1

ので，差分計算に必要なデータをプロセス間の通信により番目の格子点が全体では１番目の格子点に対応する場合，交換しなければならない．このため，配列を正味の格子点

両者の対応は iglobal=i+my_rank_x*lx−m で与えられだけでなく，差分計算に必要な近隣プロセス上の格子点もる．例えば，my_rank_x=0 の分割領域では iglobal=i

含めて定義しておくとよい．このように差分計算のために−m となる．
次 に y 方 向 に は 周 期 境 界 で な く，全 体 格 子 上 の拡張した領域は袖領域またはのりしろ領域と呼ばれる．両側 m 点ずつ参照して差分計算を行うものとすると，各プロ

jglobal=1 と jglobal=ly_global の物理量はあらかじ

セスの格子点数は袖領域を含めて(lx+2*m,ly+2*m)となめ与えられるものとしよう．y 方向についても x 方向と同る．

様に分割領域格子と全体格子を対応させるものとすると，
そ れ ぞ れ の 格 子 番 号 j と jglobal は jglobal=j+my_

ここで計算量と通信量の関係を整理しておこう．各プロ
１２０

Lecture Note

3. Coding Techniques of MHD Simulations

H. Miura et al.

lystart=1，lyend=ly，（２）my_rank_y=nprocess_y−1 の 領 域 で は，jglobal=j+my_rank_y*ly−2*m,lystart=2*m+1，lyend=ly+2*m，（３）これら以外の領域では，jglobal=j+my_rank_y*ly−m，lystart=m+1，
lyend=ly+m とそれぞれ設定される．この設定は，全体格子上で開発された非 MPI 並列差分計算プログラムを MPI並列化する際にも役立つ．
図２に示した袖領域の通信では，ユーザーが新たに定義するコミュニケータを使用して通信の対象を限定するとよい．ここでは mpi_comm_split を用いてコミュニケータ
を定義する．x方向の通信ではmy_rank_yの値が等しいプ
ロセス間でデータを交換するので，コミュニケータ mpi_
comm_x を以下のように定義する．
icolor = my_rank_yikey = my_rank_x
call

mpi _ comm _ split ( mpi _ comm _ world,

icolor,ikey,mpi_comm_x,mpi_err)
図１

ここで，icolor はこの値が同じプロセス間で通信するこ分割領域（my_rank_x=0, my_rank_y=0）の格子（i，j）と全体格子（iglobal, jglobal）の対応．黒丸と白丸は内部領域と袖領域の格子点をそれぞれ示している．灰色で塗られた部分は各分割領域の内部領域に対応し，（ ）内の数
字は（my_rank_x, my_rank_y）を示している．

とを意味し，この通信グループ内のランクは ikey の順番に０から定められる．y 方向通信のコミュニケータ mpi_
comm_y は，icolor = my_rank_x，ikey = my_rank_yとして同様に定義される．このように新しく定義したコ

rank_y*ly−m で関係づけられる．分割領域格子と全体格

ミュニケータを図３に示す．

子の対応を図１に示す．m が２以上の場合には，jglobal

次に袖領域の通信について説明する．

=2,…,m お よ び jglobal=ly_global−m+1,…,ly_real(8)::aaa(lx+2*m,ly+2*m)

global−1 の格子点では境界に近い側で m 個の格子点を参照することができず，他の格子点とは異なる差分法を適用と宣言されている配列 aaa の通信を考える．y 方向の通信することになる．これらの全体格子点に対応する格子点をは，対象となるデータがメモリ上に連続に格納されている含む my_rank_y=0 および nprocess_y−1 のプロセスでので簡単である．この場合は上方と下方の隣接プロセスとは，境界付近の専用差分法を実行する必要がある．

の通信が発生する．my_rank_y=0 または my_rank_y=

分割領域のうち袖領域を除いた正味の格子の x 方向の最nprocess_y−1の場合は，それぞれ下方または上方の隣接

初と最後の点をそれぞれ lxstart，lxend，y 方向にはそプロセスとの通信が発生しないので，通信先として mpi_

れぞれ lystart と lyend とする．ここまでに述べた分割

proc_null を指定して通信を抑制する．非同期型通信

領域格子と全体格子の対応を用いると，lxstart=m+1，

mpi_isendとmpi_irecvを使うとデッドロックの心配がlxend=lx+m，lystart=m+1，lyend=ly+m で あ る．x

ない．

方向に隣接する領域間（my_rank_x および my_rank_x+1）

node_up = my_rank_y + 1のデータ交換は図２のように行われる．y 方向についても

node_down = my_rank_y − 1

同様である．

if(my_rank_y.eq.0) node_down =ここまでに説明した分割領域格子と全体格子の対応で

mpi_proc_null

は，my_rank_y=0 および nprocess_y−1 のプロセスにおif(my_rank_y.eq.nprocess_y−1) node_up =いて，境界付近の専用差分法を実行する必要があることをすでに述べた．ランクによる計算の分岐があるとプログラ
ムが複雑になり，バグが発生しやすくなる．MHD シミュ
レーションでは差分計算を頻繁に行うので，差分計算はランクに依存した計算の分岐を含まないことが望ましい．そこで，各プロ セ ス で j=2,…,m お よ び j=ly−m+1,…,
ly−1 の格子点では共通の専用差分法を適用することを考える．この場合には，my_rank_y=0 の領域において全体
格 子 と 分 割 領 域 格 子 の 下 端 が 一 致 し，my_rank_y=

図２nprocess_y−1 の領域では上端が一致しなければならない．
具体的には，（１）my_rank_y=0 の領 域では，jglobal=j，
１２１x方向に隣接する領域間（my_rank_xおよび my_rank_x+1）のデータ交換．黒丸と白丸は内部領域と袖領域，矢印の向きはデータの通信方向をそれぞれ示している．

Journal of Plasma and Fusion Research Vol.89, No.2 February 2013

real(8)::up_send_x(m,ly+2*m),up_recv_x(m,ly+2*m)
real(8)::down_send_x(m,ly+2*m),down_recv_x(m,ly+2*m)
隣接プロセスのランクは，周期境界条件を考慮して
node_up = mod(my_rank_x + 1,nprocess_x)
node_down = mod(my_rank_x − 1 + nprocess_x,
nprocess_x)と定義される．次に送信用バッファ配列にデータを格納する．
do j = 1,ly+2*m__up_send_x(1:m,j) = aaa(lxend−m+1:lxend,j)__down_send_x(1:m,j) = aaa(lstart:lxstart+m−1,j)end do通信は以下のようになる．
図３

新 し く 定 義 し た コ ミ ュ ニ ケ ー タ（a）
mpi_comm_x と（b）
mpi_comm_y．正方形と長方形は MPI プロセスと通信グ
ループをそれぞれ表している．各 MPI プロセスの中心の数
字 は 新 コ ミ ュ ニ ケ ー タ で の ラ ン ク，（ ）内の数字は
mpi_comm_world における（my_rank_x, my_rank_y）を示している．

isize = (ly+2*m)*m
call mpi_isend(up_send_x(1,1),isize,
mpi_double_precision,&
node_up,1,mpi_comm_x,
isend1,mpi_err)
call

mpi_proc_null

mpi_isend(down_send_x(1,1),isize,
mpi_double_precision,&

jsize =(lx+2*m)*m
call

node _ down, 1, mpi _ comm _ x,

mpi_isend(aaa(1,lyend−m+1),jsize,

isend2,mpi_err)

mpi_double_precision,&call mpi_irecv(up_recv_x(1,1),isize,

node _ up, 1, mpi _ comm _ y,

mpi_double_precision,&

jsend1,mpi_err)

node _ up, 1, mpi _ comm _ x,call mpi_isend(aaa(1,lystart),jsize,

irecv1,mpi_err)

mpi_double_precision,&call

mpi_irecv(down_recv_x(1,1),isize,

node_down,1,mpi_comm_y,

mpi_double_precision,&

jsend2,mpi_err)

node _ down, 1, mpi _ comm _ x,call mpi_irecv(aaa(1,lyend+1),jsize,

call

irecv2,mpi_err)

mpi_double_precision,&call mpi_wait(isend1,my_status,mpi_err)

node _ up, 1, mpi _ comm _ y,call mpi_wait(isend2,my_status,mpi_err)

jrecv1,mpi_err)

call mpi_wait(irecv1,my_status,mpi_err)

mpi_irecv(aaa(1,lystart−m),jsize,

call mpi_wait(irecv2,my_status,mpi_err)

mpi_double_precision,&

最後に，バッファ配列に受け取ったデータを元の配列に格

node_down,1,mpi_comm_y,

納する．

jrec2,mpi_err)
call mpi_wait(jsend1,my_status,mpi_err)

do j = 1,ly+2*m

call mpi_wait(jsend2,my_status,mpi_err)__aaa(lxend+1:lxend+m,j) = up_recv_x(1:m,call mpi_wait(jrecv1,my_status,mpi_err)j)call mpi_wait(jrecv2,my_status,mpi_err)__aaa(lxstart−m:lxstart−1,j)= down_recv_
x(1:m,j)x 方向の通信では参照されるデータがメモリ上に連続にend do

格納されていないため，バッファ配列を使用する．以下のバッファ配列を宣言しておく．

MPI 通信では毎回オーバーヘッドが発生するので，複数
１２２Lecture Note

3. Coding Techniques of MHD Simulations

H. Miura et al.

の配列について通信を行う場合は，以下で宣言されるよう散して配置されているものとする．ただし，ここでは２次なバッファ配列を用いて通信を１回にまとめるとよい．こ

元領域分割を考えるので，一つの次元については全データ

こで narray は通信する配列の個数である．

が一つのノード内で共有されている．FFT では座標線方向にすべてのデータを使って計算するため，データが複数のreal(8)::up_send_xm(m,ly+2*m,narray),up_

ノードにまたがっている場合には MPI による通信を行う

recv_xm(m,ly+2*m,narray)必要がある（MPI 通信の初期化等については３．
２と同様な

real(8)::down_send_xm(m,ly+2*m,narray),のでそちらを参考されたい．
）．FFT の煩瑣なバタフライ

down_recv_xm(m,ly+2*m,narray)

演算を複数のノードにまたがったデータに対して行うのは相当に大変であるため，次のような手順を考える．最初に，具体的な手順は配列１個の場合と同様なので割愛する．

全区間データが１つのノードに格納されている座標方向に

３．
３ ３次元 FFT

ついて FFT によるフーリエ変換を行う．次に他の方向の次に，スペクトル法で用いられる３次元FFTの並列化に

フーリエ変換を行う際にはデータの転置を行い，一旦全区ついて述べる．多くの計算機センターのスパコンにはすで間データをすべて同一ノードに格納してから FFT を実行にチューニングされたFFTが備わっている．並列計算機向する．このため，３．
２の差分法とは異なり，おおがかりな

けのFFTとして有名なものにはFFTW［１］とそれをもとに

MPI 通信が複数回発生する．この結果，大規模分散並列計

した P3DFFT［２］があり（これらについては最後に触れる

算では時として，全計算時間の５０％以上が３次元FFTの通ことにする），いまさら私家版の３次元 FFT を開発すると信時間という事態が発生し得る．以下では，このような実いうのはいささか世間の流れに抗っている感がないわけで態を踏まえて転置の部分を主眼に解説する．誌面の都合はない．しかし，FFT ではそのデータ配置こそが高速化の上，バタフライ演算など FFT の細部については割愛する．

観点からも使い勝手の観点からも最も重要なポイントの一

適宜参考書や論文を参考にされたい．

つである．レディーメイドの服はファッショナブルで安価，
多 く の 複 素 表 現 の FFT で は 関 数 値 "
'
(
$%"
$)

ですぐ手に入るが，個人の体型やサイズに必ずしもフィッ

$%"
"
!
!
!
"
! に対して変換，逆変換を返すように作られていトしているわけではなく我慢しなければならないこともあるので，波数空間での $&をかけるという数学演算は単純でる．一方，オーダーメードの服は高価だが，体型にフィッはなく K(i) = (i − 1) − [(i − 1)/(N/2)]*N のように修トし使い勝手がよくそれだけに長い期間愛用することがで正された波数をかける操作になる．ここでは数学的表現ときる（著者は誂えた服を着たことはないが）．同じことが，プログラム上での表現が同等になるように複素フーリエ変各人の必要とする３次元 FFT や他の基本ルーチンにもあ換を構成する．実空間の離散点 (
$#
!)
'
'
'
#
(
""
#"
$)での関数(（'値を "(，フーリエ係数を "(
'
'
'
&#"
&$)
&""
$， &
$
""
#"
$)

てはまると思う．

は整数）とすると３次元フーリエ級数を以下のように定義

ここで紹介する３次元 FFT は，乱流および MHD 乱流のする．

スペクトル法による直接数値シミュレーションに向けて開発されたものである．高レイノルズ数での乱流では高い空

!#
#!"

"(
$
"
%
"
&)
% #

間解像度が要求されるので，１方向あたりの格子点数 N

!#
#!"

#

!#
#!"

#

((
"
&#"
&$)
&""

&"%!!## &"%!!## &"%!!##

１０

が２ 以上の３次元 FFT を念頭に置いている．また，乱流$ &'##%'&!$ ("&'"&')"! "" ## $$理論との関連性やスペクトル法では波数空間で方程式を取り扱うことが多いので，以下では波数空間を中心にした解$説を行う．プログラムは Fortran で書かれているものと想((((実数条件より "
&#"
&$)
&""
%"
!&""!&#"!&$)であるか定する．後述の公開された FFT ライブラリーは C 言語など

ら，メモリ消費を抑えるためにこの特性を用いて波数空間の半分の領域のみを用いる（図４左図）．例えば，&$'!ので書かれているものが多いが，多くの場合は Fortran プロ((
領域のみで "&""&#"&$)を定義する（&$%!では実数条件

グラムからサブルーチンをコールできるように配慮されてを 課 し た 冗 長 な 配 列 を 用 い る）．一 方，物 理 空 間 で の

いる．

#&'
!!#
$
"
%
"
&)
$&!"に お け る 実 数 関 数 配 列 #(２と同様に各ノードに分
関数 "(
'
'
"
'
""
#"
$)のデータは３．

図４

３次元 FFT における変数の格納の仕方．

１２３

Journal of Plasma and Fusion Research Vol.89, No.2 February 2013%%
は，複 素 数 宣 言 さ れ た 配 列 "$!
%!
&&に 対 し て
%%
#!"）
#%
$
!
%
!
&!! &（ !$&$!#
"
$
!
%
!
&&
$
!
%
!
&&
##%
"'!"
#
という具合に，#%$!
%!
&&のうち &が負の領域をそのまま%%!##だけ第３変数の正の方向に平行移動して， "
$
!
%
!
&&
の虚数部に格納する（図４右図）．３番目の変数についての変換は後で述べる転置を行って第２変数の位置に移動させることにより，長さ ! の領域を確保したうえで ("!の領域に展開して変換を行うようにする．以下，順にみていく．
１．関数は倍精度複素変数とする．配列サイズが大きいので領域を分割する．大事なのは，フーリエ変換（三角関数

図５と変数との積と和によるバタフライ演算が主要部）を行う

３次元 FFT における空間分割．

変数の方向にはプロセス間通信を入れないようにすることは以下のようになる．

である．第１変数については疑似ベクトル処理機能を有効に働かせるようデータの連続アクセスをはかり，キャッ
シュの大きさに応じてベクトル長を最適なものにするよう!$OMP PARALLEL DO private(i,j,k,....)

分割数 Npx を決める．第３変数については Npz 個のプロセ<−− k

ス並列に加えてスレッド並列を行う．格子点数が少ない場

軸方向スレッド並列

do k = 0,(N/2)/Npz−1

合には第１変数の大きさを N にとってフーリエ変換を行う

do j=1,n<−− j 軸方向のフーリエ変換におことも可能である．しかし，６４kB のキャッシュでは格子点

ける n 段のバタフライ演算ルー

サイズが１０２４を超えるとキャッシュに入りきらなくなるのプ（!#$' に注意）

で，第２変数の方向に ! 個の格子点をとって FFT を行う

j1=..; j2=..; j3=..; j4=..

ようにする．したがって，第１，第３変数について２次元分<−− バタフライ演算で用いられる割を行う（図５）．各変数の動く範囲は，0≦i≦ N/Npx−1，

４個の点の配列インデックス−N/2≦j≦ N/2−1,0≦k≦ N/Npz−1 となる．直感的には，

trig2=..; trig3=..; trig4=.. <−− 三角関数

第１，３変数の組で構成される(N/Npx)×(N/Npz)の面に

do i=0,(N/Npx)−1<−− i 軸方向連続アクセス

ついて第２変数方向に FFT をかけるというイメージにな

w1 = f(i,j1,k)+f(i,j3,k)る．分割にあたっては，２種類のコミュニケータを必要と

w2 = f(i,j2,k)+f(i,j4,k)するので，mpi_comm_split を用いて，

w3 = f(i,j1,k)−f(i,j3,k)
w4 = f(i,j2,k)−f(i,j4,k)f(i,j1,k) = w1+ w2!!!−−− Split MPI Rank (Type A)

f(i,j2,k) = (w3−ai*w4)*trig2

color = mod(rank,Npx) ; key = rank/Npx<−−ai=(0.0,1.0)虚数単位

call MPI_comm_free(Ncomm,err)
call

f(i,j3,k) = (w1− w2)*trig3

MPI _ comm _ split ( comm, color, key,f(i,j4,k) = (w3+ai*w4)*trig4

Ncomm,err)

end do; end docall MPI_comm_rank(Ncomm,Nrank,err)

end docall MPI_comm_size(Ncomm,Nnpe ,err)!$OMP end parallel do'のように行う．変換の基数が４で !#$ の時には，FFT

図６

３次元 FFT におけるブロック間転置とブロック内転置（左図），およびバタフライ演算の方向と２回の転置（右図）．

１２４

Lecture Note

3. Coding Techniques of MHD Simulations

キャッシュを有効に使うため，同じ配列 f(i,j1,k)につ

H. Miura et al.

←---- 荷ほどき!!!−−− Unpack

いてロードとストアを行うようにする．このため４，５行

do iz=0,Nklen1

目の配列変数と回転因子（三角関数表）をうまくならべる必_do kx=0,Npx−1; kxg = k1+kx*Nilen

要があるが，この詳細は１次元 FFT の解説書にゆだねる．

___do iy=0,Nilen1____do ix=0,Nilen1;

２．第１，
３変数についての変換には，第１，３変数を第２______ixp=kxg+ix; nl=Nilen*Nilen*iz

変数の位置に持ってくるように転置が必要になる．転置は+Nilen*iy+ix

ブロック転置とブロック内転置からなり，通信が発生する______c(iy,ixp,iz)=c2(nl,kx)←---- 局所転置も行うのはブロック転置においてである．転送するデータの塊を大きくかつ変数が連続になるようにして，できるだけ MPI

end do; end do; end do ; end do

通信の回数を減らすことが望ましい．そのために，例えdeallocate(c2)ば，２次 元 配 列 に 詰 め 替 え た 後（packing），全 対 全（MPI_All_to_Allv）通信を用いてデータを交換する（図

データを密にパックして送り出すやり方でかなり計算時間６）．ブロック転置が済んだら，２次元配列をまた３次元配が変わってくる．上のプログラム例では，わかりやすくす

列に 戻 す（unpacking）と 同 時 に ブ ロ ック 内 転 置 を 行るために２次元配列を用いたが，mpi_type_vector など

う．プログラムは以下のようになる．

を使って新しい変数を定義してパックせず直接送る仕方もある．パッキングをするかしないかは送るデータのサイズ!!!−−− packNklen1=(N/2)/Npz−1;←---- 詰め替えにも依存する．もしパッキングするのであれば，バタフラ

Nilen=N/Npx;

イ演算（３行目）の最終段（j=n）の時に，ストアと同時にNilen1=Nilen−1

パックしてしまう方法もあるが，ストアのコストも考慮すbtsize2 =Nilen*Nilen*Nklen

る必要がある．

←---- １ブロックあたりの３．配列の転置にかかる全対全の通信に多くの時間がかか

転送データサイズ
do i=0,Npx−1

るので，転置の回数は最小限（２回）にとどめる．図７に___msize2(i) = btsize2示した例では，波数空間では !
!!!
!"!
!#"の順に並んでいた____disp2(i) = i*(btsize2+1)

ものが，実空間では !
"
"
"
"!
#!
!"の順に並ぶので注意が必要

end doallocate( b2(0:btsize2,0:Npx−1) )do iz=0,Nklen1___do ky=0,Npx−1; kyg=−N/2+ky*Nilen____do iy=0,Nilen1; iyp=kyg+iy_______do

ix=0,Nilen1;

nl=Nilen*

Nilen*iz+Nilen*iy+ix__________b2(nl,ky)=c(ix,iyp,iz)←---- ２次元配列にまとめる
end do; end do; end do ; end do
deallocate(c)allocate( c2(0:btsize2,0:Npx−1) )図７

３次元 FFT における２回の転置による独立変数の順序．

!!!−−− Block Transpose in x and y Directions←---- ブロック転置
call

全対全通信

MPI_alltoallv(b2(0,0),msize2,
disp 2, MPI _ complex16,&
c2(0,0),msize2,
disp 2, MPI _ complex16,&
Ncomm,err )deallocate(b2)allocate(c(0:Nilen,k1:k2,0:Nklen1))
図８

１２５

３次元 FFT のスケーラビリティ．

Journal of Plasma and Fusion Research Vol.89, No.2 February 2013

である．順序は転置の順番をどうとるかによって変化する

触れたい．最初に述べたように，有名な FFT ライブラリとので，都合の良いように入れ替えればよい．たとえば，実しては FFTW［１］と，これを利用した P３DFFT［２］があ空間での方程式を扱うために "
"
"
"
!!
"!
##の並びとしたけれ

る．FFTW は C 言語で記述されているが，Fortran のためば，波数空間で "!#!
!!!
!"#の並びになるように入れ替えれのインターフェースが用意されている．また，スレッド並ばよい．

列および MPI を用いた並列化のための関数が用意されてい る．P３DFFT は FFTW な ど を 利 用 し た３次 元 FFT

４．レイテンシーの問題を回避しメモリアクセス回数を逓である．本章で紹介した３次元FFTと同様に２次元の領域減するためFFTの基数を大きいものにとる．上のプログラ

分割が行われているため，メニーコアのスーパーコン

ム例では基数を４とした．キャッシュサイズをみながら，

ピュータで高効率の３次元 FFT が実行可能となっており，

領域分割数（Npx）を多めにとり基数８の FFT を導入する

乱流の大規模３次元シミュレーションなどに利用されていことも考えられよう．

る．この他の２次元領域分割型 FFT については［３］が参考になる．このウェブサイトには FFTW 以外の公開 FFT

２次元分割における分割パターン（Npx，Npz）によってライブラリや，２次元領域分割に関する説明もあり，公開

も FFT の性能に違いが出てくる．図８は，Fujitsu FX１にライブラリの利用者に親切な作りになっている．ここに紹おいて計算時間がプロセス数とともにどれほど減少するか介されている多くの３次元 FFT は２次元領域分割を行うを各格子点数についてみたものであり，曲線の傾きが−１

が，FFTW-3D［４］のように３次元領域分割を行うものもあに近いほど理想的なスケーラビリティとなる．プロセス数

る．（著者はまだこの３次元領域分割ライブラリを使って（Npx×Npz）の増加とともに計算時間が減少していること

いないが．）ウェブサイト［３］では紹介されていない公開ラが見て取れる．しかし，プロセス数がある程度大きくなる

イブラリとしては，たとえば CREST 事業で開発された国と通信のオーバーヘッドが大きくなり性能は低下する．た産 ラ イ ブ ラ リ FFTSS［５］が 上 げ ら れ る．FFTSS で はだし，このデータを取った時点では１プロセス１スレッド

FFTW３と イ ン タ ー フ ェ イ ス レ ベ ル で 互 換 性 が あ り，（Flat MPI）の計算を行っていた．これは，MPI 並列と自動FFTW からの移植が容易であるとされている．この他にも並列によるハイブリッド並列を行うと，前者の場合より約多数のFFTが公開されているので，興味を持った読者は検１０％から２０％程度，計算時間 が 長 く な る た め で あった．

索してみるとよい．

MPI 並列と Open MP によるハイブリッド並列を用いると，トーラスプラズマの MHD シミュレーションではトロイ

Flat MPI 並列と同等の計算時間となった．

ダル方向，ポロイダル方向にはフーリエ変換を行うが，動３次元FFTのチューニングは未だ発展途上であり，上に径方向には差分法などを使う場合，あるいはトロイダル方

述べた以外にも高速化，高効率化の手法が考えられる．た向にはフーリエ変換を使い，ポロイダル断面上では有限要とえば（a）物理変数が実数の時には，複素数表現の実数部素法を使う場合が多い．これらの場合には，FFTW などのと虚数部を使って２つの関数を同時に変換する方法，（b）ライブラリから一次元 FFT を用いて必要な方向にフーリ

ある !#面で FFT の計算を行っているときにすでに FFT

エ変換を行うと良い．フーリエ変換を行う次元についてはが完了した !#!!の面内の計算結果を転送できれば，通信

共有メモリー内でのスレッド並列化を行う一方で，動径方時間を少なくすることができる（通信の隠ぺい）．

向には差分法による離散化を行い，MPI による並列化を行

この FFT は３次元の流体力学乱流や MHD乱流のシミュ

うのが一番簡単であろう．差分法を用いた場合の MPI によレーション向けに開発されたものであるが，乱流の普遍的る並列化については前節を参考にされたい．

統計法則は高いレイノルズ数において現れると考えられており，この状態を計算機上で実現するためにはより高解像

参考文献

度の計算が必要である．流体力学の世界では，格子点数［１］http://www.fftw.org/［２］http://www.sdsc.edu/us/resources/p3dfft/;
http://code.google.com/p/p3dfft/［３］http://www.2decomp.org/［４］http://www.gromacs.org/Developer̲Zone/Programming̲Guide/FFTW-3D/［５］http://www.ssisc.org/fftss/index.html.ja.euc-jp

２０４８３での計算は今や当たり前になってきており，格子点数４０９６３での計算も特定のグループで行われ，８
１９２３での計算にも手が届きつつある．このような大規模な計算が可能な３次元 FFT の開発は乱流研究だけでなくあらゆる分野
で必須である．
最後に，公開されているFFTライブラリの利用について

１２６

3. Coding Techniques of MHD Simulations

ひであき

三 浦 英 昭核融合科学研究所ヘリカル研究部所属．流体力学・一様等方性流体の乱流シミュレーション，渦同定研究などを経て，現在は高温プラズマ・不安定性シミュレーションの研究に従事する．

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

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

ご

とう

とし

ゆき

後 藤 俊 幸名古屋工業大学大学院，研究分野：乱流理論と乱流の計算科学．福井県出身．荷風にならって，街中をぶらぶら歩くことが好きです．それと，近くの日帰り温泉につかってのほほんとしています．

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

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

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

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

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

うら

どう

やすし

藤 堂

とう

泰

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

み!!!!!!!!!!!!

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

H. Miura et al.

総合研究大学院大学物理科学研究科核融合科学専攻および核融合科学研究所・教授．
流体シミュレーションと粒子シミュレー
ションを連結して，核融合プラズマにおいて高エネルギー粒子が不安定化するアルフ
ベン固有モードとアルフベン固有モードによる高エネルギー
粒子輸送を研究している．この研究を物理として発展させた
い．学生募集中．

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

Lecture Note!!!!!!!!!!!!!!!!!!!!!!!!!!

１２７