From a251e2effac6106b40b94facdf60d21d389a1eca Mon Sep 17 00:00:00 2001 From: keyhell Date: Mon, 5 Oct 2026 09:34:47 +0200 Subject: [PATCH] Normalise the parabolic ohmic current density to the plasma current Without a gEQDSK, _setup_profs normalised the current-density shape to 1 on axis while _P_OH_prof used J = Itot*j_prof, so J peaked at Ip amperes per square metre whatever the plasma cross-section and int J dA_pol != Ip. Pohmic and vloop were off by (int j_prof dA_pol)^-2: x0.25 on the SPARC example, x29 on ITER, x0.009 on NSF_NT_CMOD. The gEQDSK branch already divides by the cross-section integral (Jrms_norm = jrms / Jrmsint). Divide the parabolic shape by its poloidal cross-section integral, dA_pol = dV / (2 pi R), so int Itot*j_prof dA_pol = Itot. Add a regression test on the SPARC settings. MANTA golden values unchanged; 160 tests pass. Co-Authored-By: Claude Fable 5.1 --- openpopcon/core.py | 7 +++++++ tests/test_invariants.py | 18 ++++++++++++++++++ 2 files changed, 25 insertions(+) diff --git a/openpopcon/core.py b/openpopcon/core.py index 448487d..5da2639 100644 --- a/openpopcon/core.py +++ b/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) * ( diff --git a/tests/test_invariants.py b/tests/test_invariants.py index c11eb7b..25a29a9 100644 --- a/tests/test_invariants.py +++ b/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)