Post-processing

Summary

• While some postprocessing can be done without ascotpy, the Python interface to ASCOT5, most advanced features require it.

  • ascotpy is enabled when the C library libascot.so has been compiled.

• Input data can be interpolated and visualized using the input_* methods of Ascot.

• The run group, RunGroup, has methods for accessing and visualizing the output.

  • a5.data.active.getstate(*quantities, mode=..., ids=..., endcond=...) fetch marker ini/endstate.

  • a5.data.active.getorbit(*quantities, ids=..., endcond=...) fetch orbit data.

  • a5.data.active.getdist(dist) fetch distribution data and a5.data.active.getmoments(dist) calculates the moments.

  • a5.data.active.getwall evaluates wall loads.

• What output is present depends on what diagnostics were enabled in the simulation options.

• Simulations can be run directly from Python using the simulate_* methods of Ascot.

• Supporting tools BioSaw and AFSI are run via a5.biosaw and a5.afsi objects, respectively, while BBNBI is a separate program that must be compiled and run similarly as ASCOT5.

Simulation output is accessed via RunGroup object, i.e. a5.data.active is by default the object containing the recent simulation results. Its contents depend on what diagnostics were enabled in the simulation options, except marker ini- and endstate (marker phase space position and associated quantities at the beginning and at the end of the simulation) which are always present.

Many post-processing (and also some pre-processing) tools make use of the Python interface to ASCOT5, referred as ascotpy. Using ascotpy requires that libascot.so has been compiled as per the installation instructions.

Note

Most plotting routines accept arguments for Axes where the plot is created (and axes for the colorbar if applicable). This can be used to create figures for publishing as the axes can be customized and reused. The axes argument is optional; a new figure is created if the argument is not supplied.

Using the Python interface

Instances of Ascot automatically establish the interface if libascot.so has been compiled. The interface provides direct access to the C routines and it can be used to interpolate inputs exactly as they are interpolated during the simulation. It can also be used to simulate markers directly from Python.

To use the interface, start by initializing the input data with input_init():

a5 = Ascot("ascot.h5")

# Initialize currently active bfield
a5.input_init(bfield=True)

Note

If the VERBOSE level is not zero, ASCOT5 prints information on the initialized inputs to stdout. Same is true for libascot.so but in the Python interface these messages are supressed. These messages can be made visible again with a5 = Ascot("ascot.h5", mute="no") which can be helpful for investigating when there are problems with inputs.

This function reads the input data from the HDF5 file, and passes it to C-routines where it is initialized and stored in the memory. After initialization, one can use input_eval() to interpolate the data using the very same C-routines which ASCOT5 uses when it traces markers:

import unyt # units are not necessary but otherwise warnings are issued
a5 = Ascot("ascot.h5")
rgrid = np.linspace(4,8,100)*unyt.m
zgrid = np.linspace(-2, 2, 100)*unyt.m
phi   = 0*unyt.deg
t     = 0*unyt.s

# Evaluate rho and psi at given (R,phi,z,t) locations
rho, psi = a5.input_eval(rgrid[0], phi, zgrid[0], t, "rho", "psi")

# Evaluate rho at given (R,phi,z,t) grid (returns 4D matrix)
psirz = a5.input_eval(rgrid, phi, zgrid, t, "rho", grid=True)

Complete list of quantities that can be evaluated can be found with input_eval_list(). Note however, that input_eval() is just one of the many functions that use ascotpy. Once you are done, free the resources with input_free():

a5.input_free()

An exception is raised if one tries to call a routine that expects specific inputs to be initialized. It is in user’s control to decide what input is used and it doesn’t have to be the same as was used in the simulation. Below is an example on different ways to specify what input is initialized:

# Init active input
a5.input_init(bfield=True)

# Init input with a given QID
qid = a5.data.bfield.TAG.get_qid()
a5.input_init(bfield=qid)

# Init all inputs used by a given run
a5.input_init(run=True)

# Init just the bfield of a run with a given QID
qid = a5.data.TAG.get_qid()
a5.input_init(run=qid, bfield=True)

Note

Initializing inputs, especially the magnetic field, can be memory-intensive since this is usually done on a desktop/laptop and not in a dedicated computing node. However, while simulations require magnetic field with fine resolution, using one with a coarser grid, and thus less memory-intensive, for post-processing is acceptable.

Classes and methods for input processing

Ascot	Primary tool for processing ASCOT5 data.
Ascot.input_init	Initialize input so it can be accessed via Python interface.
Ascot.input_free	Free input used by the Python interface.
Ascot.input_eval	Evaluate input quantities at given coordinates.
Ascot.input_plotrz	Plot input quantity on a (R, z) plane at given coordinates.
Ascot.input_eval_list	Return quantities that can be evaluated from inputs.
Ascot.input_rhotheta2rz	Convert (rho, theta, phi) coordinates to (R,z) positions.
Ascot.input_rhovolume	Evaluate flux surface volumes.
Ascot.input_eval_safetyfactor	Evaluate safety factor and associated quantities.
Ascot.input_eval_ripple	Evaluate various ripple quantities.
Ascot.input_plotwallcontour	Plot intersection of the wall and the poloidal plane at the given toroidal angle.
Ascot.input_getplasmaspecies	Get species present in plasma input (electrons first).
Ascot.input_eval_collcoefs	Evaluate Coulomb collision coefficients for a given test particle.
Ascot.input_eval_atomicsigma	Deprecated.

End conditions

Marker simulation ends when one or more of the active end conditions is met. The marker end condition is stored as a bit array, where the active bits indicate what were the end conditions resulting for that marker’s termination. The post-processing routines use more human-readable approach where end conditions are expressed with strings. This string can either be a single end condition, e.g. endcond="WALL", in which case all markers with at least that end condition are selected. Providing multiple end conditions, e.g. endcond="[POLMAX, TORMAX]", selects all markers that have at least all of those end conditions. Finally, selecting markers that don’t have a given end condition can be done with NOT keyword as endcond="NOT WALL".

The following list contains all possible end conditions.

ABORTED	Marker simulation terminated in an error.
NONE	No active end condition meaning the marker hasn't been simulated yet.
TLIM	Simulation time limit or maximum mileage reached.
EMIN	Minimum energy reached.
THERMAL	Local thermal energy reached.
WALL	Marker intersected a wall element.
RHOMIN	Minimum radial coordinate (rho) reached.
RHOMAX	Maximum radial coordinate (rho) reached.
POLMAX	Maximum poloidal turns reached.
TORMAX	Maximum toroidal turns reached.
CPUMAX	Simulation for this marker exceeded the set CPU time.
HYBRID	Marker guiding center simulation terminated and the simulation continue in gyro-orbit.
NEUTRAL	Ion marker neutralized.
IONIZED	Neutral marker ionized.

Warning

If marker was aborted, something went wrong in the simulation. The cause should always be investigated as the underlying error could have also affected the markers that finished normally. Most common causes of aborted markers is that the marker input contains non-physical data, the input data grids don’t cover the whole plasma or that the magnetic field rho is not properly normalized.

To list all end conditions present in the simulation output, and possible errors, use getstate_summary(). The listed errors shows the error message and the line and file in which it originated.

Ini- and endstate

The data in marker ini- and endstate can be accessed with getstate() method which returns the queried quantity as an array sorted by marker ID. The marker state contains both particle and guiding-center phase space coordinates, and the method also automatically evaluates some derived quantities. Whether a given quantity is evaluated in particle or guiding-center phase space is dictated by mode parameter.

Name of the queried quantity must be preceded by “ini” or “end” to indicate whether the quantity is evaluated at the ini- or endstate. Multiple quantities can be evaluated simultaneously and it is more efficient to do so than to call getstate() separately for each quantity.

Full list of available quantities is given by getstate_list(). Note that evaluating some of the quantities requires input data to be initialized.

Here are examples on how to query marker states:

# Final (R,z) coordinates of marker guiding center and particle positions
a5.data.active.getstate("end r", "end z", mode="gc")
a5.data.active.getstate("end r", "end z", mode="prt")

# Initial guiding center psi coordinate (requires magnetic field to be initialized)
a5.data.active.getstate("ini psi")

# Markers can be pruned by their ID and end condition, and the returned quantity
# can be a difference/relative difference between ini- and ednstate.
a5.data.active.getstate("diff ekin", endcond="TLIM", ids=[1, 2])

Marker states can be visualized conveniently with plotstate_scatter() and plotstate_histogram():

# Quering quantities acts in a similar way as in getstate()
a5.data.active.plotstate_scatter("ini x", "end y", "diff z", cc="log reldiff ekin",
                                 xmode="prt", endcond="TLIM", ids=[1, 2])

# 2D scatter plot with a single color
a5.data.active.plotstate_scatter("end x", "end y", cc="black")

# 2D scatter plot with color as a coordinate
a5.data.active.plotstate_scatter("end x", "end y", cc="end ekin")

# 3D scatter plot with color as a coordinate
a5.data.active.plotstate_scatter("end x", "end y", "end z", cc="end ekin")

# 1D stacked histogram where markers with different end conditions are separated
a5.data.active.plotstate_histogram("ini rho")

# 2D histogram
a5.data.active.plotstate_histogram("end rho", "end phimod")

Marker evaluation and plotting

RunGroup.getstate	Evaluate a marker quantity based on its ini/endstate.
RunGroup.getstate_list	List quantities that can be evaluated with getstate().
RunGroup.getstate_markersummary	Return a summary of marker end conditions and errors present in the data.
RunGroup.plotstate_scatter	Make a scatter plot of marker state coordinates.
RunGroup.plotstate_histogram	Make a histogram plot of marker state coordinates.

Marker orbits

The orbit output consists of marker phase-space coordinates collected on several points along the marker trajectory. The data is either i) gathered at given intervals or ii) collected when specific toroidal, poloidal or radial surfaces are crossed. The latter case allows one to plot Poincaré plots from the orbit data.

The orbit output differs from marker state in that only the phase-space coordinates corresponding to the simulation mode are stored, i.e., GO simulations store only particle coordinates and GC simulations only the guiding-center coordinates. Also the returned values are first sorted by marker ID and then by the mileage. Other than that, the orbit output is accessed and plotted in same manner as the state output.

# Return rho value and the change in energy for markers that were not
# lost to the wall. Orbit data can also be parsed by marker ID
rho, dekin = a5.data.active.getorbit("rho", "diff ekin", endcond="NOT WALL", ids=None)

# Plot the trajectories of all markers in 2D+1D (color)
a5.data.active.plotorbit_trajectory("r", "z", c="pitch")

Marker evaluation and plotting

RunGroup.getorbit	Return orbit data.
RunGroup.getorbit_list	List quantities that can be evaluated with getorbit().
RunGroup.plotorbit_trajectory	Plot orbit trajectories in arbitrary coordinates.
RunGroup.plotorbit_poincare	Create a Poincaré plot where the color separates different markers or shows the connection length.

Distributions

ASCOT5 can be used to collect N-dimensional histograms that represent distribution of the simulated particle population in phase-space. There are in total of five distributions that can be collected:

• 5d Distribution in (R, phi, z, ppar, pperp).

• 6d Distribution in (R, phi, z, pR, pphi, pz).

• rho5d Distribution in (rho, theta, phi, ppar, pperp).

• rho6d Distribution in (rho, theta, pR, phi, pz).

• com Distribution in (mu, ekin, ptor), i.e. in constants of motion space (first adiabatic invariant, energy, canonical toroidal angluar momentum).

While all distributions can be collected irrespective of the simulation mode, there are few considerations. If you are only interested in the distribution function itself, use either 5d/rho5d or 6d/rho6d depending on whether the simulation is in guiding -center (5D) or gyro-orbit (6D) mode and in what real-space basis you need the distribution in. com is for special cases. However, if you are also interested in moments of the distribution, choose one of the 5D distributions: rho5d can be used to calculate 1D profiles of the moments while 5d can be used for (R, z) profiles.

Note

In addition to the phase-space coordinates stated here, all distributions (except com) can also have charge state and time as additional coordinates. Other than that there is no way to separate contributions from different particles. Therefore, it is advisable to run separate simulations for separate sources, e.g. when there are several NBIs, so that their contributions can be separated. Finally, all markers in a simulation should have same mass when distributions are collected as the mass is not a separate coordinate.

Distributions and their moments are represented with DistData and DistMoment objects, respectively. The distribution object is created when getdist() is called, and then the object can be processed (which has no effect on the data on the disk). For 5D distributions it is possible to change the momentum space from parallel and perpendicular momentum to energy and pitch.

Example illustrating how distributions can be processed:

# Return 5d distribution
dist = a5.data.active.getdist("5d")

# Same 5D distribution but where momentum space is converted to energy-pitch
exidist = a5.data.active.getdist("5d", exi=True)

# Distribution object has both methods and attributes that can be used
dist.abscissae           # List containing names of the abscissae
dist.abscissa("r")       # Abscissa bin center values
dist.abscissa_edges("r") # Abscissa bin edges

dist.distribution() # The physical distribution function
dist.histogram()    # Distribution histogram (particles per bin)

The resulting distribution data can be easily sliced, integrated and interpolated using the methods of DistData:

dist = a5.data.active.getdist("5d")

# Integrate over time dimension. Also integrate over the charge dimension
# but use charge as a weight.
dist.integrate(time=np.s_[:], charge=dist.abscissa("charge"))

# Take slices and interpolate. These operations won't reduce dimensionality but
# they reduce the size of the dimensions to one.
dist.slice(ppar=np.s_[0], pperp=np.s_[0])
dist.interpolate(phi=2*unyt.deg)

One can use these methods to process the distribution until it has only one or two non-singleton dimensions left so that it can be visualized:

dist.plot()

Moments are calculated with getdist_moments(). Note that some of the moments require that the input data such as the magnetic field is initialized first.

Example on how to calculate and visualize moments:

# Calculate and plot toroidally averaged density on (R,z) plane
dist = a5.data.active.getdist("5d")
mom = a5.data.active.getdist_moments(dist, "density", "chargedensity")
mom.plot("density")

# Calculate and plot 1D profile of a toroidally and poloidally averaged density
dist = a5.data.active.getdist("rho5d")
mom = a5.data.active.getdist_moments(dist, "density", "chargedensity")
mom.plot("density")

Retrieving distributions and moments

RunGroup.getdist	Return distribution function.
RunGroup.getdist_moments	Calculate moments of distribution.

Distribution and moment data

DistData	Distribution data object.
DistData.abscissa	Return abscissa values.
DistData.abscissa_edges	Return abscissa edges.
DistData.distribution	Return the distribution function.
DistData.histogram	Return the distribution as a histogram (particles per bin).
DistData.phasespacevolume	Calculate phase-space volume of each bin.
DistData.plot	Plot distribution in 1D or 2D.
DistData.integrate	Integrate distribution along the given dimension.
DistData.slice	Take slices of distribution.
DistData.interpolate	Perform (N-dim) linear interpolation on the distribution.
DistMoment	Class that stores moments calculated from a distribution.
DistMoment.ordinate	Return stored moment.
DistMoment.add_ordinates	Add moments.
DistMoment.list_ordinates	List all moments
DistMoment.plot	Plot radial or (R,z) profile of a distribution moment.

Interactive simulations

Simulations can also be run directly from Python via the interface provided by ascotpy. These simulations are equivalent to running ascot5_main with the exception that the output is not stored in ascot.h5. These interactive simulations are intended for simulating a handful of markers at most which is why the distribution output is not implemented. The interface is useful for visualizing marker trajectories or Poincaré plots and evaluating quantities that require markers to be traced only for few orbits.

Before an interactive simulation can be launched, one has to initialize and pack inputs, and set options and markers.

a5 = Ascot("ascot.h5")

# You may use any options and markers but in this example we take them from
# the input file.
mrk = a5.data.marker.active.read()
opt = a5.data.options.active.read()

a5.simulation_initinputs()
a5.simulation_initmarkers(**mrk)
a5.simulation_initoptions(**opt)

The marker and options data are dictionaries with same format as required by their corresponding write_hdf5 routines. The Ascot.simulation_initinputs() method requires that no inputs are initialized before calling it. It differs from the Ascot.input_init() method in that it initializes all inputs required for the simulation (except options and markers) and packs them. Packing means that input data is allocated in a single array which is what ASCOT5 does internally when the data is offloaded for simulation. Once the inputs are packed, the input data cannot be deallocated or altered but otherwise it is possible to use methods like Ascot.input_eval(), that evaluate and plot input data, normally.

Running a simulation returns a VirtualRun object that acts similarly as an ordinary run object RunGroup would.

# Running simulations returns a virtual run that acts the same way as normal runs
vrun = a5.simulation_run()
r, z = vrun.getorbit("r", "z")

Since the results are now kept in the memory and not stored in the disk, rerunning a simulation requires that the previous results are free’d first. Note that the virtual run object cannot be used once its results have been free’d.

# Options can be changed when results are freed
a5.simulation_free(diagnostics=True)
a5.simulation_initoptions(**opt)
vrun = a5.simulation_run()

# Simulations can be rerun with new markers if markers are freed as well
a5.simulation_free(markers=True, diagnostics=True)
a5.simulation_initmarkers(**mrk)
vrun = a5.simulation_run()

# Free all resources
a5.simulation_free()

Tools for the interactive simulations

Ascot.simulation_initinputs	Prepare input fields for the interactive simulation.
Ascot.simulation_initmarkers	Create markers for the interactive simulations.
Ascot.simulation_initoptions	Set simulation options for the interactive simulation.
Ascot.simulation_run	Run the interactive simulation using inputs, options and markers that were set.
Ascot.simulation_free	Free resources used by the interactive simulation.
VirtualRun	Virtual RunGroup whose data exists solely in the memory.

Marker generation

Markers can be generated in three ways. The first way is to create markers explicitly by manually filling the marker dictionary. Here generate() can be used to create an dictionary prefilled with the data of a given marker species:

from a5py.ascot5io.marker import Marker

# Returns a dictionary where charge, mass, anum and znum are pre-filled (other
# quantities are initialized as zeros)
mrk = Marker.generate("gc", n=100, species="alpha")
mrk["energy"][:] = 3.5e6
mrk["pitch"][:]  = 0.99 - 1.98 * np.random.rand(100,)
mrk["r"][:]      = np.linspace(6.2, 8.2, 100)

The second way is to initialize markers from the end state of an existing simulation. This is usually done if the markers hit CPU time limit or if the initial simulation traced markers to separatrix in guiding center mode and the follow-up simulation uses gyro-orbit mode to trace markers until they hit the wall. Use method getstate_markers() to create marker input from the endstate:

a5 = Ascot("ascot.h5")
# Get IDs of markers to be used in the input
ids = a5.data.MYRUN.getstate("ids", endcond="RHOMAX")
mrk = a5.data.MYRUN.getstate_markers("gc", ids=ids)

In actual simulations markers are usually sampled from a source distribution e.g. fusion birth profile provided by AFSI. For this purpose MarkerGenerator can be used to sample markers from any 5D distribution represented by DistData object. Instance of this class can be found in markergen attribute in Ascot:

alphadist = a5.data.AFSIRUN.getdist("prod1")

# "Extra" dimensions (time and charge) must be removed first
alphadist.integrate(time=np.s_[:], charge=np.s_[:])

# Generate markers
nmrk   = 10**6
anum   = 4
znum   = 2
mass   = 4.014*unyt.amu
charge = 2.0*unyt.e
mrk = a5.markergen.generate(nmrk, mass, charge, anum, znum, alphadist)

However, generating markers like this means that all markers have same weights. Marker re-weighting and importance sampling should be considered always to increase the computational efficiency in simulation. For example, when simulating alpha particle losses there is little point in simulating well-confined alpha particles in the core. Instead, one should initialize more markers on the edge and re-weight them so that they still accurately represent the total alpha particle distribution. MarkerGenerator can do this for you if one provides it with and additional 5D distribution that represents the (probability) distribution from which markers are sampled.

# Construct a 5D distribution that is uniform
markerdist = alphadist._copy()
markerdist._distribution[:] = 1

# Generate markers
nmrk   = 10**6
anum   = 4
znum   = 2
mass   = 4.014*unyt.amu
charge = 2.0*unyt.e
mrk, mrkdist, prtdist = a5.markergen.generate(
    nmrk, mass, charge, anum, znum, alphadist, markerdist=markerdist,
   return_dists=True, minweight=1)

Finally, MarkerGenerator has tools to create the marker distribution e.g. from a given \(\rho\) profile.

rho  = np.linspace(0, 1, 100)
prob = np.ones((100,))
prob = (1.0+rho)**3

a5.input_init(bfield=True)
markerdist = a5.markergen.rhoto5d(
    rho, prob, alphadist.abscissa_edges("r"),
   alphadist.abscissa_edges("phi"), alphadist.abscissa_edges("z"),
   alphadist.abscissa_edges("ekin"), alphadist.abscissa_edges("pitch"))
a5.input_free()

MarkerGenerator	Tool for processing distributions into markers.
MarkerGenerator.generate	Generate weighted markers from marker and particle distributions.
MarkerGenerator.rhoto5d	Maps a 1D profile to a normalized 5D dist.

AFSI

AFSI calculates distribution of fusion products from the interaction of two arbitrary populations. It can be used to create either a fast particle source for ASCOT5 simulations or a neutron source e.g. for codes such as Serpent. AFSI has three modes of operation: thermal where two Maxwellian populations are interacting, beam-thermal where one population is Maxwellian and the other is given by arbitrary distribution (e.g. beam ion slowing-down distribution), or beam-beam where both populations are given by arbitrary distributions. The distributions must have the format of a 5D distribution.

AFSI uses libascot.so to perform the calculations efficiently, but the interface is completely in Python. Each AFSI5 run is stored in HDF5 and its Python counterpart, , operates similarly to . The main difference is that the group contains only distribution data: prod1 and prod2 which are DistData objects representing the fusion products.

Running AFSI requires that the inputs contain magnetic field and plasma data. AFSI can automatically detect the reactants from the plasma input when running thermal or beam-thermal simulations.

To run AFSI, use Ascot attribute afsi which is an instance of Afsi. Here’s a demonstration of running AFSI withh different modes:

a5 = Ascot("ascot.h5")
# Lists all available reactions
a5.afsi.reactions()

# Thermal mode requires defining the abscissae for the product distributions.
# For the other modes the output distributions have the same size and dimensions
# as the beam input distribution.
a5.afsi.thermal(
   "DT_He4n",
   rmin, rmax, nr, zmin, zmax, nz,
   minphi=0, maxphi=2*np.pi, nphi=1, nmc=1000,
   minppara=-1.3e-19, maxppara=1.3e-19, nppara=80,
   minpperp=0, maxpperp=1.3e-19, npperp=40)

# In beam-thermal, the beam is the first reactant unless swap=True
a5.afsi.beamthermal("DD_Tp", beamdist, swap=False)

# Beam-beam can either be two separate beams or the beam interacting with itself
a5.afsi.beambeam("DD_Tp", beamdist, beam2=None)

# Results are accessed as usual
dist = a5.data.active.getdist("prod1")

Afsi	ASCOT Fusion Source Integrator AFSI.
Afsi.thermal	Calculate thermonuclear fusion between two thermal (Maxwellian) species.
Afsi.beamthermal	Calculate beam-thermal fusion.
Afsi.beambeam	Calculate beam-beam fusion.
Afsi.reactions	Return reaction data for a given reaction.

BBNBI

BBNBI implements a beamlet based model for generating neutrals from an injector geometry. The code follows generated neutrals until they are either ionized or they hit the wall (and become shinethrough).

BBNBI is a separate program that has to be compiled separately: make bbnbi5 in the main folder. This creates binary ./build/bbnbi5 that uses the same input file and shares some inputs with the main program. A bbnbi5 run requires following inputs: nbi, bfield, plasma, wall, and options. Options are only used partially as only the distribution options are used to bin ionized markers to create particle source distribution. The nbi input refers to NBI which consists of a single or a bundle of injectors represented by Injector.

Note

It is important to note that there are now way to separate results by injector if the injectors were bundled. Only bundle injectors if they have identical geometries or otherwise there is no reason to separate the results.

For testing purposes one can use to create a simple injector and there are scripts to generate injectors for various machines, so generating the injector data explicitly is rarely necessary. To generate injectors explicitly, refer to Injector for the required data. In short, one has to provide the location of the beamlets, their direction and beam divergence, injected species, the energy fractions and the beam power.

Note

Scripts to generate injector geometries exist for ITER, JET, ASDEX Upgrade, MAST-U, JT-60SA, NSTX, and TCV. These are stored in a private repository so contact the ASCOT team in order to gain access.

Once all the required inputs are present, BBNBI is run in a same way as ASCOT5 with the exception that BBNBI is not MPI parallelized:

./bbnbi5 --in=ascot.h5 --d="MYTAG My description" --n=100000 --t0=0.0 --t1=1.0

Here the command line arguments specify the number of markers to be injected, n, and the time-interval [t0,t1] when the beam is enabled. The time interval is only relevant for time-dependent simulations and all it does is to randomly pick a birth time for an injected marker from the interval.

Each BBNBI run is stored in HDF5 and the data is accessed via BBNBIGroup which in many ways is similar to RunGroup. Instead of storing separate ini- and endstates, only the final state of the markers is stored and it can be used to evaluate wall loads in post-processing or to initialize markers for the follow-up ASCOT5 simulation.

# File which has BBNBI5 run tagged "BBNBI"
a5 = Ascot("ascot.h5")
# List all inputs and outputs
a5.data.BBNBI.ls()
# Marker summary and data (note that there is only one "state")
a5.data.BBNBI.getstate_markersummary()
r, z = a5.data.BBNBI.getstate("r", "z")
# Shinethrough
a5.input_init(bfield=True)
a5.data.active.plotwall_torpol()
a5.input_free()
# Create marker input from ionized markers
ids = a5.data.active.getstate("ids", endcond="IONIZED")
mrk = a5.data.active.getstate_markers("gc", ids=ids)
a5.data.create_input("gc", **mrk)
# Read the particle source distribution
dist = a5.data.getdist("5d")

BBNBIGroup	Node containing BBNBI results and methods to process them.
BBNBIGroup.getstate	Evaluate a marker quantity based on its state.
BBNBIGroup.getstate_markers	Convert state to marker input.
BBNBIGroup.getstate_markersummary	Return a summary of marker end conditions and errors present in the data.
BBNBIGroup.getdist	Return distribution function.
BBNBIGroup.getwall_figuresofmerit	Get peak power load and other 0D quantities related to wall loads.
BBNBIGroup.getwall_loads	Get wall loads and associated quantities.
BBNBIGroup.getwall_3dmesh	Return 3D mesh representation of 3D wall and associated loads.
BBNBIGroup.plotwall_3dstill	Take a still shot of the mesh and display it using matplotlib backend.
BBNBIGroup.plotwall_3dinteractive	Open vtk window to display interactive view of the wall mesh.
BBNBIGroup.plotwall_torpol	Plot the toroidal-poloidal projection of the wall loads.

BioSaw

BioSaw is a solver using Biot-Savart law to conveniently include perturbation from external coils to ASCOT5 simulation. It computes the magnetic field from a coil geometry and turns it into a B_3DS input. The produced data is divergence-free although the level of divergence in simulation, when the data is interpolated with splines, depends on the grid resolution. BioSaw uses libascot.so to perform the calculations efficiently, but the interface is completely in Python.

To run BioSaw, you need a file that has either B_2DS or B_3DS magnetic field present or set as active. One starts by defining the coil geometry:

rad   = 3.0
ang   = np.linspace(0, 2*np.pi, 40)
coilx = rad * np.cos(ang) + 6.2
coily = np.zeros(ang.shape)
coilz = rad * np.sin(ang)
coilxyz = np.array([coilx, coily, coilz])

Note

The coil geometry does not have to form a closed loop as BioSaw simply calculates the field arising from each line segment.

BioSaw is represented by class BioSaw whose instance exists as an attribute named biosaw in Ascot. The method BioSaw.calculate() can be used to calculate the magnetic field explicitly, but the methods BioSaw.addto2d() and BioSaw.addto3d() are more convenient to use when one wants to include the perturbation to an existing 2D or 3D field.

To include perturbation from a coil to field that is already 3D:

# Current is in Ampéres
b3d = a5.biosaw.addto3d(coilxyz.T, current=30e4, revolve=None)

# Multiple coils can be given as lists
b3d = a5.biosaw.addto3d([coilxyz.T, ...], current=[30e4, ...], revolve=None)

In 2D, the most common application is the inclusion of TF ripple (which overrides the 2D toroidal field already present in the data). Here it is not necessary to know the coil current, as the current in TF coils can be scaled so that the average toroidal magnetic field on axis has the given value. Furthermore when multiple identical coils has a toroidally symmetrical arrangement, the field can be calculated with just a single coil and the total field computed by revolving and summing this result. This is achieved with revolve which defines the number of toroidal grid points (not the toroidal angle) the field is revolved in each rotation until we have gone whole turn. So if there are e.g. 18 TF coils, the number of grid points must be divisible by 18 and revolve is the number of grid points divided by the number of coils.

# Since the original field is 2D, one can now define the toroidal resolution for the 3D field
b3d = a5.biosaw.addto2d(
      coilxyz.T, phimin=0, phimax=360, nphi=180,
      revolve=10, b0=5.3)

# Write the data to HDF5
a5.data.create_input("B_3DS", **b3d, desc="TFRIPPLE")

BioSaw	Biot-Savart solver for computing magnetic field from a coil geometry.
BioSaw.addto2d	Generate 3D tokamak magnetic field input from coil geometry and
BioSaw.addto3d	Generate new 3D tokamak magnetic field input where contribution from
BioSaw.calculate	Calculate magnetic field at given points due to given coil(s).