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

Open
  • Is something wrong?
#497 0 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 3 days

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
48/100
Issue type
Bug
Clarity
Mostly clear
Activity status
Active
Tech stack
julia
Domain
hpc

Research direction

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.

Written by the indexing model from the issue text.

Description

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

Dominant language
Julia
Stars
11
Forks
7
Avg merge
4d 3h
Merged PRs (30d)
70

Getting set up

  • No Dockerfile or Docker Compose file
  • Has a pull request template
  • No contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OpenFUSIONToolkit/GPEC

All issues in OpenFUSIONToolkit/GPEC

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.