νＧＣに関するN体計算の解説と展望
矢作 日出樹Introduction
N-body Methods
Adaptive Mesh Refinement (AMR)
Optimization of the AMR code
Why should I optimize the code?
How was the code changed?

N-body Methods
Open Boundary Problems (Isolated system, monolithic collapse)
Direct Sum Method
Tree method

Periodic Boundary Problems(statistics of objects)
Particle-Mesh (PM) Method
Particle-Particle Particle-Mesh (P3M) Method
Tree-PM Method
Adaptive Mesh Refinement (AMR)

N-body Methods
Direct sum method
Calculate force from each particle
O(N2)N-body Methods
Tree method
Distant particles are bundled upO(N log N)

N-body Methods
Particle-mesh (PM) method
Mass assignment
Poisson solver
Force interpolation
O(N log N)

N-body Methods
Particle-Mesh Method
Mass assignment (Cloud-in-Cell Scheme)

N-body Methods
Particle-Mesh Method
Mass assignment (Cloud-in-Cell Scheme)for (i=0; i<np; i++){ix=(int) pos[ip][0]; dx=pos[ip][0] – (double)ix;
iy=(int) pos[ip][1]; dy=pos[ip][1] – (double)iy;
iz=(int) pos[ip][2]; dz=pos[ip][2] – (double)iz;
rho[ix ][iy ][iz ] += mass[ip] * (1-dx) * (1-dy) * (1-dz);
rho[ix ][iy ][iz+1] += mass[ip] * (1-dx) * (1-dy) * ( dz);rho[ix ][iy+1][iz ] += mass[ip] * (1-dx) * ( dy) * (1-dz);
rho[ix ][iy+1][iz+1] += mass[ip] * (1-dx) * ( dy) * ( dz);rho[ix+1][iy ][iz ] += mass[ip] * ( dx) * (1-dy) * (1-dz);
rho[ix+1][iy ][iz+1] += mass[ip] * ( dx) * (1-dy) * ( dz);rho[ix+1][iy+1][iz ] += mass[ip] * ( dx) * ( dy) * (1-dz);
rho[ix+1][iy+1][iz+1] += mass[ip] * ( dx) * ( dy) * ( dz);}

N-body Methods
Particle-Mesh Method
Poisson Solver (Convolution method)• Δφ=4πGρ
• k2φk=4πGρk• ρー(FFT)→ρkー(conv.)→φkー(FFT)→φ
fft3d_real_forward (n, rho);for (kx=0; kx<n; kx++){for (ky=0; ky<n; ky++){for (kz=0; kz<n; kz++){phi[kx][ky][kz]=4.0 * PI * G * rho[kx][ky][kz]/ (kx*kx+ky*ky+kz*kz);}}}
fft3d_real_backward(n, phi);N-body Methods
Particle-Mesh Method
Force interpolation (Cloud-in-Cell Scheme)for (i=0; i<np; i++){ix=(int) pos[ip][0]; dx=pos[ip][0] – (double)ix;
iy=(int) pos[ip][1]; dy=pos[ip][1] – (double)iy;
iz=(int) pos[ip][2]; dz=pos[ip][2] – (double)iz;
acc[ip] =grv[ix ][iy ][iz ] * (1-dx) * (1-dy) * (1-dz)+
grv[ix ][iy ][iz+1] * (1-dx) * (1-dy) * ( dz)+
grv[ix ][iy+1][iz ] * (1-dx) * ( dy) * (1-dz)+
grv[ix ][iy+1][iz+1] * (1-dx) * ( dy) * ( dz)+
grv[ix+1][iy ][iz ] * ( dx) * (1-dy) * (1-dz)+
grv[ix+1][iy ][iz+1] * ( dx) * (1-dy) * ( dz)+
grv[ix+1][iy+1][iz ] * ( dx) * ( dy) * (1-dz)+
grv[ix+1][iy+1][iz+1] * ( dx) * ( dy) * ( dz);}

N-body Methods
Particle-particle particlemesh (P3M) method
Distant particles• PMNearby particles• Dierct sum

N-body Methods
Adaptive P3M Method

Couchman (1991)

N-body Methods
TPM (Xu 1995; Bode+ 2000;Bode&Ostriker2003)

P3MのPPをTreeにする

TreePM (Bagla 2002; Dubinski+ 2004;
Yoshikawa&Fukushige 2005; Springel 2005)

P3Mとは異なる重力の分解をする(Bagla 2002)

N-body Methods
Adaptive Mesh Refinement (AMR)
Divide cells only
where higher
resolution required

N-body Methods
AMRAMRPMYahagi & Yoshii (2001)

AMR の歴史
Berger & Collela (1989)
数値流体法で分解能が必要な個所に、より細かい構造格子を再帰的に配置
BC(89)では斜めに格子を配置している。
• この手法は、その後余り受け入れられていないようである。

Moving Mesh との違い
• AMR:
決まった分解能を最小の計算機資源で実現• Moving Mesh: 決まった計算機資源で最大の分解能を実現

AMRの宇宙論業界での歴史
Villumsen (1989)
Anninos, Norman, & Clarke (1994)
Kravtsov, Klypin, & Khokhlov (1997; ART)
Norman & Bryan (1999; Enzo)
Knebe, Green, & Binney (2001)
Yahagi & Yoshii (2001)
Teyssier (2002; RAMSES)

Villumsen (1989)

Anninos, Norman, & Clarke (1994)

Norman & Bryan (1999)Bryan, Abel, & Norman (astro-ph/0112089)

Kravtsov, Klypin, & Khokhlov (1997)

Yahagi & Yoshii (2001)

Adaptive Mesh Refinement

Adaptive Mesh Refinement
Data structurechildneighborchildcell cell
cell octet
parentneighbor(Khokhlov 1997)

Adaptive Mesh Refinement
Data structurechildneighborchildcell cell
cell octet
parentneighbor(Khokhlov 1997)

Adaptive Mesh Refinement
Data structureAdaptive Mesh Refinement
Data structureVectorization
Particle loops
Integration of the equation of motion

Mesh loops
Poisson equation

Particle-Mesh loops
Mass assignment
force interpolation
particle sieve

Vectorization
Loop order exchange• Tree (Hernquist 1990)• Hash (Yahagi 2002)

Vectorization
Loop split
Vectorized mass assignmentAMRコードの速度

AMRコードの速度高速化の動機
地球シミュレータ
並列化効率が50%になる台数まで利用可能
VPP5000/32PEで約50%
VPP5000とSX-6ではプロセッサあたりのFlops値はほぼ同じ
プロセッサあたりのメモリは地球シミュレータの方が小さい
つまり、天文台でできる計算より大きい計算が出来ない

階層格子のデータ分割
Morton ordering
Warren & Salmon (1995), Fryxell et al. (2000)

Warren & Salmon (1995)

階層格子のデータ分割
Morton ordering階層格子のデータ分割

階層格子のデータ分割

階層格子のデータ分割

階層格子のデータ分割

階層格子のデータ分割

粒子のデータ分割(旧コード)

粒子のデータ分割(旧コード)

粒子のデータ分割(旧コード)
Force
Massinterpolation
assignment

process 0

process 1

process 2

粒子のデータ分割(新コード)

process 0

process 1

process 2

粒子データの分割(新コード)

粒子データの分割(新コード)

粒子データの分割(新コード)

粒子データの分割(新コード)

この先危険私が遠くに言ってしまったら、この世に戻ってくるよう、遠慮なく呼び戻してください。

粒子データの分割
粒子の速度を更新
粒子の位置を更新
階層格子を更新
粒子の速度を更新

階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ
階層格子の消去
階層格子の追加
粒子リンクリストの局所化
粒子データの転送

階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ

階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ
階層格子の消去
階層格子の追加
粒子リンクリストの局所化
粒子データの転送

階層格子の更新と粒子データの転送
階層格子の消去階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ
階層格子の消去
階層格子の追加
粒子リンクリストの局所化
粒子データの転送

階層格子の更新と粒子データの転送
階層格子の追加

階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ
階層格子の消去
階層格子の追加
粒子リンクリストの局所化
粒子データの転送

階層格子の更新と粒子データの転送
粒子リンクリストの局所化
階層格子の追加が行われた後に、階層格子はMorton順序に従って整列させられる。
その結果、階層格子とその配下にある粒子のリンクリストが異なるプロセッサ上にくることがある。
そのため、階層格子と同じプロセッサ上へ粒子のリンクリストを転送する必要が生じる。

階層格子の更新と粒子データの転送
現在の階層格子に合わせて粒子をふるいわけ
階層格子の消去
階層格子の追加
粒子リンクリストの局所化
粒子データの転送

階層格子の更新と粒子データの転送
粒子データの転送
今まで転送されているのは粒子のリンクリストとインデックスだけである。
そこで、仕上げに、粒子の位置及び速度情報を送りあう。

AMR N体コードの改造#define LEV_NPTCL 23#define MAX_NPTCL (1<<LEV_NPTCL)#define MASK_PTCL_PROC (MAX_NPTCL-1)…
for (ip=0; ip<np; ip++){proc[ip] = ptcl[ip] >> LEV_NPTCL;
local[ip] = ptcl[ip] & MASK_PTCL_PROC;}

32ビットの壁

VPP5000(32PE) 5123体

プロセッサID
25=3210243体

27=128ES(128PE)

プロセッサ内の粒子ID
2*222=8M
2*223=16M

ES(128PE)10243体はぎりぎり入るが、粒子数分布の平均からのずれが1を超えると破綻• そこで、二つの4バイト整数に分割する

変更箇所
粒子のデータ構造(配列から延びるリンクリスト)の変更
粒子の移動
セルは変更せずに、粒子のリンクリストのつなぎ替えをする
オープンセルにフラグを立てる
フラグの無いセルに入っている粒子を上の階層に渡していく
セルを上のレベルから消したり足したりしていく、と同時に粒子を下の階層へ渡していくセルは変更せずに、粒子のリンクリストのつなぎ替えをする
下のレベルから上のレベルへ向けて以下の操作を行う
セル内に留まる粒子はそのままにして、セルの壁を越える粒子をリンクリストから抜き出す
抜き出した粒子を含むセルが異なるPEにある場合は、その粒子を転送する。
抜き出した粒子を含むセルがバッファセルの場合、粒子を上の階層へ転送する。

セルは変更せずに、粒子のリンクリストのつなぎ替えをする
下のレベルから上のレベルへ向けて以下の操作を行う
セル内に留まる粒子はそのままにして、セルの壁を越える粒子をリンクリストから抜き出す
抜き出した粒子を含むセルが異なるPEにある場合は、その粒子を転送する。
抜き出した粒子を含むセルがバッファセルの場合、粒子を上の階層へ転送する。

セルは変更せずに、粒子のリンクリストのつなぎ替えをする
上のレベルから下のレベルへ向けて以下の操作を行う
粒子を含むセルに子供がいない、又はいてもバッファセルである場合は、粒子をリンクリストに残す粒子を含むセルの子供がオープンセルなら下の階層へ粒子を転送

変更箇所
粒子のデータ構造(配列から延びるリンクリスト)の変更
粒子の移動
セルは変更せずに、粒子のリンクリストのつなぎ替えをする
オープンセルにフラグを立てる
フラグの無いセルに入っている粒子を上の階層に渡していく
セルを上のレベルから消したり足したりしていく、と同時に粒子を下の階層へ渡していくフラグの無いセルに入っている粒子を上の階層に渡していく
下の階層から上の階層へ以下の操作を行う
オープンセルだったが次のステップで消える、又はバッファセルになるセルオクテット配下のセル
中の粒子のリンクリストを結合する
結合した粒子のリンクリストを上の階層へ転送する
オープンセルのままでいつづける粒子のデータを詰める

上の階層から下の階層への操作は次にまとめて行う

変更箇所
粒子のデータ構造(配列から延びるリンクリスト)の変更
粒子の移動
セルは変更せずに、粒子のリンクリストのつなぎ替えをする
オープンセルにフラグを立てる
フラグの無いセルに入っている粒子を上の階層に渡していく
セルを上のレベルから消したり足したりしていく、と同時に粒子を下の階層へ渡していく

セルを上のレベルから消したり足したりしていく、と同時に粒子を下の階層へ渡していく
上の階層から下の階層へ以下の操作を行う
子セルオクテットがオープンセルオクテットである
セル中の粒子をリンクリストから取り出す
それら粒子が新しい子セルオクテット中のどのセルに入るかを調べる
粒子を新しいセルに振り分ける

まとめ
32ビット整数を使う場合、粒子のインデックスにローカルのインデックスとPEのインデックス両方入れるのは限界である。
ローカルのインデックスとPEのインデックスを分けることにした。
これで、いよいよES向けチューニングを始められる。

HPC用語の基礎知識
高速化率(speed-up)
台数を増やして、何処まで速くなったか
sn=T1/Tn

並列化効率
pn=sn/n

並列化率
経過時間のうち並列実行された個所の割合
αHPC用語の基礎知識
Amdahlの法則(並列版)

並列化できない個所

並列化できる個所

T1(1-α)T1αT1Tn
T1(1-α)T1α/n

並列化率の求め方
T1が計算できない場合は良くある異なるプロセッサ数で二回計測することにより、並列化率(α)を求めることが出来る。

並列化効率の求め方

ある並列化効率を満たす台数

要求されている並列化率

白色雑音問題での計測
N=2563

計時解析 IΛCDM simulations
Ω=0.3, λ=0.7, h=0.7, n=0.9, Ωb=0.048

使用機種
VPP5000@ADAC/NAOJ

粒子数
N=1283, 2563, 5123

CPU数
Nproc=2, 4, 8, 16, 32

計時解析 IIΛCDMでの計測

Parallel fraction (α)

10.9

0.8
1e+06

z=zstart
z=01e+07
Number of particles

1e+08ΛCDMでの計測
1

Parallel efficiency (p32)

0.8

0.6
z=zstart
z=00.40.2

01e+06

1e+07
Number of particles

1e+08

各ルーチン毎の計時解析
AMR construction
及び Assignment
and Interpolation では粒子データを取り扱うが、これら二つのケースでも理想的な状況に近いスピードアップが実現できている。

ΛCDMシミュレーションの経過時間を以下の三つの関数型でフィッティングした。

N=5123

Nproc

T[sec]

T[sec]

N=1283T[sec]

N=2563

Nproc

Nproc

33
N=256
N=256

Parallel Efficiency

前項のフィッティングからT1を見積もり、並列化効率pn=（T1/(nTn))を見積もった。
また、T64, T128 も前項のフィッティングから見積り、p64, p128を算出した。その結果から、n<=128まで、どのような見積りをしても、 50%以上の並列化効率を維持できるものと予測される。
N=1283

Nproc

Parallel Efficiency

Parallel Efficiency

N=5123

Nproc

Nproc

56 10粒子分割(a) L=Lbase=95

Number of cells

5 1054 10

Nparticle=1
Nparticle=253 10

Nparticle=352 105

1 100(b) L=12

Nparticle=1
Nparticle=25Nparticle=3

Number of cells

4 1053 1052 105

1 100(c) L=Lmax=154

6 104

Number of cells

5 10

Nparticle=1
Nparticle=244 10

Nparticle=3
Nparticle>=44

3 104

2 104

1 10001 0

2 0

Index of processors

3 0

上段は一番粗い階層の格子数、下段は一番細かい階層の格子数、中段はその中間の階層の格子数を表している。
赤、緑、青、水色は、粒子を一つ含む格子数、二つ含む格子数、三つ含む格子数、四つ以上含む格子数をそれぞれ表している。
粒子分割は深い階層の方が浅い階層より均等に分割されている。

6

Number of particles

1.5 10L=9L=10L=11L=12L=13L=14L=156

1.0 1055.0 10001 02 03 0

Index of processors
67×106

Number of particles

6×10

65×10

64×10

63×10

62×10

6

1×10001020

Index of processors

30

粒子分割
上段は各階層毎の各CPUの担当する粒子数、下段は各CPUの全粒子数を表している。
各格子に入る粒子数に制限のない一番細かい格子が均等に分割できないと思われるが、実際には、逆に最も良く均等分割されている。
均等分割から、最大で約50%の不均等が生じている。

誤差解析 I共同研究者: 牧野淳一郎
Hernquist, Hut, & Makino (1993; HHM93)
定常解では、全体の力学的エネルギーだけでなく、個々の粒子の力学的エネルギーが保存する
N体シミュレーションで Plummer 解を計算し、個々の粒子の力学的エネルギーの変化を調べた

誤差解析 II
G=1, M=1, E=-1/4
ε=1/32, δt=2
t=0 から t=8 までの4区間について平均を取る

スケーリングAMRN体コードでは、最も細かい格子間隔を1とする単位系を採用。
以下のスケーリングをした。
G=1,E=-1/4,,.
T ∝M5/2 |E|-3/2.

HHM93誤差解析 III
HHM93の図の上に、今回の結果を重ね合わせた(上図)。およそ近い値の誤差になっている。
をHHM93の図の上にプロットしたものである(下図)。ＡＭＲの結果はツリー法のθ=1.0の結果に近い結果となった。
NHHM93

Potential energy

ポテンシャル
大局的にはGRAPE-6の結果と、AMRの結果に差は無い。N=65536.

Potential energyrr

中心付近を拡大すると、ＧＲＡＰＥ-6よりAMRの結果の方が解析解に近い。
GRAPE-6:1.91367e-05
AMR:1.51096e-05

Hernquist Model

NEMOのmkhernquistを使って初期データを生成し、外側の粒子(M(r)/M>0.99)を取り除き、スケーリングを行う。(a=16, M=1)
AMR
LBASE=8, LDYNR=15

N=2048φr

N=65536

φr

N=65536

φrN=2048arN=65536arN=65536

ar

σφ2N

σa2N

まとめ
Hernquist モデルの力とポテンシャルの誤差を求めた。
力、ポテンシャルの両方で、中心付近ではAMRコードの方が解析解に近いふるまいをする。
ポテンシャルについては、粒子が少ない場合はGRAPEの方が誤差が小さい。
力については、どの粒子数でもAMRコードの方が誤差が小さい。
何故、力とポテンシャルで結果が変わるのか?

この先危険(再)地球シミュレータ固有の話です
飛ばしますか、そうですか並列化チューニング
メモリによる制限のため、10243体計算には64ノード(512AP)が必要
しかし、並列化でつまずき相談員に相談する
そうしたら、ESでは通信中に計算を実行することが出来ていないことが判明。I系非同期通信を全てSendrecvに置き換えることに。

標準の集団通信を使いましょう
MPIに含まれる集団通信で置き換えられる個所は置き換える。
for (jproc=0; jproc<N_PROC; jproc++)
MPI_Irecv(nrip+jproc,1,MPI_INT,jproc
TAG_POST_MOD_SIEVE_HM_NI,
MPI_COMM_WORLD, rnireq+jproc);for (jproc=0; jproc<N_PROC; jproc++)
MPI_Issend(tnsip+jproc,1,MPI_INT,jproc,
TAG_POST_MOD_SIEVE_HM_NI,
MPI_COMM_WORLD, snireq+jproc);
MPI_Waitall (N_PROC, rnireq, MPI_STATUSES_IGNORE);
MPI_Waitall (N_PROC, snireq, MPI_STATUSES_IGNORE);を以下のように置き換える
MPI_Alltoall (tnsip, 1, MPI_INT, nrip,1,MPI_INT,MPI_COMM_WORLD);

I系ではなくSendrecvを使いましょう
集団通信で置き換えられない、MPI_Issend, MPI_Irecv, MPI_Wait の組み合わせ
は MPI_Sendrecv で置き換えるfor (jproc=0; jproc<nsiproc; jproc++){kproc = lsiproc[jproc];
MPI_Issend (slp+hsip[kproc], nsip[kproc], MPI_INT, kproc,
TAG_POST_MOD_SIEVE_HM_I, MPI_COMM_WORLD, ireq->sreq+jproc);}for (jproc=0; jproc<nriproc; jproc++){kproc = lriproc[jproc];
MPI_Irecv (rlp0+hrip[kproc], nrip[kproc], MPI_INT, kproc,
TAG_POST_MOD_SIEVE_HM_I, MPI_COMM_WORLD, ireq->rreq+jproc);}を以下のように置き換える
nsrproc = 0;
for (jproc=0; jproc<N_PROC; jproc++){kproc = lkproc2[jproc];
if (nsip[kproc] > 0 || nrip[kproc] > 0) lsrproc[nsrproc++] = kproc;}
for (jproc=0; jproc<nsrproc; jproc++){kproc = lsrproc[jproc];
MPI_Sendrecv (slp+hsip[kproc], nsip[kproc], MPI_INT, kproc,
TAG_POST_MOD_SIEVE_HM_I, rlp0+hrip[kproc], nrip[kproc],
MPI_INT, kproc, TAG_POST_MOD_SIEVE_HM_I, MPI_COMM_WORLD, MPI_STATUS_IGNORE);}

一見なんてことないように見えるが、実は、lkproc2 を決めるのは自明ではない。

Sendrecvペアの作り方
前と通信するのか、後ろと通信するのか、それが問題だ
Np=4

Sendrecvペアの作り方
lkproc2の構築個所 (!(N_PROC & 1))
lkproc2[0] = ikproc;
for (jproc=1; jproc<N_PROC/2; jproc++){kproc = 1;
while (!(jproc & kproc)) kproc <<= 1;
if ( (ikproc / kproc) & 1){
lkproc2[jproc*2-1] = lkproc[N_PROC-jproc];
lkproc2[jproc*2] = lkproc[jproc];} else {lkproc2[jproc*2-1] = lkproc[jproc];
lkproc2[jproc*2] = lkproc[N_PROC-jproc];}}
lkproc2[N_PROC-1] = lkproc[N_PROC/2];

Np=8

jproc=1

jproc=２

jproc=3

jproc=4

もっとエレガントなペアの作り方
似鳥君の提案for (jproc=0; jproc<N_PROC; jproc++)
lkproc2[jproc] = ikproc ^ jproc;
A ^ J = B なら B ^ J = A なのでペアになる。
ペア間の距離は一定ではない• PE間の通信量がMorton順序の距離に応じて減少する

ならば、通信待ちが生ずることになるが、実際に試していないので何ともいえない。

時間計測
日時

8/04:12

9/04:14

9/05:15

9/05:15

9/05:15

粒子数

2563

5123

5123

5123

5123

ステップ

20

2020

20

80

32AP(sec)

65.189274.364271.311300.8481116.764AP(sec)

58.355

183.577205.549

193.414679.557

α0.8946

0.98430.96790.9877

0.9914

MAX-AP 10.49

64.62

32.11

81.97

117.4

mpi_ssend_init mpi_send_init mpi_sendrecv

地球シミュレータ
前回までのあらすじ
通信を隠蔽するMPIの登場• 以前のMPIはI系コマンドを使っても、計算によって通信を隠蔽することが出来なかった• その為、Sendrecv を使うことが推奨されたいた

8ノード64CPUを使い並列化率を計測• 20ステップの経過時間から2ステップの経過時間を引くと、並列化率0.99733(375.9CPU)が達成• 暫定利用申請へ地球シミュレータ
並列化率的には376CPUにしか達していな
かったが、512CPUを使う場合、計算規模が大きくなるので、512CPUを利用するだけの並列化率に達する見込みは大きい、と申請したところ、無事通過。
暫定利用64ノード(512CPU)をした結果、128ノード(1024CPU)でも並列化効率50%を達することが判明

地球シミュレータ
暫定128ノードで計算したところ、今度は1024ノードでも並列化効率が50%を維持することが判明。
暫定1024ノードで計算をしようとしたところ、計算が止まる。原因を究明しようにもジョブはなかなか回らない。そうこうしている間に暫定利用期間が終了。

地球シミュレータ
128ノード本申請が許可された段階で、10243体計算を64ノード使って開始。
現在は z=3 辺りまで計算が進んでいるが、計算が止まっている。原因は調査中。といっても、5月から中断中なので、この研究会から帰ったら、この問題に取り掛かる。

展望
無衝突系
開放境界条件• Tree– GRAPEで加速させることもできる

周期境界条件• Tree-PM– メモリ使用量が少ない– SPH• AMR– 高速(少なくとも私のは)– Euler法

展望
SuprimeCAM
100Mpc立法の計算の視野と同程度

Hyper-SuprimeCAM
11月の研究会では、嶋作さんが、zがなんぼかで、一回の撮像で300Mpcの領域が撮像できるようになると言っていた
質量分解能を上げるのではなく、領域を大きくする必要あり

Springel+(2005)