Fix inverted q95 to qbar transform in Peng TART scalings

[dallonby]

Closes #4532

Overview

The q95 → q̄ transform in the two i_plasma_current=2 (Peng TART/STAR) call sites multiplied by 1.3*(1-eps)^0.6 where inverting the relation q95 = 1.3*q̄*(1-eps)^0.6 requires division.

Source for the direction of the relation: Muldrew et al., "PROCESS": Systems studies of spherical tokamaks, Fusion Eng. Des. 154 (2020) 111530, Eq. (19): "q̄ is the 'edge' safety factor and is related to q95 through q₉₅ = 1.3 q̄ (1.0 − ε)^0.6", with the Peng current in Eq. (18) written in terms of q̄. The pre-#3320 code implemented exactly this direction: the input safety factor was taken as q̄ and q95 was derived as q * 1.3*(1-eps)**0.6 (v2.5.0 process/physics.py:52-55; Fortran physics_variables.f90 documents q as "mean edge safety factor qbar" for icurr=2). #3320 (Oct 2024) removed that branch and introduced qbar = q95 * 1.3*(1-eps)**0.6 at both call sites, and rewrote the docs section to match — so the current docs page describes the inverted code rather than Muldrew Eq. (19). Detail in the issue.

Changes

• process/models/physics/plasma_current.py (calculate_plasma_current_peng) and process/models/physics/plasma_fields.py (calculate_surface_averaged_poloidal_field, TART branch): divide by 1.3e0 * (1 - eps)**0.6e0; comments cite Muldrew Eq. (19). Also removes the stale "(= q-bar for i_plasma_current=2)" from the q95 parameter docstring (a vestige of the two-variable scheme; wrong under either form).
• documentation/source/physics-models/plasma_current/plasma_current.md (STAR/Peng section): the printed relation corrected to Muldrew Eq. (19) with the existing [^10] citation; the embedded interactive figure's callback corrected. The same formula in documentation/scripts/plotting_scripts/profile_peng_qbar.py was also mis-parenthesised ((1 − ε^0.6) for (1−ε)^0.6) — fixed; the figure's baked-in initial data array should be regenerated from the script when the docs are next built.
• New tests/unit/models/physics/test_plasma_current.py (first dedicated test file for this module — Peng/bpol goldens already existed in test_physics.py): reference values for the corrected current and bpol at ST parameters, an algebraic round-trip on the transform, and a consistency check that both call sites share the same q̄.
• tests/unit/models/physics/test_physics.py: the four existing Peng/bpol expected values at A=2.7 and A=3.0 are updated. Justification per CONTRIBUTING: they were generated from the inverted code and pinned it; fixed/old ratios are 0.970038 (−3.00%) and 1.038908 (+3.89%) respectively, i.e. (1.3*(1-eps)^0.6)^2.

Behavioural impact

Selected-model outputs change only for i_plasma_current=2 runs; none of the eight shipped regression inputs selects it. At fixed inputs the corrected Peng current and bpol scale by (1.3*(1-eps)^0.6)^2 relative to the old code: −62% at A=1.4, −36% at A=1.8, −3.0% at A=2.7, +3.9% at A=3.0. A re-optimised design need not show that drop directly — on the shipped ST case switched to icurr=2 the optimiser holds Ip nearly constant and moves q95 instead (Ip 23.2 → 22.9 MA; q95 9.08 → 5.89; P_net 132 → 102 MW), i.e. the old code was reporting ~55% more edge-q margin than the machine had.

Cross-check at fixed inputs (shipped ST geometry, q95 = 6.41): corrected Peng 21.0 MA vs the independent icurr=9 FIESTA scaling 22.9 MA (−8%), IPDG89 23.0 MA, Peng-analytic 22.1 MA; the old code gave 32.9 MA (+43% vs FIESTA). The −8% is consistent with Muldrew's own Peng-vs-FIESTA comparison (Fig. 2, up to 20% under at high κ).

Regression CI: the always-on diagnostic c_plasma_peng_double_null (Peng current reported for every tokamak run regardless of the selected model) changes by the fix's factor on the five tokamak inputs (e.g. st_regression 3.291e7 → 2.102e7 A; +3.9% on the A=3 machines). In the CI logs each failing regression comparison shows exactly that one variable; on the local st_regression run it is the only non-metadata MFILE difference. Those references need re-baselining on merge. improve_quality_check also fails on this branch: calculate_plasma_current C901 complexity 13 > 12 — pre-existing on main, the function is untouched.

Verification

• Full unit suite passes locally (Python 3.12.7): 850 collected, 846 passed / 4 skipped; CI unit-test (ubuntu, macos) green.
• Out of scope but noted for a follow-up: constraint 45 compares q95 against q95_min computed from Muldrew's Eq. (17) q̄ᵐⁱⁿ = 3(1+2.6ε^2.8) — a q̄ quantity — since :3263 update the plasma current section of the docs to show and explain all models #3320 merged q into q95; and the icurr=2 OUT.DAT label "Mean edge safety factor (q95)" (plasma_current.py:158) is another vestige of the same scheme.

Correction to an earlier version of this description

An earlier version of this text, and the test docstring, justified the direction with "q̄ > q95, which the forward relation requires for any ε > 0". That was wrong as stated — the relation gives q̄ > q95 only for A < 2.82 (1.3(1−ε)^0.6 < 1), and q̄ < q95 at A = 3 — and was in any case a property of the fitted relation, not physics. It has been removed from the docstring in this branch's second commit. The fix rests on Muldrew Eq. (19) and the pre-#3320 implementation.

Found during an independent audit of v3.4.2.

[dallonby]

Correction to the 'Behavioural impact' section above: the regression MFILEs are not fully unchanged — the diagnostic comparison entry c_plasma_peng_double_null (the Peng-scaling current reported alongside every run's chosen scaling) changes by the fix's factor even on inputs that use other i_plasma_current models. Verified on st_regression: every design variable is identical to the last digit; only this one diagnostic row moves (3.291e7 → 2.102e7 A, −36%). So the 0.2%/5% regression jobs will flag exactly that one variable per tokamak input. All actual design outputs are unaffected.

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 49.32%. Comparing base (c0ae5b2) to head (39e09e1).

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #4533   +/-   ##
=======================================
  Coverage   49.32%   49.32%           
=======================================
  Files         150      150           
  Lines       29802    29802           
=======================================
  Hits        14700    14700           
  Misses      15102    15102           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mkovari]

@chris-ashe @ajpearcey
This is an ST issue.