Utilities - BUGFIX! - Compute the trapped-particle fraction from the pitch-angle integral

オープン
  • 何か問題がありますか？
#497 コメント 0 件 リアクション 0 件 担当者 0 名 GitHub で見る

メンテナーはふだん 3 日以内に返信

まだ誰も着手していません。

評価

難易度
4/5
見積もり時間
3〜5日
初心者へのやさしさ
48/100
issue の種類
バグ
明瞭さ
おおむね明確
活発さ
活発
技術スタック
julia
領域
hpc

調査の方向性

Read Utilities.NeoclassicalResistivity.trapped_fraction and trapped_fraction_eps, then follow resist_geometry and its θ loop to understand the available field values and flux-surface weights. Implement the stated pitch-angle integral and check the proposed circular-field comparisons, shaped-field bounds, and ε → 0 behavior; also correct the p_local docstring if it is in scope.

索引モデルが issue の本文から書いたものです。

説明

Summary

Utilities.NeoclassicalResistivity.trapped_fraction(avg_B, avg_Bsq, B_min, B_max) does not return the trapped-particle
fraction it documents. It overestimates f_t by 11–44 % (worst at small inverse aspect ratio), which makes the Sauter and
Redl neoclassical resistivity too large. Every consumer of SauterNeoModel / RedlNeoModel is affected: build_ggj_inputs
and build_slayer_inputs (through ResistGeometry.f_trap), so τ_R, S and tearing growth rates shift.

What the code does

h   = B_min / B_max
f_t = 1 − (⟨B⟩² / ⟨B²⟩) · (1 − √(1 − h)·(1 + h/2))

The docstring attributes this to Lin-Liu & Miller, Phys. Plasmas 2, 1666 (1995). Their upper bound uses
h = ⟨B⟩/B_max with prefactor ⟨B²⟩/⟨B⟩² (not B_min/B_max and ⟨B⟩²/⟨B²⟩), and their recommended estimate is
0.75·f_tu + 0.25·f_tl with a separate lower bound f_tl. Using B_min/B_max makes 1 − h ≈ 2ε instead of ≈ ε, so the
large-aspect-ratio limit is f_t ≈ 2.1√ε instead of 1.46√ε.

Reference: the definition

f_t = 1 − (3/4)·⟨B²⟩ ∫₀^{1/B_max} λ dλ / ⟨√(1 − λB)⟩ (Hirshman & Sigmar 1981; the f_t that the Sauter 1999 F₃₃ fit uses)

For a circular large-aspect-ratio field B ∝ 1/(1 + ε cos θ), flux-surface averages weighted by dl/B:

ε exact Lin-Liu–Miller (0.75 f_tu + 0.25 f_tl) trapped_fraction trapped_fraction_eps code / exact
0.01 0.1460 0.1459 0.2097 0.1462 1.44
0.05 0.3228 0.3226 0.4489 0.3239 1.39
0.10 0.4492 0.4488 0.6028 0.4506 1.34
0.30 0.7230 0.7224 0.8686 0.7238 1.20
0.50 0.8645 0.8641 0.9589 0.8646 1.11

The ε-only fit already in the module (trapped_fraction_eps) agrees with the exact integral to 0.3 %, so the two
functions in the same file disagree with each other by up to 44 %.

Impact

On a circular ε = 0.2 equilibrium with rational surfaces at local ε ≈ 0.13 (q = 2) and 0.18 (q = 3):
f_t = 0.659 / 0.751 from the code against 0.487 / 0.563 exact, and Sauter η is 19 % / 9.5 % too high. Resistive
growth rates scale roughly as η^(3/5), so a 5–10 % shift in γ. Comparisons between layer models inside the code are
unaffected (they share η); comparisons with other codes and absolute numbers are.

Proposed fix

Evaluate the definition directly in resist_geometry, which already has B(θ) and the jac/V′ weights in its θ loop:
Gauss–Legendre in λ on [0, 1/B_max] with the existing periodic flux-surface average for ⟨√(1 − λB)⟩. Convergence of
f_t on 259 θ points: 2.9e-5 (8 nodes), 1.8e-6 (16), 1.2e-7 (32). The exact integral is preferred over the Lin-Liu–Miller
estimate: it costs about the same, and the 0.75/0.25 weighting is off by ~7 % on a narrow deep magnetic well
(bounds 0.11 / 0.26 around an exact 0.205), although it is within 0.1 % for ordinary tokamak shaping.

Also: the ResistGeometry docstring lists p_local as pressure in Pa; it holds μ₀P.

Tests to add

  • circular large-aspect-ratio field: f_t within 0.5 % of trapped_fraction_eps for ε ∈ {0.01, 0.05, 0.1, 0.3, 0.5}
  • the Lin-Liu–Miller bounds bracket the result on a shaped field (f_tl ≤ f_t ≤ f_tu)
  • f_t → 0 as ε → 0

Results move for any run using Sauter or Redl resistivity, hence the !.

主要言語
Julia
スター
11
フォーク
7
平均マージ
4日 1時間
マージ済み PR（30日）
68

環境構築

  • Dockerfile・Docker Compose ファイルなし
  • プルリクエストのテンプレートあり
  • コントリビューションガイドなし

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

OpenFUSIONToolkit/GPEC のほかの issue

OpenFUSIONToolkit/GPEC の issue をすべて見る

似ている issue

Julia の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。