Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 0 additions & 5 deletions .flake8

This file was deleted.

5 changes: 0 additions & 5 deletions .pre-commit-config.yaml
Original file line number Diff line number Diff line change
@@ -1,9 +1,4 @@
repos:
- repo: https://github.com/PyCQA/flake8
rev: 7.1.1
hooks:
- id: flake8

- repo: https://github.com/PyCQA/isort
rev: 5.13.2
hooks:
Expand Down
35 changes: 21 additions & 14 deletions 01-freeboundary.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,39 +8,45 @@

tokamak = freegs4e.machine.TestTokamak()

eq = freegs4e.Equilibrium(tokamak=tokamak,
Rmin=0.1, Rmax=2.0, # Radial domain
Zmin=-1.0, Zmax=1.0, # Height range
nx=65, ny=65, # Number of grid points
boundary=freegs4e.boundary.freeBoundaryHagenow) # Boundary condition
eq = freegs4e.Equilibrium(
tokamak=tokamak,
Rmin=0.1,
Rmax=2.0, # Radial domain
Zmin=-1.0,
Zmax=1.0, # Height range
nx=65,
ny=65, # Number of grid points
boundary=freegs4e.boundary.freeBoundaryHagenow,
) # Boundary condition


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

profiles = freegs4e.jtor.ConstrainPaxisIp(1e3, # Plasma pressure on axis [Pascals]
2e5, # Plasma current [Amps]
2.0) # Vacuum f=R*Bt
profiles = freegs4e.jtor.ConstrainPaxisIp(
1e3, 2e5, 2.0 # Plasma pressure on axis [Pascals] # Plasma current [Amps]
) # Vacuum f=R*Bt

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

xpoints = [(1.1, -0.6), # (R,Z) locations of X-points
(1.1, 0.8)]
xpoints = [(1.1, -0.6), (1.1, 0.8)] # (R,Z) locations of X-points

isoflux = [(1.1,-0.6, 1.1,0.6)] # (R1,Z1, R2,Z2) pair of locations
isoflux = [(1.1, -0.6, 1.1, 0.6)] # (R1,Z1, R2,Z2) pair of locations

constrain = freegs4e.control.constrain(xpoints=xpoints, isoflux=isoflux)

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

freegs4e.solve(eq, # The equilibrium to adjust
profiles, # The toroidal current profile function
constrain) # Constraint function to set coil currents
freegs4e.solve(
eq, # The equilibrium to adjust
profiles, # The toroidal current profile function
constrain,
) # Constraint function to set coil currents

# eq now contains the solution

Expand Down Expand Up @@ -78,6 +84,7 @@
# Safety factor

import matplotlib.pyplot as plt

plt.plot(*eq.q())
plt.xlabel(r"Normalised $\psi$")
plt.ylabel("Safety factor")
Expand Down
7 changes: 3 additions & 4 deletions 02-read-geqdsk.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
from freegs4e import geqdsk
from freegs4e import machine
from freegs4e import geqdsk, machine
from freegs4e.plotting import plotEquilibrium

#tokamak = machine.MAST_sym()
# tokamak = machine.MAST_sym()
tokamak = machine.TestTokamak()

#with open("g014220.00200") as f:
# with open("g014220.00200") as f:
with open("lsn.geqdsk") as f:
eq = geqdsk.read(f, tokamak, show=True)

Expand Down
48 changes: 29 additions & 19 deletions 03-mast.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,44 +8,53 @@

tokamak = freegs4e.machine.MAST()

eq = freegs4e.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
eq = freegs4e.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

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

profiles = freegs4e.jtor.ConstrainPaxisIp(3e3, # Plasma pressure on axis [Pascals]
7e5, # Plasma current [Amps]
0.4) # vacuum f = R*Bt
profiles = freegs4e.jtor.ConstrainPaxisIp(
3e3, 7e5, 0.4 # Plasma pressure on axis [Pascals] # Plasma current [Amps]
) # vacuum f = R*Bt

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

xpoints = [(0.7, -1.1), # (R,Z) locations of X-points
(0.7, 1.1)]
xpoints = [(0.7, -1.1), (0.7, 1.1)] # (R,Z) locations of X-points

isoflux = [(0.7,-1.1, 1.45, 0.0) # Outboard midplane, lower X-point
,(0.7,1.1, 1.45, 0.0) # Outboard midplane, upper X-point
# ,(0.7,-1.1, 1.5,-1.9) # Lower X-point, lower outer leg
# ,(0.7,1.1, 1.5, 1.9) # Upper X-point, upper outer leg
]
isoflux = [
(0.7, -1.1, 1.45, 0.0), # Outboard midplane, lower X-point
(0.7, 1.1, 1.45, 0.0), # Outboard midplane, upper X-point
# ,(0.7,-1.1, 1.5,-1.9) # Lower X-point, lower outer leg
# ,(0.7,1.1, 1.5, 1.9) # Upper X-point, upper outer leg
]

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

constrain(eq)

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

freegs4e.solve(eq, # The equilibrium to adjust
profiles, # The plasma profiles
constrain, # Plasma control constraints
show=True) # Shows results at each nonlinear iteration
freegs4e.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

Expand All @@ -68,4 +77,5 @@

# Call matplotlib show so plot pauses
import matplotlib.pyplot as plt

plt.show()
8 changes: 3 additions & 5 deletions 04-read-mast-geqdsk.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,11 @@
from freegs4e import geqdsk
from freegs4e import machine
from freegs4e import geqdsk, machine
from freegs4e.plotting import plotEquilibrium

# Reading MAST equilibrium, up-down symmetric coils
tokamak = machine.MAST()

#with open("g014220.00200") as f:
# with open("g014220.00200") as f:
with open("mast.geqdsk") as f:
eq = geqdsk.read(f, tokamak, show=True)
eq = geqdsk.read(f, tokamak, show=True)

plotEquilibrium(eq)

29 changes: 17 additions & 12 deletions 05-fixed-boundary.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,18 +10,24 @@
# Boundary conditions
import freegs4e.boundary as boundary

profiles = freegs4e.jtor.ConstrainPaxisIp(1e3, # Plasma pressure on axis [Pascals]
1e5, # Plasma current [Amps]
1.0) # fvac = R*Bt

eq = freegs4e.Equilibrium(Rmin=0.1, Rmax=2.0,
Zmin=-1.0, Zmax=1.0,
nx=65, ny=65,
boundary=boundary.fixedBoundary)
profiles = freegs4e.jtor.ConstrainPaxisIp(
1e3, 1e5, 1.0 # Plasma pressure on axis [Pascals] # Plasma current [Amps]
) # fvac = R*Bt

eq = freegs4e.Equilibrium(
Rmin=0.1,
Rmax=2.0,
Zmin=-1.0,
Zmax=1.0,
nx=65,
ny=65,
boundary=boundary.fixedBoundary,
)

# Nonlinear solver for Grad-Shafranov equation
freegs4e.solve(eq, # The equilibrium to adjust
profiles) # The toroidal current profile function
freegs4e.solve(
eq, profiles # The equilibrium to adjust
) # The toroidal current profile function

print("Done!")

Expand All @@ -31,6 +37,5 @@

# Plot equilibrium
from freegs4e.plotting import plotEquilibrium
plotEquilibrium(eq)


plotEquilibrium(eq)
32 changes: 16 additions & 16 deletions 06-xpoints.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,22 @@
#
# Example demonstrating functions for creating and finding X-points

import matplotlib.pyplot as plt

import freegs4e

# Plotting routines
from freegs4e.plotting import plotEquilibrium, plotCoils, plotConstraints

import matplotlib.pyplot as plt
from freegs4e.plotting import plotCoils, plotConstraints, plotEquilibrium

tokamak = freegs4e.machine.TestTokamak()
eq = freegs4e.Equilibrium(tokamak=tokamak, nx=256,ny=256)
eq = freegs4e.Equilibrium(tokamak=tokamak, nx=256, ny=256)

##########################################################
# Calculate currents in coils to create X-points
# in specified locations
#
#

xpoints = [(1.1, -0.8), # (R,Z) locations of X-points
(1.1, 0.8)]
xpoints = [(1.1, -0.8), (1.1, 0.8)] # (R,Z) locations of X-points

control = freegs4e.control.constrain(xpoints=xpoints)
control(eq) # Apply control to Equilibrium eq
Expand All @@ -34,27 +33,28 @@

##########################################################
# Find critical points (O- and X-points)
#
#
#
#

import freegs4e.critical as critical

opt, xpt = critical.find_critical(eq.R, eq.Z, psi)

print("=> Found O- and X-points")

ax = plotEquilibrium(eq, show=False, oxpoints=False)
for r,z,_ in xpt:
ax.plot(r,z,'ro')
for r,z,_ in opt:
ax.plot(r,z,'go')
for r, z, _ in xpt:
ax.plot(r, z, "ro")
for r, z, _ in opt:
ax.plot(r, z, "go")
psi_bndry = xpt[0][2]
sep_contour=ax.contour(eq.R, eq.Z,psi, levels=[psi_bndry], colors='r')
sep_contour = ax.contour(eq.R, eq.Z, psi, levels=[psi_bndry], colors="r")
plt.show()

##########################################################
# Create a mask array, 1 in the core and 0 outside
#
#
#
#

mask = critical.core_mask(eq.R, eq.Z, psi, opt, xpt)

Expand Down
8 changes: 3 additions & 5 deletions 07-increase-resolution.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

import freegs4e
from freegs4e import geqdsk
from freegs4e.equilibrium import refine
Expand All @@ -7,9 +6,9 @@
# Reading MAST equilibrium, up-down symmetric coils
tokamak = freegs4e.machine.MAST()

#with open("g014220.00200") as f:
# with open("g014220.00200") as f:
with open("mast.geqdsk") as f:
eq = geqdsk.read(f, tokamak, show=True)
eq = geqdsk.read(f, tokamak, show=True)

# Increase resolution by a factor of 2
eq2 = refine(eq)
Expand All @@ -19,7 +18,6 @@

# Save to G-EQDSK
with open("mast-highres.geqdsk", "w") as f:
geqdsk.write(eq2, f)
geqdsk.write(eq2, f)

plotEquilibrium(eq2)

Loading