相対論的流体力学における安定・高精度な数値手法
Based on arXiv:2306.12696

Speaker: 村瀬功一 (京大基研)
Nathan Touroux (阪大理、SUBATECH, 京大基研)
北沢正清 (阪大理、京大基研)
Marlene Nahrgang (SUBATECH)
熱場の量子論とその応用 (TQFT2023)、KEK、つくば
2023/8/29

1大自由度(場)・非線形・実時間発展

相対論的流体力学における安定・高精度な数値手法
Based on arXiv:2306.12696

Speaker: 村瀬功一 (京大基研)
Nathan Touroux (阪大理、SUBATECH, 京大基研)
北沢正清 (阪大理、京大基研)
Marlene Nahrgang (SUBATECH)
熱場の量子論とその応用 (TQFT2023)、KEK、つくば
2023/8/29

2

相対論的流体力学の数値計算

PartX (1/1)

高エネルギー重イオン衝突反応
QGP

nuclei

Lorentzcontraction

流体の発展

Hadrons

STAR Collaboration

QCD相図(予想図)

中性子星連星合体
想像図

数値相対論= 計量gμνの発展+ 流体Tμνの発展2023/8/29

http://www.astronomy.com/news/2017/10/

3数値流体力学(CFD)の世界にようこそ(?)

2023/8/29

4計算効率と数値精度

PartX (1/1)

数値計算のため時間・空間を離散化 → Δt, Δx
空間メッシュ

CFL条件
Δt < cΔx

計算量 ~ O(1/Δx3 Δt) > O(1/Δx4)→ 離散化幅 Δx は計算効率に強い影響
動機・目標: できるだけ大きな Δx で計算量を抑えつ
つ (Δx, Δt) に対する近似の次数を上げて計算したい+ 計算の安定性。
2023/8/29

5空間高次化 & 時間高次化

PartX (1/1)

空間高次化(空間再構築)
周囲のセルからセル内の流体場のプロファイルを再構築(有限体積法)MUSCL, PPM, PRM, ENO, WENO, …
線形補間+minmod

愚直な高次化

時間高次化

2次補間 有理関数+minmod+minmod

複数のn次補間+選択 or 重み付け和→ 振動 (衝撃波や急激な物理量の変動で✕)→ Flux limiter (～あぶないとき自動で次数を下げる)
時間1ステップでステップ幅 Δt の高次まで一致させる

Euler法, TVDルンゲクッタ, 陰的ルンゲクッタ, …
高次化: 解が時間的・空間的になめらかでない時 → ??

keywords: 線形安定性、全変動減少(TVD)、…
2023/8/29

6空間高次化 & 時間高次化

PartX (1/1)

空間高次化(空間再構築)
周囲のセルからセル内の流体場のプロファイルを再構築(有限体積法)

MUSCL, PPM, PRM, ENO, WENO, …
線形補間+minmod

愚直な高次化

時間高次化

2次補間 有理関数+minmod+minmod複数のn次補間+選択 or 重み付け和→ 振動 (衝撃波や急激な物理量の変動で✕)→ Flux limiter (～あぶないとき自動で次数を下げる)
時間1ステップでステップ幅 Δt の高次まで一致させる

Euler法, TVDルンゲクッタ, 陰的ルンゲクッタ, …
高次化: 解が時間的・空間的になめらかでない時 → ??

keywords: 線形安定性、全変動減少(TVD)、…
2023/8/29

本研究7陰解法と線形安定性

PartX (1/1)

常微分方程式→線形化&対角化

Euler法

2023/8/29(Dahlquist eq)

陰的Euler法

8陰解法と線形安定性

PartX (1/1)

常微分方程式→線形化&対角化

Euler法

(Dahlquist eq)

陰的Euler法

硬い方程式(stiff equation) = λmax 大な方程式

陽解法は硬い方程式で☓2023/8/29

陰解法は線形安定9アイディア

2023/8/2910我々の手法(アイディア)

PartX (1/1)

Implicit (陰解法)
1段Gauss-Legendre法(2次)vs比較用: Explicit (陽解法)
Heun法(2次)

遅い? → 2つのポイント① k について解く必要がある (f は”流体方程式”)→ Fixed point 法: 収束するまで逐次 k を代入し直す(最初の k の guess として前ステップの k を再利用)② 空間の領域によって収束の速さが著しく異なる→ 早く収束した領域は逐次更新をスキップ※収束の速さは方程式の局所的な硬さに関係
収束判定 = 自動硬さ検出
2023/8/29

11数値精度&計算コストの検証

PartX (1/1)★比較:数値計算結果 vs 既知の解析解 ←定性的な確認～ Δx → 0, Δt → 0

今回の解析解 = Riemann問題 / Gubser flow↓今回は時間離散化に着目するのでこちらを定量的に調べる★比較:数値計算結果 vs 小さな Δt の数値計算結果～ Δx finite, Δt → 0“Δt → 0 の収束の具合”
数値精度の指標
相対差

セル xj 毎に求まる→ max or 平均

計算コストの指標

nKT = (時間発展を通したセルあたりの時間微分評価回数の平均)= (時間発展右辺評価回数) / (セル数)2023/8/2912

解析解1: Riemann 問題

PartX (1/1)

Lora-Clavijo, et al (2013)= 一点で不連続な初期条件からの時間発展問題の総称 (衝撃波)
初期条件 (t = 0)

Δt 経過後?

方程式系、問題設定、近似などにより色々な種類あり

相対論的完全流体における解析解 (模式図) (状態方程式 p=e/3)(解析的に式で書けるが省略)
縦軸 エネルギー密度

横軸 x ~ = x/t※時間が立つと時刻に比例して横に広がっていく。
2023/8/29

13解析解2: Gubser flow S. S. Gubser, PRD 82, 085027 (2010)

PartX (1/1)

完全流体軸対称解: 座標 r = √(x2+y2), τ = √(t2–z2)
エネルギー密度動径方向の flow rapidity(u0 = cosh k)

エネルギー密度 vs r の時間発展

2023/8/29

流速 u0 (Lorentz factor) vs r の時間発展

14

結果

2023/8/29

15解析解1: Riemann 問題で比較

PartX (1/1)

数値精度&計算コスト←数値精度定性チェック(解析解と比較)

t = 15 fm

計算コスト →

衝撃波・不連続点近くで誤差大
2023/8/29

16解析解1: Riemann 問題で比較

PartX (1/1)

数値精度&計算コスト←数値精度定性チェック(解析解と比較)

t = 15 fm

計算コスト →

衝撃波・不連続点近くで誤差大
2023/8/29

17解析解1: Riemann 問題で比較

PartX (1/1)

数値精度&計算コスト←数値精度定性チェック(解析解と比較)

t = 15 fm
衝撃波・不連続点近くで誤差大
2023/8/29

計算コスト →

★陰解法の方がより速く正確に解くことができた(※) 18

解析解1: Riemann 問題で比較

PartX (1/1)

Fixed point 法の繰り返し回数空間分布①多くの場合、繰り返し回数 1 回で十分。
②局所的に硬い部分で繰り返し回数が多くなる。
2023/8/29

19解析解2: Gubser flowで比較

PartX (1/1)

時間※空間離散化 = KT/MUSCL スキームなので精度は良くない

Fixed point
繰り返し回数厳密解との相対差エネルギー密度x-y平面内2023/8/29→ ほとんどの所: fixed point 繰り返し数は 1 or 2. 危ない所で 3-4 回

20解析解2: Gubser flowで比較

PartX (1/1)

Riemann問題のときと同様の結果

2023/8/29

21

まとめ

2023/8/2922

まとめ

PartX (1/1)

背景: 相対論的流体力学 in 高エネルギー重イオン衝突• 安定・高効率・高精度な数値計算が重要
空間の離散化と時間の離散化• 時間陰解法 = 線形安定
手法: 陰解法(2次Gauss-Legendre法)と陽解法(Heun法)を相対論的完全流体に対して実装・比較• Fixed point 法 (初期guessに前ステップの値を利用)• 自動局所硬さ検出 (局所的にfixed point収束を判定)→陰解法が陽解法より高精度・高速
今後: 粘性流体・測定量の解析
2023/8/29

23

BACKUP

2023/8/2924

方程式系のまとめ(粘性)PartX (1/1)

保存則

エネルギー運動量保存
電荷保存
共変微分

状態方程式
構成方程式 (例)

2次 緩和

Navier-Stokes
物質微分

空間微分
2023/8/29

25

離散化と曲がった座標における保存則

2023/8/2926保存則・拘束条件

PartX (1/1)

経験上、既知の拘束条件・etc.を離散化誤差で破ってしまうとすぐに計算が破綻する→基本的に既知の拘束条件はすべて保つようにする
相対論的流体が満たすべき条件
12

保存則: 保存密度の空間積分が時間変化しない
エネルギー密度は正3固有値 (内部エネルギー+主応力3つ) が正

4

流体は時間的 (静止質量)5

散逸流は空間的2023/8/29

27

数値流体力学

PartX (1/1)

離散化
連続の場は直接計算機で取り扱えないので有限自由度で表す。
空間を小さなマス目に分割してマス目上で場の値を定義。

保存スキーム= 保存則に離散化誤差がないもの∂μTμν = 0, ∂μNμ = 0→ 有限体積法 (FVM) を使うことが多い※各セルの中の保存電荷の総量

動的変数 = 保存量の分布
時間発展 = 保存量の移動
2023/8/29

28

膨張座標系

PartX (1/1)

高エネルギー原子核衝突反応における縦膨張
初期の細かい構造が最後の大きな構造に反映される※原子核ローレンツ収縮
γ ～ 2700-3700
2023/8/2929

膨張座標系

PartX (1/1)

τ-η座標系 (R1,1Milne座標系⊕R2)
単なる膨張系ではなくこのように取る
物理的理由もあるが略2023/8/2930

一般座標上の保存則

PartX (1/1)

Minkowski座標での保存則 (flatな時空)

bar 添字 = Minkowski座標基底
bar なし添字 = 曲線座標の自然基底

divergence 型なのでFVMで厳密に保存させられる
任意の座標における保存則変な形!!

テンソル基底の変化分だけおつりの項が出てくる。
(縮約を取っている添字についてはヤコビアン √(-g) に吸収させられる)
2023/8/2931一般座標上の保存則

PartX (1/1)

テンソル密度 (一般の座標変換性を考慮した共変量)→ FVM OK!
→ ???
何がいけないのか?
2023/8/2932一般座標上の保存則

PartX (1/1)

保存するエネルギー運動量
そもそもの保存量はMinkowski座標での並進対称性に由来するエネルギー運動量

Minkowski座標系での保存量密度の各成分について、一般座標系での保存流を考える

OK?
☓ 原子核衝突の対称性 (現象のBoost invariance) を破る☓ 計算格子の端の方で γ ～ 3000 の分だけ保存量が無駄に大きくなり誤差が蓄積する
2023/8/2933

一般座標上の有限体積法

PartX (1/1)

作戦(1) この段階で離散化をしてから、(2) その後で一般座標のテンソル密度で書き表し直す

一般座標の上での保存するFVMの公式

y→x のベクトルの平行移動

時空が歪んでいる場合にも自然に適用可能
2023/8/2934

一般座標上の有限体積法

PartX (1/1)

一般座標の上での保存するFVMの公式③①②
y→x のベクトルの平行移動

時空が歪んでいる場合にも自然に適用可能★因みに他のところ(構成方程式)に現れる共変微分もすべて「平行移動」を使って実装している③t2023/8/29②x①35

散逸流の直交性

2023/8/2936

散逸流と流速の直交性

PartX (1/1)

拘束条件
微分形の方程式はもちろんこれを保持する。
しかし時間方向に離散化すると崩れてしまう。

もっと簡単な例: 相対論的粒子の運動

拘束条件
Euler 法で有限の時間幅だけ進むと確実に曲面の外に出てしまう(A) はみ出た後に射影する作戦。
射影 = 一番近い点に移す→ どうやって”距離”を定義する?
(B) 初めから曲面上のパラメータについての運動方程式を考える
2023/8/2937

散逸流の表現

PartX (1/1)

局所静止系で見ると局所静止系(ガリレイ系)のテンソル基底で表現すれば良い!
※局所静止系の基底をとると言っても局所静止系におけるSO(3)回転の自由度は残っている

2023/8/2938

局所静止系の基底の構成

PartX (1/1)①座標の自然基底を正規直交化②流速の方向にBoost (回転せずに一回のBoostで移る)

2023/8/2939

接続の計算PartX (1/1)

変なテンソル基底を取っているので共変微分が滅茶苦茶になる。
→接続係数を計算する。

2023/8/2940

構成方程式PartX (1/1)

構成方程式の本質的構造

回転速度 a に伴うコリオリ力的な力

2023/8/2941

高エネルギー重イオン衝突実験

PartX (1/1)

目的⚫ QGP などのQCD物質の平衡の性質⚫ 状態方程式、粘性・拡散係数などの輸送係数⚫ 臨界点、一次相転移、…
⚫ ジェット(高エネルギー粒子)⚫ エネルギー損失の仕方⚫ ハドロンの性質⚫ ハドロン間相互作用、物質中での性質の変化⚫ …
⚫ 衝突反応自体のダイナミクスの理解⚫ 初期状態⚫ ハドロン化
2023/8/29

42

PartX (1/1)
Xe+Xe衝突(5.44 TeV) の様子 (3D)

研究: 原子核の変形の効果と熱ゆらぎ

Deformed
Xe+Xe

T [MeV kB]

各設定4000事象の計算• 変形/球形 Xe+Xe• 流体ゆらぎあり/なし
SphericalXe+Xe

アニメーションは変形&流体ゆらぎあり

2023/8/29

43