多次元ボルツマン輻射流体コードによる 超新星計算 長倉 洋樹(京大基礎物理学研究所) 共同研究者 岩上わかな(早稲田、基研)、住吉光介(沼津高専)、山田章一(早稲田) 古澤峻(国立天文台)、松古栄夫(KEK)、今倉暁(筑波) HPCI分野5全体シンポジウム@ 紀尾井フォーラム 2015/3/11-‐12 (分野5 課題3) 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 Diffusion 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, iffusion 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 Diffusion 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-‐Diffusion (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軸対称超新星計算を京で実行予定 ポスト京へ:3D+数値相対論+曲がった時空上でのニュートリノ輸送 多次元一般相対論的ボルツマンニュートリノ輻射流体計算 現状報告 ベースコードの概要 チューニング状況 京での小規模テストラン (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計算(1Dチェック) 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太陽質量の2モデル 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) 注:以下、一種類ニュートリノに対しての演算数の見積もり(実際はこの3倍) 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" 2048ノード: 実行効率 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. 今年度までの成果として ベースコード開発の完了 (1Dテスト計算で先行研究とconsistent) チューニングを進めた (2048 nodeで実行性能11%、Strong Scale 80%) 小規模テストランを京で実行 (prompt convecUonを確認) 3. 今後は4月から本格計算に向けてコードの整備 (解析コードの構築、小規模ランの詳細解析、データ管理 etc..) 以下補足 多次元ボルツマン流体計算の困難 次元が多い(空間3次元+運動量空間3次元+時間1次元) => 計算コストが大 (解像度チェックを行うにも大変) 新しい数値計算アルゴリズムの開発が必要 これまでの1次元球対称計算とは全く違った手法が必要 Lagrangian Code (流体comoving系をベースに解く) Eulerian Code (実験室系をベースに解く) ニュートリノ輸送(ボルマン)ソルバー 重力ソルバー 流体ソルバー 滝脇超新星計算とBoltzmann-‐Hydro計算との比較 滝脇 – 固武 − 諏訪 グループの超新星計算 来年度行う Boltzmann-‐Hydro 計算 3次元流体 + 1次元重力(Newton) 2次元流体 + 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 ニュートリノのフラックスが、 背景流体の非一様性に引きずられる エントロピー分布(3D滝脇計算) 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
© Copyright 2025 ExpyDoc