Normalise the parabolic ohmic current density to the plasma current

openpopcon/core.py

@@ -1530,6 +1530,13 @@ def _setup_profs(self) -> None:
             self.j_prof = (1 - self.j_offset) * (
                 1 - rho**self.j_alpha1
             ) ** self.j_alpha2 + self.j_offset
+            # Scale so that the poloidal cross-section integral of Itot*j_prof
+            # is Itot (dA_pol = dV / (2 pi R)), as the gEQDSK branch already
+            # does with Jrms_norm = jrms / int jrms dA. Without this, J peaks
+            # at Ip amperes per square metre whatever the plasma size and
+            # P_OH = int eta J^2 dV is low by (int j_prof dA_pol)^-2.
+            I_shape = np.trapezoid(self.j_prof, self.volgrid) / (2 * np.pi * self.R)
+            self.j_prof = self.j_prof / I_shape
             self._jdefined = True
         if not self._nedefined:
             self.ne_prof = (1 - self.ne_offset) * (

tests/test_invariants.py

@@ -67,3 +67,21 @@ def test_output_roundtrip(solved_manta, tmp_path):
             err_msg=f"roundtrip changed {field}",
         )
     plt.close("all")
+
+
+def test_parabolic_current_profile_integrates_to_ip():
+    """
+    Without a gEQDSK the current density used for the ohmic power is
+    Itot * j_prof. Its poloidal cross-section integral must be the plasma
+    current, as it already is on the gEQDSK branch (Jrms_norm).
+    """
+    from helpers import SPARC_SETTINGS, SPARC_PLOTSETTINGS
+
+    pc = op.POPCON(settingsfile=SPARC_SETTINGS, plotsettingsfile=SPARC_PLOTSETTINGS)
+    pc.settings.verbosity = 0
+    pc.run_POPCON(setuponly=True)
+    al = pc.algorithms
+    rho = al.sqrtpsin
+    J = al.state.Itot * 1e6 * al.get_profile(rho, 0)  # A/m^2
+    ip_from_j = al.volume_integral(rho, J) / (2 * np.pi * al.state.R) / 1e6  # MA
+    np.testing.assert_allclose(ip_from_j, al.state.Ip, rtol=1e-3)