サイエンティフィック・システム研究会(SS研) 2016年10月27日@ANAクラウンプラザホテル神戸

スーパーコンピューターを用いた太陽のシミュレーション
千葉大学

堀田英之
共同研究者：
Matthias Rempel [HAO/NCAR]
横山央明 [東京大学]

SS研

1

自己紹介
名前：堀田英之
所属：千葉大学
性別：男
年齢：31歳
身分：特任助教(テニュアトラック)
単位系：cgs guass
経歴
2005-2009年 東京大学
2009-2014年 東京大学大学院
横山央明准教授指導2014-2015年 HAO/NCAR
海外学振2015年現職

SS研2

目次・黒点数１１年周期の問題
・太陽ダイナモ計算の現状と問題・数値計算法・結果・これからの研究課題とコード開発の展望

SS研

3

太陽黒点

太陽黒点：
太陽表面に見える強磁場・低温度の領域
普通の領域：10 G・6000 K
黒点：3000 G・4000 K
黒点領域では放射層からの熱対流による熱輸送を抑制することで低温になっている。

SS研

4太陽黒点とその周期(NASA/JAXA/岡本丈典)(SDO/HMI)

地球から見える太陽半球における黒点の面積(％)

SS研

5

蝶形図

Latitude11年周期は規則正しいルールを持っているü 太陽黒点の出現緯度は周期の間に赤道に移動する。
ü 全球双極子磁場は11年に一度反転する。
SS研

6

ヘールニコルソンの法則
1996-2008

E-++-2008-

W

E+--+

W

90%以上の黒点ペアが従う極性のルール.
この法則は太陽内部の大規模磁場の存在を示唆している。
SS研

7

太陽ダイナモのざっくりとした考え方
1. 対流層の底で温め、表面で冷やすことで熱対流が発生
対流層
2. 乱流的な熱対流が、回転の影響を受け、非等方的になる
3. 非等方な乱流が角運動量を運び、差動回転や子午面還流を作る。

放射層

4. 差動回転、子午面還流、非等方な乱流のコンビネーションにより大規模な磁場をつくる。
熱対流をよく理解することが重要SS研

8

解かなければいけない乱流
エネルギー注入領域
運動エネルギー

慣性領域

⌫v

拡散領域

⇠ 1 cm2 s 1⇠ 104 cm s 1

L ⇠ 2 ⇥ 1010 cmvL
Re ⇠⇠ 1014⌫O(0.1 cm)

440万 km 20万 km
一番大きな見積もり

慣性領域をより多く分解したいが、太陽が大きくてエネルギー注入領域を分解するだけで精一杯。
SS研

9

スケール間ギャップ(Hinode/SOT)(Rempel, 2014, Numerical Calculation)

SS研

10

これまでのダイナモ計算
The Astrophysical Journal, 731:69 (19pp), 2011 April 10

disrupt the convection at the mid-latitudes. In co
dynamo case D5 the differential rotation is substan
in both radial and latitudinal angular velocity coevident in Figure 3(a), the convective cells fill iquite completely.
In our prior hydrodynamic simulations of cyounger suns we reported on localized nests o(Brown et al. 2008), with those most prominent arotation rates. Though there is some modulation w
Figure 1.in
Magneticwreaths achieved(a) case D3 (PaperI) andthis(b) case
the equatorial
rollin amplitudes
here,modu
D5 (this paper). Shown are time–latitude plots of mean(Brown+2010)
toroidal (longitudinal)
in either caseD5scalingor H5thanin ourmagnetic extreme
field ⟨Bφ ⟩ at mid-convection
zone, withvaluesindicated.
Case preD3 buildsrotating
persistent, time-independent
wreathsbutthewreaths
achievedinsimulations of stellar convectioncase(e.g.,
D5
undergoquasi-regularreversals of polarity (three shown here, with roughly
最近の3次元計算にでは、コヒーレントな磁場やその周期変化が再現でき
Brown
et al.
2008).
difference
appears to bea 1500 dayperiod). The
dynamic
range Thisof the color
bars is indicated.
backgroundstratificationfrom ther(A color versionof this figureis Racine+2011,
available in theand
onlinefeedbacks
journal.)
るようになっている(Ghizaru+2010, Brown+2010,
near the top of the domain. Here, we have attemp
Käpylä+2012, Masada+2013, Fan+2014,
Warnecke+2015,
regionof ofinfluence
the unresolved
SGS
the twothe
strongwreaths
oppositelyofdirected
polarities havewhichcarries fluxconnectivity
out the topof thesubstantial
cross-equatorial(Figure2(a)).domain
This (Karak+2015)。
cross-equatorial
flux
appears
toplayanimportantroleinthe2008). This thinner thermal boundary has largergpolarityareversals
thatareobserved.
steeper profile of the background entropy gradien
Unlike the wreaths of case D3, here magnetism fills the entire
thuszoneslightlyhigherRayleighnumbers
and主な生成機構はα効果(乱流誘導起電力)とΩ効果(差動回転によるひきのば
convection
includingthe polar
regions (Figure2(b)).
On radiIn our broader
studytheofwreaths
rapidlynearrotating
dynamtheir high-latitude(polar) edges,
the equatorし)である。
are connectedmagnetic
structuresof weaker
amplitude
found tothat
strongly
localized
activenests ofandconve
oppositepossible
polarity atpolar caps.
These polarin thedynamo
simulations
at structures
the mostare
rapid
relic wreaths from the previous cycle that propagate toward(Ω ! 10 Ω⊙ ).
SS研
11solar-lik
The convection establishes a prominent高解像度ダイナモ計算で弱る大規模磁場
一方で、(発表者が確認している限り)全てのグローバルダイナモ計算はとても低い解像度でおこなわれている。解像度をあげる、すなわち粘性・磁気拡散を下げると、小スケール乱流によって大局的磁場は弱くなる。 磁気拡散：η=2.6x1012 cm2 s-1
磁気拡散：η=1.2x10 12 cm2 s-1(Nelson+2014)

ローカルダイナモでは、引き延ばしが最も重要なプロセスであり、最も小さなスケールが、最も大きな線形成長率を持つ。
とても小さな粘性(O(1 cm2 s-1))・磁気拡散(〜O(104 cm2 s-1))を持つ実
際の太陽はどのように大規模磁場を作っているかは謎。
SS研

12理想的乱流シミュレーションからの示唆
高磁気レイノルズ数(Rm〜2500)での2次元運動学的ダイナモ(Tobias+2013, Cattaneo+2014, Nature)
exponentialな成長成分を規格化した
平均磁場
helical
重要な提言は、大ス
no shear
ケールダイナモ達成のためには、小スケールダイナモを抑non-helical制することが必要でwith shearあるということ
この論文では強いシアを加えて小スケールダイナモを殺している。

helicalwith shear
SS研13

小スケールダイナモの非線形効果
Energy spectra

小スケールダイナモ
✓◆
D B
B
=·r v
Dt ⇢
⇢

Kinetic
Magnetic

大スケールダイナモ

kinematic phase

Rempel, 2014@hBi
= r ⇥ (hvi ⇥ hBi) + r ⇥ (hv0 ⇥ B0 i)@t

大規模流れや乱流の非等方性がなくても、乱流場そのものによって小スケールの磁場をつくることができ、解像度が確保されれば小スケールで磁場エネルギーが運動エネルギーを上回り、小スケール流れ場を抑制
慣性領域を一桁以上分解することが必須の様子。
SS研

14

全球大規模計算の解像度-5/3-6

これまでの世界最高の解像度はMiesch+2008の計算(磁場なし)。
慣性領域のようなものはせいぜいエネルギー注入領域の半分くらいまでしか解けていない。
磁場があったとしても小スケールダイナモは不活発。

Spherical harmonic degree l
SS研

当然他の計算はもっと解けてない
15おおざっぱな動機
大きな問題は、解像度をあげればあげるほど大規模磁場は弱くなっていき現実からずれていくこと。
太陽ダイナモ計算の暗黙の了解は解像度をあげすぎないこと。
当然、これでは現実の太陽を理解できたことにならないので、なんとか打開策を探る。
現在のところ最高解像度の計算では乱流の慣性領域はほとんど解けておらず、小スケールのダイナモ効果が効率的ではない。
太陽全球を扱っても、小スケールのダイナモ効果が効率的になるような高解像度を達成することを目標とする。

SS研

16

Anelastic近似
太陽対流層内部では音速がとても速いので、陽解法で解くとCFL条件がきつい∆x∆t <cs + vc + ca

108cs: 音速 vc: 熱対流速度, ca: アルフベン速度cm s-1107

そこで採用されるのが、アネラスティッ
ク近似。そこでは連続の式を以下のように変形する。

0 = r · (⇢0 v) $ 0 = r · v

音波の速度は無限大になり音波の伝搬は解かれない。CFLは熱対流速度だけで決まる。

SS研

Velocities

音速

106
105
104

対流速度

103
0.70 0.75 0.80 0.85 0.90 0.95 1.00
r/RO•
17Anelastic近似の問題
音速無限大の仮定：音速が無限大で情報を運ぶので毎時間ステップすべてのCPUの情報をすべてのCPUが必要とする。通信負荷大高解像度化の限界(ASH限定)：球面調和関数を用いたスペクトル法を用いた手法が一般的であるが、フーリエ変換(O(NlogN))と違いルジャンドル陪関数を用いた展開には高速なものがないので、負荷はO(N2)となる。
また、MHDでよく使われるslope-limiter的な手法が使いにくいので、実質的に低解像度になる。(詳細は後述)
近似の限界：太陽表面(>0.99Rsun)では、音速と熱対流速度が同程度になるので、Anelastic近似の仮定自体が破綻。表面が取り扱えない。

SS研

18音速抑制法(RSST)
密度の擾乱ρ1をξ2倍する(ρ0は時間変化しない背景密度場)

ρ = ρ0 + ρ1

ρ = ρ0 + ξ 2 ρ1

この操作によって実効的な音速は1/ξになる。連続の式を以下のように変形したことと同等になる。

@⇢1=@t1r·(⇢v)⇠2

この手法ならば、局所通信のみが必要となる。この方法の妥当性は我々の過去研究で確認されている (Hotta+2012)。
抑制された音速で測ったマッハ数が0.7よりも小さかったら熱対流に対して影響がないことが確かめられた。
非一様なξを使ってもよく。対流層の底で大きなξ、太陽表面近くではξ=1を用いることで物理を失うことなく太陽表面まで解ける。
SS研

19

解く方程式@⇢1= r · (⇢v)
@t@v⇢
= ⇢(v · r)v@t

rp1⇢1 ger

1+2⇢v ⇥ ⌦0 +(r ⇥ B) ⇥ B
4⇡@s1⇢T= ⇢T (v · r)s +@t✓ ◆✓ ◆@p@p
p1 =⇢1 +s1@⇢ s@s ⇢@B= r ⇥ (v ⇥ B)@t
⇢ = ⇢0 + ⇠ 2 ⇢1

球対称な静水圧平衡モデルρ0,P0, T0を用いる。
球対称場は、太陽標準モデル(Model S)を使う。
そこからの擾乱場のみを解く。
対流層の中だと、

⇢1
p1
T1
vc2
⇠⇠
⇠ 2 ⇠ 10 6
⇢0p0T0cs
状態方程式は、電離の効果を入れるために背景場より計算しておいた値を使う(線形近似)。
エントロピーの方程式を解くModel Sは以下からダウンロードできる。
http://astro.phys.au.dk/~jcd/solar̲models/
SS研

20

解く状況
恣意的に与えるエネルギーフラックス
総量は太陽放射

太陽対流層
計算領域

現状では、計算資源の制約から太陽表面まで一貫して解くことは難しいのである程度のところに人工的な上部境界条件をおいて、そこを冷やして熱対流を駆動している。基本的には、太陽放射が注入されて、太陽放射が抜かれるシステム。
対流層の底

太陽表面

空間スケール60,000 km

300 km

速度

100 m/s

1,000 m/s

1週間5分

時間スケール

太陽標準モデルで与えられるエネルギーフラックス

空間スケールは密度スケール長

Fe = cp ⇢hvr0 T 0 i
SS研

21数値アルゴリズム
時間積分：メモリ節約の４段Runge-Kutta@U= R(U)@tt
U1/4 = U0 +R(U0 )4t
U1/3 = U0 +R(U1/4 )3U1/2 = U0 +U1 = U0 +t
R(U1/3 )2tR(U1/2 )空間微分：４次中央差分✓@u@x◆

1=( ui+2 + 8ui+112 xi

この手法だとクーラン数がお得。

t = CCFL
SS研

8ui 1 + ui 2 )x, CCFL < 2
cs + vc + cA
22

人工粘性
4次精度の中央差分では、ほとんど数値粘性が入らないために乱流問題を解くと不安定化(カスケードの先で散逸させる機構は必要)
そこで人工粘性を採用。風上法で１次精度と２次精度を切り替えるのと原
理的には同じ(Rempel+2009, Rempel, 2014参照)

@Qi
=@t

Fi+1/2Fi 1/2x

Fi+1/2 =1ci+1 x (ur2ul , ui+1ui )(urul = ui + 0.5 ui
ur = ui+1

0.5 ui+1(ui+1ui = minmod(h= max 0, 1 + h=0ui 1 )2⇣, ✏(ui+1ur ului+1 uiui ), ✏(ui

⌘i1(ur(urSS研ui 1 )ul ) · (ui+1ul ) · (ui+1ui ) > 0ui )  023ul )人工粘性の理念 (1/2)
Fi+1/2 =1
ci+1 x (ur2ul , ui+1ui )(urul )

@Q
@F=@t
@x
Qni 1 + Qni+1Fi+1 Fi 1+ t22 xnnnQ2Q+QFi+1 Fi 1i1ii+1nQi ++ t
22 xnn2 Qn
2Q+Q(x)Fi+1 Fi 1i 1ii+1n
Qi +t+
22 t( x)2 x

Lax-Friedrich法(空間時間１次精度)を基本に改訂Qn+1i===粘性係数(Δx)2/(2Δt)の粘性をかけているのと同じ。CFLの限界までかける( x)2( x)2t<! ⌫<2⌫
2 t
SS研

24

人工粘性の理念 (2/2)
Fi+1/2 =1ci+1 x (ur2ul , ui+1ui )(urul )x( x)21t=! ⌫<= c x
c2 t
22 Qn(x)i 1nQn+1=Q+ii
2 t

時間4次精度に変える

2Qni + Qni+1Fi+1 Fi 1t+x2 x

ローカルに粘性係数を決めてminmod関数で最適化

空間4次精度に変える

4次精度を用いて、最大でLax-Friedrichの粘性フラックスを用いることで、安定で高解像度なスキームを達成する。
SS研

25

粘性フラックスの最適化 (1/2)
ul = ui + 0.5 ui
0.5 ui+1(ui+1ui = minmod

Fi+1/2 =ur = ui+1ui 1 )2

1ci+1 x (ur2, ✏(ui+1ui ), ✏(ui

ul , ui+1ui )(urul )ui 1 )

Generalized minmod limiter: 英語wikipedia flux-limitter参照

ulururuli-1ii+1i+2i-1ur ul=1ui+1 uiii+1i+2ur ul<1ui+1 ui

この状況のときには、Lax-Friedrichの粘性を最大までかけて安定化

この状況のときは、非常に小さい粘性粘性がまったくかからないこともある

SS研

26

粘性フラックスの最適化 (2/2)

1ci+1 x (ur2Fi+1/2 =(h= max 0, 1 + h=0

⇣

ur ul
ui+1 ui

ul , ui+1⌘i1(ur(urui )(urul ) · (ui+1ul ) · (ui+1ul )ui ) > 0ui )  0

h>1のとき、非常になめらかな構造に対しては、人工粘性を一切かけない。ダイナモ計算ではh=2を使っている。最大値はなく、磁気プラントル数がとても小さい計算では、h=200程度を用いても問題はない。

SS研27

人工粘性の成果
102
Eφ [km3 s-2]

10010-2MURaM
Stagger

10-4
10-6

Ekin [erg cm-3]

10-8
106105
104

101-1.5-6103
102
101
100
0.001

0.010
0.100
k/2π [Mm-1]

1.000

拡散領域が4-5点で分解できる。

102103104
spherical hamonic degree (l)105Hyperdiffusivityというのはたくさん種類がある。6次精度中央差分を用いている
Staggerコードの計算と比較すると、一方向の格子点が1/4程度でも同じようなスペクトルが再現できている。つまり3次元計算ならば256倍得をする。
SS研28

Yin-Yang grid
(a)

(b)

片方の座標の境界で3次精度補完
オリジナルと違う(Kageyama & Sato, 2004)
SS研

29

データ取り扱い

とりあえずの結果を見るために、毎出力ごとに3次元分割のデータを動径方向(重力方向)の1次元分割にリマッピング。一つの水平データを一つのファイルに出力。
また、エネルギーフラックス、平均二乗速度などの統計量を出力とともに計算して、一つのファイルに出力。ダウンロードすべき情報を減らす。
SS研

30

並列化
・三次元領域分割でのMPI・京を使うときはノード内はコンパイラによる自動並列化
・京を使うときはFJMPIでノード配置を調べて配置(Fortranで使えない？)
京でないときはPeano-Hilbert曲線で配置。

SS研

31

音速抑制法を用いたコード

192x128x128 grid/node (weak)

Performance to the peak(%)

1012

cell updates/sec

96x64x64 grid/node (weak)
96x64x64 grid/node

1011@82944 node (strong)

1010

109
103

451010No. of Cores

10630

20100103

104105
No. of Cores106

ほぼ理想的なスケーリングを京全システム(663,552 core)まで示している。
理論ピーク性能に対して24%の実行性能を示している。
SS研

32

これまでの計算との違い

Hotta+2014

Miesch+2008
SS研

33

これまででもっとも大きい計算
4096 32400 32352=4.3兆格子点確認した限り世界でもっとも大きい流体力学計算
京全システムを8時間使用(82944node/663552core)
1.14 PBのメモリを使用(京は1.26 PB)
2.54 PFLOPSを達成！

SS研

34

成果１：表面勾配層

1960年頃に発見された太陽表面付近にある角速度の勾配層を世界で初めて数値計算で実現し、その生成・維持機構を明らかにした。
SS研35

成果２：効率的な小スケールダイナモ
vr in MHD

Br in MHD

Emagnetic/Ekinetic

3.0

1.02.5High resolution

0.8

2.0
Ekin/Emag

Beq, Brms [104 G]

Superhigh resolution

Typical Global

vr in HD1.50.60.41.0

Low resolution0.20.50.00.00.75

0.80

0.85
x/Rsun0.900.95

0.75

0.80

0.85
x/Rsun0.900.95(Hotta+2015b)

これまでのダイナモ計算では低解像度ゆえに、熱対流の運動エネルギーに対して、数％の磁場エネルギーしか達成できていなかったが、領域を制限して、非常に高解像度にすると、運動エネルギーに対して90%程度まで磁場が増強された。
SS研36

今日の話の動機の復習
大きな問題は、解像度をあげればあげるほど大規模磁場は弱くなっていき現実からずれていくこと。
太陽ダイナモ計算の暗黙の了解は解像度をあげすぎないこと。
当然、これでは現実の太陽を理解できたことにならないので、なんとか打開策を探る。
現在のところ最高解像度の計算では乱流の慣性領域はほとんど解けておらず、小スケールのダイナモ効果が効率的ではない。
太陽全球を扱っても、小スケールのダイナモ効果が効率的になるような高解像度を達成することを目標とする。

SS研

37

解像度ごとの比較
粘性が低いLow̲D
ν=κ=1 1012 cm2 s-1

Medium64x192x384

High
256x768x1536

vr at
0.95R

Bφ at
0.72R
Low̲D (Fan+2014)では対流層の底付近でコヒーレントな磁場が生成されている。すこし実効的な解像度をあげると(Medium)大規模な磁場はあまり形成されていないように見える。
さらに解像度をあげると、再度大規模な磁場が生成されているように見える。より詳細に調べるためには蝶形図を見る必要あり。
SS研38

京で挑んだ高解像度ダイナモ計算<Emag>: 平均磁場のエネルギー<Bφ> at r=0.72Rsun

Less diffusiveLow̲D
ν=η=1 10 12 cm2 s-1<Emag>=2.9 10 4 erg cm-3

Medium64x192x384<Emag>=1.3 10 4 erg cm-3

High
256x768x1536<Emag>=2.4 10 4 erg cm-3高解像度では高いレイノルズ数にもかかわらず、コヒーレントな磁場、そしてその周期が観測される。
SS研39

なぜ大規模磁場が再度つくられるのか
高解像度にして、小スケールダイナモが効率的になることになることによって、小スケール磁場エネルギーが小スケール運動エネルギーを超えるほどになる。

Ekin and Emag [erg cm-3]

107106105104103
1MediumHigh
10100
spherical harmonic degree l

実線：運動エネルギー
破線：磁場エネルギー

その結果、ローレンツ力が小スケールの運動を抑制する。
もともと小スケールの運動は大規模な磁場を壊す向きに働きがちなので、その破壊運動を抑制することにより、α効果・Ω効果がもう一度効くようになり大規模磁場が構築される。
SS研40

Ekin and Emag [erg cm-3]

さらなる高解像度化
107106105
104 High-S
High103
1
10100
1000
spherical harmonic degree l

最高解像度の途中結果を利用して、補完して短時間のみ計算。小スケールダイナモはさらに活発になり、小スケール流れ場はさらに抑えられている→限界はどこに？
SS研

41

さらなる未来への高解像度へ
http://solarscience.msfc.nasa.gov/images/magbfly.jpg

実際の太陽周期は規則正しく綺麗
現在のところのシミュレーション
結果は、周期の長さもまちまちで対称性も悪い。
太陽のような規則正しい周期は、非常に高粘性・高磁気拡散の時に達成される。太陽では小スケールダイナモが活発であることの証拠か？→さらなう高解像度計算
SS研42

さらなる未知の領域へ
太陽で観測される磁場周期は黒点数の周期であるので、黒点を再現できないことには解決できない。
(例えば、太陽全放射の変動は11年周期で0.1%程度であり磁場エネルギーが11年周期でどう変動しているかはよくわからない)
黒点を再現するためには、太陽表面である光球まで一貫して解く必要がある。
光球の放射を取り入れる必要があるが、Abbett+2012で簡単ないい方法が提案されているので、線形近似している状態方程式をテーブルで読み込むのとのコンビネーションにすることで、一貫して取り扱える。
2016年中に対流層全部200 Mm 300 Mm 300 Mmの計算を格子間隔100 km 200 km 200 kmで実行予定(したい)。
100日分の計算に現在のコードだと100万ノード時間＠Kほど必要。
SS研

43

まとめ
これまでの太陽ダイナモ計算の一般的な手法であったアネラスティック近似＋スペクトル法を抜本的に見直し音速抑制法＋有限差分法を用いた方法に変えることで、現在のところ世界最高解像度の太陽熱対流・ダイナモ計算を実行することに成功している。
その結果、これまでの謎であった表面勾配層を解決するとともに、小スケールダイナモを低く見積もっていたことを指摘した。
小スケールダイナモは大スケールダイナモを理解する上で重要であることも新たに発見した。
今後はさらなる高解像度化に挑戦するとともに、光球を目指した計算を実行予定である。
SS研

44