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

未关闭
  • 哪里不对？
#497 0 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 3 天内回复

还没有人认领这个 Issue。

评估

难度
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
平均合并
3 天 14 小时
30 天内合并 PR
69

环境准备

  • 没有 Dockerfile 或 Docker Compose 文件
  • 有 Pull Request 模板
  • 没有贡献指南

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

OpenFUSIONToolkit/GPEC 的其他 Issue

查看 OpenFUSIONToolkit/GPEC 的全部 Issue

相似的 Issue

更多 Julia Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。