diff --git a/process/physics.py b/process/physics.py index b09fff0af0..1de78bf5f3 100644 --- a/process/physics.py +++ b/process/physics.py @@ -84,11 +84,9 @@ def nui(j, zmain, ni, tempi, amain): @nb.jit(nopython=True, cache=True) -def _plascar_bpol(aspect, kappa, delta): +def _plascar_bpol(aspect, eps, kappa, delta): # Original coding, only suitable for TARTs [STAR Code] - eps = 1.0 / aspect - c1 = kappa**2 / (1.0 + delta) + delta c2 = kappa**2 / (1.0 - delta) - delta @@ -117,7 +115,7 @@ def _plascar_bpol(aspect, kappa, delta): @nb.jit(nopython=True, cache=True) -def bpol(icurr, ip, qbar, aspect, bt, kappa, delta, perim, rmu0): +def bpol(icurr, ip, qbar, aspect, eps, bt, kappa, delta, perim, rmu0): """Function to calculate poloidal field author: J Galambos, FEDC/ORNL author: P J Knight, CCFE, Culham Science Centre @@ -141,13 +139,13 @@ def bpol(icurr, ip, qbar, aspect, bt, kappa, delta, perim, rmu0): if icurr != 2: return rmu0 * ip / perim - ff1, ff2, _, _ = _plascar_bpol(aspect, kappa, delta) + ff1, ff2, _, _ = _plascar_bpol(aspect, eps, kappa, delta) return bt * (ff1 + ff2) / (2.0 * np.pi * qbar) @nb.jit(nopython=True, cache=True) -def plasc(qbar, aspect, rminor, bt, kappa, delta): +def plasc(qbar, aspect, eps, rminor, bt, kappa, delta): """Function to calculate plasma current (Peng scaling) author: J Galambos, FEDC/ORNL author: P J Knight, CCFE, Culham Science Centre @@ -167,7 +165,7 @@ def plasc(qbar, aspect, rminor, bt, kappa, delta): Fusion Technology, 21, 1729 """ - ff1, ff2, d1, d2 = _plascar_bpol(aspect, kappa, delta) + ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) e1 = 2.0 * kappa / (d1 * (1.0 + delta)) e2 = 2.0 * kappa / (d2 * (1.0 - delta)) @@ -2094,7 +2092,7 @@ def culcur( if icurr == 1: # Peng analytical fit fq = (1.22 - 0.68 * eps) / ((1.0 - eps * eps) ** 2) * sf**2 elif icurr == 2: # Peng scaling for double null divertor; TARTs [STAR Code] - curhat = 1.0e6 * plasc(qpsi, asp, rminor, bt, kappa, triang) / bt + curhat = 1.0e6 * plasc(qpsi, asp, eps, rminor, bt, kappa, triang) / bt elif icurr == 3: # Simple ITER scaling (simply the cylindrical case) fq = 1.0 elif icurr == 4: # ITER formula (IPDG89) @@ -2161,7 +2159,9 @@ def culcur( physics_variables.normalised_total_beta = ( 1.0e8 * physics_variables.beta * rminor * bt / plascur ) - bp = bpol(icurr, plascur, qpsi, asp, bt, kappa, triang, pperim, constants.rmu0) + bp = bpol( + icurr, plascur, qpsi, asp, eps, bt, kappa, triang, pperim, constants.rmu0 + ) # Ensure current profile consistency, if required # This is as described in Hartmann and Zohm only if icurr = 4 as well... diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index dfc8b50251..d569b594c3 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -642,6 +642,7 @@ def test_culcur(culcurparam, monkeypatch, physics): { "qbar": 2.5, "aspect": 2.7, + "eps": 0.37037037, "rminor": 1.5, "bt": 12, "kappa": 1.85, @@ -653,6 +654,7 @@ def test_culcur(culcurparam, monkeypatch, physics): { "qbar": 2.5, "aspect": 3.0, + "eps": 0.33333333, "rminor": 1.5, "bt": 12, "kappa": 1.85, @@ -675,6 +677,7 @@ def test_plasc(arguments, expected): "ip": 1.6e7, "qbar": 2.5, "aspect": 2.7, + "eps": 0.37037037, "bt": 12, "kappa": 1.85, "delta": 0.5, @@ -689,6 +692,7 @@ def test_plasc(arguments, expected): "ip": 1.6e7, "qbar": 2.5, "aspect": 3.0, + "eps": 0.33333333, "bt": 12, "kappa": 1.85, "delta": 0.5, @@ -703,6 +707,7 @@ def test_plasc(arguments, expected): "ip": 1.6e7, "qbar": 2.5, "aspect": 3.0, + "eps": 0.33333333, "bt": 12, "kappa": 1.85, "delta": 0.5,