Heat load calculation (Wendelstein 7-X)

This example shows how to create heat-load distribution plots for Wendelstein 7-X. First, we need to perform some setup to make sure W7-X data are available.

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

Now we need to setup the geometry and magnetic field configuration.

field = w7x.standard()
geometry = w7x.op21Geometry()

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

field = field.compute(grid)
geometry = geometry.index(w7x.defaultGeometryGrid())

await field
await geometry

The tracing calculations underlying the anisotropic diffusion process are very expensive and dominate the computation time. To reduce them, we can pre-compute a global mapping of the field lines, which reduces field line tracing to a series of interpolations. Here, we use the triangular planes of W7-X, where both the plasma core and the islands have favorable aspect ratios, as our remapping planes.

mapping = fsc.flt.computeMapping(
    field,
    np.radians([72]), # Mapping planes
    r = np.linspace(4, 7, 200), # R values for mapping grid
    z = np.linspace(-1.5, 1.5, 200), # Z values for mapping grid
    stepSize = 0.01, # Step size for tracing computation
    distanceLimit = 7 * 2 * np.pi / 8, # Distance limit for tracing computation (no field line needed to be longer than 1/10th of device)
    toroidalSymmetry = 5
)
await mapping

In addition, we also map the geometry into straight field-line coordinates to be able to use extremely long step sizes for the intersection tests.

geoMapping = mapping.mapGeometry(geometry.triangulate(0.03), toroidalSymmetry = 5, nPhi = 32, nU = 100, nV = 100)

From now on, we can accelerate all computations with the mapping by passing it as an additional argument. Before we can run the calculation, we need to determine the seeding points for it. Our strategy here is to find the LCFS, then move a bit inward.

pInside = [5.5, 0, 0]
pOutside = [6.5, 0, 0]

lcfs = fsc.flt.findLCFS(field, geometry, pOutside, pInside, mapping = geoMapping, stepSize = 13)

pStart = lcfs + [-0.0012, 0, 0]
print(lcfs, pStart)

[6.20675202 0.         0.        ] [6.20555202 0.         0.        ]

startPlanes = np.linspace(0 , 2 / 5 * np.pi, 100, endpoint = False)
nTurns = 1000
xSurf, ySurf, zSurf, _, _ = fsc.flt.poincareInPhiPlanes(pStart, field, startPlanes, nTurns, mapping = mapping, stepSize = 1, distanceLimit = 1e5)

Let’s take a look at our starting surface

rSurf = np.sqrt(xSurf**2 + ySurf**2)
plt.scatter(rSurf[0], zSurf[0], s = 0.1)
geometry.plotCut(startPlanes[0], color='k')
plt.axis("equal")
xSurf.shape

(100, 1000)

Since the surface looks OK we can now start to do our tracing on the nTurns * nPlanes = 220 * 100 = 22000 points we obtained. Since we want to do diffusive tracing, we enable the parallel diffusion model. Because we use field line mapping for the parallel tracing calculation, the only limitation for our step size is how much we will tolerate the geometry intersection tests (which are based on line-casts between adjacent field line positions) to deviate from the exact field line.

import time
start = time.time()

traceResult = fsc.flt.trace(
    [xSurf, ySurf, zSurf], field, geometry,
    distanceLimit = 0, # Disable distance limit
    collisionLimit = 1, # Terminate field line on first collision
    meanFreePath = 1, meanFreePathGrowth = 0.1,
    parallelDiffusionCoefficient = 1,
    # parallelConvectionVelocity = 1,
    isotropicDiffusionCoefficient = 1e-7,
    mapping = geoMapping, stepSize = 13 # The step size should be about 1/3rd of machine circumference (limiting is errors in angle unwrapping)
)

end = time.time()
print(end - start)

237.8221664428711

The result of flt.trace is a superset of the Poincare and connection-length functions (in fact these functions just call trace and return part of the result). Of interest for us is the dict entry ‘endPoints’, which contains 4 pieces of information about the end points: Their x, y, and z coordinates as well as the length of the field line (which we don’t care about here). Since we know the final distribution to be Stellarator-symmetric, we symmetrize the heat loads between our divertors to get better statistics.

print(traceResult.keys())
print(traceResult['endPoints'].shape)
xyz = traceResult['endPoints'][0:3]
xyz = fsc.flt.symmetrize(xyz, nSym = 5, stellaratorSymmetric = True)
print(xyz.shape, xyz[0].size)

dict_keys(['endPoints', 'poincareHits', 'stopReasons', 'fieldLines', 'fieldStrengths', 'endTags', 'numSteps', 'responseSize'])
(4, 100, 1000)
(3, 10, 100, 1000) 1000000

To obtain the heat loads from the point cloud, we use the ‘hfcam’ module. To do that, we need to follow the following steps:

• Define a projection view

• Build a camera that contains a rasterized geometry mesh and infos about the rasterization.

• Render the points into the camera image.

camGeo = w7x.divertor('OP21')# + w7x.baffles('OP21')
camProj = fsc.hfcam.toroidalProjection(
    w = 2500, h = 1200,
    #w = 625, h = 300,
    phi = np.radians(7), rTarget = 5.8, zTarget = -0.5,
    verticalInclination = np.radians(30), horizontalInclination = 0,
    distance = 0.5, viewportHeight = 2.5,
    fieldOfView = np.radians(0)
)
cam = fsc.hfcam.make(camProj, camGeo, depthTolerance = 0)

camProj

width: 2500
height: 1200
transform: [-0.058497284834470767, -0.49627307582066105, -0.85957018174986644, -0, 0.47642215278783445, -0.060934671702573738, -0.10554194733138723, -0, 6.6613381477509375e-18, 0.86602540378443837, -0.49999999999999989, -0, 1.6653345369377344e-18, 3.3330127018922191, 5.2729473419497426, 1.25]

matplotlib.rc('font', size = 22)

camData = cam.getData()
detBuf = np.asarray(camData.determinantBuffer)
cam.clear()
cam.addPoints(xyz, r = 0.015)
result = cam.get()

plt.figure(figsize = (24, 14))

# Show mesh in background
plt.imshow(detBuf.T, origin ='lower', cmap = 'gray')

plt.imshow(np.asarray(result.image).T / xyz.size * 3, origin = 'lower', cmap = 'magma', alpha=0.7)
plt.axis('off')
plt.title("Normalized divertor heat loads")

ax = plt.gca()
cax = plt.gcf().add_axes([ax.get_position().x1+0.01,ax.get_position().y0,0.02,ax.get_position().height])

plt.colorbar(label = "Heat-load normalized against total power ($m^{-2}$)", cax=cax)

<matplotlib.colorbar.Colorbar at 0x13cab4b2950>