Describe the bug
PlasmaCurrent.calculate_plasma_current_peng and the TART branch of PlasmaFields.calculate_surface_averaged_poloidal_field derive qbar as:
# Transform q95 to qbar
qbar = q95 * 1.3e0 * (1.0e0 - (1.0 / aspect)) ** 0.6e0
The Peng/STAR relation between the two quantities is q95 = 1.3 * qbar * (1 - eps)^0.6, so recovering qbar from q95 requires division by 1.3*(1-eps)^0.6, not multiplication. Before the pure-Python refactor of #3320, the input safety factor for icurr=2 was interpreted as qbar directly and q95 was derived from it with exactly the forward relation above (v2.5.0 process/physics.py):
if physics_variables.icurr == 2:
physics_variables.q95 = (
physics_variables.q * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0
)so the pre-refactor pair was self-consistent. The current code applies the forward relation in the direction its own comment calls the inverse.
Impact
qbar is under-computed by (1.3*(1-eps)^0.6)^2, and current/bpol are overestimated by the same factor. At conventional aspect ratio the factor is close to 1 (−3% at A=2.7, +4% at A=3.0 — which is why the existing golden tests, both at conventional A, could not see it), but at the ST aspect ratios this scaling targets it is large: +57% in plasma current at A=1.8.
Cross-check: at the st_regression geometry (q95=6, A=1.8, a=2.5 m, B_T=3 T, kappa=2.8), the corrected Peng scaling gives 22.5 MA, consistent with the 22.9 MA that the independent i_plasma_current=9 FIESTA scaling produces for the same machine; the current code gives 35.2 MA.
No shipped regression input uses i_plasma_current=2 and plasma_current.py had no unit tests, so nothing caught the inversion.
Proposed fix
Divide instead of multiply at both call sites, update the two conventional-A golden values accordingly, and add unit tests for the Peng scaling (hand-derived reference values + a consistency check between the two call sites). PR ready (found during an independent audit of v3.4.2).
Describe the bug
PlasmaCurrent.calculate_plasma_current_pengand the TART branch ofPlasmaFields.calculate_surface_averaged_poloidal_fieldderive qbar as:The Peng/STAR relation between the two quantities is
q95 = 1.3 * qbar * (1 - eps)^0.6, so recovering qbar from q95 requires division by1.3*(1-eps)^0.6, not multiplication. Before the pure-Python refactor of #3320, the input safety factor foricurr=2was interpreted as qbar directly and q95 was derived from it with exactly the forward relation above (v2.5.0process/physics.py):so the pre-refactor pair was self-consistent. The current code applies the forward relation in the direction its own comment calls the inverse.
Impact
qbar is under-computed by
(1.3*(1-eps)^0.6)^2, and current/bpol are overestimated by the same factor. At conventional aspect ratio the factor is close to 1 (−3% at A=2.7, +4% at A=3.0 — which is why the existing golden tests, both at conventional A, could not see it), but at the ST aspect ratios this scaling targets it is large: +57% in plasma current at A=1.8.Cross-check: at the
st_regressiongeometry (q95=6, A=1.8, a=2.5 m, B_T=3 T, kappa=2.8), the corrected Peng scaling gives 22.5 MA, consistent with the 22.9 MA that the independenti_plasma_current=9FIESTA scaling produces for the same machine; the current code gives 35.2 MA.No shipped regression input uses
i_plasma_current=2andplasma_current.pyhad no unit tests, so nothing caught the inversion.Proposed fix
Divide instead of multiply at both call sites, update the two conventional-A golden values accordingly, and add unit tests for the Peng scaling (hand-derived reference values + a consistency check between the two call sites). PR ready (found during an independent audit of v3.4.2).