#392 Sauter's formulae for plasma current and geometric parameters.

process/physics.py

@@ -2015,7 +2015,7 @@ def culcur(
         p0,
         pperim,
         q0,
-        qpsi,
+        q95,
         rli,
         rmajor,
         rminor,
@@ -2037,6 +2037,8 @@ def culcur(
         5 = Todd empirical scaling I
         6 = Todd empirical scaling II
         7 = Connor-Hastie model
+        8 = Sauter scaling (allowing negative triangularity) Issue #392
+            'Geometric formulas for system codes including the effect of negative triangularity'
         iprofile : input integer : switch for current profile consistency
         0 = use input alphaj, rli
         1 = make these consistent with q, q0
@@ -2045,7 +2047,7 @@ def culcur(
         p0       : input real :  central plasma pressure (Pa)
         pperim   : input real :  plasma perimeter length (m)
         q0       : input real :  plasma safety factor on axis
-        qpsi     : input real :  plasma edge safety factor (= q-bar for icurr=2)
+        q95      : input real :  plasma safety factor at 95% flux (= q-bar for icurr=2)
         rli      : input/output real : plasma normalised internal inductance
         rmajor   : input real :  major radius (m)
         rminor   : input real :  minor radius (m)
@@ -2055,11 +2057,9 @@ def culcur(
         bp       : output real : poloidal field (T)
         qstar    : output real : equivalent cylindrical safety factor (shaped)
         plascur  : output real : plasma current (A)
-        This routine performs the calculation of the
-        plasma current, with a choice of formula for the edge
-        safety factor. It will also make the current profile parameters
+        This routine calculates the plasma current based on the edge
+        safety factor q95. It will also make the current profile parameters
         consistent with the q-profile if required.
-        AEA FUS 251: A User's Guide to the PROCESS Systems Code
         J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code,
         unpublished internal Oak Ridge document
         Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992,
@@ -2082,7 +2082,7 @@ def culcur(
         # Calculate the function Fq that scales the edge q from the
         # circular cross-section cylindrical case
 
-        # First check for negative triangularity using unsuitable current scaling
+        # Only the Sauter scaling (icurr=8) is suitable for negative triangularity:
 
         if icurr != 8 and triang < 0.0:
             raise ValueError(
@@ -2092,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, eps, rminor, bt, kappa, triang) / bt
+            plascur = 1.0e6 * plasc(q95, asp, eps, rminor, bt, kappa, triang)
         elif icurr == 3:  # Simple ITER scaling (simply the cylindrical case)
             fq = 1.0
         elif icurr == 4:  # ITER formula (IPDG89)
@@ -2130,7 +2130,8 @@ def culcur(
             w07 = 1.0  # zero squareness - can be modified later if required
 
             fq = (
-                (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2)
+                (4.1e6 / 5.0e6)
+                * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2)
                 * (1.0 + 0.09 * triang + 0.16 * triang**2)
                 * (1.0 + 0.45 * triang * eps)
                 / (1.0 - 0.74 * eps)
@@ -2141,26 +2142,23 @@ def culcur(
         else:
             raise ValueError(f"Invalid value {icurr=}")
 
-        if icurr == 8:
-            curhat = 4.1e6 * rminor**2 / (rmajor * qpsi) * fq
-        elif icurr != 2:
-            curhat = 5.0e6 * rminor**2 / (rmajor * qpsi) * fq
+        if icurr != 2:
+            plascur = 5.0e6 * rminor**2 / (rmajor * q95) * fq * bt
         # == 2 case covered above
 
         qstar = (
             5.0e6
             * rminor**2
-            / (rmajor * curhat)
+            / (rmajor * plascur / bt)
             * 0.5
             * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3))
         )
 
-        plascur = curhat * bt
         physics_variables.normalised_total_beta = (
             1.0e8 * physics_variables.beta * rminor * bt / plascur
         )
         bp = bpol(
-            icurr, plascur, qpsi, asp, eps, bt, kappa, triang, pperim, constants.rmu0
+            icurr, plascur, q95, asp, eps, bt, kappa, triang, pperim, constants.rmu0
         )
 
         # Ensure current profile consistency, if required

process/plasma_geometry.py

@@ -222,24 +222,8 @@ def geomty(self):
             physics_variables.kappa,
             physics_variables.triang,
         )
-
-        #  Poloidal perimeter
-        physics_variables.pperim = 2.0e0 * (xo * thetao + xi * thetai)
-        physics_variables.sf = physics_variables.pperim / (
-            2.0e0 * numpy.pi * physics_variables.rminor
-        )
-
-        #  Volume
-        physics_variables.vol = physics_variables.cvol * self.xvol(
-            physics_variables.rmajor,
-            physics_variables.rminor,
-            xi,
-            thetai,
-            xo,
-            thetao,
-        )
-
-        #  Surface area
+        #  Surface area - inboard and outboard.  These are not given by Sauter but
+        #  the outboard area is required by DCLL and divertor
         xsi, xso = self.xsurf(
             physics_variables.rmajor,
             physics_variables.rminor,
@@ -249,10 +233,45 @@ def geomty(self):
             thetao,
         )
         physics_variables.sareao = xso
-        physics_variables.sarea = xsi + xso
 
-        #  Cross-sectional area
-        physics_variables.xarea = self.xsecta(xi, thetai, xo, thetao)
+        # icurr = 8 specifies use of the Sauter geometry as well as plasma current.
+        if physics_variables.icurr == 8:
+            (
+                physics_variables.pperim,
+                physics_variables.sf,
+                physics_variables.sarea,
+                physics_variables.xarea,
+                physics_variables.vol,
+            ) = self.Sauter_geometry(
+                physics_variables.rminor,
+                physics_variables.rmajor,
+                physics_variables.kappa,
+                physics_variables.triang,
+            )
+
+        else:
+
+            #  Poloidal perimeter
+            physics_variables.pperim = 2.0e0 * (xo * thetao + xi * thetai)
+            physics_variables.sf = physics_variables.pperim / (
+                2.0e0 * numpy.pi * physics_variables.rminor
+            )
+
+            #  Volume
+            physics_variables.vol = physics_variables.cvol * self.xvol(
+                physics_variables.rmajor,
+                physics_variables.rminor,
+                xi,
+                thetai,
+                xo,
+                thetao,
+            )
+
+            #  Cross-sectional area
+            physics_variables.xarea = self.xsecta(xi, thetai, xo, thetao)
+
+            #  Surface area - sum of inboard and outboard.
+            physics_variables.sarea = xsi + xso
 
     def xparam(self, a, kap, tri):
         """
@@ -534,3 +553,40 @@ def xsect0(self, a, kap, tri):
         xsect0 = xlo**2 * (thetao - cto * sto) + xli**2 * (thetai - cti * sti)
 
         return xsect0
+
+    def sauter_geometry(self, a, r0, kap, tri):
+        """
+        Plasma geometry based on equations (36) in O. Sauter, Fusion Engineering and Design 112 (2016) 633–645
+        'Geometric formulas for system codes including the effect of negative triangularity'
+        Author: Michael Kovari, issue #392
+        a      : input real :  plasma minor radius (m)
+        r0     : input real :  plasma major radius (m)
+        kap    : input real :  plasma separatrix elongation
+        tri    : input real :  plasma separatrix triangularity
+        """
+        w07 = 1
+        eps = a / r0
+
+        # Poloidal perimeter (named Lp in Sauter)
+        pperim = (
+            2.0e0
+            * numpy.pi
+            * a
+            * (1 + 0.55 * (kap - 1))
+            * (1 + 0.08 * tri**2)
+            * (1 + 0.2 * (w07 - 1))
+        )
+
+        # A geometric factor
+        sf = pperim / (2.0e0 * numpy.pi * a)
+
+        # Surface area (named Ap in Sauter)
+        sarea = 2.0e0 * numpy.pi * r0 * (1 - 0.32 * tri * eps) * pperim
+
+        # Cross-section area (named S_phi in Sauter)
+        xarea = numpy.pi * a**2 * kap * (1 + 0.52 * (w07 - 1))
+
+        # Volume
+        vol = 2.0e0 * numpy.pi * r0 * (1 - 0.25 * tri * eps) * xarea
+
+        return pperim, sf, sarea, xarea, vol

process/sctfcoil.py

@@ -804,7 +804,7 @@ def supercon(
             jtol = 1.0e4
 
             for lap in range(100):
-                if ttest <= 0.0e0:
+                if ttest <= delt:
                     error_handling.idiags[0] = lap
                     error_handling.fdiags[0] = ttest
                     error_handling.report_error(157)

tests/unit/test_physics.py

@@ -498,7 +498,7 @@ class CulcurParam(NamedTuple):
 
     q0: Any = None
 
-    qpsi: Any = None
+    q95: Any = None
 
     rmajor: Any = None
 
@@ -541,7 +541,7 @@ class CulcurParam(NamedTuple):
             p0=0,
             pperim=24.081367139525412,
             q0=1,
-            qpsi=3.5,
+            q95=3.5,
             rmajor=8,
             rminor=2.6666666666666665,
             sf=1.4372507312498271,
@@ -569,7 +569,7 @@ class CulcurParam(NamedTuple):
             p0=626431.90482713911,
             pperim=24.081367139525412,
             q0=1,
-            qpsi=3.5,
+            q95=3.5,
             rmajor=8,
             rminor=2.6666666666666665,
             sf=1.4372507312498271,
@@ -616,7 +616,7 @@ def test_culcur(culcurparam, monkeypatch, physics):
         p0=culcurparam.p0,
         pperim=culcurparam.pperim,
         q0=culcurparam.q0,
-        qpsi=culcurparam.qpsi,
+        q95=culcurparam.q95,
         rmajor=culcurparam.rmajor,
         rminor=culcurparam.rminor,
         sf=culcurparam.sf,