An example of a dead start (from scratch) for the MAST tokamak. A solution to the Grad-Shafranov equation requires
In [1]:
# Step 1: Specify the locations of the coils, and the domain you want to solve over
from freegs import machine
from freegs.equilibrium import Equilibrium
# Define the poloidal field coil set
tokamak = machine.MAST()
# Define the domain to solve over
eq = Equilibrium(tokamak=tokamak,
Rmin=0.1, Rmax=2.0, # Radial domain
Zmin=-2.0, Zmax=2.0, # Height range
nx=65, ny=65) # Number of grid points
Define toroidal current function. The RHS of the Grad-Shafranov equation is determined by the toroidal current density
$j_{tor} = -R\frac{\partial p}{\partial \psi} - \frac{1}{\mu_0 R}f\frac{\partial p}{\partial \psi}$
The profiles can be adjusted to match global constrants, such as the pressure on axis and the total plasma current $I_p$
The toroidal current function jtorfunc calculates the toroidal current density $j{tor}$ given $\psi$ on an R-Z mesh
In [2]:
# Step 2: Specify the profiles of pressure and f=R*Bt.
# Currently quite simple functions are supported
from freegs.jtor import ConstrainPaxisIp
profiles = ConstrainPaxisIp(3e3, # Plasma pressure on axis [Pascals]
7e5, # Plasma current [Amps]
0.4) # vacuum f = R*Bt
In [3]:
# Step 3: Specify the control system and feedback variables.
# The control system adjusts the currents in the poloidal field coils
# to produce X-points in the desired locations, and ensure that the desired
# pairs of locations have the same poloidal flux.
from freegs import control
xpoints = [(0.7, -1.1), # (R,Z) locations of X-points
(0.7, 1.1)]
# Contstrain these pairs of (R,Z, R,Z) locations to have the same poloidal flux
# This is needed for radial and vertical position control of the plasma.
isoflux = [(0.7,-1.1, 1.45, 0.0) # Lower X-point, Outboard midplane
,(0.7,1.1, 1.45, 0.0) # Upper X-point, Outboard midplane
]
constrain = control.constrain(xpoints=xpoints, gamma=1e-12, isoflux=isoflux)
constrain(eq)
In [4]:
# With these three components (coils, profiles and constraints), solve the nonlinear
# system with a Picard iteration. This modifies the "eq" object.
from freegs import picard
picard.solve(eq, # The equilibrium to adjust
profiles, # The toroidal current profile function
constrain) # Constraint function to set coil currents
In [5]:
print("Plasma current: %e Amps" % (eq.plasmaCurrent()))
In [6]:
tokamak.printCurrents()
In [7]:
%matplotlib inline
from freegs.plotting import plotEquilibrium
plotEquilibrium(eq)
Out[7]:
In [8]:
xpoints = [(0.7, -1.0), # (R,Z) locations of X-points
(0.7, 1.0)]
isoflux = [(0.7,-1.0, 1.4, 0.0),(0.7,1.0, 1.4, 0.0), (0.7,-1.0, 0.3, 0.0)]
constrain = control.constrain(xpoints=xpoints, gamma=1e-12, isoflux=isoflux)
constrain(eq)
plotEquilibrium(eq)
Out[8]:
In [9]:
picard.solve(eq, # The equilibrium to adjust
profiles, # The toroidal current profile function
constrain) # Constraint function to set coil currents
In [10]:
plotEquilibrium(eq)
Out[10]:
In [ ]: