freegs.equilibrium.Equilibrium

class freegs.equilibrium.Equilibrium(tokamak=Machine(coils=[], wall=None), Rmin=0.1, Rmax=2.0, Zmin=-1.0, Zmax=1.0, nx=65, ny=65, boundary=<function freeBoundary>, psi=None, current=0.0, order=4, check_limited=False)

Bases: object

Represents the equilibrium state, including plasma and coil currents

R[nx, ny]

Read-only

Z[nx, ny]

Read-only

Rmin, Rmax

Read-only

Zmin, Zmax

Read-only

tokamak

The coils and circuits

_applyBoundary

_solver

Grad-Shafranov elliptic solver

_profiles

An object which calculates the toroidal current

_constraints

Control system which adjusts coil currents to meet constraints e.g. X-point location and flux values

__init__(tokamak=Machine(coils=[], wall=None), Rmin=0.1, Rmax=2.0, Zmin=-1.0, Zmax=1.0, nx=65, ny=65, boundary=<function freeBoundary>, psi=None, current=0.0, order=4, check_limited=False)

Initialises a plasma equilibrium

Parameters:

• tokamak (Machine) – The set of coils and currents

• Rmin (float) – Range of major radius R [m]

• Rmax (float) – Range of major radius R [m]

• Zmin (float) – Range of height Z [m]

• Zmax (float) – Range of height Z [m]

• nx (int) – Resolution in R. This must be 2^n + 1

• ny (int) – Resolution in Z. This must be 2^m + 1

• boundary – The boundary condition, either freeBoundary or fixedBoundary

• psi (Optional[float]) – Magnetic flux. If None, use concentric circular flux surfaces as starting guess

• current (float) – Plasma current

• order (int) – The order of the differential operators to use. Valid values are 2 or 4

• check_limited (bool) – Checks if the plasma is limited

Methods

Bpol(R, Z)	Total poloidal magnetic field
Br(R, Z)	Total radial magnetic field
Btor(R, Z)	Toroidal magnetic field
Btot(R, Z)	Total magnetic field
Bz(R, Z)	Total vertical magnetic field
Rgeometric([npoints])	Locates major radius R of the geometric major radius.
Rmagnetic()	The major radius R of magnetic major radius
Zgeometric([npoints])	Locates the height z of the geometric axis.
Zmagnetic()	The height Z of magnetic axis
__init__([tokamak, Rmin, Rmax, Zmin, Zmax, ...])	Initialises a plasma equilibrium
aspectRatio([npoints])	Calculates the plasma aspect ratio.
betaN([npoints])	Calculate normalised plasma beta
callSolver(psi, rhs)	Calls the psi solver, passing the initial guess and RHS arrays
effectiveElongation([npoints])	Calculates plasma effective elongation using the plasma volume
elongation([npoints])	Calculates the elongation, kappa, of the plasma.
elongationLower([npoints])	Calculates the lower elongation, kappa_l, of the plasma.
elongationUpper([npoints])	Calculates the upper elongation, kappa_u, of the plasma.
ffprime(psinorm)	Return ff' at given normalised psi
fpol(psinorm)	Return f = R*Bt at specified values of normalised psi
fvac()	Return vacuum f = R*Bt
geometricAxis([npoints])	Locates geometric axis, returning [R,Z].
getForces()	Calculate forces on the coils
getMachine()	Returns the handle of the machine, including coils
innerOuterSeparatrix([Z])	Locate R co ordinates of separatrix at both inboard and outboard poloidal midplane (Z = 0)
internalInductance([npoints])	Calculates plasma internal inductance li
internalInductance1([npoints])	Calculates li1 plasma internal inductance
internalInductance2()	Calculates li2 plasma internal inductance
internalInductance3([npoints])	Calculates li3 plasma internal inductance
intersectsWall()	Assess whether or not the core plasma touches the vessel walls.
inverseAspectRatio([npoints])	Calculates inverse of the plasma aspect ratio.
magneticAxis()	Returns the location of the magnetic axis as a list [R,Z,psi]
minorRadius([npoints])	Calculates minor radius of the plasma, a.
plasmaBr(R, Z)	Radial magnetic field due to plasma Br = -1/R dpsi/dZ
plasmaBz(R, Z)	Vertical magnetic field due to plasma Bz = (1/R) dpsi/dR
plasmaCurrent()	Plasma current [Amps]
plasmaVolume()	Calculate the volume of the plasma in m^3
plot([axis, show, oxpoints])	Plot the equilibrium flux surfaces
poloidalBeta()	Calculate plasma poloidal beta by integrating the thermal pressure and poloidal magnetic field pressure over the plasma volume.
poloidalBeta2()	Return the poloidal beta betap = (8pi/mu0) * int(p)dRdZ / Ip^2
poloidalBeta3()	Calculates alterantive poloidal beta definition.
pprime(psinorm)	Return p' at given normalised psi
pressure(psinorm)	Returns plasma pressure at specified values of normalised psi
pressure_ave()	Calculate average pressure, Pa.
printForces()	Prints a table of forces on coils
psi()	Total poloidal flux (psi), including contribution from plasma and external coils.
psiN()	Total poloidal flux (psi), including contribution from plasma and external coils.
psiNRZ(R, Z)	Return poloidal flux psi at given (R,Z) location.
psiRZ(R, Z)	Return poloidal flux psi at given (R,Z) location
q([psinorm, npsi])	Returns safety factor q at specified values of normalised psi
qcyl()	Cylindrical safety factor.
rhotor([psi])	Calculates normalised toroidal flux at specified values of poloidal flux.
separatrix([npoints])	Returns an array of npoints (R, Z) coordinates of the separatrix, equally spaced in geometric poloidal angle.
setSolver(solver)	Sets the linear solver to use.
setSolverVcycle([nlevels, ncycle, niter, direct])	Creates a new linear solver, based on the multigrid code
shafranovShift([npoints])	Calculates the plasma shafranov shift [delta_shafR,delta_shafZ] where
solve(profiles[, Jtor, psi, psi_bndry])	Calculate the plasma equilibrium given new profiles replacing the current equilibrium.
tor_flux([psi])	Calculates toroidal flux at specified values of poloidal flux.
toroidalBeta()	Calculate plasma toroidal beta by integrating the thermal pressure and toroidal magnetic field pressure over the plasma volume.
totalBeta()	Calculate plasma total beta
triangularity([npoints])	Calculates plasma triangularity, delta.
triangularityLower([npoints])	Calculates plasma upper triangularity, delta_u.
triangularityUpper([npoints])	Calculates plasma upper triangularity, delta_u.
w_th()	Stored thermal energy in plasma, J.

Bpol(R, Z)

Total poloidal magnetic field

Br(R, Z)

Total radial magnetic field

Btor(R, Z)

Toroidal magnetic field

Btot(R, Z)

Total magnetic field

Bz(R, Z)

Total vertical magnetic field

Rgeometric(npoints=360)

Locates major radius R of the geometric major radius.

Rmagnetic()

The major radius R of magnetic major radius

Zgeometric(npoints=360)

Locates the height z of the geometric axis.

Zmagnetic()

The height Z of magnetic axis

aspectRatio(npoints=360)

Calculates the plasma aspect ratio.

A = R0/a where R0 = major radius, a = minor radius.

betaN(npoints=360)

Calculate normalised plasma beta

callSolver(psi, rhs)

Calls the psi solver, passing the initial guess and RHS arrays

Parameters:

• psi – Initial guess for the solution (used if iterative)

• rhs

Returns:

Solution psi

Return type:

float

effectiveElongation(npoints=360)

Calculates plasma effective elongation using the plasma volume

elongation(npoints=360)

Calculates the elongation, kappa, of the plasma. A large number of points should be supplied such that any primary xpoint(s) on the LCFS are captured. The R,Z of the primary x-point is NOT itself included in the R,Z of the LCFS as the plasma may be limited. P2 is the point at the upper extent of the plasma, and P4 is the point at the lower extent of the plasma.

kappa = (Z(P2) - Z(P4))/a

elongationLower(npoints=360)

Calculates the lower elongation, kappa_l, of the plasma. A large number of points should be supplied such that any primary xpoint(s) on the LCFS are captured. The R,Z of the primary x-point is NOT itself included in the R,Z of the LCFS as the plasma may be limited. P2 is the point at the upper extent of the plasma.

kappa_u = (z0 - Z(P4))/a

elongationUpper(npoints=360)

Calculates the upper elongation, kappa_u, of the plasma. A large number of points should be supplied such that any primary xpoint(s) on the LCFS are captured. The R,Z of the primary x-point is NOT itself included in the R,Z of the LCFS as the plasma may be limited. P2 is the point at the upper extent of the plasma.

kappa_u = (Z(P2) - z0)/a

ffprime(psinorm)

Return ff’ at given normalised psi

fpol(psinorm)

Return f = R*Bt at specified values of normalised psi

fvac()

Return vacuum f = R*Bt

geometricAxis(npoints=360)

Locates geometric axis, returning [R,Z]. First locates the extrema points in R of the LCFS, wherein P3 is at the IMP and P1 is at the OMP.

R0 = R(P3) + 0.5*(R(P1)-R(P3)) z0 = 0.5*(Z(P1)+Z(P3))

getForces()

Calculate forces on the coils

Returns a dictionary of coil label -> force

getMachine()

Returns the handle of the machine, including coils

innerOuterSeparatrix(Z=0.0)

Locate R co ordinates of separatrix at both inboard and outboard poloidal midplane (Z = 0)

internalInductance(npoints=360)

Calculates plasma internal inductance li

li = 4/(mu0*R0*Ip^2) * int(2piR*(Bp^2/2mu0)*dR*dZ)

= 2/(mu0^2*R0*Ip^2)*int(2piR*Bp^2*dR*dZ)

internalInductance1(npoints=360)

Calculates li1 plasma internal inductance

internalInductance2()

Calculates li2 plasma internal inductance

internalInductance3(npoints=360)

Calculates li3 plasma internal inductance

intersectsWall()

Assess whether or not the core plasma touches the vessel walls. Returns True if it does intersect.

inverseAspectRatio(npoints=360)

Calculates inverse of the plasma aspect ratio.

epsilon = 1/A A = R0/a where R0 = major radius, a = minor radius.

magneticAxis()

Returns the location of the magnetic axis as a list [R,Z,psi]

minorRadius(npoints=360)

Calculates minor radius of the plasma, a. First locates the extrema points in R of the LCFS, wherein P3 is at the IMP and P1 is at the OMP.

a = 0.5*(R(P1) - R(P3))

plasmaBr(R, Z)

Radial magnetic field due to plasma Br = -1/R dpsi/dZ

plasmaBz(R, Z)

Vertical magnetic field due to plasma Bz = (1/R) dpsi/dR

plasmaCurrent()

Plasma current [Amps]

plasmaVolume()

Calculate the volume of the plasma in m^3

plot(axis=None, show=True, oxpoints=True)

Plot the equilibrium flux surfaces

Parameters:

• axis – Specify the axis on which to plot

• show – Call matplotlib.pyplot.show() before returning

• oxpoints – Plot X points as red circles, O points as green circles

Return type:

Axes

poloidalBeta()

Calculate plasma poloidal beta by integrating the thermal pressure and poloidal magnetic field pressure over the plasma volume.

poloidalBeta2()

Return the poloidal beta betap = (8pi/mu0) * int(p)dRdZ / Ip^2

poloidalBeta3()

Calculates alterantive poloidal beta definition.

pprime(psinorm)

Return p’ at given normalised psi

pressure(psinorm)

Returns plasma pressure at specified values of normalised psi

pressure_ave()

Calculate average pressure, Pa.

printForces()

Prints a table of forces on coils

psi()

Total poloidal flux (psi), including contribution from plasma and external coils.

psiN()

Total poloidal flux (psi), including contribution from plasma and external coils. Normalised such that psiN = 0 on the magnetic axis and 1 on the LCFS.

psiNRZ(R, Z)

Return poloidal flux psi at given (R,Z) location. Normalised such that psiN = 0 on the magnetic axis and 1 on the LCFS.

psiRZ(R, Z)

Return poloidal flux psi at given (R,Z) location

q(psinorm=None, npsi=100)

Returns safety factor q at specified values of normalised psi

psinorm is a scalar, list or array of floats betweem 0 and 1.

>>> safety_factor = eq.q([0.2, 0.5, 0.9])

If psinorm is None, then q on a uniform psi grid will be returned, along with the psi values

>>> psinorm, q = eq.q()

Note: psinorm = 0 is the magnetic axis, and psinorm = 1 is the separatrix.

If you request a value close to either of these limits, an extrapolation based on a 1D grid of values from 0.01 to 0.99 will be used. This gives smooth and continuous q-profiles, but comes at an increased computational cost.

qcyl()

Cylindrical safety factor.

rhotor(psi=None)

Calculates normalised toroidal flux at specified values of poloidal flux.

>>> rhotor = sqrt ( tor_flux/max(tor_flux)).

Maximum toroidal flux shoud be at LCFS.

separatrix(npoints=360)

Returns an array of npoints (R, Z) coordinates of the separatrix, equally spaced in geometric poloidal angle.

setSolver(solver)

Sets the linear solver to use. The given object/function must have a __call__ method which takes two inputs

solver(x, b)

where x is the initial guess. This should solve Ax = b, returning the result.

setSolverVcycle(nlevels=1, ncycle=1, niter=1, direct=True)

Creates a new linear solver, based on the multigrid code

nlevels - Number of resolution levels, including original ncycle - The number of V cycles niter - Number of linear solver (Jacobi) iterations per level direct - Use a direct solver at the coarsest level?

shafranovShift(npoints=360)

Calculates the plasma shafranov shift [delta_shafR,delta_shafZ] where

delta_shafR = Rmagnetic - Rgeo delta_shafR = Zmagnetic - z0

solve(profiles, Jtor=None, psi=None, psi_bndry=None)

Calculate the plasma equilibrium given new profiles replacing the current equilibrium.

This performs the linear Grad-Shafranov solve

Parameters:

• profiles –

  An object describing the plasma profiles. At minimum this must have methods:

  • .Jtor(R, Z, psi, psi_bndry) -> [nx, ny]

  • .pprime(psinorm)

  • .ffprime(psinorm)

  • .pressure(psinorm)

  • .fpol(psinorm)

• Jtor – If supplied, a 2D array specifying the toroidal current at each (R,Z) point. If not supplied, Jtor is calculated from profiles by finding O,X-points

• psi_bndry – Poloidal flux to use as the separatrix (plasma boundary). If not given then X-point locations are used.

tor_flux(psi=None)

Calculates toroidal flux at specified values of poloidal flux. >>> q = drho/dpsi

toroidalBeta()

Calculate plasma toroidal beta by integrating the thermal pressure and toroidal magnetic field pressure over the plasma volume.

totalBeta()

Calculate plasma total beta

triangularity(npoints=360)

Calculates plasma triangularity, delta.

Here delta is defined as the average of the upper and lower triangularities.

triangularityLower(npoints=360)

Calculates plasma upper triangularity, delta_u. P4 is the point at the lower extent of the plasma.

tri_l = (R0 - R(P4))/a

triangularityUpper(npoints=360)

Calculates plasma upper triangularity, delta_u. P2 is the point at the upper extent of the plasma.

tri_u = (R0 - R(P2))/a

w_th()

Stored thermal energy in plasma, J.