MAST-U example, improvements

Centre column coil modelled as a short solenoid

Getting a super-X configuration requires some manual changing
of coil currents at the moment; the control system isn't really very good.

Changed tolerance when removing X-points, as the wrong
primary X-point was being chosen in many cases.

08-mast-upgrade.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+
+import freegs
+from freegs.plotting import plotConstraints
+
+#########################################
+# Create the machine, which specifies coil locations
+# and equilibrium, specifying the domain to solve over
+
+tokamak = freegs.machine.MASTU()
+
+
+eq = freegs.Equilibrium(tokamak=tokamak,
+                        Rmin=0.1, Rmax=2.0,    # Radial domain
+                        Zmin=-2.1, Zmax=2.1,   # Height range
+                        nx=65, ny=65)          # Number of grid points
+
+#########################################
+# Plasma profiles
+
+profiles = freegs.jtor.ConstrainPaxisIp(2e4, # Plasma pressure on axis [Pascals]
+                                        6e5, # Plasma current [Amps]
+                                        1.0) # vacuum f = R*Bt
+
+#########################################
+# Coil current constraints
+#
+# Specify locations of the X-points
+# to use to constrain coil currents
+
+Rx = 0.509
+Zx = 1.291
+
+Rmid = 1.34   # Outboard midplane
+Rin = 0.3581  # Inboard midplane
+
+xpoints = [(Rx, -Zx),   # (R,Z) locations of X-points
+           (Rx,  Zx)]
+
+isoflux = [(Rx,-Zx, Rmid, 0.0)   # Outboard midplane, lower X-point
+           ,(Rx,Zx, Rmid, 0.0)   # Outboard midplane, upper X-point
+
+           ,(Rmid, 0.0, Rin, 0.0)
+
+           # ,(Rx,-Zx, 0.73, -1.6)
+           # ,(Rx, Zx, 0.73,  1.6)
+
+           # ,(Rmid, 0.0, 0.73, -1.6)
+           # ,(Rmid, 0.0, 0.73,  1.6)
+
+           # ,(Rx,-Zx, 0.45, -0.95)
+           # ,(Rx, Zx, 0.45,  0.95)
+
+           # ,(Rin, 0.0, 0.45, -0.95)
+           # ,(Rin, 0.0, 0.45,  0.95)
+           # ,(Rin, 0.0, 0.45, -0.95)
+           # ,(Rin, 0.0, 0.45,  0.95)
+           # ,(Rin, 0.0, 0.45, -0.95)
+           # ,(Rin, 0.0, 0.45,  0.95)
+
+
+            ,(Rx,-Zx, 0.95, -1.77)
+            ,(Rx, Zx, 0.95,  1.77)
+
+           # ,(Rx,-Zx, 1.25, -1.9)
+           # ,(Rx, Zx, 1.25,  1.9)
+
+           ]
+
+constrain = freegs.control.constrain(xpoints=xpoints, gamma=1e-5, isoflux=isoflux)
+
+constrain(eq)
+
+#########################################
+# Nonlinear solve
+
+freegs.solve(eq,          # The equilibrium to adjust
+             profiles,    # The plasma profiles
+             constrain,   # Plasma control constraints
+             show=False)   # Shows results at each nonlinear iteration
+
+isoflux = [(Rx,-Zx, Rmid, 0.0)   # Outboard midplane, lower X-point
+           ,(Rx,Zx, Rmid, 0.0)   # Outboard midplane, upper X-point
+
+           ,(Rmid, 0.0, Rin, 0.0)
+
+           ,(Rx,-Zx, 0.95, -1.77)
+           ,(Rx, Zx, 0.95,  1.77)
+
+           ,(Rx,-Zx, 0.76, -1.58)
+           ,(Rx, Zx, 0.76,  1.58)
+
+           ,(Rx,-Zx, 1.25, -1.8)
+           ,(Rx, Zx, 1.25,  1.8)
+
+           ]
+
+constrain = freegs.control.constrain(xpoints=xpoints, gamma=1e-12, isoflux=isoflux)
+
+# Turn off feedback control for all coils
+for label, coil in tokamak.coils:
+    coil.control = False
+
+# Centre column coil
+tokamak["Pc"].current = -3e4
+
+# Turn on vertical feedback control
+tokamak["P61"].control = True
+
+# Coil in the "nose" of the divertor
+tokamak["Dp"].current = -1000.0
+
+# At top of divertor chamber
+tokamak["D6"].current = -500.0
+
+# X-point location
+tokamak["Px"].current = 2000.0
+
+tokamak["D1"].current = 1000.0
+
+# Coil in outer corner 
+tokamak["D5"].current = 1500.
+
+# Coil at bottom centre
+tokamak["D3"].current = 2800
+
+freegs.solve(eq,          # The equilibrium to adjust
+             profiles,    # The plasma profiles
+             constrain,   # Plasma control constraints
+             show=True)   # Shows results at each nonlinear iteration
+
+# eq now contains the solution
+
+print("Done!")
+
+print("Plasma current: %e Amps" % (eq.plasmaCurrent()))
+print("Pressure on axis: %e Pascals" % (eq.pressure(0.0)))
+print("Poloidal beta: %e" % (eq.poloidalBeta()))
+
+eq.tokamak.printCurrents()
+
+axis = eq.plot(show=False)
+constrain.plot(axis=axis)
+
+##############################################
+# Save to geqdsk file
+
+from freegs import geqdsk
+
+with open("mast-upgrade.geqdsk", "w") as f:
+    geqdsk.write(eq, f)
+
+# Call matplotlib show so plot pauses
+import matplotlib.pyplot as plt
+plt.show()

freegs/critical.py

@@ -218,8 +218,8 @@ def remove_dup(points):
             # rather than the distance in space
 
             maxp = amax(pline)
-            if (maxp - pline[-1])/(maxp - pline[0]) > 0.05:
-                # More than 5% drop in psi from maximum to X-point
+            if (maxp - pline[-1])/(maxp - pline[0]) > 0.001:
+                # More than 0.1% drop in psi from maximum to X-point
                 # -> Discard
                 continue
 
@@ -229,7 +229,7 @@ def remove_dup(points):
                 continue
             xpt_keep.append(xpt)
         xpoint = xpt_keep
-
+        
     # Sort X-points by distance to primary O-point in psi space
     psi_axis = opoint[0][2]
     xpoint.sort(key=lambda x: (x[2] - psi_axis)**2)

freegs/machine.py

@@ -874,7 +874,8 @@ def TCV():
 def MASTU():
     coils = [
         ("Solenoid", Solenoid(0.19475, -1.581, 1.581, 324)),
-        ("Pc", Coil(0.067, 0.0, turns=142)),
+        #("Pc", Coil(0.067, 0.0, turns=142)),
+        ("Pc", Solenoid(0.067, -0.6, 0.6, 142)),
         ("Px", Circuit([("PxU", Coil(0.2405,  1.2285, turns=44), 1.0),
                         ("PxL", Coil(0.2405, -1.2285, turns=44), 1.0)])),
         ("D1", Circuit([("D1U", Coil(0.381,  1.555, turns=35), 1.0),
@@ -904,7 +905,7 @@ def MASTU():
                          ("P62L", Coil(1.2575, -1.0575, turns=2), -1.0)]))
         ]
 
-    rwall = [1.2503, 1.3483, 1.47, 1.47, 1.45, 1.45, 1.3214, 1.1904,
+    rwall = [1.6, 1.2503, 1.3483, 1.47, 1.47, 1.45, 1.45, 1.3214, 1.1904,
              0.89296, 0.86938, 0.83981, 0.82229, 0.81974, 0.81974,
              0.82734, 0.8548, 0.89017, 0.91974, 0.94066, 1.555, 1.85,
              2, 2, 2, 2, 1.3188, 1.7689, 1.7301, 1.35, 1.09, 1.09,
@@ -913,7 +914,7 @@ def MASTU():
              0.4788, 0.333, 0.333, 0.275, 0.334, 0.261, 0.261, 0.244,
              0.261, 0.261, 0.244, 0.261, 0.261]
 
-    zwall = [1, 0.86, 0.86, 0.81, 0.81, 0.82, 0.82, 1.007, 1.304,
+    zwall = [1, 1, 0.86, 0.86, 0.81, 0.81, 0.82, 0.82, 1.007, 1.304,
              1.3312, 1.3826, 1.4451, 1.4812, 1.4936, 1.5318, 1.5696,
              1.5891, 1.5936, 1.5936, 1.567, 1.08, 1.08, 1.7, 2.035,
              2.169, 2.169, 1.7189, 1.68, 2.06, 2.06, 2.06, 1.8786,