Iota Calculation, Fourier Surface Extraction, and Perturbations (Wendelstein 7-X)

This example shows how to calculate the rotational transform for Wendelstein 7-X.

import fusionsc as fsc
from fusionsc.devices import w7x
import numpy as np
import matplotlib.pyplot as plt

First, we need to perform some setup to make sure W7-X data are available.

Note: The W7-X geometry is currently protected by the W7-X data access agreement. Therefore, the file referenced below is available to W7-X team members (including external collaborators who signed this agreement) upon request.

field = w7x.standard()

grid = w7x.defaultGrid()
grid.nR = 32
grid.nZ = 32
grid.nPhi = 32

field = field.compute(grid)
await field

Finally, we need to decide in which phi planes we want to evaluate our iota and start our surfaces from.

#rStart = np.concatenate([np.linspace(6.1, 6.18, 5), np.linspace(6.225, 6.245, 5)])[None, :]
#phiStart = np.linspace(0, 2 * np.pi, 5, endpoint = False)[:,None]
rStart = np.asarray([6.19] + [6.225, 6.235, 6.245] * 5)
phiStart = np.asarray(
    [0] +
    [0 * 2 * np.pi / 5] * 3 + [1 * 2 * np.pi / 5] * 3 + [2 * 2 * np.pi / 5] * 3 + [3 * 2 * np.pi / 5] * 3 + [4 * 2 * np.pi / 5] * 3
)


xStart = rStart * np.cos(phiStart)
yStart = rStart * np.sin(phiStart)
zStart = 0 * xStart

Now it’s time to run our calculation.

pos, ax = fsc.flt.findAxis(field, startPoint = [xStart, yStart, zStart], islandM = 5, targetError = 1e-4, minStepSize = 1e-3, nTurns = 50)

iotas = fsc.flt.calculateIota(
    field, [xStart, yStart, zStart],
    200, # Turn count
    unwrapEvery = 10, distanceLimit = 1e4, 
    targetError = 1e-4, minStepSize = 1e-3,
    islandM = 5, axis = ax
)

The iota profile is returned in a shape matching the start point shape, and can be easily plotted.

print(iotas.shape, rStart.shape)

(16,) (16,)

plt.scatter(rStart + 0 * phiStart, iotas)

<matplotlib.collections.PathCollection at 0x2631eaadfc0>

Additionally, the field line tracer can also extract the Fourier decomposition of the magnetic surfaces from the field line.

modes = fsc.flt.calculateFourierModes(
    field, [xStart, yStart, zStart],
    30, # Turn count
    nMax = 15, mMax = 5, toroidalSymmetry = 5,
    unwrapEvery = 10, recordEvery = 1,
    targetError = 1e-4, distanceLimit = 1e5, maxStepSize = 0.01,
    stellaratorSymmetric = False, aliasThreshold = 0.05, # The island surfaces are not Stellarator symmetric
    islandM = 5
)

The modes are returned with a dict containing rotational transform and the Fourier expansion of the modes, as well as an encoding of the result as a FourierSurfaces object that can also be used in VMEC inputs.

modes.keys()

dict_keys(['surfaces', 'iota', 'theta', 'rCos', 'zSin', 'mPol', 'nTor', 'rSin', 'zCos'])

Of particular interest are two components: The first one is the ‘iota’ array, which returns the rotational transforms.

modes["iota"]

array([0.94893875, 1.08870193, 1.09037153, 1.09210265, 1.08833397,
       1.09001014, 1.09174102, 1.087973  , 1.08964504, 1.09137735,
       1.0876151 , 1.08928435, 1.09101242, 1.0872471 , 1.08891643,
       1.09065221])

Secondarily relevant is the ‘surfaces’ element, which contains the magnetic surfaces. This is an instance of the class fusionsc.magnetics.SurfaceArray, which can be sliced, added, and multiplied similar to a regular NumPy array.

modes["theta"]

array([-0.01311462,  3.08886114,  3.08667328, -0.01027982, -1.86227137,
       -1.8612186 ,  1.32774406, -0.53052922, -0.52729869,  2.66455488,
        0.80054195,  0.80604205, -2.28296756,  2.1296124 ,  2.13764443,
       -0.94786024])

To plot the modes, we need to multiply the Fourier coefficients with the appropriate angles.

surfaces = modes["surfaces"]

def plotCut(phi):
    plt.figure()
    surfaces.asGeometry().plotCut(np.radians(phi))
    plt.axis('equal')
    plt.xlim(5, 6.5)
    plt.ylim(-1.2, 1.2)

plotCut(0)
plotCut(18)
plotCut(36)

surf = surfaces[0]

fieldModes = field.calculateRadialModes(surf, field, nSym = 5, mMax = 15, nMax = 15, nTheta = 400, nPhi = 400)
cc = fieldModes["cosCoeffs"]
sc = fieldModes["sinCoeffs"]

m = fieldModes["mPol"]
n = fieldModes["nTor"]

phi = fieldModes["phi"]
theta = fieldModes["theta"]

tot = np.sqrt(cc**2 + sc**2)

plt.figure(figsize=(16,9))
plt.imshow(fieldModes["radialValues"].T, origin='lower', extent = [0, 360, 0, 360])
plt.colorbar()
plt.title("Brad / B")
plt.xlabel('Phi [deg]')
plt.ylabel('Theta [deg]')

plt.figure(figsize=(16,9))
plt.imshow(tot.T, origin = "lower", norm = "log")
#plt.xlabel(f"n[0, {n[1]:.0f} ... {max(n):.0f}, {min(n):.0f}, ..., {n[-1]:.0f}]")
plt.xlabel("n")
plt.ylabel("m")
plt.xticks(range(len(n)), n)
plt.title("Fourier components of Brad / B")
plt.colorbar()
plt.show()

C:\Users\Knieps\Documents\repos\fsc\src\python\fusionsc\_api_markers.py:15: UserWarning: The function fusionsc.magnetics.MagneticConfig.calculateRadialModes is part of the unstable API. It might change or get removed in the near future. While unlikely, it might also not be compatible across client/server versions.
  warnings.warn(f"The function {f.__module__}.{f.__qualname__} is part of the unstable API. It might change or get removed in the near future. While unlikely, it might also not be compatible across client/server versions.")
C:\Users\Knieps\Documents\repos\fsc\src\python\fusionsc\magnetics.py:402: UserWarning: calculateRadialModes can only use the FFT fast path if nPhi == 2 * nMax + 1 and nTheta == 2 * mMax + 1. Other values are not recommended
  warnings.warn("calculateRadialModes can only use the FFT fast path if nPhi == 2 * nMax + 1 and nTheta == 2 * mMax + 1. Other values are not recommended")

The calculated surfaces can also be converted into geometries to plot.

asPv = surfaces[0].asGeometry(nPhi = 360, nTheta = 360).asPyvista()
asPv.plot(jupyter_backend = 'static')