Fix inverted q95 to qbar transform in Peng TART scalings

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;

process/models/physics/plasma_current.py

@@ -245,7 +245,7 @@ def calculate_plasma_current(
         len_plasma_poloidal : float
             Plasma perimeter length (m).
         q95 : float
-            Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2).
+            Plasma safety factor at 95% flux.
         rmajor : float
             Major radius (m).
         rminor : float
@@ -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
+        # (Muldrew et al., Fusion Eng. Des. 154 (2020) 111530, Eq. 19). 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

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
+        # (Muldrew et al., Fusion Eng. Des. 154 (2020) 111530, Eq. 19),
+        # 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)
 

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,
         ),
         (
             {

tests/unit/models/physics/test_plasma_current.py

@@ -0,0 +1,113 @@
+"""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.
+    """
+    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)
+
+
+def test_calculate_surface_averaged_poloidal_field_peng(plasma_fields):
+    """TART branch of <Bp(a)>: 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 <Bp> = B_T (ff1+ff2) / (2 pi qbar)
+    imply I / <Bp> 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)
+    )