2023.10.27 科学技術計算分科会 2023年度会合 富岳スペシャル4.0 〜拡がる応用〜 @神戸国際会議場

富岳で迫る太陽風の起源
飯島 陽久 名古屋大学 宇宙地球環境研究所/高等研究院

太陽風とは？
太陽から吹き出す太陽風■ 太陽から吹き出す超音速のプラズマ流
20世紀中ごろに発見
以降、その形成機構が長年の謎
太陽の方向
プラズマの尾(太陽風の方向)
ダストの尾(彗星の軌道方向)©ESA/NASA
SOHO/LASCO Coronagraph

彗星が持つ2つの尾と太陽風 (Biermann & Lust, 1958)

なぜ太陽を研究するのか
磁気嵐
太陽宇宙線

オーロラ

コロナ質量放出

惑星宇宙天気電離圏擾乱

原始惑星系円盤

惑星大気
ダイナモ

フレア
放射光

太陽太陽風

(1) 最も近い恒星「太陽」の理解(2) 太陽圏内の環境変動への理解(3) 宇宙天気予報の精度向上©NASA宇宙天気とは？
■ 宇宙天気: 太陽活動が生み出す地球磁気圏・電離圏擾乱の総称

電波通信の障害や宇宙線被爆など社会システムへ影響; 予測精度の向上が研究課題

本研究の概要■ これまでの研究: 太陽風のみ計算

太陽風の磁気乱流

太陽内部(~100 km)から太陽風(~1千万 km)まで大きな時間・空間スケールのギャップが存在
密度も約20桁変化、物理量のギャップも大きい
磁気乱流を分解するため格子点数も必要

エネルギー注入⇒ コスト、安定性のため一貫した計算が困難■ 今回の研究(1) エネルギー整合性の高い計算手法を開発と(2) スーパーコンピュータ「富岳」の利用により太陽内部から太陽風までの一貫した計算が実現
恣意的な仮定なく観測と整合的な太陽風を再現

太陽表面の磁気乱流

目次● はじめに● (1) 計算手法面での工夫: エネルギー整合磁気流体スキーム● (2) 実装面の詳細: 並列化、高速化● (3) 科学成果: 「富岳」による太陽風の輻射磁気流体計算● まとめ(1) 計算手法面での工夫:

エネルギー整合磁気流体スキーム

エネルギー整合磁気流体スキームの概要■ 格子を使った微分方程式の数値解法

fj−2 fj−1

fj

fj+1 fj+2( f′)j ≈ ( fj+1 − fj−1)/2■ 従来手法の問題点( f ⋅ g)′ ≠ f′ ⋅ g + f ⋅ g′
⇒ 計算が破綻(温度が負に)異なる場所への輸送
熱エネルギー
エネルギーの交換
運動エネルギー

磁気エネルギー■ 提案手法 (Iijima, 2021, JCP)

( f ⋅ g)′ = f′ ⋅ g + f ⋅ g′
⇒ 整合性・安定性が向上

プラズマ-磁場間のエネルギー交換

太陽は巨大なプラズマのかたまり■ 通常時の原子■ プラズマ(電離気体)

正の電荷を持った陽イオンと負の電荷を持った電子が結合

クーロン力

+

粒子の持つエネルギーが大きく、陽イオンと電子がバラバラに運動-

高温・低圧+

プラズマの数値モデリング: 粒子 v.s. 流体■ プラズマを粒子として近似
短周期のプラズマ波動まで考慮
計算量: 大きい

主な用途: 波の分散性・散逸過程
トカマク型核融合装置©Max Planck Institute for Plasma Physics

粒子スケール ≈ 空間スケール■ プラズマを流体として近似
長周期のプラズマ波動のみ考慮
計算量: 小さい

本研究

主な用途: 時空間構造、エネルギー論
太陽大気の観測画像

Cranmer & Weinbarger (2019)

粒子スケール ≪ 空間スケール

磁気流体(MHD)方程式: 流体プラズマの基礎方程式∂ρ
質量保存則:+ ∇ ⋅ (ρV) = 0∂t

黒: 通常の流体
赤: 磁気流体

ローレンツ力 粘性力 アンペールの法則∂
運動方程式: ρ+ V ⋅ ∇ V = − ∇P + J × B + ∇ ⋅ (ν ∇V), J = ∇ × B( ∂t)
オームの法則∂B
電磁誘導の法則:+ ∇ × E = 0, E + V × B = ηJ
∂t

粘性加熱

ジュール加熱∂eint22
熱エネルギー方程式:+ ∇ ⋅ (eintV) = − P∇ ⋅ V + ν ∇V + η ∇ × B∂t
磁気流体は通常の流体とよく似ている ⇒ (圧縮性)流体の数値解法が利用可能

レイノルズ数が非常に高く境界条件が単純な太陽の場合、構造格子上の圧縮性乱流の数値解法が向いている素直に各保存則を有限差分法・有限体積法等で解くと…

何もしてないのに計算が壊れる

負の温度による計算の破綻■ 多くのMHDコードでは全(熱+運動+磁気)エネルギー保存則を解く∂ (eint + ρV /2 + B /2)
22∂t+ ∇ ⋅ [(eint + P + ρV /2) V + E × B] = 02

■ 負の温度が発生する状況: 熱E ≪ 運動E+磁気E
eint = (eint + ρV /2 + B /2)−(ρV /2 + B /2)2222= (小さい熱E + 大きい運動・磁気E) − (大きい運動・磁気E)= 熱Eに対する数値誤差の発生 (負の温度の発生)

太陽大気では「磁気E ≈ 105 × 熱E」。容易に計算が破綻する。

どのような離散化手法が必要なのか■ エネルギーを厳密に保存

太陽風の磁気乱流

太陽表面の磁気乱流が持つエネルギーのうちほんの10万分の1が太陽風を加速
残りは太陽大気を加熱■ 極限状況でも安定に計算可能

エネルギー注入- 低マッハ数 (音速 ≫ 流速) の太陽内部- 強磁場 (磁気E ≈ 105 熱E) の太陽大気- 超音速 (音速 < 流速) の太陽風

■ 超高レイノルズ数の乱流を解像
太陽はレイノルズ数 10

10以上、非常に乱流的表面から太陽風まで伝播する波を散逸せず追跡

太陽表面の磁気乱流

MHD方程式のエネルギー整合性
磁気流体方程式には、3つの熱・運動・磁気エネルギーが存在。
これらはそれぞれ輸送項と相互作用項で時間発展。

∂eint+ ∇ ⋅ Fent = −WP + Qvis + Qres∂t2∂ρV /2+ ∇ ⋅ Fkin = WP + WL − Qvis
∂t2
∂B /2+ ∇ ⋅ Fmag = −WL − Qres∂tQvis2

ρV /2Fenteint
WP

FkinWL

計算が壊れる原因がエネルギー整合性の破綻に起因するなら、エネルギー整合性を保てる離散化であれば安定に計算出来るはず。
物理的にも良い性質だし、不要な数値拡散も抑えられる。

Qres
2B /2Fmag

離散的な積の微分法則とエネルギー保存■ 解析的な積の微分法則:∂V∂U ∂UV
U+V=∂x∂x∂x(UV)j+1/2 の空間2次精度近似■ 離散的な積の微分法則(2次精度中央差分法版):

Uj

Vj+1 − Vj−12Δx+ Vj

Uj+1 − Uj−12ΔxUjVj+1 + Uj+1Vj Uj−1Vj + UjVj−11=−Δx (22)

■ 線形スカラー移流方程式の離散的なエネルギー保存則:∂Uj+aUj+1 − Uj−1=0⇒2∂Uj /2Δx2(閉じた境界条件なら)総エネルギーが保存: Σj(Uj /2)Δx = cnst .
∂t2Δx∂t+(aUjUj+1 /2) − (aUj−1Uj /2)=0エネルギー整合MHDスキーム (Iijima, 2021, JCP)
以下の離散化で熱・運動・磁気エネルギー間の整合性を保持可能:δρVj
δ1Dj(ρ)∂ρ+=∑j δ1xj∂t ∑j δxj
δBj δBiδ1Dj(ρVi)∂ρVi
1 δρViVj Vi δρVj ρVj δViδP+++=−−Bj−+∑j ( 2 δxj∂t2 δxj
2 δxj )
δxi ∑j ( δxiδxj ) ∑j δ1xj
δ (eint + P) Vj
δ1Dj(eint)∂eint
δP+=Vj+ Qvis + Qres +∑j∑j δxj∑j δ1xj∂tδxj
δ(Vj Bi − ViBj)
δ1Dj(Bi)∂Bi+=∑j∑j δ1xj∂tδxj

粘性・ジュール加熱率 Qvis, Qres は
拡散フラックスがある条件を満たすと
δu/δx: 中央差分演算子(コロケート格子上、任意精度) 正値性を保つ(熱力学第二法則)

δ1u/δ1x: 隣り合うセル境界面で定義された物理量の2次精度中央差分
D(U): セル境界面で定義された物理的・数値的な任意の拡散フラックス

エネルギー整合MHDスキーム: 数値実験結果■ MHD点源爆発問題■ MHD渦問題

磁気エネルギー・熱エネルギーの比が1千万倍以上という極限状況でも安定

熱・運動・磁気エネルギーの輸送・交換を離散化後も厳密に表現

熱エネルギー/磁気エネルギー

運動エネルギーの収支関連する離散化手法■ 移流項の離散化を工夫することによるエネルギー保存解法は非圧縮性乱流計算 (Harlow & Welch 1965; Morinishi et al. 1998)で提案:∂V+ (V ⋅ ∇) V = − ∇P, ∇ ⋅ V = 0∂tこの方程式系ではエネルギー方程式がなく、離散的なエネルギー保存が重要■ 類似する構造保存数値解法:↓ KEEPスキームとして「富岳」で利用

圧縮性流体 (Kok et al., 2009; Morinishi et al., 2010; Kuya et al., 2018),
プラズマ運動論 (Idomura et al., 2007),

非圧縮性MHD (Ni et al., 2007), 二温度流体プラズマ (Shiroto., 2018).

提案手法は、これらの圧縮性MHD版 (Iijima, 2021)と位置づけられる。

(2) 実装面の詳細: 並列化、高速化

RAMENS (Iijima & Yokoyama, 2015/2017; Iijima et al., 2023)
RAdiation Magnetohydrodynamics Extensive Numerical SolverFortran 2003 + MPI + OpenMP (依存ライブラリ: FFTW)
圧縮性磁気流体(MHD)方程式 on 直交曲線座標系
現実的な状態方程式・不透明度(2次元テーブル化)
多次元輻射輸送(低層大気)+局所輻射損失(高層大気)
磁力線に沿った硬い熱伝導

エネルギー整合MHDスキームによるロバスト性向上 (Iijima, 2021)
6次精度中央差分法 + 人工粘性
4段2次精度Jameson-Baker Runge-Kutta法
彩層・コロナ界面の数値的処方 (Iijima & Imada, 2021)
Super-Time-Stepping法による硬い熱伝導の陽的高速化 (Meyer et al., 2014)
準相対論近似による時間刻みの緩和 (Gombosi et al., 2002; Rempel, 2017)
短い特性線法による空間局所性の高い輻射輸送解法 (Auer & Mihalas, 1973)

プロセス並列化: 空間領域の分割■ プロセスごとに空間領域を3次元分割、袖領域をMPI通信
並列化効率を高く計算したい ⇒ 可能な限り隣接通信のみのアルゴリズムを選択

全計算領域

1プロセス1プロセス

1プロセス

1プロセス

1プロセス

1プロセス

1プロセス

1プロセス

1プロセス

yx

広いパラメータレンジに由来する位相速度のばらつき
流速

典型速度 [km/s] 陽解法の時間刻み [s]プラズマ波動流速

1-10010-0.1

プラズマ波動1000-10000

0.01-0.001

熱伝導-<10-5

熱伝導■ 解決策1: 陰解法の利用■ 解決策2: 位相速度の人工的な抑制

プラズマ波動や熱伝導は陰解法

方程式や離散化を変更し位相速度を抑制

流速に合わせて時間刻みを設定⇒ グローバル通信の発生(並列化効率の悪化)

物理的性質を損なわないように⇒ 陽的に積分可能、隣接通信のみでOK

プラズマ波動の位相速度抑制∂
運動方程式: ρ+ V ⋅ ∇ V = − ∇P + ( ∇ × B) × B + ∇ ⋅ (ν ∇V)( ∂t)∂+ V ⋅ ∇ V = − ∇P + ( ∇ × B) × B + ∇ ⋅ (ν ∇V)位相速度の制限: ξρ( ∂t)ξ: 慣性の増幅率(実際にはテンソル形式)
流速
本来のプラズマ波動
抑制後のプラズマ波動■ 慣性の増幅による位相速度抑制
定常解を保ったまま時間刻みを伸ばすことが可能
非定常解への適用可否は問題依存

「硬い」熱伝導方程式の陽的高速化
離散的性質 計算コスト
良大■ 陰解法: 情報伝達速度 → ∞ (グローバル通信が必要)
Tn+1n−Tn+1= ∇ ⋅ (κbb ⋅ ∇T )Δt■ Super-Time-Stepping法: N段陽解法でCFLを緩和ΔtΛ=∇ ⋅ (κbb ⋅ ∇T),T(計算量: 𝒪(1/Nargmax
μj s.t. Πj∑
μ1,...,μN j

I − μjΛ ≤ 1NN)倍)本研究で利用■ 双曲型方程式への変形: 最大位相速度 c で制限
悪小2∂T∂q c= ∇ ⋅ (qb),=(κb ⋅ ∇T − q)∂t∂tκ熱伝導計算におけるテーブル型状態方程式の回避■ 状態方程式は2次元テーブルを補間

ek+1 Tj−1,k+1
ek

Tj,kek−1 Tj−1,k−1
ρj−1

Tj,k+1
Tj,kTj,k−1ρjTj+1,k+1
Tj,kTj+1,k−1
ρj+1■ Super-Time-Stepping法での熱伝導計算
熱伝導式を約30億ステップ計算(ρ = const.)∂e= ∇ ⋅ F(T) ⇒ 状態方程式がボトルネック∂t■ テーブル状態方程式の回避
方程式を連立方程式に変形
状態方程式の評価は初期値の1回のみ初期値: T′(t = 0) = T(t = 0)

STS法:

∂e
∂T′
= ∇ ⋅ F(T′)( ∂T ) ∂tρ
∂e= ∇ ⋅ F(T′) 注: 右辺は同一∂t「富岳」に向けたコードの最適化■ 定番の手法で試行錯誤「富岳」関連の公開資料、教科書、コード読み等で勉強・実験の繰り返し■ 基本形:
xループをベクトル化、yzループをスレッド並列■ SPILLの回避
必要に応じて一時変数を導入、最内ループを分割■ 指示行・コンパイラオプション・環境変数
意図どおり解釈されない場合などで利用■ 最終的な性能
FLOPSピーク性能比: 6-7%並列化効率: 97.1% (192-768ノード間で弱スケーリング)(3) 科学成果:「富岳」による太陽風の輻射磁気流体計算

太陽風とは？

再掲

太陽から吹き出す太陽風■ 太陽から吹き出す超音速のプラズマ流
20世紀中ごろに発見
以降、その形成機構が長年の謎
太陽の方向
プラズマの尾(太陽風の方向)
ダストの尾(彗星の軌道方向)©ESA/NASA
SOHO/LASCO Coronagraph

彗星が持つ2つの尾と太陽風 (Biermann & Lust, 1958)

太陽風のエネルギー源 = 太陽内部の熱対流■ 太陽内部の熱対流が持つ2つの役割:(1) ダイナモ作用による磁場構造の形成(2) 各種プラズマ波動の励起
可視光放射

Cranmer & van Ballegooijen (2005)

表面磁場強度

プラズマ波動の励起・伝播
太陽表層のダイナモ作用の数値シミュレーション (Rempel, 2014)これまでの研究: 太陽風のみを計算領域とする数値計算■ 太陽内部から太陽風までのパラメータギャップ

太陽内部から太陽風まで大きな時間・空間スケールのギャップが存在
密度も約20桁変化、物理量のギャップも大きい
高レイノルズ数(Re > 1010)の乱流を分解するため格子点数も必要⇒ 計算コスト、安定性のため一貫した計算が困難
太陽内部
太陽表面
太陽大気
太陽風

パラメータギャップ

時間スケール [分]

空間スケール [km]典型密度 [1/cc]

104104102111021017

103

1041010104107102~104~105~104

本研究: エネルギー源である太陽内部からの包括計算■ 今回の研究

太陽風の磁気乱流

太陽内部から太陽風までの一貫した計算
熱対流の性質は太陽の明るさ・重力・元素比が決定
恣意的な仮定を排除して太陽風を再現出来るか検証■ 利点

エネルギー注入◯ 磁場構造・プラズマ波動が自発的に形成◯ 上部境界の悪影響がない(計算領域外への流れ)■ 欠点✕ 数値計算コストの増大 ⇒ 「富岳」で解決✕ パラメータギャップ ⇒ 計算手法の改良で解決

太陽表面の磁気乱流

計算設定: 水平方向に局所的、太陽内部から太陽風まで計算領域太陽表面の典型的な水平構造を包含しつつ、対流層上部から太陽風加速の完了までを一貫して追跡■ 格子点数: 256×256×3072 (~2億)■ 時間ステップ数(MHD): ~2千万29R⊙

太陽風磁力線■ 時間ステップ数(熱伝導): ~30億
■ 計算資源: 富岳, 不老 (名大; 最適化), ATERUI II (NAOJ; 初期条件)■ 利用ノード数(富岳): 768ノード (3072プロセス)

0.07R⊙
対流層(太陽内部)

太陽風速度@r = 29R⊙

強い時間変動黒: 水平平均太陽風密度

灰: 水平方向の最大/最小

太陽風鉛直磁場

磁場の極性反転
乱流的・間欠的な太陽風

大気からのX線強度

太陽風の速度構造

ループ構造の形成、突発的発光

表面磁場

ダイナモで作られた微細磁場構造

太陽各層における典型統計量をほぼ全て再現
太陽風速度@r = 29R⊙

■ 表面磁場強度: ~50 G■ 大気X線強度: 観測平均値程度

太陽風密度■ 太陽風速度: 300-600 km/s
■ 太陽風磁極反転の再現-

太陽風鉛直磁場

X線強度

太陽風速度

表面磁場反転の頻度は1桁不足(おそらく格子点数不足)■ 太陽風質量損失率: 1-4x平均値

恣意的な仮定なく自発的に太陽風を再現(世界初の第一原理的計算)新たな科学的知見: 磁気リコネクションの重要性■ 磁気リコネクションとは(電気抵抗がほぼゼロの状態で)反対向きの磁力線が急激につなぎ変わる現象
磁場から熱・運動量への効率的な変換方法
太陽フレア(観測)数値モデル

磁気リコネクションの概念図■ 太陽物理における磁気リコネクション
太陽フレア(爆発的発光)をはじめとして磁気リコネクションが重要とされる現象は多い
Tsuneta+ (1996)

Yokoyama+ (2001)

しかし、太陽風に寄与するかは未解明だった本数値モデル内では磁気リコネクションが頻繁に発生
プラズマ温度磁力線(色: プラズマ温度)

磁気リコネクション由来の上昇流
鉛直速度

黒線: 太陽風へ寄与する領域の境界面

ダイナモによる磁気ループ構造形成

磁気ループと開いた磁力線の間でリコネクション

数値モデル全体でのエネルギー収支
15% (太陽風へ)
太陽風加速

25% (X線放射)

15%

20%
(磁気リコネクションによる輸送)
60% (X線放射)

開いた領域
20%

大気加熱(X線放射)
25%

20%

閉じたループ80%

大気へのエネルギー注入量: 100%
開いた磁力線

閉じた磁気ループ
コロナ底部

20% (開いた領域へ)

80% (磁気ループへ)■ 太陽風が持つエネルギーの約半分が磁気リコネクション由来

今後の展開: 太陽全体としての太陽風形成過程の理解■ 現在進行中のプロジェクト

高解像度化と小スケール構造の寄与
全球スケール磁場構造の影響

黒点(超強磁場領域)からの太陽風
太陽の全球スケール磁場構造Cranmer et al. (2017)(1) 太陽全球での太陽風形成過程の理解⇒ (2) 過去の太陽、他の恒星への応用(3) 宇宙天気予報の精度向上

太陽の超強磁場領域(黒点)「富岳」成果創出プログラム「シミュレーションとAIで解き明かす太陽地球環境変動」■ 大規模シミュレーションと機械学習を連携させて計算を高速化・パラメータ空間充填
期間: 2023〜2025年度
飯島はサブ課題B代表

https://sites.google.com/view/fugaku-heliophysまとめ: 富岳で迫る太陽風の起源

SOHO/LASCO Coronagraph©ESA/NASA● 計算手法面での工夫: エネルギー整合磁気流体スキーム

離散的な「積の微分法則」を利用して方程式間の整合性向上
磁気エネルギーが熱エネルギーの1千万倍という極限状況で安定● 実装面の詳細: 並列化、高速化

パラメータギャップにより「硬い」方程式が出現
方程式レベル・離散化レベルでの実効位相速度抑制による陽解法での高速化● 科学成果: 「富岳」による太陽風の輻射磁気流体計算

恣意的な仮定を置かずに観測と整合的な太陽風を再現(世界初の第一原理的計算)
太陽風が持つエネルギーの約半分が磁気リコネクション由来であることを発見