多次元ボルツマン輻射流体コードによる超新星計算

多次元ボルツマン輻射流体コードによる超新星計算長倉　洋樹（京大基礎物理学研究所）
共同研究者
岩上わかな（早稲田、基研）、住吉光介（沼津高専）、山田章一（早稲田）古澤峻（国立天文台）、松古栄夫（KEK）、今倉暁（筑波）
HPCI分野５全体シンポジウム@ 紀尾井フォーラム 2015/3/11-‐12
（分野５　課題３）
H23
H24
H25
H26
H27
課題： 3+1+1
3+3+1ボルツマン流体コードを開発し、 3
1
H25
超新星爆発に対する輻射流体計算を行う (r,
384*128*256*20
京では2次元軸対称の科学的計算をH27年度に実行予定現状： 70m
3
ベースコードの開発を完了 FX10及び京を用いて、コードのチューニングを進めた r
小規模テストランの実行 H26
3+3+1
3+3+1
H27
課題代表者：柴田さんのスライド(去年)を拝借
重力崩壊型超新星爆発
多次元的な流体不安定性＋ニュートリノによるエネルギー輸送が爆発の鍵（右図は滝脇さんスライドより）
Catoon From Iwakami D thesis
ν-加熱は中間領域で起こる：!
超新星爆発におけるボルツマン計算の必要性
ニュートリノ輻射の詳細が必要：�
ex.
!
D
ν-加熱率�
Janka A&A (1996)
C
(
MeV % Lν Eν2
Q ≈ 110
X*
'
s⋅ N & R72 < µ > i )
i
ν
A
freestreaming
B
€
R
€
heating region
領域C
< µ >=< cosθν >= 0 ~ 1
Eddington factor:!
1
f =< cos 2 θν >= ~ 1
3
€
shock
ν"
領域D
Flux factor:!
10~20 %
outside
外層（低密度）
shockwave
Average energy, flux: Eν, Lν!
€
f (Eν ,θν )
!
領域B
ボルツマン計算が必要な領域 θν
（近似的な取り扱いが難しい）
ν-sphere
領域A
diffusion
ν"
ν"
center
中心（高密度）
4�
Cartoon fro Janka and Sumiyoshi
近似的ニュートリノ輸送スキーム
The Astrophysical Journal, 747:73 (12pp), 2012 March 1
Ray-‐by-‐Ray Approach (MPA, Oak Ridge, Kotake-‐Takiwaki-‐Suwa)
1
2
Neutrino-‐AdvecUon is essenUally considered under spherical symmetry.
Isotropic Diﬀusion Source ApproximaUon (IDSA) 1
2
(Basel, Kotake-‐Takiwaki-‐Suwa)
constant.
kb denotes
the Boltzmann
Neutrinos are decomposed into trapped and streaming parts. constant.
Figure 9. Illustration of the “ray-by-ray” transport approximation.
and the solid lines represent two in
§2. tMoment
formalism
of Thorne
Two reduced equaUons are coupled by each source erm, which is represents the neutrinosphere
SchemaUc picture for approximately described under dFirst,
iﬀusion treatment. Ray-‐by-‐Ray we review the Thorne’s moment formalism.
In theapproach first step, h
(Lentz e
t a
l. 2012)med
(See e.g., Berninger et al. 2013)
an unprojected moment of massless particles associated with a moving
Moment method “rays” in the RbR approximation. The dashed lines are tange
neutrinosphere and indicate the regions that contribute to the neu
at points 1 and 2. The “blob” on the neutrinosphere below point 1 is a
where the temperature is2)
higher than the rest of the neutrinosphere. F
the RbR method will compute the neutrino field as if the entire neut
has the properties of the hot spot, overestimating the neutrino flux an
For point 2, the RbR misses the contribution of the hot spot by assu
the neutrinosphere properties are only those of the cooler region dire
′ the neutrino flux and heating.
it, underestimating
′α1 ′α2
′αk
′
f (p′α , xβ )δ(ν − ν )
=
p p · · · p dVp ,
ν ′k−2
reduce computational
costs and
code devel
Shibata et simplify
al. 2011
(MPA, Kyoto, Caltech, Basel (Kuroda))
CHIMERA, Vertex, and Zeus+IDSA
break the ′µ
no
′
where f is the distribution function of the(lateral,
relevant
radiation,
ν through
= −utheµ p“ray
or angular)
spatial coupling
(RbR) approximation,
and Vulcan/2D breaks the c
Neutrino angular direcUon is iquency
ntegrated. he so-‐called “closure relaUon” is groups
imposed of theTradiation
in the
rest-frame
of theenergy
medium
the restbetween
and(i.e,
neutrinoin
species.
µ
µ
In the
RbR approximation,
is co
in the higher moment.
the fiducial observer) with u being medium’s
four
velocity, the
p neutrino
the transport
four-mo
as a number of independent, spherically symmetric pr
of the radiation, and dVp the invariant integration
element
on for
thethelight
referred to as “rays,”
which allows
reuse ofc
1D neutrino transport codes. (See2)
Figure 9 for a sc
here, is positive integer, 1, 2, · · · . As pointed
out
byRbRThorne,
choic
illustration
of the
approximation.)the
RbR methods
good
parallel
scaling
for
large
numbers
of
independe
fiducial observer is crucial when deriving arays,
good
truncated formalism from
(Oak Ridge, Princeton, Caltech)
which can be evolved without communicatio
computing
the neutrino
Typically, in RbR
ment formalism. In the following, the fluid,
coupled
withtransport.
the radiation,
i
the neutrinos in opaqueα
regions
arekadvected laterally
1 α2 ···α
2),
9),
10)
always
de
as the
Namely, the
frequency,
ν, in
M(ν) to the pressure.
Neutrino Transports are treated as medium.
the Energy-‐Dependent Diﬀusion quaUon.
fluid E
motions
and contribute
The indep
of the rays artificially sharpens the lateral variation
frequency measured in the rest-frame of theneutrino
fluid luminosity
throughout
thisabove
paper.
Th
and heating
the proto-NS
results terms
in some regions
of the
hot mantle being
ov
is crucially helpful when computing the source
of the
radiation
equat
α1 α2 ···αk β
M(ν)
(x )
!
MulU-‐Group Flux-‐Limited-‐Diﬀusion (MGFLD) Boltzmann-‐Hydro Code 開発
1次元球対称一般相対論的ボルツマンニュートリノ輻射流体計算 (Yamada 1997, 1999 and Sumiyoshi et al. 2005)
多
次
元
化
数値アルゴリズムの大幅な変更（流体、重力、ニュートリノ輸送の解法、これら全てを変更）数値コストの拡大 6次元ボルツマンコード開発
Sumiyoshi and Yamada, 2012 Sumiyoshi et al., 2014
6次元ボルツマン + 流体コード開発
Nagakura et al. 2014
来年度に2D軸対称超新星計算を京で実行予定ポスト京へ：３D＋数値相対論＋曲がった時空上でのニュートリノ輸送
多次元一般相対論的ボルツマンニュートリノ輻射流体計算
現状報告
ベースコードの概要チューニング状況京での小規模テストラン (preliminary)
where εfr (≡0pt ≡ −uµ pµ ) d
n-shell condition pµ pµ = −m2ν , in which mν is a neutrino
ical are
schemes
(e.g.,
tube
problems)
√
fluid
restframe.
Here
uµwer
is th
, only three of the four components
independent
andshock
r
⎟
⎜ −Lorentz
gGtransformati
The
Nagakura
al. (2011).
s why only spatial components appear
in the et
second
term
⎜ of√neutrino
θ⎟
relation
energies
i
e left-hand side; τ stands for the affine
parameter
of
the
− gG
⎟
⎜
Although our Boltzmann
solver
is
fully
SR,
th
frame
⎜as √ φ ⎟fr. lb
ino trajectory. The left-hand side of Equation (1) expresses Wi =
− √gG
ics solver
is right-hand
Newtonian. As a⎜matter
of fact,
⎟ε =itεcaγ
odesic motion in the phase space,
while the
ニュートリノ輸送（ボルマン）ソルバー
⎝
t ⎠ for it
code
(Nagakura
2008)
except
−
gG
symbolically denotes the so-called
collision
terms, i.e.,&
theYamada
where v, √
γ (≡ (1 − v 2 )−1/
s that give the rate of changes inwhich
f due tois
neutrino–matter
Newtonian and corresponding
based−on the
MICCG
Lorentz
factor
gΓ
actions.
流体ソルバー
重力ソルバー
that indicates the flight
direc
gakura
et
al.
2011).
The
implementation
of
an
E
lb
n the spherical coordinates in flat spacetime, which are the
frame. The
factor
D ≡ γ (1o
solver
is
currently
underway,
the
perspective
dinates we employ for the laboratory
frame
in
our
Eulerian
Note that Wi correspondsthetoDoppler
the interactions
b
shift of neutrino
oach, Equation (1) is cast into the
following conservation
mentioned
in Section
8.expressions
we can obtain the
relation
bet
ボルツマン方程式 (
Sn m
ethod)
and
matter
(the
explicit
will
be
pre
流体、レプトン保存 :
√
frames as
The basic
equations
Newtonian
hyd
2
(shock cof
apturing scheme) facto
# and
g(=
r
sin
θ)
denotes
the
volume
"lb
are written in the !
followin
1 − µ2ν spherical
cos φν ∂ coordinates
∂f µν ∂ 2
（各コンポーネントは補足参照）
δf
j
+ 2
(r f ) +
(sinOther
θf )
,
v
, an=
coordinates.
variables,
ρ,
p,
e,
Y
e
∂t #r ∂r
r sin θ
∂θ
δt col
j
pressure,
internal
∂t Q + energy
∂j U =density,
Wh + Welect
1 − µ2ν sin φν ∂f density,
1 ∂
i,
2
+
[(1 − µν )f ]
+
The Lorentz transformatio
and Newtonian gravitational
poten
r sin θ
∂φ velocity,
r ∂µν
the flight directions in the flu
#
"
!
lb
2
where
each
is given重力 asis (solved
Newtonian
self-gravity
with the
1 − µν cos θ ∂ The
δf term
as
MICCG)
−
(sin φν f ) =
,
(2)
%
$⎞
r
sin θ ∂φν
δt col
⎛
√
frgρ
fr
lb
lb
ε
n
=
ε
+
n
∆ψ
=
4πρ.
衝突項移流項 √
e r, θ, and φ denote
the spatial variables.
For the three
⎜ √gρvr ⎟
(Collision Term)we do
(AdvecUon erm) four-momentum,
pendent components
of Tneutrino
⎜ fr
⎟
gρv
⎜
⎟syste
se the spacial components but adopt
the energy
and two
vect
Here n
θ the unit
In our
central
scheme,
the
above
√denotes
計算コードの概要
8
1e+09
7
0
6
-1e+09
5
-2e+09
entropy
2e+09
速度
-3e+09
4
2
-5e+09
-6e+09
エントロピー
HN
Sumi
3
-4e+09
1
HN
Sumi
0
0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
mass coordinate (Mso)
2
0
0.2
0.4
0.6 0.8 1 1.2 1.4
mass coordinate (Mso)
1.6
1.8
2
0.5
14
HN
Sumi
12
0.45
10
温度
8
Electron fracUon
0.4
Ye
T (MeV)
velocity (cm/s)
Boltzmann-‐Hydro計算（１Dチェック）
6
0.35
4
0.3
2
0
0
0.2 0.4 0.6 0.8
1
1.2 1.4 1.6 1.8
mass coordinate (Mso)
2
0.25
HN
Sumi
0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
mass coordinate (Mso)
2
現状報告
ベースコードの概要チューニング状況京での小規模テストラン (preliminary)
モデルと計算サイズ
親星モデル 11.2、15太陽質量の２モデル The Astrophysical Journal, 771:27 (17pp), 2013 July 1
Yamamoto
3
3
1D
計算領域中心から4000km 2D
1
計算時間 Post bounce 後 400〜500 ms
51
[10
0.3
0.1
（400 ms 以内に爆発する必要あり）
See Yamamoto et al. 2013 →
E
E
[10
51
erg]
erg]
1
Eexp
Enuc
Eν
0.3
0.1
Eexp
Enuc
Eν
0.6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0.6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
tpb
tpb
[s]
[s]
Figure 21. Individual contributions of neutrino heating and nuclear reactions to the explosion energy for all the 1D and 2D models.
(A color version of this figure is available in the online journal.)
グリッド数及び解像度
qualitatively different from the 2D SASI we have studi
this paper (Iwakami et al. 2008). It should also be recalle
Real Space: 384 (Nr) × 128 (Nth) (2010)
that shock revival is even e
ΔR min
〜 100 Nordhaus
m, Δet cal.os θ claimed
= 2/Nth
in 3D than in 2D, although controversies are still contin
Momentu Space：　20 (Ne) × 10 (Na) × 6 (Nb)
(Hanke et al. 2012). If the critical luminosity is much
0.3
[M ]
0.25
1D
2D
0.2
0.15
in 3D than in 2D, the yield of 56 Ni may be reduced furth
3D. The complex flow patterns also have some influence o
nickel yield. We are particularly concerned with how the all
region in the shock-relaunch time that is opened in 2D is m
fied in 3D. The relative importance of nuclear reactions fo
explosion energy compared with neutrino heating is the hi
in 1D. We are certainly interested in how these proportions
or may not change in 3D.
One of the greatest uncertainties in the present study
effect of the inner boundary condition that is imposed by
The artificial treatment employed in this paper results in a
mass injection from the inner boundary of about 7 × 10−
in the 1D fiducial model. This injection contributes 2%
to the explosion energies and the 56 Ni masses. Although
may be a slight underestimate (Arcones et al. 2007), w
MNi56
ノード数： 2048 (64×32) : 6 (ir) × 4 (ith) 1 node あたり
1
024 (
64×18) :
6 (ir) × 8 (ith) Real Spaceを並列化
0.1
0.05
0
0.1
0.2
0.3
0.4
0.5
tpb
ステップ数
0.6
0.7
0.8
0.9
[s]
200〜300万ステップ　（Δt 〜 10 -‐7
s）
Figure 22. Comparison of 56 Ni masses in the ejecta between the 1D and 2D
models.
(A color version of this figure is available in the online journal.)
energy and nickel mass in the appropriate range. This hap-
Boltzmannパートのチューニング
Ns （spaUal mesh size）, Ne（neutrino energy mesh size） Na ( neutrino angular mesh size ), Nb ( neutrino angular mesh size ) Nite (number of matrix iteraUon), Np (preconditoner factor) 注：以下、一種類ニュートリノに対しての演算数の見積もり（実際はこの３倍）
Step1: 反応レートの計算 Step2: 衝突項の行列要素計算
Step3: OpUcal Depthの計算
Step4: 移流項の行列要素計算
Step5: 行列前処理計算
Step6: 行列計算
（反応率標準セット） Ns × Ne × (Na × Nb)^2 （電子散乱） Ns × (Ne × Na × Nb)^2
(Ns)^α × Ne × Na × Nb （with 通信） α=2 (1D), α=3/2 (2D), α=4/3 (3D)
Ns × Ne × Na × Nb （with通信）（反応率標準セットケース） Ns × Ne × (Na × Nb)^2 × (Nite + Np) (with通信) （電子散乱込みケース） Ns × (Ne × Na × Nb)^2 × (Nite + Np) (with通信)
チューニング状況
-‐7
ストロングスケール (main loop)
Δ t 〜 10 sで時間発展
Average MFLOPS/PEAK (%) 200 Ume stepのElapsed Time (s)
14"
1600"
12"
1400"
1200"
10"
1000"
２０４８ノード：実行効率 11 % を達成 (Strong Scale 〜80％) 8"
6"
800"
1"
600"
ProducUve runではI/Oを含め
て8〜9％前後の性能
4"
2"
0"
0"
500"
1000"
1500"
ノード数
2000"
2500"
400"
Current Code
100％スケール
200"
0"
0"
500"
1000"
1500"
2000"
ノード数
単ノードチューニング：If 分岐の削除とメモリアクセスの改善 → SIMD化を徹底！
通信チューニング：通信の呼び出し回数を削減（Packing）。毎ステップ行う必要のない領域を洗いだし、通信量を削減。
マトリックスチューニング : Matrix をブロック毎に分割。全体のElapsed Timeを削減前処理を単精度化。
2500"
現状報告
ベースコードの概要チューニング状況京での小規模テストラン (preliminary)
Post–bounce 約 50 ms 計算 on 京 (preliminary) 岩上さん（早稲田、京大）による計算
Entropy
Lateral Velocity
まとめ
1. 開発した多次元ボルツマン流体コードを用いて、来年度に京でプロダクティブランを実行予定 2. 今年度までの成果としてベースコード開発の完了（１Dテスト計算で先行研究とconsistent）チューニングを進めた（2048 nodeで実行性能11％、Strong Scale 80％）小規模テストランを京で実行 (prompt convecUonを確認) 3. 今後は４月から本格計算に向けてコードの整備（解析コードの構築、小規模ランの詳細解析、データ管理 etc..）
以下補足
多次元ボルツマン流体計算の困難
次元が多い（空間３次元＋運動量空間３次元＋時間１次元） => 計算コストが大　（解像度チェックを行うにも大変）
新しい数値計算アルゴリズムの開発が必要これまでの１次元球対称計算とは全く違った手法が必要
Lagrangian Code (流体comoving系をベースに解く)
Eulerian Code (実験室系をベースに解く)
ニュートリノ輸送（ボルマン）ソルバー
重力ソルバー
流体ソルバー
滝脇超新星計算とBoltzmann-‐Hydro計算との比較
滝脇 – 固武 − 諏訪グループの超新星計算来年度行う Boltzmann-‐Hydro 計算
３次元流体＋１次元重力（Newton）
2次元流体＋２次元重力（Newton）
ニュートリノ輸送（IDSA）
ニュートリノ輸送 5次元Boltzmann with full v/c order
ニュートリノ反応標準セット EOS Lasmer & Swesty ニュートリノ反応標準セット +
電子散乱原子核の電子捕獲反応 (GSI Table) EOS Furusawa EOS （NSE with light Nuclei）
– 39 –
3次元超新星プロファイル中でのBoltzmann計算 – 35 –
Sumiyoshi et al. 2014
t = 100 ms
t = 150 ms
t = 200 ms
ニュートリノのフラックスが、背景流体の非一様性に引きずられる
エントロピー分布（３D滝脇計算）
1.— Profiles of the 3D supernova core adopted in the current study. Entropy iso-surfaces
Fig.
shown for the snapshot at 100, 150 and 200 ms after the bounce in the 3D supernova
ution by Takiwaki et al. (2012) from top, middle and bottom panel, respectively.
ニュートリノ密度のコンターマップ
5.— Iso-surface of density of electron-type anti-neutrinos (ν̄e ) for the 3D supernova core
at 150 ms after the bounce. Arrows represent the flux vector of neutrinos.
number of iterations
odesic motion in the phase space, while the right-hand
symbolically denotes the so-called collision terms, i.e., the
where v, γ (≡ (1 − v 2 )−1/
s that give the rate of changes in f due to neutrino–matter
corresponding Lorentz factor
actions.
Boltzmannパートの概要
that indicates the flight direc
lb
n the spherical coordinates in flat spacetime, which are the
frame.
The
factor
D
≡ γ (1
PTEP
2012, 01A301
dinates we employ for the laboratory frame in our
Eulerian
the Doppler shift of neutrino
oach, Equation (1) is cast into the following conservation
we can obtain the relation bet
ボルツマン方程式
MxN
:
M
frames as
#
!x 3 "lb
x2
x1
1 − µ2ν cos φν ∂
∂f µν ∂ 2
δf
+ 2
(r f ) +
(sin θf )
2=
∂t #r ∂r
r sin θ
∂θ
δt col
1 − µ2ν sin φν ∂f 1 ∂
1
+
[(1 − µ2ν )f ]
+
The
Lorentz
transformatio
r sin θ
∂φ r ∂µν
the flight directions in the flu1
#
! "lb
2
1 − µν cos θ ∂
δf
as
−
(sin φν f ) =
,
(2)
%
$
r
sin θ ∂φν
δt col
fr fr
lb
lb
ε
n
=
ε
+
n
衝突項移流項 e r, θ, and φ denote
the spatial variables.
For the three
(Collision Term)we do
(AdvecUon erm) four-momentum,
pendent components
of Tneutrino
se the spacial components but adopt the energy and two
Here nfr denotes the unit vect
s, θν and φν (see Figure 3). µν is defined as µν ≡ cos θν .
of a neutrino in the fluid res
Fig. 21. Left: The pattern of the sparse matrix appear
cretization of the Boltzmann equations. N and M denote
(Semi) Implicit Methodを用いるので、大規模行列計算が必要
4
grids (Nθν , Nφν , Nε ), respectively. For the studies on c
size of the diagonal black matrices (gray) is Nθν Nφν . Rig
step for different pre-conditioners, i.e., the point Jacob
t
h
j
i
where each term is given as
⎛
⎞
√
√ gρ
⎜ √gρvr
⎟
⎜
⎟
⎜ √ gρvθ ⎟
Q=⎜
⎟,
⎜√ gρvφ ⎟
⎝ g(e + 1 ρv 2 )⎠
√ 2
gρYe
(13)
√
⎞
gρv j
⎜ √g(ρvr v j + pδrj ) ⎟
⎟
⎜ √
⎜ g(ρvθ v j + pδ j ) ⎟
θ
⎟
Uj =⎜
⎜ √g(ρv v j + pδ j ) ⎟ ,
φ
⎜√
φ ⎟
⎝ g(e + p + 1 ρv 2 )v j ⎠
2
√
gρYe v j
(14)
⎛
流体及びレプトン数保存式
⎛
with
ρ, Y
6
⎞
0
(
⎜√gρ −ψ,r + r(v θ )2 + r sin2 θ (v φ )2 + 2p ⎟
⎜
rρ ( ⎟
⎜√ '
⎟
p
cos
θ ⎟
2
φ
2
⎜
Wh = ⎜ gρ −ψ,θ r + sin θ cos θ (v ) + ρ sin θ ⎟ , (15)
√
⎜
⎟
− gρψ,φ
⎜
⎟
√
⎝
⎠
l
− gρv ψ,l
0
'
⎞
0
√
⎜ − gGr ⎟
√
⎟
⎜
⎜ − gGθ ⎟
Wi = ⎜ √ φ ⎟ .
⎜− √gG ⎟
⎝ − gGt ⎠
√
− gΓ
⎛
(16)
Note that Wi corresponds to the interactions between neutrinos
and matter (the explicit expressions will be presented in Step 5)
In
ent e
Sect
a rat
of th
putin
take
stead
large
grids
inter
In
the f
1.
2.
Th
hype
of ex
The
fact,
just
parti
the v
in th
Th
Figu
for e
A bu
disti
grid
large
not t

Download Report