斎藤貴之
東京工業大学地球生命研究所

1

Contents• 研究背景• ASURA: 並列N体/SPHコード• Density Independent SPH• まとめ2

3© 内海洋輔(広島大学)

NASA/JPL-Caltech/ESO/R. Hurt

4

Credit: S. Beckwith & the HUDF Working Group (STScI), HST, ESA, NASA

我々の宇宙我々の宇宙我々の宇宙NASA / WMAP Science Team

銀河形成に関連するスケール
宇宙の果て

銀河中心
巨大BH

バルジ

銀河円盤

~106M◎~1010M◎~10-4pc~1kpc

銀河-銀河相互作用~1013-15M◎

大規模構造~1011M◎

銀河の力学的スケール~1012M◎~10kpc~100kpc~Mpc10Mpc~100Mpc-Gpc

星形成のスケール
分子雲のサイズ～10pc
分子雲の質量～106M◎
星形成領域のサイズ～1pc
恒星質量～0.1-102M◎

超新星爆発のスケール
超新星爆発のサイズ～10pc
はき出される金属量～1M◎

銀河形成シミュレーションモデル• 暗黒物質– 位相空間からサンプリングした粒子に置き換えて重力相互作用を計算• バリオン(通常の物質)– 圧縮性流体として扱う– ガスの放射冷却加熱

サブグリッドモデル• 通常事前計算したテーブルを利用– 低温高密度ガスから星形成• 質量の一部を重力相互作用だけする恒星集団粒子へ– 超新星爆発• 大質量星がエネルギーと金属を放出
10Simulation:Takayuki Saitoh, Visualization: Takaaki Takeda, Sorahiko Nukatani

Contents• 研究背景• ASURA: 並列N体/SPHコード• Density Independent SPH• まとめ12並列 N体/SPH コード

ASURA• C (C99) + MPI• 領域分割:Orthogonal Recursive Bisection• 重力:Parallel Tree+GRAPE– Hardware accelerators : GRAPE series– Software accelerator : Phantom-GRAPE(Tanikawa+2012)• SIMD命令で書かれた重力計算ライブラリ– Symmetrized Plummer Potential (Saitoh&Makino 2012)• 流体:Density Independent SPH(Saitoh&Makino 2013)• 時間積分:Leap-frog+Individual time steps+Time-step limiter (Saitoh&Makino 2009)+FAST (Saitoh&Makino 2010)• サブグリッドモデルライブラリ:– CELib(Saitoh 2017), ASRCH, ASRFLX並列 N体/SPH コード

ASURA• C (C99) + MPI• 領域分割:Orthogonal Recursive Bisection• 重力:Parallel Tree+GRAPE– Hardware accelerators : GRAPE series– Software accelerator : Phantom-GRAPE(Tanikawa+2012)• Assembler tuned software library!!

– Symmetrized Plummer Potential (Saitoh&Makino 2012)• 流体:Density Independent SPH(Saitoh&Makino 2013)• 時間積分:Leap-frog+Individual time steps+Time-step limiter (Saitoh&Makino 2009)+FAST (Saitoh&Makino 2010)• サブグリッドモデルライブラリ:– CELib(Saitoh 2017), ASRCH, ASRFLX

独立時間刻み法• 粒子ごとに異なる時間刻みを持たせて時間積分する
方法(McMillan 1986)– dtは粒子のローカルな物理状態から決め、explicit に積分• E.g., クーラン条件(~λ/cs)、～sqrt(e/a)、～(e/v)– dt=2-nT (nは正の整数)に切りつめる (Makino 1991)• 幅広い時間スケールを扱う銀河形成では標準的な手法• 作用反作用は満たさない：通常はその影響は小さい
粒子0
粒子1
粒子2
時間

対称性の破れが問題になる例Saitoh & Makino 2010
• 銀河形成sim.だと、explicit に決めた dt の間に周囲で SN ショックが発生したとき– TISM=10107K(=TSN):~1000

周りのガスが反応するまでに~1000step 進む• 相互作用相手に短い時間刻みがあれば explicit に決めた自分のdtも短くするLimiter導入:dt=min(dt,f*dtnb)、fは任意定数

密度点源爆発問題

半径Saitoh & Makino 2010

並列 N体/SPH コード

ASURA• C (C99) + MPI• 領域分割:Orthogonal Recursive Bisection• 重力:Parallel Tree+GRAPE– Hardware accelerators : GRAPE series– Software accelerator : Phantom-GRAPE(Tanikawa+2012)• Assembler tuned software library!!

– Symmetrized Plummer Potential (Saitoh&Makino 2012)• 流体:Density Independent SPH(Saitoh&Makino 2013)• 時間積分:Leap-frog+Individual time steps+Time-step limiter (Saitoh&Makino 2009)+FAST (Saitoh&Makino 2010)• サブグリッドモデルライブラリ:– CELib(Saitoh 2017), ASRCH, ASRFLX

“銀河形成シミュレーション”• 暗黒物質とバリオン– 重力、流体力学、放射冷却、星形成、超新星爆発• シミュレーション方法– Tree SPH 法＋独立時間刻み法• シミュレーション中で最も厳しいのは小粒子の積分の時のツリー構築(毎ステップフルスクラッチすると)– 何とかする方法として、1.Dynamic update (McMillan &Aarseth 1993)、2.ツリー構築をサボる(Wetzstein+2009,Nelson +2009)

Wadsley + 2004

もっとも短い時間刻み幅はどこから来る？？

23

FAST: A Fully Asynchronous
Split Time-Integrator fora Self-Gravitating Fluid

Saitoh & Makino 2010

• 一つの粒子の重力相互作用と流体相互作用に異なる時間刻み幅を与えて別々に積分する

シンプレクティック積分法
ハミルトニアンHに対する正準変換は次のように書ける
q,pをfと表して、ポアッソン括弧を用いて書くと次のようなオペレータを定義：
形式解：

ハミルトニアンHをq,pの項に分離する
形式解はBCH公式を用いて、“FAST構築

Saitoh & Makino 2010

自己重力流体のハミルトニアンは以下のように書ける
断熱を仮定しHを分離する
二次精度の形式解は次のようになる
Hhydroを二つに分ける

Hhydroの積分にも二次精度積分を用いる
Δtgrav = l Δthydro l≧1とすると

Saitoh & Makino 2010

Mergersimulation test

並列 N体/SPH コード

ASURA• C (C99) + MPI• 領域分割:Orthogonal Recursive Bisection• 重力:Parallel Tree+GRAPE– Hardware accelerators : GRAPE series– Software accelerator : Phantom-GRAPE(Tanikawa+2012)• Assembler tuned software library!!

– Symmetrized Plummer Potential (Saitoh&Makino 2012)• 流体:Density Independent SPH(Saitoh&Makino 2013)• 時間積分:Leap-frog+Individual time steps+Time-step limiter (Saitoh&Makino 2009)+FAST (Saitoh&Makino 2010)• サブグリッドモデルライブラリ:– CELib(Saitoh 2017), ASRCH, ASRFLXChemical Evolution libraryfor Galaxy Formation: CELib

C言語(C99)による実装、C++ からも利用可能Saitoh 201730

CELib distribution site

https://bitbucket.org/tsaitoh/celib
MIT license

31

Contents• 研究背景• ASURA: 並列N体/SPHコード• Density Independent SPH• まとめ32

並列 N体/SPH コードASURA• C (C99) + MPI• 領域分割:Orthogonal Recursive Bisection• 重力:Parallel Tree+GRAPE– Hardware accelerators : GRAPE series– Software accelerator : Phantom-GRAPE(Tanikawa+2012)• Assembler tuned software library!!

– Symmetrized Plummer Potential (Saitoh&Makino 2012)• 流体:Density Independent SPH(Saitoh&Makino 2013)• 時間積分:Leap-frog+Individual time steps+Time-step limiter (Saitoh&Makino 2009)+FAST (Saitoh&Makino 2010)• サブグリッドモデルライブラリ:– CELib(Saitoh 2017), ASRCH, ASRFLXSmoothed Particle
Hydrodynamics (SPH) とは• SPH 法は、Lucy (1977)、Gingold &Monaghan (1977) により開発された圧縮性流体の解放– ラグランジュ法の一種– 流体物理量は粒子からの寄与の畳み込みで与えられる

Muller+’03 SIGGRAPH~cm scale

Saitoh et al.
~1022cm scale従来の SPH 法の問題点• Agertz+2007– SPH 法と Euler 法の比較
– SPH は接触不連続面の扱いが苦手
不安定性成長の抑制
fundamental difference(see also Okamoto+2003)– 原因はSPHの定式化に密度の微分可能性を用いているから
GridSPH

Agertz+2007SPHGrid

メッシュの問題：!=ガリレイ不変
数値拡散

エントロピー↗
メッシュコアが溶けている

Springel 2010

Tasker+2008

Formulation ofStandard SPH• 位置 r の物理量 f を次のように評価する• 体積要素 ΔV を用いて離散化• 体積要素 ΔV = m/ρ と f = ρ から

• 密度がスムーズであるという仮定が要請されている– 接触不連続面で破綻

37

接触不連続面における圧力の評価

Saitoh & Makino 2013

密度8:1の接触不連続面、圧力=1(=一定)過大評価

過小評価• 接触不連続面の密度gap圧力エラー斥力となる• 密度から他の物理量を求めるのに無理がある– 離散化からやり直す必要がある

38

Density Independent
Saitoh & Makino 2013
SPH

• 次の体積要素を定式化に用いる：• 物理量 f は次のようになる• f に q を入れる：

39
理想気体では q は圧力に比例する (P=(γ-1)q)運動方程式とエネルギー方程式• 運動方程式

Saitoh & Makino 2013
See also Hopkins 2012

• エネルギー方程式40接触不連続面における圧力の評価 (DISPH)

Saitoh & Makino 2013

Pressure

圧力を基本量にしているので、接触不連続面でなめらかに分布Hydrostatic equilibrium
tests

Saitoh & Makino 2013

ρ=1
ρ=4

Initial condition

DISPH

SSPH

(q ∝ P)
See Hopkins 2013, Hosono TRS et al. 2013, Yamamoto TRS et al. 2015,
43
Read et al. 2010, Ritchie & Thomas 2001; See also Inutsuka 2002, Price 2008

Saitoh & Makino 201345

Kelvin-Helmholtz
instability tests• シアー起源の流体
不安定性の試験• 初期条件：密度比
1:2、圧力2.5、相対速度1• 境界面にy方向速度の摂動を入れる
λ=1/6、A=0.025-0.5

Saitoh & Makino 2013

密度1、圧力 2.5

密度2、圧力 2.5
N=131072-0.5

密度1、圧力 2.5
N=655220.5

Kelvin-Helmholtz
instability tests

Saitoh & Makino 2013

Rayleigh-Taylor instability
tests
Saitoh & Makino 2013
• 静水圧平衡にある流体中で、重力により引き起こされる不安定性• 初期条件：密度比1:2• 境界面にy方向速度の摂動を入れる

重力Rayleigh-Taylor instability
tests
Saitoh & Makino 2013

Blob tests• Agertz+2007– SPH/Eular codes の比較– 大体tkh=1ぐらいから表面に発生した不安定性で壊れる

Blob tests

54

Saitoh & Makino 2013

Two phase fluid mixing

Saitoh & Makino 2013

Point like Explosion tests

Saitoh & Makino 2013

57

他の方法

58

Moving mesh: AREPO

Springel 2010

• 移動メッシュでガリレイ不変性の回復
AREPO オリジナルの実装は実質一次精度
時間積分、勾配評価の更新 Pakmor et al. 2016

Hopkins 2015

Meshfree method• 形状関数を用いて局所的な物理量分布を再構築– SPH の問題の一つであった
P.U. の回復
空間精度の向上• GIZMO:

http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html

60

スキームの影響:
サンタバーバラクラスタ Saitoh & Makino 2016
• 重力+断熱ガスによるクラスタ形
成 : Mesh/SPH
でエントロピープ
ロファイルに系統的な違い• DISPH/Moving
mesh/Meshfreeはほぼ一致• SSPH は非物理的表面張力の影響、メッシュは数値粘性
61Contents• 研究背景• ASURA: 並列N体/SPHコード• Density Independent SPH• まとめ62

まとめ• 並列 N 体/SPH コード”ASURA”を開発• 銀河形成研究のためのアルゴリズムを開発して ASURA に実装– Time-step limiter (Saitoh & Makino 2009)– FAST (Saitoh & Makino 2010)– Symmetrized Plummer Softening Tree(Saitoh & Makino 2012)– DISPH (Saitoh & Mkaino 2013)• 銀河形成及び広く天体物理現象に応用している

63