From 39e09e10b92af4d99530761a8558d6c54a2aa1f4 Mon Sep 17 00:00:00 2001 From: DavidA Date: Sat, 15 Aug 2026 14:34:38 +0100 Subject: [PATCH 1/2] Fix inverted q95 to qbar transform in Peng TART scalings calculate_plasma_current_peng (i_plasma_current=2) and the TART branch of calculate_surface_averaged_poloidal_field derived qbar as q95 * 1.3 * (1 - eps)^0.6. The documented Peng/STAR relation is q95 = 1.3 * qbar * (1 - eps)^0.6, so recovering qbar from q95 requires division; before the pure-Python refactor (#3320) the input safety factor was interpreted as qbar directly and q95 was derived from it with exactly this forward relation, keeping the pair self-consistent. The multiply/divide inversion under-computed qbar by (1.3 * (1 - eps)^0.6)^2, overestimating plasma current by ~57% at A=1.8. With the fix, the Peng scaling at the ST regression geometry (q95=6, A=1.8, a=2.5 m, B_T=3 T, kappa=2.8) gives 22.5 MA, consistent with the 22.9 MA the independent i_plasma_current=9 scaling produces for the same machine; the inverted transform gave 35.2 MA. Adds the first unit tests for plasma_current.py: hand-derived reference values for the current and bpol, a round-trip check on the transform, and a consistency check that both call sites share the same qbar. Co-Authored-By: Claude Fable 5 --- process/models/physics/plasma_current.py | 7 +- process/models/physics/plasma_fields.py | 6 +- tests/unit/models/physics/test_physics.py | 16 ++- .../models/physics/test_plasma_current.py | 116 ++++++++++++++++++ 4 files changed, 137 insertions(+), 8 deletions(-) create mode 100644 tests/unit/models/physics/test_plasma_current.py diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py index cc1c12da37..a1972219a8 100644 --- a/process/models/physics/plasma_current.py +++ b/process/models/physics/plasma_current.py @@ -758,8 +758,11 @@ def calculate_plasma_current_peng( 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), 1729-1738. https://doi.org/10.13182/FST92-A29971 """ - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - (1.0 / aspect)) ** 0.6e0 + # Transform q95 to qbar by inverting q95 = 1.3 * qbar * (1 - eps)^0.6 + # (Peng, Galambos & Shipe 1992). Prior to #3320 the input safety factor + # was interpreted as qbar directly and q95 was derived from it with the + # forward relation above. + qbar = q95 / (1.3e0 * (1.0e0 - (1.0 / aspect)) ** 0.6e0) ff1, ff2, d1, d2 = self.plascar_bpol( aspect=aspect, eps=(1.0 / aspect), kappa=kappa, triang=triang diff --git a/process/models/physics/plasma_fields.py b/process/models/physics/plasma_fields.py index e178ebf2e3..9d6774cfd0 100644 --- a/process/models/physics/plasma_fields.py +++ b/process/models/physics/plasma_fields.py @@ -87,8 +87,10 @@ def calculate_surface_averaged_poloidal_field( aspect=aspect, eps=(1 / aspect), kappa=kappa, triang=triang ) - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - (1 / aspect)) ** 0.6e0 + # Transform q95 to qbar by inverting q95 = 1.3 * qbar * (1 - eps)^0.6 + # (Peng, Galambos & Shipe 1992), consistent with + # PlasmaCurrent.calculate_plasma_current_peng + qbar = q95 / (1.3e0 * (1.0e0 - (1 / aspect)) ** 0.6e0) return b_plasma_toroidal_on_axis * (ff1 + ff2) / (2.0 * np.pi * qbar) diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index 28deeeff83..601cd4ec6f 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -1116,6 +1116,10 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): @pytest.mark.parametrize( ("arguments", "expected"), [ + # Expected values updated when the q95 -> qbar transform was + # corrected from multiplication to division by 1.3*(1-eps)^0.6: + # relative to the old values this scales the current by + # (1.3*(1-eps)^0.6)^2 (-3.0% at A=2.7, +3.9% at A=3.0). ( { "q95": 2.5, @@ -1125,7 +1129,7 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): "kappa": 1.85, "triang": 0.5, }, - 38.00677211030666, + 36.86802557693355, ), ( { @@ -1136,7 +1140,7 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): "kappa": 1.85, "triang": 0.5, }, - 31.290463480745593, + 32.50792056059886, ), ], ) @@ -1149,6 +1153,10 @@ def test_calculate_plasma_current_peng(arguments, expected, physics): @pytest.mark.parametrize( ("arguments", "expected"), [ + # Expected values for the i_plasma_current=2 (Peng TART) branch + # updated when the q95 -> qbar transform was corrected from + # multiplication to division by 1.3*(1-eps)^0.6 (scales bpol by + # (1.3*(1-eps)^0.6)^2). ( { "i_plasma_current": 2, @@ -1160,7 +1168,7 @@ def test_calculate_plasma_current_peng(arguments, expected, physics): "triang": 0.5, "len_plasma_poloidal": 24, }, - 3.5258772213675047, + 3.4202360358630237, ), ( { @@ -1173,7 +1181,7 @@ def test_calculate_plasma_current_peng(arguments, expected, physics): "triang": 0.5, "len_plasma_poloidal": 24, }, - 2.902807218476584, + 3.0157503585409073, ), ( { diff --git a/tests/unit/models/physics/test_plasma_current.py b/tests/unit/models/physics/test_plasma_current.py new file mode 100644 index 0000000000..681326fb7b --- /dev/null +++ b/tests/unit/models/physics/test_plasma_current.py @@ -0,0 +1,116 @@ +"""Unit tests for the Peng TART plasma current scaling (i_plasma_current=2).""" + +import numpy as np +import pytest + +from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel +from process.models.physics.plasma_fields import PlasmaFields + + +@pytest.fixture +def plasma_current(): + return PlasmaCurrent() + + +@pytest.fixture +def plasma_fields(): + return PlasmaFields() + + +def test_calculate_plasma_current_peng_qbar_transform(plasma_current): + """The q95 -> qbar transform must invert q95 = 1.3 * qbar * (1-eps)^0.6. + + Expected value derived by hand for q95=6, A=1.8, a=2.5 m, B_T=3 T, + kappa=2.8, triang=0.5: + qbar = q95 / (1.3 * (1 - 1/A)^0.6) = 6 / 0.79916... = 7.507881... + and the remaining geometry factors evaluated from plascar_bpol as in + the function body. Note qbar > q95, as required by the forward + relation q95 = 1.3 * qbar * (1 - eps)^0.6 < qbar for any eps > 0. + """ + current = plasma_current.calculate_plasma_current_peng( + q95=6.0, + aspect=1.8, + rminor=2.5, + b_plasma_toroidal_on_axis=3.0, + kappa=2.8, + triang=0.5, + ) + + assert current == pytest.approx(22.45707000674425) + + +def test_peng_qbar_roundtrip(): + """Applying the documented forward relation to the derived qbar returns q95.""" + q95 = 6.0 + aspect = 1.8 + factor = 1.3e0 * (1.0e0 - (1.0 / aspect)) ** 0.6e0 + qbar = q95 / factor + assert qbar * factor == pytest.approx(q95) + # qbar exceeds q95 for the Peng relation at any finite aspect ratio + assert qbar > q95 + + +def test_calculate_surface_averaged_poloidal_field_peng(plasma_fields): + """TART branch of : expected value derived by hand. + + bpol = B_T * (ff1 + ff2) / (2 * pi * qbar) with + qbar = q95 / (1.3 * (1 - 1/A)^0.6) and ff1/ff2 from plascar_bpol at + A=1.8, kappa=2.8, triang=0.5. + """ + bpol = plasma_fields.calculate_surface_averaged_poloidal_field( + i_plasma_current=PlasmaCurrentModel.PENG_DIVERTOR_SCALING, + cur_plasma=0.0, # unused by the Peng branch + q95=6.0, + aspect=1.8, + b_plasma_toroidal_on_axis=3.0, + kappa=2.8, + triang=0.5, + len_plasma_poloidal=1.0, # unused by the Peng branch + ) + + assert bpol == pytest.approx(0.915149073471596) + + +def test_peng_current_and_bpol_share_qbar(plasma_current, plasma_fields): + """The two Peng-scaling call sites must use the same qbar transform. + + I = 5 * kappa * rminor * B_T / (2 * pi^2 * qbar) * (asin(e1)/e1 + + asin(e2)/e2) * (ff1 + ff2) and = B_T (ff1+ff2) / (2 pi qbar) + imply I / is independent of qbar; verify both functions agree on + the qbar-dependent part instead by reconstructing qbar from each. + """ + q95, aspect, rminor, bt, kappa, triang = 4.0, 2.0, 1.5, 2.5, 2.2, 0.4 + + current = plasma_current.calculate_plasma_current_peng( + q95=q95, + aspect=aspect, + rminor=rminor, + b_plasma_toroidal_on_axis=bt, + kappa=kappa, + triang=triang, + ) + bpol = plasma_fields.calculate_surface_averaged_poloidal_field( + i_plasma_current=PlasmaCurrentModel.PENG_DIVERTOR_SCALING, + cur_plasma=0.0, + q95=q95, + aspect=aspect, + b_plasma_toroidal_on_axis=bt, + kappa=kappa, + triang=triang, + len_plasma_poloidal=1.0, + ) + + ff1, ff2, d1, d2 = plasma_current.plascar_bpol( + aspect=aspect, eps=(1.0 / aspect), kappa=kappa, triang=triang + ) + e1 = (2.0 * kappa) / (d1 * (1.0 + triang)) + e2 = (2.0 * kappa) / (d2 * (1.0 - triang)) + shape = (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) * (ff1 + ff2) + + qbar_from_current = rminor * bt * 5.0 * kappa * shape / (2.0 * np.pi**2 * current) + qbar_from_bpol = bt * (ff1 + ff2) / (2.0 * np.pi * bpol) + + assert qbar_from_current == pytest.approx(qbar_from_bpol) + assert qbar_from_current == pytest.approx( + q95 / (1.3e0 * (1.0e0 - (1.0 / aspect)) ** 0.6e0) + ) From 3b515493f4ff4502f5f7eb4f1cf76f6f279e74da Mon Sep 17 00:00:00 2001 From: DavidA Date: Tue, 18 Aug 2026 10:07:02 +0100 Subject: [PATCH 2/2] Cite Muldrew (2020) Eq. 19 for the q95-qbar relation; fix docs and test rationale - Attribute the relation q95 = 1.3*qbar*(1-eps)^0.6 to Muldrew et al., Fusion Eng. Des. 154 (2020) 111530, Eq. (19), in both call-site comments (previously attributed to Peng, Galambos & Shipe 1992, which was not verified). - Correct the STAR/Peng section of plasma_current.md, which stated the inverse (qbar = q95 * 1.3(1-eps)^0.6, describing the post-#3320 code), and cite [^10]. Fix the mis-parenthesised formula in the interactive figure's callback and in documentation/scripts/plotting_scripts/ profile_peng_qbar.py ((1 - eps^0.6) -> (1-eps)^0.6, and the direction); the figure's baked-in initial data array should be regenerated from the script at the next docs build. - Remove the physical-ordering rationale ('qbar > q95 for any eps > 0') from the test docstring and the corresponding assert: the relation gives qbar > q95 only for A < 2.82, and the ordering is a property of the fit, not a first-principles requirement. The fix rests on Muldrew Eq. (19) and the pre-#3320 implementation. - Remove the stale '(= q-bar for i_plasma_current=2)' from the q95 parameter docstring in calculate_plasma_current. Co-Authored-By: Claude Fable 5 --- .../scripts/plotting_scripts/profile_peng_qbar.py | 4 ++-- .../physics-models/plasma_current/plasma_current.md | 8 +++++--- process/models/physics/plasma_current.py | 8 ++++---- process/models/physics/plasma_fields.py | 4 ++-- tests/unit/models/physics/test_plasma_current.py | 5 +---- 5 files changed, 14 insertions(+), 15 deletions(-) diff --git a/documentation/scripts/plotting_scripts/profile_peng_qbar.py b/documentation/scripts/plotting_scripts/profile_peng_qbar.py index 1e5f24e818..4c769050a7 100644 --- a/documentation/scripts/plotting_scripts/profile_peng_qbar.py +++ b/documentation/scripts/plotting_scripts/profile_peng_qbar.py @@ -6,7 +6,7 @@ from bokeh.plotting import figure, output_file, save x = np.linspace(1.0, 3.0, 500) -y1 = 5.0 * 1.3 * (1.0 - (1.0 / x) ** 0.6) +y1 = 5.0 / (1.3 * (1.0 - 1.0 / x) ** 0.6) y2 = 5.0 * (1.0 + 2.6 * (1.0 / x) ** 2.8) # Initial data for the second line source = ColumnDataSource(data={"x": x, "y1": y1, "y2": y2}) @@ -50,7 +50,7 @@ const B = qbar.value; const x = source.data['x']; - const y1 = x.map(xi => A * 1.3 * (1.0 - (1.0 / xi) ** 0.6)); + const y1 = x.map(xi => A / (1.3 * (1.0 - 1.0 / xi) ** 0.6)); const y2 = x.map(xi => B * (1.0 + 2.6 * (1.0 / xi) ** 2.8)); // Example transformation for the second line source.data['y1'] = y1; diff --git a/documentation/source/physics-models/plasma_current/plasma_current.md b/documentation/source/physics-models/plasma_current/plasma_current.md index f7d31f95c9..ab384ac294 100644 --- a/documentation/source/physics-models/plasma_current/plasma_current.md +++ b/documentation/source/physics-models/plasma_current/plasma_current.md @@ -154,10 +154,12 @@ $$ $$ The STAR code document states that the $\bar{q}_0$ is normally equal to 3.0 -In PROCESS this is implemented differently as the function: +In PROCESS this is implemented differently: $\bar{q}$ is obtained from the available $q_{95}$ parameter using the relation given in [^10] (Eq. 19 therein), $$ -\bar{q} = q_{95} \times 1.3(1-\epsilon)^{0.6} +q_{95} = 1.3\,\bar{q}\,(1-\epsilon)^{0.6} +\quad\Rightarrow\quad +\bar{q} = \frac{q_{95}}{1.3(1-\epsilon)^{0.6}} $$ This is to allow the use of the available $q_{95}$ parameter as **the origin and definition of $\bar{q}_0$ is not fully known.** @@ -205,7 +207,7 @@ The coefficient values $q_{95}$ and $\bar{q}_0$ for the PROCESS and STAR code im