富岳における無衝突自己重力系のVlasovシミュレーション
神戸大学

惑星科学研究センター CPSセミナー
2021年5月11日 16：00~

筑波大学 計算科学研究センター
吉川 耕司

共同研究者
田中 賢 (京都大学 基礎物理学研究所)

斎藤 俊 (Missouri S&T)

吉田 直紀 (東京大学 Kavli IPMU)

簑島 敬 (JAMSTEC)

1

Contents
 自己重力系とVlasov 方程式
 VlasovシミュレーションとN体シミュレーション
 Vlasovシミュレーションの数値スキーム
 宇宙大規模構造形成におけるニュートリノのVlasovシミュレーション
 富岳における実装・最適化
 更なる応用2

無衝突自己重力系とVlasov方程式
 無衝突自己重力系

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

Vlasov 方程式

 ある Hamiltonian に従う多数の粒子の分布関数の時間発展: 分布関数 （位相空間での粒子の数密度）
 無衝突 Boltzmann 方程式（Vlasov方程式）

位相空間での移流速度
 Hamilton 方程式

 自己重力系のHamiltonian

4

N体シミュレーション
 6次元位相空間の物質分布を超粒子でサンプリング サンプリングした点の運動をVlasov方程式の特性方程式に沿って数値的に解く
 物理量（密度場・速度場）などにショットノイズが内在する。
 分布関数のテイルの部分は相対的にサンプリングされにくい
 テイル部分が重要な役割を果たす運動論的現象を表現するのが不得手
 無衝突減衰、二流不安定
 磁気プラズマの粒子シミュレーション(PIC シミュレーション）についても基本的に同じ。
5

Vlasov シミュレーション 分布関数を有限体積法で離散化することで、Vlasov方程式を直接数値的に解いて分布関数の時間発展を得る

 ショットノイズの影響を避けられる 分布関数のテイルを再現できるので、運動論的効果も計算できる
 6次元位相空間を離散化
かなり膨大なメモリ容量が必要。計算コストも大きい。
空間分解能はあまり期待できない 偏微分方程式の数値解法

数値拡散の影響は避けられない
高次精度スキームの必要性

 N体シミュレーションに替わる独立なシミュレーション手法
長所・短所は互いに相補的6

無衝突自己重力系のVlasovシミュレーション
 Vlasov-Poisson 方程式系を直接時間積分

 6次元位相空間でのVlasovシミュレーションを世界で初めて実現
 かつては低次元の位相空間で行われていた。 Fujiwara (1981, 1983)
 核融合プラズマのシミュレーションでも5次元 (gyro-kinetic Vlasov simulation)
 筑波大 CCSのT2K-Tsukubaで世界で最初に実現

T2K-Tsukuba (2008-2014)

Yoshikawa, Yoshida, Umemura (2013)7

無衝突自己重力系のVlasovシミュレーション
 6次元位相空間を直交座標系にそってレギュラーメッシュに離散化 移流方程式を有限体積法で時間積分
 MPI並列のため、位置空間を領域分割
 運動量(速度)空間は領域分割せず
 メッシュ数は最大で Nx=643、Nv=643

 directional splitting method

運動量メッシュ 時間積分

6本の1次元移流方程式に分割: ℓ 方向へのΔt だけの時間発展8Merging of two King spheres
 等質量の King 球 を正面からぶつける

質量密度

位相空間密度

9

3D gravitational instability andcollisionless damping ホワイトノイズを持つ密度揺らぎの無衝突減衰と重力不安定
 初期条件 密度揺らぎのパワースペクトルはホワイトノイズ
 Jeans wavenumber

無衝突減衰
重力不安定

10

3D gravitational instability andcollisionless damping11高次精度移流スキーム
 可能な限り数値拡散を抑え、空間分解能を上げたい 一番手っ取り早い方法は、「メッシュ数を増やす」。
 メモリ容量が膨大なVlasovシミュレーションでは、そうもいかない。

 空間高次精度スキームを採用する
 メッシュ数・メモリ容量はそのままで実効的な解像度を向上させる
数値解の単調性・正値性を保証する必要性
1次元移流方程式の空間高次精度な数値解法は数値流体の教科書にいろいろ書いてある 空間精度を向上させると、時間積分の次数も上げないと数値解が不安定に
空間5次精度時間3次精度のTVD Runge-Kutta スキーム
数値流速を1時間ステップに3回計算する必要あり

 2013年の論文では空間3次精度のPFCスキーム
宇宙大規模構造形成のシミュレーションで精度が足りず。

12

単段の時間高次精度スキーム
 保存型のSemi-Lagrange法による数値流速を使うと単段の時間積分でOK

 ポスト「京」重点課題9に参加していたJAMSTECの簑島さんが発見 単調性・正値性を保証し空間高次精度な単段時間積分スキーム
SL-MPP5スキーム：空間5次精度
SL-MPP7スキーム：空間7次精度
Tanaka, Yoshikawa, Minoshima, Yoshida (2017)13

高精度移流スキーム
 1次元自己重力系
 初期条件空間3次精度

時間3次＋空間5次精度

時間単段＋空間5次精度

時間単段＋空間7次精度

 空間3次精度では数値拡散の影響が大きい。
 空間5次精度、7次精度を用いた結果では同じメッシュ数でも細かい構造まで分解出来ている。
 Semi-Lagrange法を用いた時間単段スキームは時間3次精度のスキームと同じ精度が得られている。（計算コストは1/3。）

14

宇宙大規模構造形成におけるニュートリノに力学的影響

Yoshikawa, Tanaka, Yoshida, Saito (2021)
15

宇宙におけるニュートリノ
 ニュートリノ質量に対する制限

宇宙マイクロ波背景放射

ニュートリノ振動

バリオン音響振動 初期宇宙で photon から脱結合した宇宙論的なニュートリノの数密度

 宇宙大規模構造に対する影響
 非相対論的になる赤方偏移 非相対論的になったニュートリノがCDMと重力相互作用
 宇宙大規模構造形成に対する力学的な影響

16ニュートリノ質量のプローブとしての大規模構造
 無衝突減衰 (collisionless damping / free streaming) 臨界波数よりも小スケールの密度揺らぎが減衰 ニュートリノの速度分散： 減衰の強さがニュートリノ質量に依存
Hannestad et al. (2006)

 宇宙大規模構造の観測からニュートリノ質量の制限が可能
 ニュートリノの質量階層も観測で決定できるかも ニュートリノ密度揺らぎの非線形効果を取り入れる
ために数値シミュレーションが必要
17

N-body + Vlasov Simulation
 N-body シミュレーションとVlasov-Poissonシミュレーションのハイブリッド
CDMは速度空間ではほとんど広がりを持たない
CDMについてはN体シミュレーションを採用
 CDM：Particle-Mesh法 / TreePM法によるN体シミュレーション

 ニュートリノ：Vlasov-Poisson シミュレーション(SL-MPP5 or SL-MPP7スキーム)
宇宙論的共動座標系におけるVlasov方程式：スケール因子：成分 i の質量割合

 Poisson 方程式18

Simulation Settings cosmological parameters (Planck 2015 result) neutrino mass

 3つの質量固有値は等しいと仮定。Mν < 0.2 eV では本当は良くない。

 simulation volume

 # of particles and mesh grids
N-body simulation for CDM
19

N-body + Vlasov Simulation of
Cosmological Relic Neutrinos
CDMの密度場

ニュートリノの密度場

 ニュートリノの運動を Vlasov シミュレーションで計算した世界で初めてのシミュレーション

Yoshikawa, Tanaka, Yoshida, Saito (2020)20

N体シミュレーションとの比較
Vlasov シミュレーション
Np=1923, Nv=643
NCDM=7683

必要なメモリ容量はほぼ同程度

N体 シミュレーションNν = 15363

NCDM=7683

 N体シミュレーションではニュートリノの細かい構造がショットノイズに埋もれてしまっている。(特に Void 領域)

21

速度場

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

速度分散

22

富岳における最適化

23

A64FX プロセッサ Core Memory Group (CMG) 12 compute cores + HBM2(8GiB)
 4CMGs per socket 512 bit-wide SIMD instruction set (SVE) 32 SIMD registers / core
 16 predicate registers / core
 2FMA units / core

 Peak Performance 0.75Tflops (DP) / 1.5Tflops (SP) per CMG

24

Scalable Vector Extension (SVE) A64FXプロセッサに実装されているSIMD拡張命令
 512bit 幅のSIMDレジスタが32本
 レジスタ上のデータに対して、半精度32個、単精度16個、倍精度8個を並列実行
 コンパイラがちゃんとコードからSIMD実行可能な部分を判断してSVE命令を吐いてくれると良いのだが、、、

 ACLE (ARM C-Language Extention)
 C/C++からSVE拡張命令を直接操作するAPIメモリからレジスタにデータをload

メモリ・レジスタ間のデータの load / store
レジスタ間のデータの並べ替え

レジスタ上でいろいろ計算

各種の整数・浮動小数点・論理演算
 x86 プロセッサのSIMD intrinsic と同じ要領

結果をメモリに store
25

SIMD命令を用いた移流計算の高速化
x軸方向への移流
 y 軸に沿ってデータをSIMDレジスタに load
 複数の x 軸方向への移流計算がそのままSIMD実行可能

26

SIMD命令を用いた移流計算の高速化
ｙ軸方向への移流
 SIMD実行するためには不連続なデータを各SIMDレジスタに納める必要がある。
 データの並べ替えを行ってからSIMDレジスタに格納するのはコストが大きい

27

SIMD命令を用いた移流計算の高速化
 y軸に沿ってデータをSIMDレジスタにload SIMD レジスタ上でデータを転置
16×16のデータの転置は64回のSIMD命令で実行可能。

 ｙ軸に沿った移流計算のSIMD実行ができるデータ配置になっている。

28

SIMD命令を用いた移流計算の高速化
CMG あたりの性能 [ Gflops ]方向w/o SIMD inst.

w/ SIMD inst.

w/ transposition

vy

7.14

233.3

N/Avxvz
xyz

4.84
7.44

5.51

6.88

6.50

176.7

N/A

17.9

224.2

154.1

N/A

150.0

149.2

N/A

N/A

単精度理論ピーク性能 : 1.5 Tflops/CMG

 uz方向以外の方向に沿った移流計算は、uz方向の loop を6重loopの内側から2番目の loopとすることで効率的な SIMD 実行が可能。
 uz方向の移流は、uz軸に沿った向きにデータを SIMD レジスタに load して転置することで、高性能を達成。
 理論ピーク性能の12～15%の実行性能を達成

29

Strong & Weak Scaling Efficiencies

 Weak scaling: 69~94%
 Strong scaling: 83-86%
 最大ノード数付近で、N体シミュレーション部分が足を引っ張ってしまった

30

富岳における世界最大規模のVlasovシミュレーション
CDM

NCDM=69123
neutrino

Np=7683, Nv=643
 73728ノードを用いて3.56hで終了

31

磁気プラズマへの応用

32

磁気プラズマのVlasovシミュレーション
 PIC (Particle-In-Cell) simulation
 無衝突プラズマの粒子シミュレーション

3D PIC simulation of collisionless shock 電子・イオンの粒子の運動方程式

 電子・イオンの運動による電磁場

Matsumoto et al. (2017)

磁気リコネクションのPICシミュレーション

 無衝突衝撃波における粒子加速 磁気リコネクション・磁気回転不安定性の素過程(pCANS web page より)

33

Vlasov-Maxwell シミュレーション
 PICシミュレーションの超粒子近似によるショットノイズ
粒子加速、磁気回転不安定性

 PICシミュレーションでは電磁場を計算するためのメッシュ幅はDebye長以下でなくてはならない。
MHDスケールまでの一貫したシミュレーションが困難

 Vlasov-Maxwell 方程式を直接解くVlasov-Maxwell シミュレーション

34

Gyro Motion 磁場によるジャイロ運動

速度空間内での剛体回転移流

rigid-body rotation of a 2D gaussian profile

 剛体回転移流を数値的に解くのは難しい。
 数値拡散が著しく、100回転もするとピークが半分以下になってしまう。
 100回転ぐらいはマトモに解けないと使い物にならない35

我々の移流スキームで剛体回転移流
初期条件

SL-MPP7スキーム

空間3次精度スキーム

100回転

gaussian profile の幅の時間変化

 rigid rotation of the 2D gaussian profile 空間3次精度スキームでは数値拡散が進行し、100回転でprofileの幅が50%増しに。
 SL-MPP5スキームとSL-MPP7スキームは、1~2%の増加にとどまる。
36

Vlasov-Maxwell シミュレーション
color : Jz

: Joule heating

dissipation of magnetic field: momentum transport

momentum transport to upstream region

37

まとめ
 世界で初めての6次元位相空間でのVlasovシミュレーション 空間高次精度移流スキームの開発
 計算コストを抑えた単段時間積分スキーム

 宇宙大規模構造形成のおけるニュートリノの力学的効果
 N体シミュレーションとのハイブリッドシミュレーション 従来のN体シミュレーションによる結果の検証

 富岳における最適化
 効率的なSIMD命令の利用による高速化

 磁気プラズマへの応用 (Vlasov-Maxwell シミュレーション) 他にもいろいろ、、
38