重力多体系のVlasovシミュレーションとその先筑波大学 計算科学研究センター
吉川耕司
2023年度 第2回 HPCIC計算科学フォーラム
2024年3月22日
共同研究者：

田中 賢 (京都大学 基礎物理学研究所)
吉田 直紀 (東京大学 Kavli IPMU)
斉藤 俊 (Missouri Sci. & Tech Univ.)自己紹介 筑波大学 計算科学研究センター
宇宙物理研究部門
 専門分野 宇宙大規模構造、銀河団、銀河間物質
 宇宙物理学におけるHPC(重力多体系、流体力学、輻射輸送)

 2013年ころからVlasovシミュレーションの研究を始める
 Yoshikawa, Yoshida, Umemura, ApJ, 762, 116 (2013)
 Tanaka, Yoshikawa, Minoshima, Yoshida, ApJ, 849, 76 (2017)
 Yoshikawa, Tanaka, Yoshida, Saito, ApJ, 904, 159 (2020)
 Yoshikawa, Tanaka, Yoshida, SC ’21: Proc. of the International Conference on High
Performance Computing, Networking, Storage and Analysis, Article No. 5, (2021)

本講演の内容
 宇宙物理学における重力多体系の数値シミュレーション 重力多体系のVlasovシミュレーション
 宇宙大規模構造におけるニュートリノの数値シミュレーション
 N体シミュレーションとの比較
 VlasovシミュレーションのSIMD命令を用いた高速化
 Vlasovシミュレーションの他分野への応用

重力多体系の数値シミュレーション
 無衝突自己重力系（重力多体系）

Ishiyama et al. (2013) 非常に多数の質点で構成されていて、個々の質点の運動が系の大域的な質量分布で決定される。
 重力的に束縛されたダークマター（ダークマターハロー）
銀河・銀河団・宇宙大規模構造 重力N体シミュレーション
 無衝突自己重力系の数値シミュレーションでは、これまでほぼ唯一の選択肢。
 過去数十年にわたって、さまざまな天体形成の研究に利用されてきた。
 計算手法の改良も長い歴史があり、極めて洗練された手法。
 今日のテーマは、「これとは全く別の手法に敢えて挑戦しました」というお話。

phase space distribution

N-body Simulation
 In collisionless N-body simulations, the matter distribution inphase space is statistically sampled by a set of super-particles.
 The motion of each super-particle is integrated according to the EoM.

N-body simulation

 adaptive resolution : good spatial resolution in high-density regions
 sophisticated schemes for huge number of particles and parallelization
 shot noise contamination associated with statistical sampling not good at handing kinematic effect (e.g. free streaming) due to poor
sampling at the tail of velocity distribution

phase space distribution

Vlasov-Poisson Simulation
 In Vlasov simulations, the phase space density is discretized
directly in a finite volume manner.
 The distribution function is time-integrated in according to
the collisionless Boltzmann equation (aka Vlasov equation).

Vlasov simulation
 capability to accurately handle the free streaming effect free from shot noise large amount of memory space and computational cost
World first Vlasov simulation in 6D phase space
Yoshikawa, Yoshida, Umemura (2013)

Vlasov-Poisson Simulation
 Head-on collision of two self-gravitating spheres
 Nx=1283, Nv=643

distribution function

density field宇宙大規模構造形成におけるニュートリノ
 宇宙大規模構造 宇宙の質量の大半を占める cold dark matter(CDM) の重力によって形成される構造
 ダークマター・ダークエネルギーの正体、宇宙初期の物理状態を調べる手掛かり ニュートリノ
 質量など基本的な性質が不明な素粒子
 宇宙初期において他の成分から脱結合したニュートリノ
 ニュートリノとダークマターの重力相互作用
 宇宙大規模構造の観測によるニュートリノ質量などの解明
 ニュートリノを取り入れた大規模構造形成シミュレーション これまでの研究はN体シミュレーション

Credit: 2dF Galaxy Redshift SurveyHybrid of Vlasov/N-body Simulation
 ダークマターとニュートリノの2成分系での宇宙大規模構造の数値シミュレーション
 ダークマター成分はN体シミュレーションで計算
速度分散が極めて小さいので速度空間での広がりはほぼ無い

 ニュートリノ成分はVlasovシミュレーションで計算

 ダークマターとニュートリノが感じる重力ポテンシャル
ダークマターとニュートリノの密度分布からPoisson 方程式を用いて計算
9Hybrid of N-body + Vlasov Simulation
CDM density

neutrino density

Yoshikawa, Tanaka, Yoshida, Saito (2020)

富岳の全系によるシミュレーション
ダークマター粒子数：69123

ダークマター Vlasov シミュレーションのメッシュ数：11523ｘ643
 ノード数：147456
 Vlasovシミュレーションの部分は strong scaling/ weak scaling ともに90％以上の並列化効率
 2021年の ACM Gordon-Bell Prize のファイナリストに選ばれた。
“A 400 Trillion-Grid Vlasov Simulation on Fugaku Supercomputer:
Large-Scale Distribution of Cosmic Relic Neutrinos in a Sixdimensional Phase Space”

ニュートリノ

N体シミュレーションとの比較
 neutrino density fieldN-bodyNν = Nc

Mν=0.4 eV

Mν=0.1 eVNν = 8Nc

Nν = 64Nc

VlasovNν = 180Nc

N体シミュレーションとの比較
Nν = Nc

Mν=0.4 eV

Mν=0.1 eV

N-bodyNν = 8Nc
Nν = 64Nc

Nν = 180Nc

Vlasov

N体シミュレーションとの比較
 neutrino velocity fieldN-bodyNν = Nc

Mν=0.4 eV

Mν=0.1 eVNν = 8Nc

Nν = 64Nc

VlasovNν = 180Nc

N体シミュレーションとの比較
 neutrino velocity dispersion
N-bodyNν = Nc

Mν=0.4 eV

Mν=0.1 eVNν = 8Nc

Nν = 64Nc

VlasovNν = 180Nc

Comparison of Numerical Resources
N-bodymemory

wall time

137 GiB
2.00 h420 GiB
2.28 h

Vlasov

2688 GiB
4.83 h7386 GiB
10.39 h1122 GiB
2.30 h Vlasov simulation is clearly superior to the N-body approach in the simulation of cosmological neutrinos.

Vlasovシミュレーションの高速化・最適化
 数値計算手法

1次元の移流方程式に方向分割

有限体積法で離散化 SIMD命令による高速化(single instruction multiple data) 複数のデータに対して同じ演算を並列に行う命令 プロセッサ内の専用レジスタ上で演算を実行
 レジスタの大きさは256 bit ～ 512 bit
単精度演算を8～16個倍精度演算を4～8個並列に実行

移流スキームのSIMD命令による高速化
 x軸に沿った移流 y軸に沿った移流SIMDレジスタ上での“in-place” なデータの転置18

移流スキームのSIMD命令による高速化#define DF(ix,iy,iz) df[((iz)+NMESH_Z*((iy)+NMESH_Y*(ix)))]for(int ix=0;ix<NMESH_X;ix++) {for(int iz=0;iz<NMESH_Z;iz+=vector_length) {svfloat32_t r0, r1, r2, r3, r4;
r0=svld1(pred, &DF(ix,iy,iz));
r1=svld1(pred, &DF(ix,iy+1,iz));
r2=svld1(pred, &DF(ix,iy+2,iz));
r3=svld1(pred, &DF(ix,iy+3,iz));

4x4 データの転置

svtranspose_4x4_f32(&r0, &r1, &r2, &r3);
sweep_1d_sve(&r0, &r1, &r2, &r3);
svtranspose_4x4_f32(&r0, &r1, &r2, &r3);
svst1( pred, &DF(ix, iy, iz), r0 );svst1( pred, &DF(ix, iy+1, iz), r1 );svst1( pred, &DF(ix, iy+2, iz), r2 );svst1( pred, &DF(ix, iy+3, iz), r3 );}}移流スキームの計算

SIMD命令による転置 (SVE命令)
 svzip1 命令と svzip2 命令の組み合わせで実現 sv{type} svzip1_f32(sv{type} op1, sv{type} op2)
Intel AVXなどの unpacklo 命令に相当 sv{type} svzip2_f32(sv{type} op1, sv{type} op2)

田中君のqiita記事より転載

 4x4データ(4つの4要素レジスタ) の転置の例
svzip100010203101112132021222330

3132

3300

2001

21

1030

1102

221232

svzip2

転置svzip100

102030

3101

1121

3103

2302

122232133303132333

svzip2

 16x16 データの転置は64命令で実行可能

富岳における最適化
CMGあたりのVlasovシミュレーションのパフォーマンス

direction

w/o SIMD inst.

w/ SIMD inst.

w/ transpositionvy

7.14 Gflops

233.3 Gflops

N/Avz

7.44 Gflops

17.9 Gflops

224.2 Gflops

x5.51 Gflops

150.0 Gflops

N/Ay

6.88 Gflops

154.1 Gflops

N/Az

6.50 Gflops

149.2 Gflops

N/Avx

4.84 Gflops

176.7 Gflops

N/A(cf. 単精度理論ピーク性能 : 1.5Tflops/CMG ) Vlasovシミュレーションの演算性能は理論性能の12-15%を達成
 (参考)我々の移流スキーム：B/F＝0.23
A64FXの単精度浮動小数点演算でのHBM2メモリに対するB/Fは0.18

21Further Application : Magnetized Plasma Particle-In-Cell (PIC) simulation

3D PIC simulation of collisionless shock particle-based method for collisionless plasma EoM for protons and electrons
Matsumoto et al. (2017)

 E and B obtained by solving Maxwell eq.
 Similar drawbacks to grav. N-body simulation

PIC simulation of magnetic reconnection

 shot-noise contamination
 mesh spacing in solving Maxwell eq. must be smaller than the Debye scale.
difficult to cover the MHD scales.
(pCANS web page)Vlasov-Maxwell Simulation

 Vlasov-Maxwell equation

 Mesh spacing not restricted to the Debye scale Vlasov-Maxwell simulation of solar wind andplanets with weak intrinsic B-field.
 no significant difference between the runs withΔH=λD and 10λD
 Vlasov simulation works well even for ΔH>λDUmeda, Fukazawa (2015)

Solving Gyration around B-field rigid body rotation in velocity space numerical solution of rigid-body rotationis sensitive to numerical diffusion.

 SL-MPP5 / SL-MPP7 schemes we found our advection scheme works excellently inrigid-body rotation

artificial heating of plasma
spatially 3rd order PFC scheme

SL-MPP5 scheme

velocity dispersion
PFC scheme

SL-MPP5/7 scheme# of rotation

Vlasov-Maxwell Simulation
 Vlasov-Maxwell equation

 Vlasov-Maxwell simulation of magnetic reconnection

 Vlasov-Maxwell simulation
 free from shot noise seen in the PIC simulation
 free from the constraint of the minimum scale to the Debye length consistent simulation from the kinetic scale to the MHD scale.

(simulation by T. Minoshima)

Toward Boltzmann Simulation
 collision term in the Boltzmann equation

 self-interacting dark matter (SIDM) as a remedy for the CDM crisis ingalaxy scales collisional plasma
 particle-based approach will miss the interaction where particles do
not reside in the velocity space calculation of the collision term can be significantly accelerated by GPUs
 we can now conduct the Boltzmann-Poisson simulation in 6D phase spaceスペクトル法による衝突項の計算
 フーリエ変換で波数空間へ変換した衝突項を、分布関数の増加に寄与する部分と、減少に寄与する部分で分解する
gain項

loss項

フーリエ変換 /離散化
gain項
loss項
畳み込み計算 畳み込み計算をFFTで高速化 まともに計算すると計算コスト∝ N6。FFTを用いることで計算コスト∝MN4 log N (M << N2)

スペクトル法による衝突項の計算

GPU(A100)

CPU

空間メッシュ数(mesh)

5121

衝突項の計算時間[msec/step]

593

1126

Effective wall time
[msec/step/mesh]1.151126

 複数の空間メッシュ内での衝突項の計算を並列化
 高速フーリエ変換の並列化
 GPUの利用で空間メッシュあたり数百倍の高速化
 3次元空間(6次元位相空間)でのボルツマン方程式の時間発展の数値シミュレーションが可能。

まとめ
 6次元位相空間上でのVlasov-Poissonシミュレーションを世界で初めて実現。
 宇宙大規模構造に付随するニュートリノのダイナミクスのシミュレーションに応用。
 磁気プラズマの数値シミュレーションにも応用が可能。
 SIMD命令を用いたVlasovシミュレーションの高速化。
 衝突項を取り入れたBoltzmann シミュレーションへの拡張。