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

Đang mở
  • Có gì không ổn?
#497 0 bình luận 0 reaction 0 người được giao Xem trên GitHub

Maintainer thường phản hồi trong vòng 3 ngày

Chưa có ai nhận issue này.

Đánh giá

Độ khó
4/5
Thời gian dự kiến
3-5 ngày
Mức phù hợp với người mới
48/100
Loại issue
Lỗi
Độ rõ ràng
Khá rõ ràng
Mức độ hoạt động
Sôi nổi
Công nghệ
julia
Lĩnh vực
hpc

Hướng nghiên cứu

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.

Do mô hình lập chỉ mục viết ra từ nội dung của issue.

Mô tả

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 !.

Ngôn ngữ chính
Julia
Star
11
Fork
7
Merge trung bình
3 ngày 22 giờ
Pull request đã merge (30 ngày)
70

Chuẩn bị môi trường

  • Không có Dockerfile hay tệp Docker Compose
  • Có mẫu pull request
  • Không có hướng dẫn đóng góp

Bắt đầu từ đâu

  1. Đọc hết issue, rồi đọc hướng dẫn đóng góp của dự án.
  2. Bình luận trên issue rằng bạn sẽ nhận — tránh hai người làm cùng một việc.
  3. Fork repository và làm thay đổi trên một nhánh.
  4. Mở pull request có tham chiếu số hiệu của issue.

Issue khác của OpenFUSIONToolkit/GPEC

Tất cả issue của OpenFUSIONToolkit/GPEC

Issue tương tự

Thêm issue về Julia

Nhận issue mới trong hộp thư của bạn

Bản tóm tắt ngắn những issue GitHub phù hợp với người mới.