Estimating wall loads

This example shows how to estimate and visualize wall loads.

Accurate modelling of the wall loads is one of the most computationally expensive operations one can use ASCOT5 for because it can require millions of markers. Therefore it is worthwhile to invest in marker generation and pre-selection, but those topics are discussed in a separate example.

In this example we focus just on the wall loads. Therefore we generate neutral markers which obviously are lost instantly.

import numpy as np
import unyt
from a5py import Ascot

a5 = Ascot("ascot.h5", create=True)

a5.data.create_input("bfield analytical iter circular")
a5.data.create_input("plasma flat")
a5.data.create_input("E_TC")
a5.data.create_input("N0_1D")
a5.data.create_input("Boozer")
a5.data.create_input("MHD_STAT")
a5.data.create_input("asigma_loc")

# Neutrals with a random velocity vector
power = 1.0e7 * unyt.W
from a5py.ascot5io.marker import Marker
nmrk = 10**5
mrk = Marker.generate("gc", n=nmrk, species="alpha")
mrk["charge"][:] = 0
mrk["energy"][:] = 3.5e6
mrk["pitch"][:]  = 0.999 - 1.998 * np.random.rand(nmrk,)
mrk["zeta"][:]   = 2*np.pi * np.random.rand(nmrk,)
mrk["r"][:]      = 8.5
mrk["z"][:]      = 0.0
mrk["phi"][:]    = 180
mrk["weight"][:] = (power / (nmrk * mrk["energy"])).to("particles/s")
a5.data.create_input("gc", **mrk)

print("Inputs created")

Inputs created

The other relevant input is the 3D wall as the wall loads can’t be estimated with a 2D wall (whose linear elements have no area!). Experience has shown that while the magnetic field perturbations governs how much and where particles are lost, the wall shape has a huge impact on the values of the wall loads and where exactly the hot spots form. Hence an accurate wall geometry is essential for accurate results.

But since we are simulating neutral particles born on-axis and immune to ionization, we clearly don’t care about accuracy and our 3D wall is just a 2D contour that is revolved toroidally:

from a5py.ascot5io.wall import wall_3D

rad  = 2.0
pol  = np.linspace(0, 2*np.pi, 181)[:-1]
w2d = {"nelements":180,
       "r":7.0 + rad*np.cos(pol), "z":rad*np.sin(pol)}
w3d = wall_3D.convert_wall_2D(180, **w2d)
a5.data.create_input("wall_3D", **w3d, desc="REVOLVED")

'wall_3D_0681514430'

While engineers weep, we set up the simulation options and run the code. In options the only noteworthy parameters are that the wall hit end condition is enabled, and that we use gyro-orbit simulation mode. Guiding-center simulations also produce wall hits but they might underestimate the loads or hit “wrong” spots if the Larmor radius is large.

from a5py.ascot5io.options import Opt

opt = Opt.get_default()
opt.update({
    # Simulation mode
    "SIM_MODE":1, "FIXEDSTEP_USE_USERDEFINED":1, "FIXEDSTEP_USERDEFINED":1e-8,
    # Setting max mileage above slowing-down time is a good safeguard to ensure
    # simulation finishes even with faulty inputs. Same with the CPU time limit.
    "ENDCOND_WALLHIT":1,
    # Physics
    "ENABLE_ORBIT_FOLLOWING":1,
})
a5.data.create_input("opt", **opt, desc="WALLHITS")

'opt_2343101729'

import subprocess
subprocess.run(["./../../build/ascot5_main", "--d=\"GREATESTHITS\""])
print("Simulation completed")

ASCOT5_MAIN
Tag e389e6c
Branch docs

Initialized MPI, rank 0, size 1.

Reading and initializing input.

Input file is ascot.h5.

Reading options input.
Active QID is 2343101729
Options read and initialized.

Reading magnetic field input.
Active QID is 3359284521

Analytical tokamak magnetic field (B_GS)
Psi at magnetic axis (6.618 m, -0.000 m)
-5.410 (evaluated)
-5.410 (given)
Magnetic field on axis:
B_R = 0.000 B_phi = 4.965 B_z = -0.000
Number of toroidal field coils 0
Estimated memory usage 0.0 MB
Magnetic field read and initialized.

Reading electric field input.
Active QID is 3936038155

Trivial Cartesian electric field (E_TC)
E_x = 0.000000e+00, E_y = 0.000000e+00, E_z = 0.000000e+00
Estimated memory usage 0.0 MB
Electric field read and initialized.

Reading plasma input.
Active QID is 1537485711

1D plasma profiles (P_1D)
Min rho = 0.00e+00, Max rho = 1.00e+01, Number of rho grid points = 100, Number of ion species = 1
Species Z/A  charge [e]/mass [amu] Density [m^-3] at Min/Max rho    Temperature [eV] at Min/Max rho
   1  /  1     1  /  1.000             1.00e+21/1.00e+00                1.00e+04/1.00e+04
[electrons]   -1  /  0.001             1.00e+21/1.00e+00                1.00e+04/1.00e+04
Quasi-neutrality is (electron / ion charge density) 1.00
Estimated memory usage 0.0 MB
Plasma data read and initialized.

Reading neutral input.
Active QID is 3556211514

1D neutral density and temperature (N0_1D)
Grid:  nrho =    3   rhomin = 0.000   rhomax = 2.000
 Number of neutral species = 1
Species Z/A   (Maxwellian)
        1/  1 (1)
Estimated memory usage 0.0 MB
Neutral data read and initialized.

Reading wall input.
Active QID is 0681514430

3D wall model (wall_3D)
Number of wall elements 64800
Spanning xmin = -9.100 m, xmax = 9.100 m
         ymin = -9.100 m, ymax = 9.100 m
         zmin = -2.100 m, zmax = 2.100 m
Estimated memory usage 4.4 MB
Wall data read and initialized.

Reading boozer input.
Active QID is 4068103494

Boozer input
psi grid: n =    6 min = 0.000 max = 1.000
thetageo grid: n =   18
thetabzr grid: n =   10
Boozer data read and initialized.

Reading MHD input.
Active QID is 0271665079

MHD (stationary) input
Grid: nrho =    6 rhomin = 0.000 rhomax = 1.000

Modes:
(n,m) = ( 1, 3) Amplitude = 0.1 Frequency =   1 Phase =   0
(n,m) = ( 2, 4) Amplitude =   2 Frequency = 1.5 Phase = 0.785
Estimated memory usage 0.0 MB
MHD data read and initialized.

Reading atomic reaction input.
Active QID is 0025829278

Found data for 1 atomic reactions:
Reaction number / Total number of reactions = 1 / 1
  Reactant species Z_1 / A_1, Z_2 / A_2     = 1 / 1, 1 / 1
  Min/Max energy                            = 1.00e+03 / 1.00e+04
  Min/Max density                           = 1.00e+18 / 1.00e+20
  Min/Max temperature                       = 1.00e+03 / 1.00e+04
  Number of energy grid points              = 3
  Number of density grid points             = 4
  Number of temperature grid points         = 5
Estimated memory usage 0.0 MB
Atomic reaction data read and initialized.

Reading marker input.
Active QID is 1506640806

Loaded 100000 guiding centers.
Marker data read and initialized.

All input read and initialized.

Initializing marker states.
Estimated memory usage 0.8 MB.
Marker states initialized.
Initialized diagnostics, 0.0 MB.

Preparing output.
Note: Output file ascot.h5 is already present.

The qid of this run is 0220891846

Inistate written.
Simulation begins; 4 threads.
Simulation complete.
Simulation finished in 0.496156 s
Endstate written.

Combining and writing diagnostics.

Writing diagnostics output.

Diagnostics output written.
Diagnostics written.

Summary of results:
   100000 markers had end condition Wall collision

          No markers were aborted.

Done.
Simulation completed

Now to visualize wall loads. First we want to print the 0D quantities and then plot the wall load histogram.

a5 = Ascot("ascot.h5") # Refresh data

warea, peak = a5.data.active.getwall_figuresofmerit()
print("Wetted area: " + str(warea))
print("Peak power load: " + str(peak.to("MW/m**2")))
a5.data.active.plotwall_loadvsarea()

Wetted area: 153.29035765717975 m**2
Peak power load: 3.219736892936127 MW/m**2

Note that the histogram shows the wetted area cumulatively; the value on the \(y\)-axis corresponds to the area where the load is at least the amount given by the \(x\)-axis.

Now where on the wall the particles ended and are there any hot-spots? This plot uses the magnetic axis to calculate the poloidal angle, which is why the magnetic field data has to be initialized.

a5.input_init(bfield=True)
a5.data.active.plotwall_torpol()
a5.input_free()

While previous plot is good for giving a sense of how the loads are distributed, it skews the areas. The way to properly visualize the wall loads is with a 3D plot. Plotting in 3D requires Visualization Toolkit (VTK) and this doesn’t work that well in Jupyter notebooks, which is why the command below might give a warning and the figure might appear twice (this issue only affects notebooks and not normal work with ASCOT5).

a5.data.active.plotwall_3dstill(cpos=(0,6.9,0), cfoc=(0,7.0,0), cang=(120,0,-90), data="eload", log=True)

/usr/share/miniconda/envs/ascot-dev/lib/python3.10/site-packages/pyvista/core/pointset.py:1365: PyVistaDeprecationWarning: The current behavior of `pv.PolyData.n_faces` has been deprecated.
                Use `pv.PolyData.n_cells` or `pv.PolyData.n_faces_strict` instead.
                See the documentation in '`pv.PolyData.n_faces` for more information.
  warnings.warn(
/usr/share/miniconda/envs/ascot-dev/lib/python3.10/site-packages/pyvista/jupyter/notebook.py:37: UserWarning: Failed to use notebook backend:

No module named 'trame'

Falling back to a static output.
  warnings.warn(

You can also make an interactive plot with a5.data.active.plotwall_3dinteractive(). The most convenient way to investigate a 3D plot is to use a5gui where one can record the camera position in an interactive plot, and use it in stills.