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

Aperta
  • Qualcosa non va?
#497 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 3 giorni

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
4/5
Tempo stimato
3-5 giorni
Idoneità per principianti
48/100
Tipo di issue
Bug
Chiarezza
Abbastanza chiara
Stato di attività
Attiva
Stack tecnologico
julia
Ambito
hpc

Direzione di ricerca

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.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

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

Lingua principale
Julia
Stelle
11
Fork
7
Merge medio
4g 3h
PR unite (30g)
70

Preparare l'ambiente

  • Nessun Dockerfile né file Docker Compose
  • Ha un modello di pull request
  • Nessuna guida per i contributori

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di OpenFUSIONToolkit/GPEC

Tutte le issue di OpenFUSIONToolkit/GPEC

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.