Skip to content

Commit

Permalink
Merge pull request #70 from freegs-plasma/updateDIIID
Browse files Browse the repository at this point in the history 
Update DIIID machine and add example
  • Loading branch information
@bendudson
bendudson committed Oct 26, 2022
2 parents 477be4c + 9e037b8 commit 4e233b7
Show file tree
Hide file tree
Showing 2 changed files with 136 additions and 22 deletions.
77 changes: 77 additions & 0 deletions 16-DIIID.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@
#!/usr/bin/env python

import freegs
import matplotlib.pyplot as plt

#########################################
# Create the machine, which specifies coil locations
# and equilibrium, specifying the domain to solve over

tokamak = freegs.machine.DIIID()

eq = freegs.Equilibrium(tokamak=tokamak,
Rmin=0.1, Rmax=2.8, # Radial domain
Zmin=-1.8, Zmax=1.8, # Height range
nx=129, ny=129) # Number of grid points

#########################################
# Plasma profiles

profiles = freegs.jtor.ConstrainPaxisIp(eq,
159811, # Plasma pressure on axis [Pascals]
-1533632, # Plasma current [Amps]
-3.231962138124) # vacuum f = R*Bt

#########################################
# Coil current constraints
#
# Specify locations of the X-points
# to use to constrain coil currents

xpoints = [(1.285, -1.176), # (R,Z) locations of X-points
(1.2, 1.0)]

isoflux = [(1.285, -1.176, 1.2 ,1.2)] # (R1,Z1, R2,Z2) pair of locations

constrain = freegs.control.constrain(xpoints=xpoints, gamma=1e-12, isoflux=isoflux)

constrain(eq)

#########################################
# Nonlinear solve

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("Plasma poloidal beta: %e" % (eq.poloidalBeta()))
print("Plasma volume: %e m^3" % (eq.plasmaVolume()))

eq.tokamak.printCurrents()

# plot equilibrium
axis = eq.plot(show=False)
tokamak.plot(axis=axis, show=False)
constrain.plot(axis=axis, show=True)

# Safety factor
plt.plot(*eq.q())
plt.xlabel(r"Normalised $\psi$")
plt.ylabel("Safety factor")
plt.grid()
plt.show()

##############################################
# Save to geqdsk file

from freegs import geqdsk

with open("diiid.geqdsk", "w") as f:
geqdsk.write(eq, f)
81 changes: 59 additions & 22 deletions freegs/machine.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -922,32 +922,69 @@ def TestTokamakSensor():

def DIIID():
"""
PF coil set from ef20030203.d3d
Taken from Corsica
PF coil set and boundary from g file and GA webpage
"""

coils = [
{"label": "F1A", "R": 0.8608, "Z": 0.16830, "current": 0.0},
{"label": "F2A", "R": 0.8614, "Z": 0.50810, "current": 0.0},
{"label": "F3A", "R": 0.8628, "Z": 0.84910, "current": 0.0},
{"label": "F4A", "R": 0.8611, "Z": 1.1899, "current": 0.0},
{"label": "F5A", "R": 1.0041, "Z": 1.5169, "current": 0.0},
{"label": "F6A", "R": 2.6124, "Z": 0.4376, "current": 0.0},
{"label": "F7A", "R": 2.3733, "Z": 1.1171, "current": 0.0},
{"label": "F8A", "R": 1.2518, "Z": 1.6019, "current": 0.0},
{"label": "F9A", "R": 1.6890, "Z": 1.5874, "current": 0.0},
{"label": "F1B", "R": 0.8608, "Z": -0.1737, "current": 0.0},
{"label": "F2B", "R": 0.8607, "Z": -0.5135, "current": 0.0},
{"label": "F3B", "R": 0.8611, "Z": -0.8543, "current": 0.0},
{"label": "F4B", "R": 0.8630, "Z": -1.1957, "current": 0.0},
{"label": "F5B", "R": 1.0025, "Z": -1.5169, "current": 0.0},
{"label": "F6B", "R": 2.6124, "Z": -0.44376, "current": 0.0},
{"label": "F7B", "R": 2.3834, "Z": -1.1171, "current": 0.0},
{"label": "F8B", "R": 1.2524, "Z": -1.6027, "current": 0.0},
{"label": "F9B", "R": 1.6889, "Z": -1.578, "current": 0.0},
("FC1",ShapedCoil([[1.6042, -1.64455], [1.7736, -1.64455], [1.7736, -1.51145], [1.6042, -1.51145]])),
("FC2",ShapedCoil([[0.8354, 0.00777], [0.8862, 0.00777], [0.8862, 0.32883], [0.8354, 0.32883]])),
("FC3",ShapedCoil([[0.836, 0.34757], [0.8868, 0.34757], [0.8868, 0.66863], [0.836, 0.66863]])),
("FC4",ShapedCoil([[0.8374, 0.68857], [0.8882, 0.68857], [0.8882, 1.00963], [0.8374, 1.00963]])),
("FC5",ShapedCoil([[0.8357, 1.02937], [0.8865, 1.02937], [0.8865, 1.35043], [0.8357, 1.35043]])),
("FC6",ShapedCoil([[0.9345, 1.3876], [1.0737, 1.5268], [1.0737, 1.6462], [0.9345, 1.507]])),
("FC7",ShapedCoil([[2.5298, 0.3403], [2.703, 0.3403], [2.6949, 0.5349], [2.5217, 0.5349]])),
("FC8",ShapedCoil([[2.5387, 1.0325], [2.7267, 1.0325], [2.2078, 1.2017], [2.0198, 1.2017]])),
("FC9",ShapedCoil([[1.13435, 1.55935], [1.36925, 1.55935], [1.36925, 1.64445], [1.13435, 1.64445]])),
("FC10",ShapedCoil([[1.6043, 1.52085], [1.7737, 1.52085], [1.7737, 1.65395], [1.6043, 1.65395]])),
("FC11",ShapedCoil([[0.8354, -0.33423], [0.8862, -0.33423], [0.8862, -0.01317], [0.8354, -0.01317]])),
("FC12",ShapedCoil([[0.8353, -0.67403], [0.8861, -0.67403], [0.8861, -0.35297], [0.8353, -0.35297]])),
("FC13",ShapedCoil([[0.8357, -1.01483], [0.8865, -1.01483], [0.8865, -0.69377], [0.8357, -0.69377]])),
("FC14",ShapedCoil([[0.8376, -1.35623], [0.8884, -1.35623], [0.8884, -1.03517], [0.8376, -1.03517]])),
("FC15",ShapedCoil([[0.9329, -1.507], [1.0721, -1.6462], [1.0721, -1.5268], [0.9329, -1.3876]])),
("FC16",ShapedCoil([[2.5217, -0.5349], [2.6949, -0.5349], [2.703, -0.3403], [2.5298, -0.3403]])),
("FC17",ShapedCoil([[2.0299, -1.2017], [2.2179, -1.2017], [2.7368, -1.0325], [2.5488, -1.0325]])),
("FC18",ShapedCoil([[1.13495, -1.64525], [1.36985, -1.64525], [1.36985, -1.56015], [1.13495, -1.56015]]))
]

return Machine(coils)
R = np.array([1.016,1.016,1.016,1.016,1.016,1.016,1.016,1.016,1.016,
1.016,1.016,1.012,1.001,1.029,1.042,1.046,1.056,1.097,
1.108,1.116,1.134,1.148,1.162,1.181,1.182,1.185,1.19,
1.195,1.201,1.209,1.215,1.222,1.228,1.234,1.239,1.242,
1.248,1.258,1.263,1.28,1.28,1.28,1.31,1.328,1.361,
1.38,1.419,1.419,1.372,1.37167,1.37003,1.36688,1.36719,1.37178,
1.37224,1.38662,1.38708,1.40382,1.41127,1.41857,1.421,1.48663,1.4973,
1.49762,1.49745,1.49275,1.4926,1.49261,1.49279,1.4934,1.4947,1.49622,
1.47981,1.48082,1.48149,1.48646,1.49095,1.50305,1.59697,1.6255,1.63752,
1.647,1.785,2.07,2.128,2.245,2.33956,2.34708,2.34913,2.35103,
2.35158,2.35125,2.35051,2.34965,2.3487,2.3476,2.3402,2.32942,2.134,
1.786,1.768,1.768,1.682,1.372,1.372,1.42,1.42,1.273,
1.153,1.016,1.016,1.016,1.016,1.016,1.016,1.016,1.016])
Z = np.array([ 0.00000e+00,9.64000e-01,9.68000e-01,1.00100e+00,1.01900e+00,
1.07700e+00,1.07000e+00,1.09600e+00,1.11300e+00,1.13800e+00,
1.14700e+00,1.16500e+00,1.21700e+00,1.21700e+00,1.16240e+00,
1.16238e+00,1.16260e+00,1.16450e+00,1.16594e+00,1.16591e+00,
1.16896e+00,1.17175e+00,1.17556e+00,1.18300e+00,1.18350e+00,
1.18500e+00,1.18800e+00,1.19100e+00,1.19600e+00,1.20200e+00,
1.20800e+00,1.21400e+00,1.22100e+00,1.23100e+00,1.23800e+00,
1.24400e+00,1.25400e+00,1.27800e+00,1.29000e+00,1.33100e+00,
1.34700e+00,1.34800e+00,1.34800e+00,1.34800e+00,1.34800e+00,
1.34800e+00,1.34800e+00,1.31000e+00,1.31000e+00,1.29238e+00,
1.28268e+00,1.25644e+00,1.22955e+00,1.19576e+00,1.19402e+00,
1.16487e+00,1.16421e+00,1.15696e+00,1.15730e+00,1.16132e+00,
1.16400e+00,1.24050e+00,1.23458e+00,1.23428e+00,1.23174e+00,
1.21330e+00,1.21061e+00,1.20486e+00,1.20214e+00,1.19642e+00,
1.18511e+00,1.16070e+00,1.12426e+00,1.12256e+00,1.12138e+00,
1.11692e+00,1.11439e+00,1.11244e+00,1.09489e+00,1.08530e+00,
1.07988e+00,1.07700e+00,1.07700e+00,1.04000e+00,9.93000e-01,
7.09000e-01,4.61430e-01,4.15830e-01,2.72180e-01,1.70180e-01,
7.01200e-02, -3.17900e-02, -1.44350e-01, -2.14830e-01, -3.26690e-01,
-3.86770e-01, -4.53040e-01, -4.77570e-01, -9.73000e-01, -1.17400e+00,
-1.21100e+00, -1.25000e+00, -1.25000e+00, -1.25000e+00, -1.32900e+00,
-1.32900e+00, -1.36300e+00, -1.36300e+00, -1.36300e+00, -1.22300e+00,
-1.22300e+00, -8.30000e-01, -8.00000e-01, -4.15000e-01, -4.00000e-01,
-1.00000e-03,0.00000e+00])
wall = Wall(R,Z)

return Machine(coils,wall=wall)


def MAST():
Expand Down

0 comments on commit 4e233b7

Please sign in to comment.