Welcome to SOFTv2’s documentation!

Introduction

SOFT is a synthetic radiation diagnostic designed to be used to study synchrotron radiation and bremsstrahlung from runaway electrons in tokamaks. This documentation is for the rewrite of SOFT, called SOFT2 (or SOFTv2).

Improvements in SOFTv2 from the older version include:

• Up to 100x faster simulation times
• Improved models of the radiation
• More robust implementations of key algorithms
• Wide selection of analytical magnetic fields (though support for numerical magnetic fields of course remains)
• CODE/NORSE momentum-space distributions can be taken as input directly
• Ability to “include” configuration files in other configuration files
• Simultaneous support for HDF5 and MAT files

The many new features however also comes at a cost: the input format of the old version of SOFT is no longer supported. All old config files thus have to be converted (manually) to the new format.

SOFTv2 is written entirely in C++.

Compiling

Building SOFT2 is a two-step process. First, you must download and build the SOFT support library, called softlib. After that, you can go ahead and build the actual SOFT2 executable.

Preliminaries

Before building SOFT2 (or softlib), you need to make sure that you have the following software installed on your system:

• CMake version 3.9 or later
• A C++17 compliant compiler and standard library (i.e. gcc >= 7 or equivalent)
• GNU Scientific Library version 2.0 or later
• HDF5

If you intend to compile SOFT2 with MPI support, you also need an MPI library, such as OpenMPI (not to be confused with OpenMP, which is required as part of your compiler).

Building softlib

softlib is a SOFT2 support library that contains various helper classes and routines. You should place the directory containing softlib in the directory that holds the SOFT2 root directory (note: not in the SOFT2 root directory though).

Summary

The whole build process is described in more detail below, but can be summarized using the following chain of commands:

$ cd /path/to/softlib
$ mkdir build
$ cd build
$ cmake ../
$ make

1. Create build directory

Once downloaded, navigate to the softlib directory and create a new directory called build. Next, navigate into the newly created build directory.

2. Run CMake

Now it is time to run cmake:

$ cmake ../

This tells CMake to read the file CMakeLists.txt located in the top directory and execute all commands contained therein. There are a number of options that can be given to CMake that can be used to configure softlib. The default options are the desired options in pretty much all builds.

Each option is passed using the -D flag to CMake. The available options for softlib are

Option	Description	Values	Default
BUILD_TESTS	Build unit tests	ON or OFF	OFF
DEBUG	Build a debug version of softlib	ON or OFF	OFF
OFFICIAL_MATLAB	Link against official MATLAB libraries [1]	ON or OFF	OFF
PRECISION_DOUBLE	Use double precision floating-point numbers [2]	ON or OFF	ON

[1]	For SOFT2, PRECISION_DOUBLE should be turned on. Single-precision will yield very poor results.

[2]	If this option is OFF, MAT files can still be generated using the HDF5 API. If it is ON however, HDF5 can not be generated, and this requires an activated MATLAB installation to be present.

For example, to build softlib with an official MATLAB distribution, run:

$ cmake ../ -DOFFICIAL_MATLAB=ON

instead of the cmake ../ command shown above.

Warning

When you run cmake ../, CMake will output some information about which compiler and libraries it finds on your system. It can be worth to keep an eye on this output and verify that it agrees with your expectations. For example, if there are multiple compilers installed on the system, CMake may choose the wrong version.

Note

To explicitly tell CMake which C/C++ compilers to use, run the following two commands before running CMake:

$ export CC="/path/to/c-compiler"
$ export CXX="/path/to/c++-compiler"

(Note that you may have to clean up the build directory before running cmake again).

3. Run make

In the final step of the build we run make:

$ make

Since the build can take some time, it is recommended to use the -j option of make. Passing this option to make along with a number N will cause make to run N compilations in parallel. For example, if your computer has 4 threads available, you can run:

$ make -j 4

for optimal compilation speed.

Building SOFT2

Once softlib has been built it is time to download and build SOFT2.

Summary

The whole build process is described in more detail below, but can be summarized using the following chain of commands:

$ cd /path/to/SOFT2/build
$ cmake ../
$ make

1. Run CMake

The SOFT2 code already comes with a build/ directory. Navigate to that directory and run:

$ cmake ../

As with softlib, there are a number of options that can be given to cmake to configure the SOFT2 build. The commands are specified to cmake using the -D command-line option. The following options are available:

Option	Description	Values	Default
BUILD_TESTS	Build unit tests	ON or OFF	OFF
COLOR_TERMINAL	Enable colored output	ON or OFF	ON
DEBUG	Include debug symbols in the binary	ON or OFF	OFF
OPTIMIZE_NATIVE	Apply native compiler optimizations	ON or OFF	ON
PROFILING	Compile with profiler flags	ON or OFF	OFF
WITH_MPI	Compile with MPI support	ON or OFF	OFF

For example, to disable colored terminal output (useful if you’re redirecting stdout to a text file for example), run cmake as:

$ cmake ../ -DCOLOR_TERMINAL=OFF

Warning

When you run cmake ../, CMake will output some information about which compiler and libraries it finds on your system. It can be worth to keep an eye on this output and verify that it agrees with your expectations. For example, if there are multiple compilers installed on the system, CMake may choose the wrong version.

Note

To explicitly tell CMake which C/C++ compilers to use, run the following two commands before running CMake:

$ export CC="/path/to/c-compiler"
$ export CXX="/path/to/c++-compiler"

(Note that you may have to clean up the build directory before running cmake again).

2. Run make

In the final step of the build we run make:

$ make

Since the build can take some time, it is recommended to use the -j option of make. Passing this option to make along with a number N will cause make to run N compilations in parallel. For example, if your computer has 4 threads available, you can run:

$ make -j 4

for optimal compilation speed.

Complete Ubuntu example

On a clean Ubuntu 18.10 install, the complete installation procedure would look as follows:

$ apt update
$ apt install build-essential libgsl-dev cmake libhdf5-serial-dev git
$ mkdir SOFT
$ cd SOFT
$ git clone https://github.com/hoppe93/softlib.git
$ git clone https://github.com/hoppe93/SOFT2.git
$ cd softlib
$ mkdir build
$ cd build
$ cmake ../
$ make
$ cd ../SOFT2/build
$ cmake ../
$ make

Configuration Scripts

The SOFT2 configuration scripting language differs in many ways from that of the original version of SOFT. In this chapter you will find information about the basic structure of SOFT2 pi file (aka configuration script) as well as a list of available modules and options that can be configured.

Scripting language

SOFT uses a custom language for configuring runs. The language is used for assigning values to configuration parameters, and as such, the most common construct found in a SOFT configuration file is:

parameterName = value;

This assigns the value value to the parameter named parameterName. Note the semi-colon at the end of the line, which is mandatory.

Parameters exist on many levels—they can apply globally or to specific modules used in the simulation—and are grouped using curly brackets, { and }. Everything between a pair of curly brackets is usually referred to as a block. Assignments that are outside of a block apply globally, i.e. they do not belong to any specific module, but are rather used by all modules.

Module parameters are assigned using the syntax shown below:

@ModuleName blockName (secondary-type) {
    parameterName1 = value1;
    ...
}

First we specify which type of module we’re configuring by stating its name, prefixed with an @ sign: @ModuleName. A module is a component of the simulation that provides specific functionality. Examples include @ParticleGenerator which generates particles for our simulations, @ParticlePusher which simulates guiding-center and particle orbits, and @Radiation which simulates a radiation detector.

The blockName is specific to this block of settings. It is possible to define several configuration blocks for the same module @ModuleName. Which one(s) to use is specified by a parameter in the parent module. For example, the @Equation module (which defines which equations of motion to solve) is used by the @ParticlePusher module, and hence we could define the equation to solve using the following code:

particle_pusher = PPusher;

@ParticlePusher PPusher {
    ...
    equation = gcEquation;
    ...
}
@Equation gcEquation (guiding-center) {
    tolerance = 1e-9;
}
@Equation pEquation (particle) {
    tolerance = 1e-9;
}

Here we define two separate equations but we only use one. This is due to that we assign the configuration block gcEquation to the equation parameter in the PPusher configuration block. Note also the parentheses after the block names in the configuration blocks of the equations. The parentheses surround the “secondary type” of the block, which determines exactly which parameters can be set. In the example above, the secondary types determine whether SOFT follows particles or guiding-centers. The secondary type can also be omitted, but then the name of the configuration block must be that of the secondary type, i.e. in the example above we could define gcEquation as:

@Equation guiding-center {
    tolerance = 1e-9;
}

of course making sure to change the setting equation in PPusher.

Note also that the @ParticlePusher is an example of a module that provides core functionality in SOFT. It is therefore selected using the global parameter particle_pusher, as indicated above.

Comments

SOFT configuration files support single-line comments. They are started with the # symbol and causes SOFT to ignore everything that appears until the next newline character.

Include other files

In SOFT2, it is possible to source other configuration scripts. This is done by placing the name of the file to include between < and >. The included configuration file will be loaded and all its settings applied before proceeding in the parent configuration file. This means that whichever parameter assignment comes last applies. As an example, consider the following two configuration files:

file1:

parameter1 = value1;
parameter2 = value2;

file2:

parameter1 = value3;

<file1>

parameter2 = value4;

After passing file2 to SOFT, i.e. running soft file2, then parameter1 will be value1, whereas parameter2 (which was assigned after including file1) will be value4.

Multiple values and strings

SOFT2 allows multiple values to be assigned to a single parameter. An example of when this is useful is for vector quantities, such as detector position, which needs one value for each component of its position vector. SOFT achieves this by splitting the assigned values at each whitespace and comma. For example, the following lines assigns three different values to the parameter:

parameterName = value1 value2 value3;
parameterName = value1,value2,value3;

(Note: the two lines are shown to show off the possibilities in SOFT; the two lines are equivalent). In some cases it may be desired to assign a value containing a space to parameter though, such as when assigning the name of a file. To overcome this, SOFT allows values to be surrounded by double quotes:

parameterName = "value1 value2 value3";

The double quotes causes parameterName to now instead be assigned a single value which contains spaces. Double quotes can be used everywhere to emphasize what the value being assigned is (even surrounding numbers).

FAQ

Some frequently asked questions:

Are newlines mandatory? No, newline characters are not mandatory in SOFT configuration files. The entire configuration could be written on a single line. They are however very convenient and highly recommended!

Global options

List of global options

distribution_function

Default value:A unit (=1) distribution function Allowed values:Name of distribution function configuration object

Specifies the distribution function to use. Distribution functions are configured separately in their own @DistributionFunction environments. Only one distribution function can be set per simulation.

More details can be found on the page about @DistributionFunction.

include_drifts

Default value:no Allowed values:Boolean

When enabled, keeps terms up to and including order \(\epsilon\) in all expressions. This means that guiding-center orbits drift, and that more accurate radiation models are used in the @Radiation tool.

magnetic_field

Default value:__default__ Allowed values:Name of magnetic field configuration object

Specifies the magnetic field to use. Magnetic fields are configured separately in their own @MagneticField environments. Only one magnetic field can be set per simulation.

If this parameter is left unspecified, then SOFT will try to select a magnetic field automatically. If exactly one magnetic field has been defined in the configuration file, then SOFT will select that magnetic field. Otherwise, if zero or more than one magnetic field has been defined, SOFT will throw an error and exit.

More details can be found on the page about @MagneticField.

mpi_verbose

Default value:no Allowed values:yes and no

If yes, outputs some information about the MPI communication to help debug MPI-related issues. This option is only available when SOFT2 has been compiled with MPI support.

num_threads

Default value:OMP_NUM_THREADS Allowed values:Any positive integer

Sets the number of OpenMP threads to run SOFT on. This value overrides the value specified in the environment variable OMP_NUM_THREADS. If unspecified, this parameter is set OMP_NUM_THREADS, which is typically equal to the number of threads available on your CPU, and will thus allow SOFT to utilize your CPU to the fullest.

particle_generator

Default value:__default__ Allowed values:Name of a particle generator configuration object

Specifes the particle generator (phase-space) to use. Particle generators are configured separately in their own @ParticleGenerator environments. Only one particle generator can be set per simulation.

If this parameter is left unspecified, then SOFT will try to select a particle generator automatically. If exactly one particle generator has been defined in the configuration file, then SOFT will select that particle generator. Otherwise, if zero or more than one particle generator has been defined, SOFT will throw an error and exit.

More details can be found on the page about @ParticleGenerator.

particle_pusher

Default value:__default__ Allowed values:Name of a particle pusher configuration object

Specifies the particle pusher (orbit solver) to use. Particle pushers are configured separately in their own @ParticlePusher environments. Only one particle pusher can be set per simulation.

If this parameter is left unspecified, then SOFT will try to select a particle pusher automatically. If exactly one particle pusher has been defined in the configuration file, then SOFT will select that particle pusher. Otherwise, if zero or more than one particle pusher has been defined, SOFT will throw an error and exit.

More details can be found on the page about @ParticlePusher.

The following options are global, i.e. are set outside of any environment. Click on the options below for further details.

Option	Description
global distribution_function	Specifies configuration of distribution function to use
global include_drifts	Whether or not to include guiding-center drifts
global magnetic_field	Specifies configuration of magnetic field to use
global num_threads	Number of OpenMP threads to run SOFT on
global particle_generator	Specifies particle generator (phase-space) to use
global particle_pusher	Specifies particle pusher (orbit follower) to use

Configurable modules

@Detector

The @Detector module is a submodule of the @Radiation tool. This module configures the radiation detector to use, and sets properties such as its position, viewing direction, field-of-view, spectral range and more.

When specifying detector position, the following illustration of the coordinate system assumed in SOFT can be useful:

Summary of options

Option	Description
aperture	Aperture size / square side length
direction	Viewing direction / detector plane normal vector
optics	Detector optics settings
position	Position relative to tokamak point-of-symmetry
vision_angle	Detector vision angle (field-of-view (half) opening angle)
spectrum	Spectral range configuration of detector
roll	Angle between \(\hat{e}_1\) and the horizontal plane.

Detector plane basis

The detector plane is spanned by two vectors which we refer to as \(\hat{e}_1\) and \(\hat{e}_2\). Together with the detector viewing direction (or normal vector) \(\hat{n}\) they form an orthonormal basis in space. The vectors \(\hat{e}_1\) and \(\hat{e}_2\) are used for several purposes, but one of the more important purposes is for spanning the camera images that this @RadiationOutput produces.

The vector \(\hat{e}_1\) is defined so that it always lies in the horizontal plane when the roll angle is zero. Its mathematical definition is

\[\begin{split}\hat{e}_1 = \begin{cases} \hat{y}\cos\vartheta + \hat{z}\sin\vartheta, \quad&\text{ if } \hat{n}\cdot\hat{y} = 0,\\ \left[ \hat{x}\left(\hat{n}\cdot\hat{y}\right) - \hat{y}\left( \hat{n}\cdot\hat{x} \right)\right]\cos\vartheta + \hat{z}\sin\vartheta, \quad&\text{ otherwise}. \end{cases}\end{split}\]

where \(\vartheta\) denotes the angle between \(\hat{e}_1\) and the horizontal plane (the “roll angle”). From this, \(\hat{e}_2\) is defined so that \((\hat{e}_1, \hat{e}_2, \hat{n})\) form a right-handed orthonormal basis:

\[\hat{e}_2 = \hat{n}\times\hat{e}_1.\]

The values of both of these vectors are output on stdout by SOFT during execution.

Example configuration

This example configuration of a detector corresponds to wide-angle visible camera located in the midplane of a medium-sized tokamak:

@Detector example_detector {
    aperture     = 0.006;              # in meters
    direction    = 0, 1, 0;            # x,y,z
    position     = 0, 1.7, 0;          # x,y,z (relative to point of symmetry)
    vision_angle = 1.25 fov;           # Field-of-view (half) opening angle
    spectrum     = 440e-9, 790e-9, 40; # Lower wavelength (m), Upper (m), Number of points
}

Options

aperture

Default value: None

Allowed values: Any positive real number

Size of detector aperture. All detectors are modeled as squares in SOFT, with the aperture specified here corresponding to the square side length.

direction

Default value: None

Example line: direction = 1.3, -0.25, 0;

Allowed values: Any real 3-vector except null

Detector viewing direction, i.e. normal vector of the detector plane. This vector is normalized internally by SOFT to become a unit vector, and does not have to specified as a unit vector.

optics

Default value: korger

Allowed values: korger

Specifies which model to use for the detector optics. At the moment, the detector optics in principle only specifies whether or not to apply any polarization filters, but could in principle in the future be used to apply any kind of detector transfer function. For now, the default value of this option is the only possible one, and so there is no need to set this option specifically.

position

Default value: None

Example line: position = 0, -1.069, 0;

Allowed values: Any real 3-vector except null

Detector position relative to the tokamak point-of-symmetry. Units used for the vector components are meters.

vision_angle

Default value: None

Example line: vision_angle = 0.75 image;

Allowed values: Real number, optionally followed by either fov or image

Specifies the half opening angle of the field-of-view. If only a number is given, then it is implied that the fov opening angle is specified.

A camera image will be a square inscribed in the circular field-of-view. By appending image after the real number, you indicate that the given vision angle is the minimum angle between a vector extending from the detector to the side of the image, and the detector normal. In a 2D top view, if a cone with the specified image vision angle was plotted, this would correspond exactly to the edges of the observed image.

spectrum

Default value: no

Example line: spectrum = 785e-9, 795e-9, 10;

Allowed values: (i) Vector of two positive real numbers and one positive integer, or (ii) no

Sets the limits and resolution of the detector spectral range. If the number of spectral points is 0, or if this parameter is assigned the boolean no value, the detector will have an infinite spectral range.

The numbers given to the spectral range specify (i) the lower spectral bound, (ii) the upper spectral bound, and (iii) the number of points on the interval. The units of the bounds depend on the emission model used. The synchrotron emission models will assume that the spectrum limits are wavelengths given in units of meters. The bremsstrahlung models assume that the spectrum limits are photon energies, given in normalized units (normalized to \(m_e c^2\), the electron mass times speed-of-light squared).

roll

Default value: 0

Example line: roll = 0.1 cw;

Allowed values: Real number, optionally followed by either cw or ccw

Specifies the angle in radians by which the image x axis, \(\hat{e}_1\) is rotated with respect to the horizontal plane (the x/y plane). By default, this parameter is 0, so that \(\hat{e}_1\cdot\hat{z} = 0\). If no direction is specified, a positive angle corresponds to rotation in the counter-clockwise (CCW) direction.

Optionally, the direction of rotation can be explicitly specified by appending either cw (positive angle in clockwise direction) or ccw (positive angle in counter clockwise direction) after the roll angle.

@DistributionFunction

The distribution function module, which defines a distribution function to weigh results with. A number of different types of distribution functions can be generated or loaded in SOFT2, including the analytical primary distribution [1], the analytical avalanche distribution [2], direct output from the kinetic solvers CODE, LUKE and GO+CODE, as well as regular SOFT phase-space distribution functions.

[1]	(1, 2) Connor and Hastie, 1975 “Relativistic limitations on runaway electrons”. Nuclear Fusion 15 (3), 415 doi:10.1088/0029-5515/15/3/007.

[2]	(1, 2) Fülöp et al., 2006 “Destabilization of magnetosonic-whistler waves by a relativistic runaway beam”. Physics of Plasmas 13 (6), 062506 doi:10.1063/1.2208327.

(avalanche)

In a spatially uniform fully ionized plasma with constant electric field, when the runaway generation is dominated by the avalanche mechanism—i.e. by multiplication through large-angle collisions—an analytical expression can be found for the resulting runaway electron distribution function. This distribution function, which was first derived in [1], takes the form

\[f(p, \xi) = \frac{A(p)}{2\pi m_e c\gamma_0 p^2} \frac{\exp\left[ -\gamma/\gamma_0 - A(p)(1+\xi) \right]}{1 - \exp\left[-2A(p)\right]}\]

with

\[\begin{split}A(p) &= \frac{\hat{E} + 1}{Z_{\rm eff} + 1} \gamma,\\ \gamma_0 &= \log\Lambda\sqrt{Z_{\rm eff}+5}\end{split}\]

where \(p\) is the electron momentum, \(\xi = \cos\theta_{\rm p}\) the electron pitch with respect to the magnetic field, \(\gamma = \sqrt{1 + p^2}\), \(\hat{E}\) the electric field strength normalized to the threshold electric field, \(Z_{\rm eff}\) the plasma effective charge, \(\log\Lambda\) the plasma Coloumb logarithm, \(m_e\) the electron rest mass and \(c\) the speed of light in vacuum.

[1]	Fülöp et al., 2006 “Destabilization of magnetosonic-whistler waves by a relativistic runaway beam”. Physics of Plasmas 13 (6), 062506 doi:10.1063/1.2208327.

Summary of options

The following parameters can be set on an avalanche distribution function.

Option	Description
avalanche EHat	Electric field strength, normalized to the threshold electric field.
avalanche lnLambda	Coloumb logarithm of the plasma.
avalanche Zeff	Plasma effective charge.
avalanche radprof	Name of configuration block for radial profile.

Example configuration

This example illustrates how an analytical avalanche distribution function can be defined and used in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge
}

All options

EHat

Default value:None Allowed values:Any real number

Electric field strength, normalized to the Connor-Hastie critical electric field.

lnLambda

Default value:None Allowed values:Any real numer

Value of the Coulomb logarithm.

radprof

Default value:Uniform radial profile Allowed values:Name of any defined @RadialProfile

Specifies the radial profile object to use to generate a radial profile.

Zeff

Default value:None Allowed values:Any real number

Effective plasma charge.

(code)

The numerical tool CODE (for COllisional Distribution of Electrons) [1] [2] [3], solves the spatially homogeneous Fokker–Planck equation for electrons. The code includes many effects important to runaway electron dynamics, including electric field acceleration, time-dependent plasma parameters, avalanche multiplication, bremsstrahlung and synchrotron radiation losses, as well as partial plasma ionization. SOFT provides support for loading distribution functions calculated with CODE directly through the present module.

[1]	Landreman et al., 2014 “Numerical calculation of the runaway electron distribution function and associated synchrotron emission”. Computer Physics Communications 185 (3), 847-855 doi:10.1016/j.cpc.2013.12.004.

[2]	Stahl et al., 2016 “Kinetic modelling of runaway electrons in dynamic scenarios”, Nuclear Fusion 56 (11), 112009 doi:10.1088/0029-5515/56/11/112009.

[3]	http://ft.nephy.chalmers.se/retools/code/

Summary of options

Option	Description
code interptype	Interpolation method to use when interpolating in the distribution function.
code name	Name of file containing distribution function.
code radprof	Name of configuration block defining the radial profile to use.
code time	Index of distribution function time slice to use for simulation.

Example configuration

The following example shows how to load a CODE distribution function and use the final time step of the CODE simulation. The configuration also sets a radial profile to use with the CODE momentum-space distribution function. If code radprof is not explicitly set, a uniform radial profile is used:

distribution_function = ourDistribution;

@DistributionFunction ourDistribution (code) {
    name    = "/path/to/code/distribution.mat";
    radprof = ourRadProf;
    time    = -1;         # -1 = last time step
}

@RadialProfile ourRadProf (linear) {}

File layout

The file containing the CODE distribution function must contain the variables listed in the table below. They are all contained in the CODE output struct.

Variable	Description
delta	Normalization of momentum. Defined such that \(p = y\delta\), where \(p\) is momentum normalized to the electron rest mass \(m_e c\).
f	Distribution function matrix. Must be of size \(n_t\times n_y n_\xi\).
Nxi	Number of Legendre modes used in the CODE calculation to resolve the angular dependence.
nref	Reference density (used for un-normalising the distribution function).
Tref	Reference temperature (used for un-normalising the distribution function).
y	Momentum grid, given in thermal units.

All options

interptype

Default value:cspline Allowed values:akima, akima_periodic, cspline, cspline_periodic, linear, polynomial, steffen

Determines which interpolation method to use for interpolating in the momentum dimension.

name

Default value:None Allowed values:String

Specifies the name of the file containing the CODE distribution function to load.

radprof

Default value:Uniform radial profile Allowed values:Name of any defined @RadialProfile

Specifies the radial profile object to use to generate a radial profile.

time

Default value:-1 (last timestep) Allowed values:Any integer with absolute value less than the number of time points in the distribution function

Selects the index of the time step to take the distribution function from. Negative indices count from the back of the array, so that -1 corresponds to the last timestep, -2 to the next-to-last etc.

(connorhastie)

When the runaway generation is dominated by the primary (or Dreicer) mechanism, the distribution function is predicted to take the form given by Connor & Hastie [1]:

\[f(p, \xi) = \frac{1}{p\xi}\exp\left[ -\frac{\hat{E}+1}{2(Z_{\rm eff}+1)}\frac{1-\xi^2}{\xi}p \right],\]

where \(p\) is the electron momentum, \(\xi = \cos\theta_{\rm p}\) the electron pitch with respect to the magnetic field, \(\hat{E}\) the electric field strength normalized to the threshold electric field, and \(Z_{\rm eff}\) the plasma effective charge. The threshold electric field \(E_{\rm c}\) is given by

\[E_{\rm c} = \frac{n_ee^3\ln\Lambda}{4\pi\epsilon_0^2 m_ec^2},\]

where \(n_e\) is the electron density, \(e\) the elementary charge, \(\ln\lambda\) the Coulomb logarithm, \(\epsilon_0\) the permittivity of free space, \(m_e\) the electron mass and \(c\) the speed of light in vacuum.

[1]	Connor and Hastie, 1975 “Relativistic limitations on runaway electrons”. Nuclear Fusion 15 (3), 415 doi:10.1088/0029-5515/15/3/007.

Summary of options

The following parameters can be set on a Connor-Hastie distribution function.

Option	Description
connorhastie EHat	Electric field strength, normalized to the threshold electric field.
connorhastie Zeff	Plasma effective charge.
connorhastie radprof	Name of configuration block for radial profile.

Example configuration

This example illustrates how a Connor-Hastie distribution function can be defined and used in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = connor-hastie;

@DistributionFunction connor-hastie (connorhastie) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    Zeff      = 4;            # Effective plasma charge
}

All options

EHat

Default value:None Allowed values:Any real number

Electric field strength, normalized to the Connor-Hastie critical electric field (see above).

radprof

Default value:Uniform radial profile Allowed values:Name of any defined @RadialProfile

Specifies the radial profile object to use to generate a radial profile.

Zeff

Default value:None Allowed values:Any real number

Effective plasma charge.

(dream)

This module provides a way for the user to load a distribution function from an output file generated by the kinetic solver DREAM. The module is capable of loading both hot electron- and runaway distribution functions.

Summary of options

Option	Description
dream distribution	Name of distribution function to load (f_hot or f_re).
dream flippitchsign	If true, flips the sign of the particle pitch.
dream interptype	Interpolation method to use.
dream logarithmize	Whether or not to interpolate in the logarithm of the distribution function.
dream name	Name of file containing DREAM distribution function.
dream time	Index of time step from distribution function to use.

Example configuration

The following example configures a DREAM distribution function:

@DistributionFunction ourDistributionFunction (dream) {
    name = "/path/to/distribution.h5";
    distribution = "f_hot";
}

All options

distribution

Default value:f_re, if it exists in the output file, otherwise f_hot. Allowed values:f_hot or f_re.

Specifies the name of the distribution function to load. DREAM can simulate the evolution of two connected distribution functions simultaneously and this option specifies which of the two distribution functions to use. Generally, synchrotron radiation originates in the high-energy region of the distribution which is represented by f_re. Hence, f_re is used by default unless this option is explicitly set otherwise. Some simulations keep the entire distribution function on a single grid, in which case SOFT will pick the only available distribution function.

flippitchsign

Default value:no Allowed values:yes or no

If yes, transforms \(f(p,\xi)\to f(p,-\xi)\). This is useful in case the distribution function is defined in an inconsistent way (i.e. if a positive parallel electric field accelerates particles in the anti-parallel direction).

interptype

Default value:cubic Allowed values:cubic or linear

SOFT interpolates in the given distribution function to evaluate it at arbitrary points on the phase space grid. A linear interpolation scheme is always used to interpolate in the radial coordinate, but interpolation in the momentum coordinates (\(p\) and \(\xi\)) can either be done using bi-linear or bi-cubic splines.

logarithmize

Default value:no Allowed values:yes or no

If yes, interpolates in the logarithm of the distribution function instead of in the distribution function directly. This can aid in fitting sharply decaying ditsribution functions.

name

Default value:Nothing Allowed values:Any valid file name

Name of the file containing the distribution function.

time

Default value:-1 (last timestep) Allowed values:Any integer with absolute value less than the number of time points in the distribution function

Selects the index of the time step to take the distribution function from. Negative indices count from the back of the array, so that -1 corresponds to the last timestep, -2 to the next-to-last etc.

(gocode)

GO+CODE refers to the coupled fluid-kinetic framework consisting of the fluid code GO [1] [2], which is used to solve for the toroidal electric field evolution, and the kinetic solver CODE [3] [4] [5], which solves the spatially homogeneous Fokker–Planck equation for electrons. The framework calculates the electron distribution function with one spatial dimension and two momentum dimensions.

[1]	Smith et al., 2006 “Runaway electrons and the evolution of the plasma current in tokamak disruptions”. Physics of Plasmas 13, 102502 doi:10.1063/1.2358110.

[2]	Papp et al., 2013 “The effect of ITER-like wall on runaway electron generation in JET”. Nuclear Fusion 53 (12), 123017 doi:10.1088/0029-5515/53/12/123017.

[3]	Landreman et al., 2014 “Numerical calculation of the runaway electron distribution function and associated synchrotron emission”. Computer Physics Communications 185 (3), 847-855 doi:10.1016/j.cpc.2013.12.004.

[4]	Stahl et al., 2016 “Kinetic modelling of runaway electrons in dynamic scenarios”, Nuclear Fusion 56 (11), 112009 doi:10.1088/0029-5515/56/11/112009.

[5]	http://ft.nephy.chalmers.se/retools/code/

Summary of options

Option	Description
gocode interptype	Interpolation method to use when interpolating in the distribution function.
gocode name	Name of file containing distribution function.
gocode time	Index of distribution function time slice to use for simulation.

Example configuration

The following example shows how to load a GO+CODE distribution function and use the final time step of the simulation:

distribution_function = ourDistribution;

@DistributionFunction ourDistribution (gocode) {
    name    = "/path/to/gocode/distribution.mat";
    time    = -1;         # -1 = last time step
}

File layout

In contrast to pure CODE and LUKE distribution functions, a GO+CODE cannot be directly saved to a MAT file. Instead, the output must be converted into a special format, consisting of two levels of structs.

In the root path of the file, the following variables must be set:

Variable	Description
r	Radial grid on which the distribution is defined (corresponding to the dimensionful r variable in GO).
t	List of times at which the distribution is available.
nt	Number of elements in t.
tX	Struct, containing the distribution functions at time index X (i.e. X is an integer, ranging from 0 to nt-1.

Each struct tX (t0, t1 etc.) in turn contains several variables rY, with Y an integer ranging from 0 to nr-1, where nr is the number of elements in r. Finally, each rY contains a momentum distribution which has the format described on the page for the CODE distribution function: (code). We repeat it here for convenience:

Variable	Description
delta	Normalization of momentum. Defined such that \(p = y\delta\), where \(p\) is momentum normalized to the electron rest mass \(m_e c\).
f	Distribution function matrix. Must be of size \(n_t\times n_y n_\xi\).
Nxi	Number of Legendre modes used in the CODE calculation to resolve the angular dependence.
nref	Reference density (used for un-normalising the distribution function).
Tref	Reference temperature (used for un-normalising the distribution function).
y	Momentum grid, given in thermal units.

All options

interptype

Default value:cspline Allowed values:akima, akima_periodic, cspline, cspline_periodic, linear, polynomial, steffen

Determines which interpolation method to use for interpolating in the momentum dimension.

name

Default value:None Allowed values:String

Specifies the name of the file containing the GO+CODE distribution function to load.

time

Default value:-1 (last timestep) Allowed values:Any integer with absolute value less than the number of time points in the distribution function

Selects the index of the time step to take the distribution function from. Negative indices count from the back of the array, so that -1 corresponds to the last timestep, -2 to the next-to-last etc.

(luke)

LUKE solves the bounce-averaged Fokker-Planck equation in one spatial dimension and two momentum dimensions.

Summary of options

Option	Description
luke interptype	Interpolation method to use.
dream logarithmize	Whether or not to interpolate in the logarithm of the distribution function.
luke name	Name of file containing distribution function.

Example configuration

File layout

The file containing the LUKE distribution function must contain the variables listed in the table below:

Variable	Description
betath_ref	Reference normalized thermal speed \(\beta_\mathrm{th,ref}\) (for de-normalization).
f	Distribution function \(f(\psi,\xi,p)\).
mhu	Pitch grid in the variable \(\xi = \cos\theta_{\mathrm{p}}\).
pn	Momentum grid, with momentum normalized to \(\beta_{\mathrm{th,ref}}\).
xrhoG	Radial grid in the variable \(\psi/\psi_a\), i.e. normalized poloidal flux.

All options

interptype

Default value:cubic Allowed values:cubic or linear

SOFT interpolates in the given distribution function to evaluate it at arbitrary points on the phase space grid. A linear interpolation scheme is always used to interpolate in the radial coordinate, but interpolation in the momentum coordinates (\(p\) and \(\xi\)) can either be done using bi-linear or bi-cubic splines.

logarithmize

Default value:no Allowed values:yes or no

If yes, interpolates in the logarithm of the distribution function instead of in the distribution function directly. This can aid in fitting sharply decaying ditsribution functions.

name

Default value:Nothing Allowed values:Any valid file name

Name of the file containing the distribution function.

(numerical)

This module provides a way for the user to specify an arbitrary phase space distribution function \(f(r, p, \xi)\).

Summary of options

Option	Description
numerical flippitchsign	If true, flips the sign of the particle pitch.
numerical interptype	Interpolation method to use.
numerical logarithmize	Whether or not to interpolate in the logarithm of the distribution function.
numerical name	Name of file containing numerical distribution function.

Example configuration

The following example configures a SOFT numerical distribution function:

@DistributionFunction ourDistributionFunction (numerical) {
    name = "/path/to/distribution.mat";
}

File layout

The format for numerical distribution functions in SOFT is such that a set of momentum-space distribution functions are specified at a number of radii. The file therefore contains two levels – a top-level where the radial grid is specified along with the number of radial points nr, and nr groups (HDF5) or structs (MAT), each containing the momentum-space distribution function corresponding to a particular point on the radial grid.

The following variables must/may be present in the top-level of the distribution function file:

Variable	Required	Description
r	Yes	(Dimensionful) minor radius (one-dimensional).
type	No	Type of distribution function (must be distribution/soft)

Additionally, on the top-level, nr (number of points in r) number of groups/structs named rX must be present, each containing the variables:

Variable	Required	Description
f	Yes	Distribution function \(f(r, p, \xi)\) of size \(n_\xi \times n_p\).
p	Yes	Momentum grid vector (one-dimensional).
xi	Yes	Pitch (cosine of pitch angle) grid vector (one-dimensional).

Note

In f, the “inner” dimension (i.e. last index in C/C++/Python, first index in Matlab/Fortran) must be p. In memory, the layout of the array should be:

f(p0,xi0)  f(p1,xi0)  f(p2,xi0) ... f(pN,xi0)  f(p0,xi1)  f(p1,xi1) ...

All options

flippitchsign

Default value:no Allowed values:yes or no

If yes, transforms \(f(p,\xi)\to f(p,-\xi)\). This is useful in case the distribution function is defined in an inconsistent way (i.e. if a positive parallel electric field accelerates particles in the anti-parallel direction).

name

Default value:Nothing Allowed values:Any valid file name

Name of the file containing the distribution function.

interptype

Default value:cubic Allowed values:cubic or linear

SOFT interpolates in the given distribution function to evaluate it at arbitrary points on the phase space grid. A linear interpolation scheme is always used to interpolate in the radial coordinate, but interpolation in the momentum coordinates (\(p\) and \(\xi\)) can either be done using bi-linear or bi-cubic splines.

logarithmize

Default value:no Allowed values:yes or no

If yes, interpolates in the logarithm of the distribution function instead of in the distribution function directly. This can aid in fitting sharply decaying ditsribution functions.

(numericaldepr)

This module provides a way for the user to specify an arbitrary phase space distribution function \(f(r, p, \xi)\).

Warning

This module uses an old distribution function format and is deprecated since December 2019. Users are instead urged to use the new (numerical) distribution function module .

Summary of options

Option	Description
numericaldepr name	Name of file containing numerical distribution function.
numericaldepr interptype	Interpolation method to use.
numericaldepr logarithmize	Whether or not to interpolate in the logarithm of the distribution function.

Example configuration

The following example configures a SOFT numerical distribution function:

@DistributionFunction ourDistributionFunction (numericaldepr) {
    name = "/path/to/distribution.mat";
}

File layout

The following variables must/may be present in the distribution function file.

Variable	Required	Description
f	Yes	Distribution function \(f(r, p, \xi)\) of size \(n_p n_\xi n_r\).
fp0	No	Verification vector; \(f(r_0, p, \xi_0)\).
fr0	No	Verification vector; \(f(r, p_0, \xi_0)\).
fxi0	No	Verification vector; \(f(r_0, p_0, \xi)\).
p	Yes	Momentum grid vector with \(n_p\) elements.
punits	Yes	Unit of the given momentum.
r	Yes	Radial grid vector with \(n_r\) elements.
xi	Yes	Pitch (cosine of pitch angle) grid vector with \(n_\xi\) elements.

All options

name

Default value:Nothing Allowed values:Any valid file name

Name of the file containing the distribution function.

interptype

Default value:cubic Allowed values:cubic or linear

SOFT interpolates in the given distribution function to evaluate it at arbitrary points on the phase space grid. A linear interpolation scheme is always used to interpolate in the radial coordinate, but interpolation in the momentum coordinates (\(p\) and \(\xi\)) can either be done using bi-linear or bi-cubic splines.

logarithmize

Default value:no Allowed values:yes or no

If yes, interpolates in the logarithm of the distribution function instead of in the distribution function directly. This can aid in fitting sharply decaying ditsribution functions.

(pitch)

The pitch distribution function is an exponential in the cosine of the pitch angle:

\[f(\xi) = \exp\left( C\xi \right),\]

where \(\xi = \cos\theta_{\rm p}\) and \(C\) is an input parameter. This is the pitch part of the analytical avalanche distribution function (see (avalanche)).

Summary of options

Option	Description
pitch C	Distribution function exponent.
pitch radprof	Name of configuration block defining the radial profile to use.

Example configuration

The following example defines a pitch distribution function and uses it for the simulation:

distribution_function = ourDistribution;

@DistributionFunction ourDistribution (pitch) {
    C = 100;
}

If a pitch distribution with a varying radial distribution function is desired, this can be specified as follows:

distribution_function = ourDistribution2;

@DistributionFunction ourDistribution2 (pitch) {
    radprof = ourRadialProfile;
}

@RadialProfile ourRadialProfile (type) {
    ...
}

(see the page about @RadialProfile for details about how to configure the radial profile).

All options

The pitch distribution function has no options.

C

Default value:None Allowed values:Any positive real number

Specifies the value of the exponent in the pitch distribution.

radprof

Default value:Uniform radial profile Allowed values:Name of any defined @RadialProfile

Specifies the radial profile object to use to generate a radial profile.

(unit)

The unit distribution function is one (1) everywhere in momentum space. If no radial profile is specified, it is also one at every radius.

Summary of options

Option	Description
unit radprof	Name of configuration block defining the radial profile to use.

Example configuration

The following example defines a unit distribution function and uses it for the simulation:

distribution_function = ourDistribution;

@DistributionFunction ourDistribution (unit) {}

An alternative method to specify a unit distribution function is to set the distribution_function option to the pre-defined __unit_distribution_function__:

distribution_function = __unit_distribution_function__;

If a unit momentum distribution with a varying radial distribution function is desired, this can be specified as follows:

distribution_function = ourDistribution2;

@DistributionFunction ourDistribution2 (unit) {
    radprof = ourRadialProfile;
}

@RadialProfile ourRadialProfile (type) {
    ...
}

(see the page about @RadialProfile for details about how to configure the radial profile).

All options

The unit distribution function has no options.

radprof

Default value:Uniform radial profile Allowed values:Name of any defined @RadialProfile

Specifies the radial profile object to use to generate a radial profile.

Note

When specifying distribution functions directly to SOFT, the radial coordinate of the distribution is interpreted as a drift surface label. This is in contrast to the usual radial coordinate in SOFT which is a flux surface label. This means that minor radius \(r=0\) corresponds to the centre of the drift surface when specifying a distribution function to SOFT, and NOT to the magnetic axis as is otherwise the case.

Note also that the radial coordinate in the (green) RadiationOutput is in fact a flux surface label. Hence, if the same distribution is applied after a simulation, instead of during a simulation (using the global distribution_function option), by multiplying it with a Green’s function, then the radial coordinate of the distribution is interpreted as a flux surface label.

Types

SOFT provides several different types of distribution functions. Which type is configured by a particular block is specified by the “secondary type” of the block, i.e. by giving the type in parentheses after the block name. The available distribution function types are listed in the table below.

Type	Description
(avalanche)	An analytical distribution function based on [2].
(code)	A momentum-space distribution function directly from CODE.
(connorhastie)	An analytical distribution function based on [1].
(dream)	Distribution function from the DREAM 1D2P transport solver.
(gocode)	GO-CODE distribution function.
(luke)	Distribution function from the LUKE 1D2P kinetic solver.
(numerical)	A SOFT numerical distribution function.
(pitch)	An exponential pitch distribution function.
(unit)	A distribution function that is one everywhere.

Example configuration

Please, see the pages for the various distribution function types to view examples of how to configure a distribution function.

@Equation

This module specifies which equations of motion to solve, and the details of how to solve them.

Types

This module has two secondary types. The secondary type, given in parentheses after the configuration block name, specifies which set of equations to solve. Note that the module settings listed under All options are available regardless of the secondary type.

Option	Description
guiding-center	Specifies the guiding-center equations of motion
particle	Specifies the particle (“full-orbit”) equations of motion

Summary of options

Option	Description
@Equation method	Numerical method used to solve equations of motion.
@Equation tolerance	Relative tolerance of ODE solver (if adaptive solver is used)

Example configuration

This example illustrates how to configure the guiding-center equations of motion:

@Equation gc (guiding-center) {
    method    = rkdp45;
    tolerance = 1e-9;
}

All options

method

Default value:rkdp45 Allowed values:rkdp45

Specifies the numerical method to use for solving the ODE that consititutes the guiding-center equations of motion. Currently, only one such numerical method is available, namely an adaptive Runge-Kutta solver of order 4, with an error estimator of order 5. The method uses Dormand-Prince coefficients, and is hence often known as the RKDP45 method [1].

tolerance

Default value:1e-9 Allowed values:\(\epsilon < \text{tolerance} < 1\)

Sets the relative tolerance to use if an adaptive ODE solver is used to solve the equations of motion. The tolerance must be smaller than one, and greater than the machine epsilon \(\epsilon\), which is \(\approx 2\cdot 10^{-16}\) using double-precision floating point numbers and \(\approx 10^{-7}\) using single-precision.

[1]	Runge-Kutta-Dormand-Prince (RKDP45), Chapter XX: Integration of differential equations, Numerical Recipes, 3rd edition.

@MagneticField

The magnetic field module defines a magnetic field and domain (walls) to use in the simulation. For setting up a quick simulation, the analytical magnetic field can be used to generate a magnetic field with circular flux surfaces. For more advanced applications, SOFT2 is able to take numeric 2D magnetic fields as input.

The type of the magnetic field to load is specified as a secondary type of the block (i.e. in parenthesis after the block name). At the moment, two different types of magnetic fields are available: analytical and numeric.

More details about the magnetic field in SOFT2 can be found on the page Magnetic fields.

Relating safety factor and current profile

The safety factor type current has been derived from the current density

\[\boldsymbol{j} = \sigma_I\hat{\boldsymbol{\varphi}}j(r) = \sigma_I\hat{\boldsymbol{\varphi}} j_0\left[ 1 - \left(\frac{r}{r_0}\right)^{q_{a2}} \right],\]

where \(j_0\) is the central current density and \(r_0\) is the plasma radius (see Magnetic fields for details about the notation). The corresponding safety factor is

\[q(r) = \frac{q_0}{1 - \frac{2}{n+2}\left( \frac{r}{r_0} \right)^{q_{a2}}}.\]

From this, the plasma current \(I\) can be related to the central safety factor \(q_0\equiv q(0)\):

\[I = \frac{2\pi B_0}{\mu_0 q_0 R_{\rm m}} r_0^2 \left( 1 - \frac{2}{n+2} \right),\]

or conversely

\[q_0 = \frac{2\pi B_0}{\mu_0 R_{\rm m} I} r_0^2 \left( 1 - \frac{2}{n+2} \right),\]

where \(\mu_0\approx 4\pi\times 10^{-7}\,\text{H/m}\) is the vacuum permeability. With \(B_0\), \(\mu_0\), \(R_{\rm m}\) and \(r_0\) in SI units, and the total plasma current in mega-ampere, the central safety factor can be computed from

\[q_0 = \frac{5n}{n+2} \frac{B_0 r_0^2}{R_{\rm m} I}.\]

Summary of options

Option	Type	Description
analytical B0	analytical	Magnetic field strength (on-axis) of analytical magnetic field
analytical Rm	analytical	Major radius of analytical magnetic field
analytical zm	analytical	Vertical position of magnetic axis
analytical rminor	analytical	Minor radius of analytical magnetic field
analytical sigmaB	analytical	Toroidal field direction of analytical magnetic field
analytical sigmaI	analytical	Current direction of analytical magnetic field
analytical qtype	analytical	Safety factor type: current, constant, linear, qudratic or exponential
analytical qa1	analytical	First safety factor parameter of analytical magnetic field
analytical qa2	analytical	Second safety factor parameter of analytical magnetic field
numeric filename	numeric	Name of file containing numeric magnetic field
numeric filetype	numeric	Override filetype of file containing numeric magnetic field

Example configuration

The following defines an analytical magnetic field, similar to Alcator C-Mod:

magnetic_field = analytical-field;

@MagneticField analytical-field (analytical) {
    B0     = 5.0;
    Rm     = 0.68;
    rminor = 0.22;
    qtype  = quadratic;
    qa1    = 3.0;
    qa2    = 1.0;
}

An example of a numeric magnetic field:

magnetic_field = numeric-field;

@MagneticField numeric-field (numeric) {
    filename = "/path/to/magnetic-field.mat";
    # SOFT will automatically identify the
    # above file as a 'MAT' file, due to its
    # '.mat' filename extension.
}

Common settings

analytical

B0

Default value:None Allowed values:Any positive real number

The magnetic field strength on-axis (in Tesla).

Rm

Default value:None Allowed values:Any positive real number (greater than analytical rminor)

The tokamak major radius (in meters).

zm

Default value:0 Allowed values:Any real number

The vertical position of the magnetic axis.

rminor

Default value:None Allowed values:Any positive real number (less than analytical Rm)

The tokamak minor radius (in meters).

sigmaB

sigmaI

Default value:cw Allowed values:cw / -, or ccw / +

Sign of the toroidal field component (sigmaB) and plasma current (sigmaI). The value cw corresponds to the toroidal component being oriented in the clock-wise direction, when seen from above the tokamak, while ccw corresponds to the toroidal component being oriented in the counter clock-wise direction, when seen from above.

Instead of specifying the direction, the sign may be given directly as either - (clock-wise) or + (counter clock-wise).

qa1

qa2

Default value:1.0 Allowed values:Any real number

Constants used to define the safety factor. See analytical qtype for details about how exactly they affect the different safety factor.

qtype

Default value:None Allowed values:constant, exponential, linear, quadratic

Specifies the radial dependence of the safety factor. The functional forms for the various options are given in terms of the normalized minor radius \(a\) (normalized to the value of rminor) in the table below. The constants \(q_{a1}\) and \(q_{a2}\) are specified using the analytical qa1 and analytical qa2 options.

qtype	Functional form
current	\(q(a) = q_{a1}/\left( 1-\frac{2}{n+2}a^{q_{a2}} \right)\)
constant	\(q(a) = q_{a1}\)
linear	\(q(a) = q_{a1} a + q_{a2}\)
quadratic	\(q(a) = q_{a1} a^2 + q_{a2}\)
exponential	\(q(a) = e^{q_{a1}} + q_{a2}\)

numeric

filename

Default value:None Allowed values:String

Specifies the name of the file that contains the magnetic field to load.

filetype

Default value:Auto-determine filetype based on filename Allowed values:h5, hdf5, mat, out, sdt

Overrides the filetype of the given file. Usually SOFT tries to determine which filetype a given file is of based on its filename extension. By explicitly setting this option, this check is overriden allows you to use non-standard filename extensions for the input file.

@Orbits

The orbits tool.

@ParticleGenerator

The particle generator is in charge of sampling particles from the phase space. Thus, the particle generator must know which type of particle to sample (i.e. electron, proton or some other charged particle?) as well as the bounds of the phase space from which to pick the particle properties (energy, initial pitch angle and initial position).

Summary of options

Option	Description
@ParticleGenerator charge	Charge of particle (in units of elementary charge).
@ParticleGenerator driftshifttol	Relative tolerance when computing the drift orbit shift.
@ParticleGenerator mass	Rest mass of particle (in units of the electron rest mass).
@ParticleGenerator mpi_distribute_mode	Specifies which parameter to split among MPI processes.
@ParticleGenerator position	Whether specified position refers to particle or guiding-center.
@ParticleGenerator progress	Instructs SOFT to output info on the simulation progress.
@ParticleGenerator nzeta	The number of points to use for the gyro phase integral.

In addition to these options, the following phase space parameters can also be specified. Exactly one spatial and two momentum parameters must be specified, and the two momentum parameters chosen must span all of momentum space together.

Parameter	Type	Description
@ParticleGenerator a	Spatial	Normalized minor radius (0 on magnetic axis; 1 on last closed flux surface).
@ParticleGenerator r	Spatial	Minor radius (meters).
@ParticleGenerator rho	Spatial	Major radius (meters).
@ParticleGenerator gamma	Momentum	Particle energy, normalized to the particle rest mass (= Lorentz factor).
@ParticleGenerator p	Momentum	Particle momentum, normalized to the particle rest mass.
@ParticleGenerator ppar	Momentum	Particle momentum parallel to magnetic field, normalized to the particle rest mass.
@ParticleGenerator pperp	Momentum	Particle momentum perpendicular to magnetic field, normalized to the particle rest mass.
@ParticleGenerator thetap	Momentum	Particle pitch angle \(\theta_{\rm p}\) (in radians).
@ParticleGenerator ithetap	Momentum	Pi-complement of particle pitch angle, \(\theta_{\rm p}' = \pi - \theta_{\rm p}\) (in radians).
@ParticleGenerator xi	Momentum	Cosine of particle pitch angle, \(\xi = \cos\theta_{\rm p}\).

Example configuration

Electron — The following example sets up a phase-space for an electron with 100 points on the grid in each dimension. The mass and charge default to those of an electron, and so do not have to be specified. We also instruct SOFT to output a total of 100 progress messages during the run. Since we do not set the meaning of the position explicitly, SOFT assumes that we specify the position of the guiding-center:

@ParticleGenerator PGen_electron {
    r        = 0, 0.30, 100;    # Minor radius (m)
    gamma    = 1.1, 3.0, 100;   # Energy       (mc^2)
    thetap   = 0.02, 0.8, 100;  # Pitch angle  (rad)
    progress = 100;             # Output 100 progress messages
}

Deuterium — The following example sets up a phase-space for a deuterium ion with 100 points on the grid in each dimension. The proton-electron mass ratio is approximately \(m_{\rm p} / m_{\rm e}\approx 1836\), and hence the deuterium-electron mass ratio is approximately \(m_{\rm D} / m_{\rm e}\approx 3672\). We explicitly state that we specify the particle position using the position option:

@ParticleGenerator PGen_deuterium {
    mass   = 3672;            # Electron masses
    charge = 2;               # Elementary charges
    a      = 0, 1, 100;       # Normalized minor radius
    p      = 1e-3, 1e-1, 100; # Momentum (mc)
    thetap = 0.1, 0.3, 100;   # Pitch angle (rad)
}

With MPI — The following example is near-identical to the electron example above, but explicitly instructs SOFT to split the gamma (energy) parameter among the MPI processes:

@ParticleGenerator PGen_electron {
    r                   = 0, 0.30, 100;    # Minor radius (m)
    gamma               = 1.1, 3.0, 100;   # Energy       (mc^2)
    thetap              = 0.02, 0.8, 100;  # Pitch angle  (rad)
    progress            = 100;             # Output 100 progress messages
    mpi_distribute_mode = gamma;           # Divide the energy parameter among MPI processes
}

Options

charge

Default value:-1 Allowed values:Any non-zero real number.

Charge of particle to simulate. The charge is given in units of the elementary charge so that a value of 1 corresponds to the proton charge and -1 to the electron charge.

driftshifttol

Default value:1e-4 Allowed values:\(\epsilon < \text{tolerance} < 1\)

Tolerance for determining the guiding-center drift orbit shift (which is used to determine where to launch particles and where to sample the distribution function). In general, this parameter does not need to be changed.

mass

Default value:1 Allowed values:Any positive real number.

Rest mass of particle to simulate. The mass is given in units of the electron rest mass. In these units, the proton mass is \(m_{\rm p}\approx 1836.15267389\) [1].

mpi_distribute_mode

Default value:auto. Allowed values:1, 2, auto, radius and all of the phase-space parameters listed under Phase space parameters.

When running in MPI mode (MPI = Message Passing Interface; distributed memory parallelization), this parameter can be used to specify how the phase space should be divided among the MPI processes. In contrast to regular OpenMP parallelization, which builds a queue of points in phase space, MPI requires one phase space parameter to be divided evenly among the processes.

Which parameter to divide is specified by giving the name of the parameter, as listed under Phase space parameters, or by giving one of 1, 2 and radius. The former two cause SOFT to divide the first and second momentum parameter respectively (i.e. the alphabetically first and second momentum parameter), while the latter divides the radial parameter, whichever it may be.

If auto is specified, SOFT2 chooses the phase space parameter with the largest number of grid points. This is the default setting.

position

Default value:guiding-center Allowed values:gc, guiding-center and particle.

Specifies whether the guiding-center or particle is given as input. If, for example, the particle position is specified, but guiding-center orbits are simulated, then the guiding-centers are offset from the given position by one Larmor radius, and vice versa for the opposite case.

progress

Default value:no Allowed values:yes, no or a positive integer.

If yes or a positive integer n, outputs a message reporting the progress of the simulation a total of n times during the run. The reports are split evenly accross the phase space, meaning that if the phase space consists of N total grid points, then SOFT reports progress roughly when the number of processed grid points is a multiple of N / n.

nzeta

Default value:1 Allowed values:Any positive integer.

Determines how many points to use for the integral over gyro angle. For guiding-center calculations, this value should always be set to 1 (the only effect of setting it to some other value for guiding-center calculations is that the simulation becomes slower by the factor specified). This value only makes sense for full-orbit simulations.

[1]	https://en.wikipedia.org/wiki/Proton-to-electron_mass_ratio

Phase space parameters

a

Normalized minor radius — The initial minor radial location of the particle/guiding-center, normalized to the minor radial value of the last closed flux surface. Thus, \(a = 0\) corresponds to the magnetic axis and \(a = 1\) to the maximum radius of the last closed flux surface.

r

Minor radius — The initial minor radial location of the particle/guiding-center, given in meters.

rho

Major radius — The initial major radial location of the particle/guiding-center, given in meters.

gamma

Energy — The energy of the particle/guiding-center, normalized to the particle rest mass \(mc^2\), where \(c\) denotes the speed of light in vacuum. This quantity is also known as the Lorentz factor or relativistic factor, and can also be written \(\gamma = \left( 1 - v^2/c^2 \right)^{-1/2}\), where \(v\) is the speed of the particle.

p

Momentum — The momentum of the particle/guiding-center, normalized to the particle rest mass \(mc\), where \(c\) denotes the speed of light in vacuum. This quantity is related to the particle energy/relativistic factor through \(\gamma^2 = 1 + p^2\).

ppar

Parallel momentum — Momentum of the particle parallel to the magnetic field, normalized to the particle rest mass \(mc\), where \(c\) denotes the speed of light in vacuum.

pperp

Perpendicular momentum — Momentum of the particle perpendicular to the magnetic field, normalized to the particle rest mass \(mc\), where \(c\) denotes the speed of light in vacuum.

thetap

Pitch angle — Angle between the particle velocity vector and the magnetic field vector. Given in radians. The pitch angle ranges between \(0\) and \(\pi\). A value greater than \(\pi/2\) means that the particle is moving antiparallel to the magnetic field.

ithetap

Complementary pitch angle — Same as @ParticleGenerator thetap, except that it is defined as “pi minus @ParticleGenerator thetap”, i.e. \(\theta_{\rm p}' = \pi - \theta_{\rm p}\). This is useful when simulating particles with negative parallel momentum (moving in the antiaparallel direction of the magnetic field), since instead of specifying \(\theta_{\rm p} = 3.14159265359\) we can then set \(\theta_{\rm p}' = 0\).

xi

Cosine of pitch angle — Cosine of the pitch angle \(\theta_{\rm p}\), i.e. \(\xi = \cos\theta_{\rm p}\).

Jacobians

The following is a list of all the Jacobian determinants for transformations from a Cartesian coordinate system \((p_x, p_y, p_z)\) to other coordinate systems \((p_1, p_2, \zeta)\), where \(\zeta\) is the gyro angle.

gamma / ppar — \((\gamma, p_{\parallel})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \gamma\mathrm{d}\gamma\mathrm{d}p_{\parallel}\mathrm{d}\zeta\]

gamma / pperp — Does not contain sufficient information to determine if the guiding-center is travelling in the parallel or anti-parallel direction of the magnetic field.

gamma / thetap — \((\gamma, \theta_{\rm p})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \gamma\sin\theta_{\rm p}\sqrt{\gamma^2-1}\,\mathrm{d}\gamma\mathrm{d}\theta_{\rm p}\mathrm{d}\zeta\]

gamma / xi — \((\gamma, \xi)\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \gamma\sqrt{\gamma^2-1}\,\mathrm{d}\gamma\mathrm{d}\xi\mathrm{d}\zeta\]

p / ppar — \((p, p_{\parallel})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = p\,\mathrm{d}p\mathrm{d}p_{\parallel}\mathrm{d}\zeta\]

p / pperp — Does not contain sufficient information to determine if the guiding-center is travelling in the parallel or anti-parallel direction of the magnetic field.

p / thetap — \((p, \theta_{\rm p})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = p^2\sin\theta_{\rm p}\,\mathrm{d}p\mathrm{d}\theta_{\rm p}\mathrm{d}\zeta\]

p / xi — \((p, \xi)\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = p^2\,\mathrm{d}p\mathrm{d}\theta_{\rm p}\mathrm{d}\zeta\]

ppar / pperp — \((p_{\parallel}, p_{\perp})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = p_\perp\,\mathrm{d}p_{\parallel}\mathrm{d}p_{\perp}\mathrm{d}\zeta\]

ppar / thetap — \((p_{\parallel}, \theta_{\rm p})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \frac{p_\parallel^2\sin\theta_{\rm p}}{\cos^3\theta_{\rm p}}\,\mathrm{d}p_{\parallel}\mathrm{d}\theta_{\rm p}\mathrm{d}\zeta\]

ppar / xi — \((p_{\parallel}, \xi)\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \frac{p_\parallel^2}{\xi^3}\,\mathrm{d}p_{\parallel}\mathrm{d}\xi\mathrm{d}\zeta\]

pperp / thetap — \((p_{\parallel}, \theta_{\rm p})\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \frac{p_\perp^2}{\sin^2\theta_{\rm p}}\,\mathrm{d}p_{\perp}\mathrm{d}\theta_{\rm p}\mathrm{d}\zeta\]

pperp / xi — \((p_\perp, \xi)\)

\[\mathrm{d}p_x\mathrm{d}p_y\mathrm{d}p_z = \frac{p_\perp^2}{\left( 1 - \xi^2 \right)^{3/2}}\,\mathrm{d}p_{\perp}\mathrm{d}\xi\mathrm{d}\zeta\]

@ParticlePusher

The particle pusher generates orbits—both particle and guiding-center orbits—which can be fed to the various tools of SOFT that use them to, for example, compute radiation images. To output the orbits, use the @Orbits tool.

Summary of options

Option	Description
@ParticlePusher equation	Specifies which equation to solve.
@ParticlePusher force_numerical_jacobian	Forces numerical computation of the guiding-center Jacobian.
@ParticlePusher nt	Number of time steps to take.
@ParticlePusher nudgevalue	Override \(\mathrm{d}\rho\) used to compute GC Jacobian numerically.
@ParticlePusher time	Duration to follow orbits for.
@ParticlePusher timeunit	Unit of the time specified to @ParticlePusher time.

Example configuration

The following is the default configuration of the particle pusher object. It causes the pusher to generate guiding-center orbits that are traced during one poloidal transit time. The resulting orbit consists of 1000 time steps:

@ParticlePusher __default__ {
    equation                 = guiding-center;
    timeunit                 = poloidal;
    time                     = 1;
    nt                       = 1e3;
    force_numerical_jacobian = no;
    nudgevalue               = __default__;
}

Options

equation

Default value:guiding-center Allowed values:Name of any @Equation configuration block.

Specifies which equation to use. The two equations guiding-center and particle are pre-defined by this module with default settings. The settings of these equations can be overridden by adding configuration blocks for them.

force_numerical_jacobian

Default value:no Allowed values:A boolean value; yes or no.

Force the guiding-center Jacobian to be computed numerically. If no, an analytical expression will instead be used for the Jacobian. The numerical approach is slower, prone to instabilities and in general discouraged.

nt

Default value:1e3, i.e. 1000 points Allowed values:Any positive integer.

The number of time steps per orbit. The orbit quantities (particle position and momentum) are given in a uniformly distributed set of time points between 0 and the maximum time, set by the @ParticlePusher time parameter.

nudgevalue

Default value:__default__ (see description below) Allowed values:Any real number.

To compute the guiding-center Jacobian numerically, the derivatives of the particle position with respect to the initial radial location must be taken. This parameter is the distance \(\Delta\rho\) by which each orbit is nudged in order to evaluate the derivative using a finite difference method.

time

Default value:1 Allowed values:Any positive real number.

Final time point in which to evaluate orbit.

Note: The @Radiation module expects this parameter to be set to 1, and the @ParticlePusher timeunit parameter to be set to poloidal. Note also that the @Radiation automatically discards the final time point to prevent double-counting.

timeunit

Default value:poloidal, i.e. poloidal transit time Allowed values:poloidal and seconds.

The unit of the @ParticlePusher time parameter. If this parameter is set to poloidal, the @ParticlePusher time gives the number of poloidal transits for which each particle should be followed. If this parameter is set to seconds, @ParticlePusher time gives the number of seconds to follow each orbit.

The option ‘poloidal’ also works for particle/full orbits, even though the concept of poloidal transit time is much less well-defined for in such cases. To overcome this poor definition, SOFT first solves the corresponding guiding-center orbit in order to determine the poloidal transit time, before solving the actual particle orbit.

Note: The @Radiation module expects this parameter to be set to poloidal, and the @ParticlePusher time parameter to be set to 1.

@RadiationModel

This module determines how the angular dependence of the emitted radiation should be treated. We can conveniently express the position of an observer relative to the emitting particle using a spherical coordinate system \((r, \theta, \phi)\). Since the direction in which the observer is located relative to the emitter is fully determined by the two angles \(\theta\) and \(\phi\), and since the amount of radiation emitted towards an observer is independent of how far away the observer is located, we can speak of the angular distribution of emitted radiation.

To not assume anything about the angular distribution of the emitted radiation, the (angdist) should be used. It accounts for the full angular distribution of all types of radiation, and is therefore also somewhat heavier to run.

Directed radiation, such as bremsstrahlung and synchrotron radiation, are almost exactly emitted along the velocity vector of the emitting particle, with an angular spread that is proportional to \(\gamma^{-1}\), where \(\gamma\) denotes the relativistic factor. As an approximation, we may therefore assume that all radiation is emitted exactly along the velocity vector of the particle. We refer to this approximation as the cone approximation, and it is implemented in the (cone) module. It corresponds to assuming

\[\frac{\mathrm{d} P}{\mathrm{d}\Omega} = \frac{P_0}{2\pi}\delta\left( \hat{\boldsymbol{v}}\cdot\hat{\boldsymbol{n}} - 1 \right),\]

where \(\frac{\mathrm{d} P}{\mathrm{d}\Omega}\) is the angular distribution of radiation (see the synthetic radiation diagnostic equation on the page about @Radiation), \(P_0\) is the total amount of emitted radiation (possibly in a limited wavelength range), \(\delta\) is a Dirac delta function, \(\hat{\boldsymbol{v}}\) denotes the direction of motion of the particle and \(\hat{\boldsymbol{n}}\) the unit vector pointing out the direction along which the particle is being observed.

The cone model is significantly faster to run than the (angdist) model, and can produce accurate results if the energies of the emitting particles are sufficiently high.

Finally, we may also assume that the emitted radiation is independent of the two angles \(\theta\) and \(\phi\) and that the radiation is emitted uniformly in all directions. A special radiation model has been implemented for this case, called (isotropic), which is primarily used for benchmarking purposes and producing pretty images.

(angdist)

The angular distribution radiation model takes the full angular distribution of radiation into account. In contrast to the cone model—which assumes that all radiation is emitted exactly along the particle velocity vector—this model allows for arbitrary shapes of the angular distribution of radiation.

Summary of options

The following options are available in the angdist radiation model.

Option	Description
@RadiationModel(angdist) emission	Type of radiation emission to model.
@RadiationModel(angdist) nsamples	Number of points in each dimension of the detector surface.
@RadiationModel(angdist) qrule2d	Quadrature rule to use for integrating over the detector surface.

Example configuration

The following example configures a synchrotron radiation model with the detector surface discretized using a total of 16 points:

@RadiationModel ourModel (angdist) {
    emission = synchrotron;
    nsamples = 4;
}

All options

emission

Default value:None Allowed values:bremsstrahlung or synchrotron

Specifies the type of radiation emission to model. Currently, the two available options are bremsstrahlung and synchrotron. The appropriate formulas to use are chosen automatically, depending on whether spectral dependence is considered, and if guiding-center drifts are included or not.

nsamples

Default value:1 Allowed values:Any positive integer.

Number of points in each dimension to discretize the detector surface with. The total number of points on the detector surface is therefore nsamples^2.

qagslimit

Default value:100 Allowed values:Any positive integer

Number of points in QAGS (quadrature) workspace.

qagstol

Default value:1e-3 Allowed values:Any positive real number

Relative tolerance for the QAGS integration.

qrule2d

Default value:simpson Allowed values:simpson

Specifies the quadrature rule to use for integrating over the detector surface. Currently, only Simpson’s rule has been implemented, with the exception of the special case nsamples = 1, in which case the function is merely evaluated in the detector center-point.

(cone)

The cone radiation model is based on the cone approximation, which makes the assumption that all radiation is emitted exactly along the velocity vector of the particle. This assumption can be very useful for studying the radiation from highly relativistic particles.

Summary of options

Option	Description
@RadiationModel(cone) emission	Type of radiation emission to model.
@RadiationModel(cone) n	(Bremsstrahlung emission only) List of plasma species densities.
@RadiationModel(cone) projection	Which projection method to use for the cone model.
@RadiationModel(cone) Z	(Bremsstrahlung emission only) List of plasma species atomic numbers.
@RadiationModel(cone) Z0	(Bremsstrahlung emission only) List of plasma species net charges.

Example configuration

The following example configures a cone model for bremsstrahlung in a deuterium-tritium-tungsten plasma:

@RadiationModel ourModel (cone) {
    emission = bremsstrahlung;
    Z        = 1, 1, 74;
    Z0       = 1, 1, 72;
    n        = 5e19;
}

All options

emission

Default value:None Allowed values:bremsstrahlung, bremsstrahlung_screened, synchrotron or unit

Specifies the type of radiation emission to model. The four available types of radiation are bremsstrahlung, bremsstrahlung_screened, synchrotron and unit. The first two model bremsstrahlung, the latter taking the effect of screening of the nuclei into account. The third model synchrotron radiation, while unit is a type of radiation that is one everywhere. Note that the radiation is assumed to be emitted exactly along the velocity vector of the particle so that unit is very different from the :ref:module-rm-isotropic` radiation model. The unit emission type rather has more similarities with the bremsstrahlung emission type.

n

Default value:None Allowed values:A list of real values.

Specifies the density to use for each plasma ion species when calculating the emitted electron-ion bremsstrahlung.

projection

Default value:reverse Allowed values:reverse (or original)

The cone model relies on a projection of either the guiding-center cone onto the detector plane, or the detector surface onto the guiding-center velocity orthogonal plane. The original projection does the former, while the reverse projection does the latter.

Note that the original model at the time of writing contains a bug and should be avoided.

Z

Default value:None Allowed values:A list of real values.

Specifies the effective plasma charge to use for each plasma ion species when calculating the emitted electron-ion bremsstrahlung.

Z0

Default value:None Allowed values:A list of real values.

Specifies the plasma net charge to use for each plasma ion species when calculating the emitted electron-ion bremsstrahlung.

(isotropic)

The isotropic radiation model assumes that all particles emit radiation uniformly in all directions. Mathematically, this radiation model takes the form

(1)\[\frac{\mathrm{d} P}{\mathrm{d}\Omega} = P_0,\]

where \(P_0\) is a user-defined constant. Note that the angular distribution of radiation is independent of any particle or background plasma properties.

Summary of options

Option	Description
@RadiationModel(isotropic) value	The value of \(P_0\) in (1).

Example configuration

The following example configures an isotropic model:

@RadiationModel ourModel (isotropic) {
    value = 1;
}

All options

value

Default value:

Allowed values:

Any real value.

Power per unit solid angle emitted by the particle.

Available models

Model	Description
(angdist)	No approximations. Consider full angular distribution of radiation.
(cone)	Cone approximation. All radiation emitted exactly along velocity vector.
(isotropic)	Isotropically emitted radiation.

@RadiationOutput

The @RadiationOutput modules allow one to consider various radiation quantities, such as images, spectra, Green’s functions and more. At the moment, there are five different radiation output sub-modules available: (green), (image), (sovvolume), (space3d), (spectrum) and (topview).

(green)

The type of output generated by this module is known by many names: Green’s functions, weight functions and response functions, to name just a few. We will henceforth use the former, Green’s function.

Mathematically, the Green’s function is defined as the function which, when integrated over phase space with the distribution function yields the detected radiation:

(1)\[I = \iiint G(\rho, p_\parallel, p_\perp) f(\rho, p_\parallel, p_\perp)\, p_\perp\, \mathrm{d}\rho\mathrm{d}p_\parallel\mathrm{d}p_\perp,\]

where \(G\) denotes the Green’s function, \(f\) the distribution function, and the lone factor of \(p_\perp\) at the end of the integrand is the momentum space Jacobian. More explicit expressions for \(G\) can be obtained from the synthetic radiation diagnostic integral (see the discussion on @Radiation), but for the sake of keeping the discussion we brief we will not give those details here.

An alternative definition of the Green’s function employed by some is to move the guiding-center Jacobian out of \(G\), so that \(G = J\tilde{G}\) and (1) instead reads

\[I = \iiint \tilde{G}(\rho, p_\parallel, p_\perp) f(\rho, p_\parallel, p_\perp)\, J\, p_\perp\, \mathrm{d}\rho\mathrm{d}p_\parallel\mathrm{d}p_\perp,\]

\(J\) is the Jacobian for the guiding-center transformation.

This output module allows you to consider many different types of Green’s functions, which may be integrated over any or all of phase space, the image pixels and/or the detector spectral range. It thus provides everything needed to analyse the phase space dependence of detected radiation, or attempt to solve an inverse problem.

Summary of options

Option	Description
@RadiationOutput(green) common	List of common quantities to include in the output file.
@RadiationOutput(green) f_as_linear_array	Store the output function as a linear array.
@RadiationOutput(green) format	Green’s function format (i.e. dependences).
@RadiationOutput(green) mpi_mode	How to generate the Green’s function.
@RadiationOutput(green) output	Name of output file.
@RadiationOutput(green) pixels	Number of pixels.
@RadiationOutput(green) stokesparams	Whether or not to store Stokes parameters.
@RadiationOutput(green) suboffseti	Sub-image offset in vertical direction.
@RadiationOutput(green) suboffsetj	Sub-image offset in horizontal direction.
@RadiationOutput(green) subpixels	Number of pixels in sub-image.
@RadiationOutput(green) with_f	Multiply with the distribution function.
@RadiationOutput(green) with_jacobian	Multiply with the guiding-center Jacobian.

Example configuration

The following example generates a Green’s function with one radial dimension and two pixel dimensions, allowing various radial density profiles to be applied after the simulation to generate the corresponding radiation images:

@RadiationOutput ourOutput (green) {
    format = "rij";
    output = "/path/to/output.mat";
    pixels = 300 300;
}

Output file structure

The output file always contains the following variables:

Option	Description
colpixels	(Only if image dimension j is part of the function) The number of column pixels.
func	Array containing the Green’s function.
rowpixels	(Only if image dimension i is part of the function) The number of row pixels.
stokesparams	Set to 1 if the Green’s function contains the four Stokes parameters.
type	Green’s function format string.
wavelengths	Spectrum grid (wavelength for synchrotron, photon energy normalized to electron rest mass for bremsstrahlung).

Common quantities

By default, the following “common quantities” are also included in the output file:

Name	Description
domain	Tokamak wall or separatrix (depending on which was used).
param1	First momentum parameter grid.
param2	Second momentum parameter grid.
param1name	SOFT name of first momentum parameter.
param2name	SOFT name of second momentum parameter.
r	Radial grid.

For details about which other common quantities can be included in the output, please consult the page about the @RadiationOutput class of modules.

Subset images

Green’s functions can become very large if pixel information is desired. To reduce the amount of unnecessary information stored, it is possible to only store a certain part of a camera image, i.e. a subset image.

The subset image is defined by three parameters:

• Offset in the vertical direction (@RadiationOutput(green) suboffseti)
• Offset in the horizontal direction (@RadiationOutput(green) suboffsetj)
• Number of pixels of subset image (@RadiationOutput(green) subpixels).

In the illustration above, these parameters are suboffseti = 3, suboffsetj = 4 and subpixels = 3.

Note

The sums suboffseti + subpixels(i) and suboffsetj + subpixels(j) must both be less than or equal to the value assigned to* @RadiationOutput(green) pixels.

MPI Mode

When running with MPI, Green’s functions can be generated and stored in two different ways. These two different modes of generation are referred to as contiguous mode and chunked mode, due to the way the function is stored on disk in each.

Contiguous mode

In this mode, the Green’s function will be stored in a single file on disk. Phase space is always divided among the MPI processes, but when all MPI processes have carried out all computations for their parts of phase space, the individual Green’s functions are combined into a single function in the root process, which is then written to a single file.

This mode is useful when the purpose for running with MPI is to further parallelize and speed up the computation. The total amount of memory required is equal to the size of the final Green’s function, multiplied by the number of MPI processes and the number of threads per process. This mode can therefore be very memory-intensive.

Chunked mode

In chunked mode, each MPI process generates its own output file, containing a part of the full Green’s function. The various parts of the Green’s function can then be processed individually and combined to form a full function.

What part of the Green’s function each chunk corresponds to depends on how phase space was divided among the MPI processes, i.e. how the @ParticleGenerator mpi_distribute_mode parameter was set. If, for example, mpi_distribute_mode was set to radius, then each chunk will correspond to the radial interval processed by the MPI process that generated the chunk. Note that this means that the parameter which is divided among the MPI processes must be a part of the Green’s function.

In this mode, the output file name is specified as usual, i.e. as if only one single file were to be created. Each MPI process will then insert an index corresponding to its rank just before the file extension. Setting @RadiationOutput(green) output to ourFile.mat in chunked mode will thus result in a number of files with names ourFile0.mat, ourFile1.mat etc. being generated. MPI ranks are zero-indexed, and thus the output files are so too.

This mode is useful when generating very large Green’s functions, as it allows you take advantage of the large amount of total memory offered by distributed memory systems. To maximize the amount of memory available for the Green’s function in a simulation, set global num_threads to 1 (i.e. one thread per MPI process). This will significantly slow down the simulation, but since each thread stores its own copy of the Green’s function, it will also significantly reduce the memory usage of SOFT.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

f_as_linear_array

Default value:no Allowed values:yes and no

If yes, then func is stored as a linear array (i.e. 1-by-many matrix) instead of the default multi-dimensional format. In this mode, the array must therefore be reshaped to have the correct number of dimensions. In Matlab, this is done by calling reshape(func, [nN, nN_1, ..., n1]), where func is the Green’s function, nN is the number of elements in the last dimension of the function, nN_1 the number of elements in the next-to-last dimension etc. Thus, if the @RadiationOutput(green) format option is set to rij, then the appropriate Matlab command would be:

reshape(func, [ni, nj, nr])

In Python, on the other hand, the order of the dimensions is reversed, so that the equivalent code reads:

import numpy as np
np.reshape(func, (nr, nj, ni))

Note

By default, @RadiationOutput(green) f_as_linear_array is set to no, meaning that the Green’s function is reshaped internally by SOFT before being written to file. Thus, by default, no reshape is required.

This default is new since 2019-04-02.

format

Default value:Nothing Allowed values:Any combination of 1, 2, i, j, r and w.

Specifies the format of the Green’s function, i.e. which parameters the function should depend on and in which order the dependences should be placed. The format is a string consisting of any number of the characters in the table below, in any order.

For example, if @RadiationOutput(green) format is set to r12, the Green’s function \(G(\rho, p_1, p_2)\) will be generated, where \(p_1\) and \(p_2\) denote the momentum space parameters used for the simulation (specified in the @ParticleGenerator module; the momentum parameters are ordered alphabetically, so that \(p_1\) is the momentum parameter which’s name comes first alphabetically). The function \(G(\rho, p_1, p_2)\) will be represented as a 3-dimensional array with the \(\rho\) dependence along the first dimension, \(p_1\) dependence along the second dimension, and \(p_2\) along the third.

Format	Description
1	(Alphabetically) first momentum parameter.
2	(Alphabetically) second momentum parameter.
i	Vertical pixel dimension.
j	Horizontal pixel dimension.
r	Radial parameter.
w	Radiation spectrum.

Note

For pixels, both i and j must specified; they may be specified in any order though.

mpi_mode

Default value:contiguous Allowed values:chunked and contiguous

When SOFT2 is compiled with MPI, specifies how the Green’s function is to be generated and stored. See the discussion above about the MPI Mode for details above the two available modes.

output

Default value:Nothing Allowed values:Any valid file name.

Name of the output file in which to store the result.

pixels

Default value:0 Allowed values:Any positive integer.

Number of pixels in image (if i and j are part of the format string). If only one value is specified, the image becomes quadratic with the same number of pixels in both the vertical and horizontal directions. If two values are given, the first value is interpreted as the number of pixels in the vertical direction and second value as the number of pixels in the horizontal direction.

stokesparams

Default value:no Allowed values:yes or no.

If yes, adds information about the Stokes parameter \((I, Q, U, V)\) to the Green’s function. Another dimension is added to the output array, and becomes the new first dimension. This effectively means that instead of storing one Green’s function, four separate Green’s function corresponding to each of the Stokes parameters is stored contiguously in memory.

suboffseti

suboffsetj

Default value:0 Allowed values:Any non-negative integer.

Specifies the vertical and horizontal offset, respectively, of the subset image.

subpixels

Default value:Same as @RadiationOutput(green) pixels Allowed values:Any positive integer.

Specifies the number of pixels of the subset image. If only one value is specified, the subset image becomes quadratic with the same number of pixels in both the vertical and horizontal directions. If two values are given, the first value is interpreted as the number of pixels in the vertical direction and the second value as the number of pixels in the horizontal direction.

with_f

Default value:no Allowed values:yes or no.

If yes, multiplies the Green’s function with the distribution function. This allows the (green) module to produce proper radiation quantities, such as camera images or Dominant particles.

Note

If this option is enabled, the momentum space Jacobian will still NOT be multiplied with the result, and must be applied manually.

with_jacobian

Default value:yes Allowed values:yes or no

If yes, includes the guiding-center in the definition of the Green’s function (i.e. generates \(G\), as defined at the top of this page).

(image)

The image radiation output generates radiation pixel images. The value of each pixel is determined by the amount of radiation within a (usually very small) pyramid-like field-of-view extending from the detector. A new feature of SOFT2 is that images can be rectangular, and not just square as in the old version of SOFT.

Summary of options

Option	Description
@RadiationOutput(image) common	List of common quantities to include in the output file.
@RadiationOutput(image) output	Sets the name of the output file.
@RadiationOutput(image) pixels	Specifies the number of pixels in each direction of the image.
@RadiationOutput(image) stokesparams	If yes, measures polarization and stores on image per Stokes parameter (synchrotron only).

Detector plane basis

The detector plane is spanned by two vectors which we refer to as \(\hat{e}_1\) and \(\hat{e}_2\). Together with the detector viewing direction (or normal vector) \(\hat{n}\) they form an orthonormal basis in space. The vectors \(\hat{e}_1\) and \(\hat{e}_2\) are used for several purposes, but one of the more important purposes is for spanning the camera images that this @RadiationOutput produces.

The vector \(\hat{e}_1\) is defined so that it always lies in the horizontal plane. Its mathematical definition is

\[\begin{split}\hat{e}_1 = \begin{cases} \hat{y}\cos\vartheta + \hat{z}\sin\vartheta, \quad&\text{ if } \hat{n}\cdot\hat{y} = 0,\\ \left[ \hat{x}\left(\hat{n}\cdot\hat{y}\right) - \hat{y}\left( \hat{n}\cdot\hat{x} \right)\right]\cos\vartheta + \hat{z}\sin\vartheta, \quad&\text{ otherwise}. \end{cases}\end{split}\]

where \(\vartheta\) denotes the angle between \(\hat{e}_1\) and the horizontal plane (the “roll angle”). From this, \(\hat{e}_2\) is defined so that \((\hat{e}_1, \hat{e}_2, \hat{n})\) form a right-handed orthonormal basis:

\[\hat{e}_2 = \hat{n}\times\hat{e}_1.\]

The values of both of these vectors are output on stdout by SOFT during execution.

Example configuration

The following generates a square radiation image that is 600x600 pixels in size:

@RadiationOutput ourImage (image) {
    output = "ourImage.h5";
    pixels = 600;
}

In SOFT2, images can also be rectangular. The following example generates a rectangular image with 600 pixels along the \(\hat{e}_1\) axis (always in the horizontal plane), and 450 pixels along the \(\hat{e}_2\) axis:

@RadiationOutput ourRectangularImage (image) {
    output = "ourRectangularImage.h5";
    pixels = 600 450;
}

One can also measure the polarization of incoming radiation. This feature is only available for spectral-dependent synchrotron radiation. To enable measurement of Stokes parameters, simply enabled the stokesparams option:

@RadiationOutput ourStokesImage (image) {
    output       = "ourStokesImage.h5";
    pixels       = 600 450;
    stokesparams = yes;
}

If stokesparams is yes, then the output file will contain four variables called StokesI, StokesQ, StokesU and StokesV, instead of the usual image variable. Each of the Stokes* variables is an image with each pixel corresponding to the value of the Stokes parameters in its field-of-view.

Output file structure

The output file contains the following variables:

Variable	Included when	Description
image	stokesparams=no	Radiation image matrix.
StokesI	stokesparams=yes	Image of Stokes \(I\) parameter.
StokesQ	stokesparams=yes	Image of Stokes \(Q\) parameter.
StokesU	stokesparams=yes	Image of Stokes \(U\) parameter.
StokesV	stokesparams=yes	Image of Stokes \(V\) parameter.

Note

It is possible to output a pixel graphic image (PNG or PPM) directly, instead of a MAT/HDF5 file as described above. Such pixel graphic images can only contain the actual image, and no additional information.

Common quantities

By default, the following “common quantities” are also included in the output file:

Name	Description
detectorDirection	Unit vector representing viewing direction of detector.
detectorPosition	Vector representing position of detector.
detectorRoll	Detector view roll angle.
detectorVisang	(Full) FOV vision angle of the detector.
wall	Domain contour used for the simulation.

For details about which other common quantities can be included in the output, please consult the page about the @RadiationOutput class of modules.

Note

The actual image is contained either in the image variable if the input parameter stokesparams=no, or in the StokesI, StokesQ, StokesU and StokesQ variables if stokesparams=yes. These variables are matrices of the same size as the number of pixels specified in the input file. Note that the first dimension corresponds to \(\hat{e}_1\), which always lies in the horizontal plane, meaning the one often desires to transpose the image before showing it.

Tip

If you are uncertain about the direction of the image, you can try to move the camera vertically upwards or downwards. When doing so, you should expect the radiation spot to move in the opposite direction in the image.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

output

Default value:Nothing Allowed values:Any valid file name.

Specifies the name of the output file to generate. The file name extension determines the type of the output file. All file types supported by the SFile interface are allowed. In addition, it is also possible to directly output PNG and PPM (i.e. pixel graphic) files directly. With such files, only the image is output and no “common quantities” are included. It also only possible to output the I Stokes parameter in PNG/PPM format.

pixels

Default value:Nothing Allowed values:One or two positive integers.

Specifies the number of pixels in the image along the \(\hat{e}_1\) and \(\hat{e}_2\) directions respectively. Either one or two numbers can be given. If only one number is given, the number of pixels will be the same along both axes (and equal to the given number), yielding a square image. If two numbers are specified, the first number gives the number of pixels along the \(\hat{e}_1\) axis and the second number gives the number of pixels along the \(\hat{e}_2\) axis.

stokesparams

Default value:no Allowed values:yes or no

If yes, measures the polarization of the radiation and produces one image for each of the four Stokes parameters \(I\), \(Q\), \(U\) and \(V\). This feature is only available for synchrotron radiation and requires the detector measure in a limited spectral range.

(sovvolume)

The @RadiationOutput module for calculating the volume of the surface-of-visibility. The volume of the surface-of-visibility depends on a number of factors, including the magnetic field geometry and the detector aperture size.

Summary of options

Option	Description
@RadiationOutput(sovvolume) common	List of common quantities to include in the output file.
@RadiationOutput(sovvolume) output	Name of the output file to generate.

Example configuration

The following example shows how to configure this module:

@RadiationOutput ourVolume (sovvolume) {
    output = "ourVolume.h5";
}

Output file structure

Output files generated by this module contain the following variables:

Variable	Description
volumearray	An \(1\times n_{p_1}n_{p_2}\) array containing the volume for each SoV.

Common quantities

By default, the following common quantities are included with the output generated by this module:

Name	Description
param1	First momentum parameter grid.
param2	Second momentum parameter grid.
param1name	SOFT name of first momentum parameter.
param2name	SOFT name of second momentum parameter.

To add further common quantities, use the @RadiationOutput(sovvolume) common option.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

output

Default value:Nothing Allowed values:Any valid file name.

Specifies the name of the output file to generate. The file name extension determines the type of the output file.

(space3d)

The @RadiationOutput module stores the 3D position of any electron emitting radiation directly at the detector, effectively creating a three-dimensional map of the origin of the radiation.

Directed radiation, such as bremsstrahlung and synchrotron radiation from relativistic electrons, can only be observed when the emitting electron is moving towards the detector. This is expressed in the equation

\[\hat{\boldsymbol{b}}(\boldsymbol{x}) \cdot \frac{\boldsymbol{x} - \boldsymbol{X}}{\left| \boldsymbol{x}-\boldsymbol{X} \right|} = \cos\theta_{\rm p},\]

where \(\hat{\boldsymbol{b}}\) denotes the magnetic field unit vector, \(\boldsymbol{x}\) is the position of the particle, \(\boldsymbol{X}\) is the detector’s position and \(\theta_{\rm p}\) denotes the particle’s pitch angle at \(\boldsymbol{x}\). The solution to this equation is the so-called surface-of-visibility, which can be studied using the Space3D module.

Summary of options

Option	Description
@RadiationOutput(space3d) common	List of common quantities to include in the output file.
@RadiationOutput(space3d) output	Name of output file to generate.
@RadiationOutput(space3d) pixels	Number of pixels in each dimension.
@RadiationOutput(space3d) point0	First edge point of pixel box.
@RadiationOutput(space3d) point1	Second edge point of pixel box.

Specification of pixel box

The Space3D output module creates a 3D histogram of where radiation that reached the detector was emitted from. The histogram is generated in a subset of space that is defined by the pixel box, illustrated below.

This box is defined by its two edge points, point0 and point1, indicated as red dots in the figure above, which specify its width, height and depth. The number of cells, or “pixels”, in the box is the same in each direction, and each cell records the radiation emitted from its region of space.

Note

One can think of the coordinates of each point as the “minimum” and “maximum” points respectively of the region to record. Thus, it is a good idea to always have point0 be element-wise smaller (or larger) than point1, i.e.

\[P^{(0)}_i < P^{(1)}_i, \quad i = x,y,z\]

Example configuration

The following example configures a Space3D output module to generate a surface-of-visibility consisting of \(100^3\) cells:

@RadiationOutput ourSpace3D (space3d) {
    output = "ourSpace3D.h5";
    pixels = 100;
    #         X      Y      Z
    point0 =  0.92, -1.06, -0.68;
    point1 =  2.21,  0.60,  0.48;
}

Output file structure

The output file contains the following variables:

Variable	Description
image	Pixel box/resulting image/histogram.
pixels	Number of pixels per dimension.
xmin	The x-component of point0.
xmax	The x-component of point1.
ymin	The y-component of point0.
ymax	The y-component of point1.
zmin	The z-component of point0.
zmax	The z-component of point1.

Common quantities

By default, no common quantities are included with output generated by this module. To add common quantities, use the @RadiationOutput(space3d) common option.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

output

Default value:Nothing Allowed values:Any valid file name.

Specifies the name of the output file to generate. The file name extension determines the type of the generated file.

pixels

Default value:Nothing Allowed values:Any positive integer.

The number of pixels per dimension in the pixel box.

point0

point1

Default value:Nothing Allowed values:Any vector in 3D space.

Coordinates of the two edge points defining the pixel box. By convention, we usually assign the lower limits of each coordinate to point0 and upper limits to point1.

(spectrum)

The @RadiationOutput module for simulating spectrometers.

Summary of options

Option	Description
@RadiationOutput(spectrum) common	List of common quantities to include in the output file.
@RadiationOutput(spectrum) output	Name of the output file to generate.
@RadiationOutput(spectrum) stokesparams	Specifies whether or not to store Stokes parameters.

Example configuration

The following example shows how to configure this module:

@RadiationOutput ourSpectrum (spectrum) {
    output = "ourSpectrum.h5";
}

Output file structure

Output files generated by this module contain the following variables:

Variable	Description
I	Full spectrum (also known as Stokes \(I\) parameter).
Q	Spectrum corresponding to Stokes \(Q\) parameter.
U	Spectrum corresponding to Stokes \(U\) parameter.
V	Spectrum corresponding to Stokes \(V\) parameter.
wavelengths	List of wavelengths/photon energies corresponding to each element in the spectra.

Note

The name wavelengths is a bit misleading, as the values have the meaning of photon energies (normalized to \(m_e c^2\), the electron rest energy) when simulating bremsstrahlung.

Note

For synchrotron radiation, the unit of the wavelengths array is meters.

Common quantities

By default, no common quantities are included with output generated by this module. To add common quantities, use the @RadiationOutput(spectrum) common option.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

output

Default value:Nothing Allowed values:Any valid file name.

Specifies the name of the output file to generate. The file name extension determines the type of the output file.

stokesparams

Default value:no Allowed values:yes or no.

If yes, adds information about the Stokes parameter \((I, Q, U, V)\) to the Green’s function. Another dimension is added to the output array, and becomes the new first dimension. This effectively means that instead of storing one Green’s function, four separate Green’s function corresponding to each of the Stokes parameters is stored contiguously in memory.

(topview)

The topview @RadiationOutput is used to generate top-view maps of where radiation is emitted from. The options and output variables for this module are very similar to that of the (image) module and they mainly differ in how the image data they generate are interpreted.

Summary of options

Option	Description
@RadiationOutput(topview) common	List of common quantities to include in the output file.
@RadiationOutput(topview) output	Sets the name of the output file.
@RadiationOutput(topview) pixels	Specifies the number of pixels in each direction of the topview.

Top view geometry

The top view geometry is independent of the detector properties. The global (tokamak) x-axis is always along the horizontal in the image, while the global (tokamak) y-axis is always along the vertical. The center of the image corresponds to the axis of symmetry of the tokamak. The physical size of the image is twice the maximum major radius of the domain used (i.e. either tokamak wall or magnetic field separatrix).

An illustration of the geometry of the topview image is shown below. The dashed square indicates the edges of the image and the dotted circle indicates the location of the inner wall, i.e. the minimum radius of the domain used.

Example configuration

The following generates a topview that is 600x600 pixels in size:

@RadiationOutput ourTopview (topview) {
    output = "ourTopview.h5";
    pixels = 600;
}

Note

In contrast to (image), topviews are always square images.

Output file structure

The output file contains the following variables:

Variable	Description
image	Radiation image matrix.

Common quantities

By default, the following “common quantities” are also included in the output file:

Name	Description
detectorDirection	Unit vector representing viewing direction of detector.
detectorPosition	Vector representing position of detector.
detectorVisang	(Full) FOV vision angle of the detector.
wall	Domain contour used for the simulation.

For details about which other common quantities can be included in the output, please consult the page about the @RadiationOutput class of modules.

All options

common

Default value:none Allowed values:See the list on @RadiationOutput.

Specifies which “common quantities” to include in the output file. A full list of possible options is given on @RadiationOutput.

output

Default value:Nothing Allowed values:Any valid file name.

Specifies the name of the output file to generate. The file name extension determines the type of the output file.

pixels

Default value:Nothing Allowed values:Any positive integer.

Specifies the number of pixels in each direction of the image. Thus, the total number of pixels in the image will be the square of this number. Note that in contrast to :ref:`module-ro-image`, topviews are always square and as such only one number can be assigned to this option.

Available sub-modules

Module name	Output description
(green)	Green’s/weight functions
(image)	Camera images
(sovvolume)	Calculate the surface-of-visibility volume.
(space3d)	3D maps of radiation
(spectrum)	Radiation spectra
(topview)	Tokamak topviews of radiation

Include common data

All @RadiationOutput modules support the common option, which allows you to specify which common/meta data to include in the output file. The value assigned to the option should be a list of names of the data objects to include in the file. The available options are shown in the table below.

Name	Description
detectorAperture	Detector size/aperture.
detectorDirection	Detector viewing direction.
detectorPosition	Detector position (relative to point-of-symmetry).
detectorVisang	Detector vision angle (FOV).
domain	Orbit solution domain. Either tokamak wall or separatrix.
param1	First momentum grid
param1name	Name of first momentum grid parameter.
param2	Second momentum grid
param2name	Name of second momentum grid parameter.
r	Radial grid (major radius).
wall	Alias for domain.

In addition to these, it is also possible to specify any of the options in the table below. Those options do not however represent single objects, such as those in the table above, but instead enables or disables bulks of objects to output.

Name	Description
all	Include all available objects in the output.
default	Include the default objects in the output.
none	Do not include any common output.

Additionally, it is possible to prefix any of the options in the first table with either a + or - to indicate whether the object should be included (+) or not (-). Options specified later overrides former ones, meaning that the line:

common = all -domain -wall

would include all the available objects in the output, except for the domain and wall objects. Similarly,

common = default none +domain

would first enable all default objects, then undo that option and remove all objects, and finally add the domain. Thus, the only common object in the output would be the domain object.

To see which common quantities are included with the default option, please consult the page of the relevant @RadiationOutput module.

Example configuration

Please, see the pages for each sub-module for examples of how to configure each of them.

Note that several @RadiationOutput modules can be used with each @Radiation module. Just specify the name of each output module to use, separated by spaces or commas, in the @Radiation configuration block:

@Radiation rad {
    ...
    output = outModule1 outModule2 outModule3;
    # ...or...
    output = outModule1,outModule2,outModule3;
    # ...or even...
    output = outModule1, outModule2, outModule3;
}

To specify which common quantities to include, the option common should be included in the appropriate @RadiationOutput configuration block:

@RadiationOutput outModule1 (XXX) {
    ...
    common = default +domain;
    ...
}

@RadialProfile

To enable combinations of momentum-space only kinetic solvers with varying radial distributions of electrons, SOFT2 provides the RadialProfile module which applies a momentum-independent radial profile to a distribution function. This module must be combined with any of the @DistributionFunction modules that support coupling.

(bessel)

This module provides a radial profile that decays according to a zeroth-order Bessel function of the first kind, from a radius \(r_{\rm min}\) to a radius \(r_{\rm max}\) and is zero outside of that radial interval. More precisely, on the interval \(r_{\rm min}\leq r \leq r_{\rm max}\) the radial profile is

\[f_r(r) = J_0\left( \frac{r-r_{\rm min}}{r_{\rm max}-r_{\rm min}} x_0 \right),\]

and it is identically zero everywhere else. Here, \(J_0(x)\) denotes a zeroth-order Bessel function of the first kind and \(x_0\) is its first zero.

This type of radial profile is expected if runaway electrons are radially transported diffusively in the plasma, with a uniform diffusion coefficient.

Summary of options

The following parameters can be set on a bessel radial profile.

Option	Description
radprof-bessel amax	Normalized minor radius of the outer edge of the electron beam.
radprof-bessel amin	Normalized minor radius of the inner edge of the electron beam.
radprof-bessel rhomax	Major radius of the outer edge of the electron beam.
radprof-bessel rhomin	Major radius of the inner edge of the electron beam.
radprof-bessel rmax	Minor radius of the outer edge of the electron beam.
radprof-bessel rmin	Minor radius of the inner edge of the electron beam.

Example configuration

This example illustrates how a bessel radial profile can be defined and used together with the (avalanche) distribution function in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge

    # Set the radial profile
    radprof   = ourradprof;
}

@RadialProfile ourradprof (bessel) {
    # Set normalized minor radius
    # These settings creates a beam with radius
    # 2/3 of the device minor radius.
    amin = 0;
    amax = 0.67;
}

All options

amax

amin

rhomax

rhomin

rmax

rmin

Default value:amin = 0 and amax = 1. Allowed values:Any radial position that is inside the plasma and on the outboard side.

Specifies the inner and outer edges of the electron beam. The prefix (a*, r*, rho*) specifies whether the edge is given in normalized minor radius, regular minor radius or major radius coordinates.

(gaussian)

This module provides a Gaussian shaped radial profile. The radial profile is defined according to

\[f_r(r) = a \exp\left[ -\left(\frac{\rho-b}{c}\right)^2 \right]\]

This radial profile has been observed to match well with hollow runaway beams.

Note

This module was kindly contributed by Giorgio Ghillardi.

Summary of options

The following parameters can be set on a gaussian radial profile.

Option	Description
radprof-gaussian amax	Normalized minor radius of the outer edge of the electron beam.
radprof-gaussian amin	Normalized minor radius of the inner edge of the electron beam.
radprof-gaussian rhomax	Major radius of the outer edge of the electron beam.
radprof-gaussian rhomin	Major radius of the inner edge of the electron beam.
radprof-gaussian rmax	Minor radius of the outer edge of the electron beam.
radprof-gaussian rmin	Minor radius of the inner edge of the electron beam.
radprof-gaussian gau_a	Gaussian scale factor.
radprof-gaussian gau_b	Gaussian shift parameter.
radprof-gaussian gau_c	Gaussian width parameter.

Example configuration

This example illustrates how a gaussian radial profile can be defined and used together with the (avalanche) distribution function in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge

    # Set the radial profile
    radprof   = ourradprof;
}

@RadialProfile ourradprof (gaussian) {
    # Set normalized minor radius
    # These settings creates a beam with radius
    # 2/3 of the device minor radius.
    amin = 0;
    amax = 0.99;

    # Create gaussian profile f(r) = a * exp(-(r-b)^2/c^2)
    gau_a = 0.05;
    gau_b = 1.80;
    gau_c = 0.04;
}

All options

amax

amin

rhomax

rhomin

rmax

rmin

Default value:amin = 0 and amax = 1. Allowed values:Any radial position that is inside the plasma and on the outboard side.

Specifies the inner and outer edges of the electron beam. The prefix (a*, r*, rho*) specifies whether the edge is given in normalized minor radius, regular minor radius or major radius coordinates.

gau_a

Default value:Nothing Allowed values:Any real number.

Scale factor in front of radial profile.

gau_b

Default value:Nothing Allowed values:Any real number.

Shift of the peak in radius of the Gaussian profile.

gau_c

Default value:Nothing Allowed values:Any real number.

Width of the Gaussian profile.

(linear)

This module provides a radial profile that decays linearly from a radius \(r_{\rm min}\) to a radius \(r_{\rm max}\) and is zero outside of that radial interval. More precisely, on the interval \(r_{\rm min}\leq r \leq r_{\rm max}\) the radial profile is

\[f_r(r) = 1 - \frac{r - r_{\rm min}}{r_{\rm max} - r_{\rm min}},\]

and it is identically zero everywhere else.

Summary of options

The following parameters can be set on a linear radial profile.

Option	Description
radprof-linear amax	Normalized minor radius of the outer edge of the electron beam.
radprof-linear amin	Normalized minor radius of the inner edge of the electron beam.
radprof-linear rhomax	Major radius of the outer edge of the electron beam.
radprof-linear rhomin	Major radius of the inner edge of the electron beam.
radprof-linear rmax	Minor radius of the outer edge of the electron beam.
radprof-linear rmin	Minor radius of the inner edge of the electron beam.

Example configuration

This example illustrates how a linear radial profile can be defined and used together with the (avalanche) distribution function in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge

    # Set the radial profile
    radprof   = ourradprof;
}

@RadialProfile ourradprof (linear) {
    # Set normalized minor radius
    # These settings creates a beam with radius
    # 2/3 of the device minor radius.
    amin = 0;
    amax = 0.67;
}

All options

amax

amin

rhomax

rhomin

rmax

rmin

Default value:amin = 0 and amax = 1. Allowed values:Any radial position that is inside the plasma and on the outboard side.

Specifies the inner and outer edges of the electron beam. The prefix (a*, r*, rho*) specifies whether the edge is given in normalized minor radius, regular minor radius or major radius coordinates.

(power)

This module provides a radial profile the decays according to a power-law from a radius \(r_{\rm min}\) to a radius \(r_{\rm max}\) and is zero outside of that radial interval. More precisely, on the interval \(r_{\rm min}\leq r \leq r_{\rm max}\) the radial profile is

(1)\[f_r(r) = 1 - \left( \frac{r - r_{\rm min}}{r_{\rm max} - r_{\rm min}} \right)^b,\]

and it is identically zero everywhere else. The parameter \(b\) must be specified by the user.

Summary of options

The following parameters can be set on a power radial profile.

Option	Description
radprof-power exponent	Exponent in radial profile.
radprof-power amax	Normalized minor radius of the outer edge of the electron beam.
radprof-power amin	Normalized minor radius of the inner edge of the electron beam.
radprof-power rhomax	Major radius of the outer edge of the electron beam.
radprof-power rhomin	Major radius of the inner edge of the electron beam.
radprof-power rmax	Minor radius of the outer edge of the electron beam.
radprof-power rmin	Minor radius of the inner edge of the electron beam.

Example configuration

This example illustrates how a power radial profile can be defined and used together with the (avalanche) distribution function in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge

    # Set the radial profile
    radprof   = ourradprof;
}

@RadialProfile ourradprof (power) {
    # Set normalized minor radius
    # These settings creates a beam with radius
    # 2/3 of the device minor radius.
    amin = 0;
    amax = 0.67;

    # Create quadratically decreasing profile
    exponent = 2;
}

All options

exponent

Default value:None (i.e. must be defined) Allowed values:Any real value.

Specifies the exponent \(b\) in (1).

amax

amin

rhomax

rhomin

rmax

rmin

Default value:amin = 0 and amax = 1. Allowed values:Any radial position that is inside the plasma and on the outboard side.

Specifies the inner and outer edges of the electron beam. The prefix (a*, r*, rho*) specifies whether the edge is given in normalized minor radius, regular minor radius or major radius coordinates.

(uniform)

This module provides a radial profile that is one (1) everywhere.

Summary of options

This module has no options.

Example configuration

This example illustrates how a uniform radial profile can be defined and used together with the (avalanche) distribution function in SOFT2 simulations:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge

    # Set the radial profile
    radprof   = ourradprof;
}

@RadialProfile ourradprof (uniform) {}

Note that if the option radprof of the distribution function is not specified, a uniform radial profile is automatically generated. Thus, the following configuration yields the same result as the above:

# Global parameter specifying which distribution function to use
distribution_function = analytical_avalanche;

@DistributionFunction analytical_avalanche (avalanche) {
    EHat     = 10;            # Electric field strength (normalized to
                              # the Connor-Hastie critical electric field)
    lnLambda = 17;            # Coulomb logarithm
    Zeff      = 4;             # Effective plasma charge
}

Types

SOFT provides several different types of radial profiles, each with its own set of parameters. Which type of radial profile to use is specified using the secondary type of the configuration block, i.e. by giving the type in parentheses after the block name. The available distribution function types are listed in the table below.

Type	Function
(bessel)	\(J_0\left(\frac{r - r_{\rm min}}{r_{\rm max} - r_{\rm min}}x_0\right)\)
(gaussian)	\(a\exp\left[-\left(\frac{\rho-b}{c}\right)^2\right]\)
(linear)	\(1 - \frac{r - r_{\rm min}}{r_{\rm max} - r_{\rm min}}\)
(power)	\(1 - \left( \frac{r - r_{\rm min}}{r_{\rm max} - r_{\rm min}} \right)^b\)
(uniform)	\(1\)

Here, \(r\) denotes the minor radial location of the particle, \(r_{\rm min}\) and \(r_{\rm max}\) are the minor radial coordinates of the inner and outer edges of the electron beam, and \(b\) is a model parameter. Note that the default values \(r_{\rm min}\) and \(r_{\rm max}\) are 0 and the separatrix radius respectively. The radii may not be negative, and if \(r_{\rm min} > 0\), the electron beam will be hollow.

Example configuration

Please, see the pages for the various radial profile types to view examples of how to configure each.

@Radiation

This tool calculates various radiation quantities, including radiation images, spectra and Green’s functions. It is the tool to use for studying bremsstrahlung and synchrotron radiation.

The basic purpose of the @Radiation tool is to evaluate various forms of the radiation diagnostic integral [2].

(1)\[I = \int \Theta\left( \frac{\boldsymbol{r}}{r} \right) \frac{\boldsymbol{r}\cdot\hat{\boldsymbol{n}}}{r^3} \frac{\mathrm{d}I(\boldsymbol{x},\boldsymbol{p},\boldsymbol{r})}{\mathrm{d}\Omega} f(\boldsymbol{x},\boldsymbol{p})\,\mathrm{d}\boldsymbol{p}\,\mathrm{d} V\,\mathrm{d} A.\]

where \(I\) denotes a general radiation quantity (e.g. radiation power, spectral power, etc.), \(\Theta\) is the field-of-view step function, \(\boldsymbol{r}\) is a vector from the particle to the detector, \(\hat{\boldsymbol{n}}\) is the detector surface normal vector, \(\mathrm{d}I/\mathrm{d}\Omega\) is the angular distribution of radiation and \(f(\boldsymbol{x},\boldsymbol{p})\) is the distribution function. The integral is taken over all of momentum space (indicated by the differential \(\mathrm{d}\boldsymbol{p}\)), real space (indicated by the volume element \(\mathrm{d}V\)), and the detector surface (indicated by \(\mathrm{d} A\)).

To simplify the computation, SOFT evaluates the above integral in guiding-center coordinates. Using these coordinates, the radiation diagnostic integral (1) can be written

\[I = \int\Theta\left( \frac{\boldsymbol{r}}{r} \right) \frac{\boldsymbol{r}\cdot\hat{\boldsymbol{n}}}{r^3} \frac{\mathrm{d}I(\boldsymbol{X},\boldsymbol{p},\boldsymbol{r})}{\mathrm{d}\Omega} f(\rho,p_\parallel,p_\perp)\,J\,\underbrace{\mathrm{d}p_\parallel\mathrm{d}p_\perp\mathrm{d}\zeta}_{\mathrm{d}\boldsymbol{p}} \,\underbrace{\mathrm{d}\rho\mathrm{d}\tau\mathrm{d}\phi}_{\mathrm{d}\boldsymbol{X}}\,\mathrm{d} A.\]

where now \(\boldsymbol{X}\) denotes the guiding-center position, \(p_\parallel\) and \(p_\perp\) are the particle momenta in the directions parallel and perpendicular respectively to the magnetic field, \(\zeta\) is the gyro angle, \(\rho\) is the maximum major radius visited by a guiding-center along its orbit, \(\tau\) is the time along the orbit of the guiding-center and \(\phi\) is the toroidal angle.

[1]	Hoppe et al., 2018, “SOFT: a synthetic synchrotron diagnostic for runaway electrons”. Nuclear Fusion 58 (2), 026032 doi:10.1088/1741-4326/aa9abb.

[2]	Hoppe, 2019, “Simulation and analysis of radiation from runaway electrons”. Licentiate thesis Available online.

Summary of options

Option	Description
@Radiation detector	Specifies which detector configuration to use
@Radiation ignore_trapped	Discards all trapped orbits
@Radiation model	Specifies which radiation model to use
@Radiation ntoroidal	Toroidal resolution parameter
@Radiation output	Specifies which output module(s) to use
@Radiation torthreshold	Parameter for Toroidal optimization
@Radiation torquad	Quadrature rule to use when evaluating toroidal integral
@Radiation wall_opacity	Specifies the wall “opacity”

Example configuration

The @Radiation is merely the parent of a set modules which together produce the desired simulation output. As such, we must specify both the detector, the radiation model and output type. An example configuration of a @Radiation module, along with its required sub-modules, is:

@Radiation rad {
    detector     = det;
    model        = cone;
    ntoroidal    = 7500;
    output       = image topview;
}

@Detector det {
    aperture     = 0.006;
    direction    = 0, 1, 0;
    position     = 0, 1.7, 0;
    vision_angle = 1.25 fov;
    spectrum     = 440e-9, 790e-9, 40;
}

@RadiationModel cone (cone) {
    emission = synchrotron;
}

@RadiationOutput image (image) {
    pixels = 600;
    output = "myimage.mat";
}

Available sub-modules

There are three types of sub-modules that must be configured for the @Radiation module. In addition to a @Detector, one radiation model must specified as well as at least one output module.

Output sub-modules

Radiation output modules are specified with the block type @RadiationOutput. The secondary type of the block (in parentheses after the block name) determines which type of output the block configures. The available secondary types of @RadiationOutput are

Module name	Output description
(green)	Green’s/weight functions
(image)	Camera images
(space3d)	3D maps of radiation
(spectrum)	Radiation spectra
(topview)	Tokamak topviews of radiation

Radiation model sub-modules

Radiation model modules are specified with the block type @RadiationModel. The secondary type of the block (in parentheses after the block name) determines which type of model the block configures. The available secondary types of @RadiationModel are

Module name	Model description
(angdist)	Full angular (and spectral) distribution of radiation
(cone)	Special model for approximating directed radiation
(isotropic)	Special model for perfectly isotropic radiation

Toroidal optimization

Options

detector

Default value:Nothing Allowed values:Name of any @Detector configuration block

Specifies the name of the configuration block to use for setting the properties of the detector.

ignore_trapped

Default value:No (include trapped orbits) Allowed values:yes and no

Since trapped runaway electrons are rare, calculations can sometimes be sped up, and numerical issues avoided, by discarding trapped orbits. In particular, SOFT has problems calculating the guiding-center Jacobian numerically for trapped orbits, and so including trapped orbits can yield unphysical results if one is not careful.

model

Default value:Nothing Allowed values:Name of any radiation model configuration block

Specifies the name of the configuration block to use for setting the radiation model to use. The radiation model basically specifies how the angular distribution of radiation is handled. SOFT can take the full angular distribution of radiation into account, but usually, for synchrotron radiation, the approximative model known as the “cone model” is often used instead. A list of available radiation models can be found above, under the section Radiation model sub-modules.

ntoroidal

Default value:3500 Allowed values:Any positive integer

Number of toroidal sections to divide the tokamak into. This is the resolution parameter for the toroidal integral in the radiation diagnostic integral evaluated by the @Radiation tool.

output

Default value:Nothing Allowed values:List of names of radiation output module configuration blocks

List of names of configuration blocks setting the properties of the output modules to use.

The @Radiation tool only facilitates the computation of various radiation quantities (such as images and spectra). The actual evaluation of these quantities, as well as subsequent generation of output files, are handled by the corresponding “radiation output” modules. A full list of available radiation output modules can be found above under the section Output sub-modules.

torthreshold

Default value:0 Allowed values:Any real value between or equal to 0 and 1

Threshold for neglecting the integrand when using the maximize quadrature to evaluate the toroidal integral. The integration stops as soon as the value of the integrand is a fraction torthreshold of the maximum integrand value seen so far.

For the cone model, this parameter can safely be set to 0. When used together with the models that take the full angular distribution into account, this parameter should be set to a value greater than 0 (yet less than 1).

torquad

Default value:maximize Allowed values:maximize, trapz

Determines which quadrature rule to use for evaluating the toroidal integral. The trapz quadrature is a simple trapezoidal rule. The maximize rule is based on the trapezoidal rule, but uses an optimization algorithm to determine which parts of space that will contribute with radiation. The maximize quadrature is often between a factor 25-100 faster than the regular trapezoidal rule.

wall_opacity

Default value:semi Allowed values:opaque, semi, transparent

Sets the “opacity” level of the wall. If opaque, all walls are fully accounted for, and radiation is not allowed to pass the wall. Conversely, when set to transparent, walls are not accounted for, and the tokamak appears to be transparent, effectively allowing radiation to pass through walls unaffected.

The setting semi is a middle-ground, where only wall segments located at a radius less than the tokamak major radius are accounted for. This means that the tokamak central column is correctly accounted for, while the camera is allowed to be located outside the tokamak wall without radiation being blocked from it. This setting is a way of emulating diagnostic ports in which the radiation diagnostic may be somewhat retracted behind the regular tokamak wall boundary level.

Module	Description
@DistributionFunction	Distribution functions
@Equation	Equations of motion to solve
@MagneticField	Magnetic field module
@ParticleGenerator	Configuration of phase-space
@ParticlePusher	Particle pusher and time-stepping
@RadialProfile	Radial profile

Tool	Description
@Orbits	Simulate and output particle/guiding-center orbits
@Radiation	Synthetic radiation diagnostic tool

Radiation sub-module	Description
@Detector	Synthetic detector configuration
@RadiationModel	Radiation model (cone/angular distribution)
@RadiationOutput	Radiation output module

Magnetic fields

To learn about how magnetic fields are configured in SOFT, check out @MagneticField.

Magnetic fields in SOFT come in two flavours: analytical and numerical. The former allows the user to specify common properties of the magnetic field, such as its on-axis strength and major radius, and generates a corresponding magnetic field with circular flux surfaces. The latter takes a numerical representation of a 2D magnetic field as input, and permits more complicated magnetic equilibria.

Analytical magnetic field

The analytical magnetic field in SOFT is defined as

\[B(r,\theta) = \frac{B_0 R_{\rm m}}{R_{\rm m} - r\cos\theta}\left( \frac{r\sigma_I}{q(r)R_{\rm m}} \hat{\boldsymbol{\theta}} + \sigma_B\hat{\boldsymbol{\varphi}} \right).\]

Here, \(r\) denotes minor radius, \(\theta\) denotes poloidal angle, and \(\hat{\boldsymbol{\theta}}\) and \(\hat{\boldsymbol{\varphi}}\) are unit vectors in the poloidal and toroidal directions respectively. The parameters of the magnetic field that must be specified by the user are the on-axis field strength \(B_0\), the major radius \(R_{\rm m}\), the safety factor \(q(r)\) and the toroidal field sign \(\sigma_B\).

The radial profile of the safety factor can take four different forms in SOFT:

\[ \begin{align}\begin{aligned}q_{\rm const}(r) &= q_{a1},\\q_{\rm lin}(r) &= q_{a1} r/r_0 + q_{a2},\\q_{\rm quad}(r) &= q_{a1} (r/r_0)^2 + q_{a2},\\q_{\rm exp}(r) &= e^{q_{a1} r/r_0} + q_{a2},\\q_{\rm curr}(r) &= \frac{q_{a1}}{1 - \frac{2}{n+2}\left(\frac{r}{r_0}\right)^{q_{a2}}},\end{aligned}\end{align} \]

where \(r_0\) is the tokamak minor radius, \(r\) is the minor radius coordinate, ranging between 0 and \(r_0\), and \(q_{a1}\) and \(q_{a2}\) are constants set by the user.

Relation to plasma current (density)

The safety factor denoted \(q_{\rm curr}(r)\) above is the safety factor resulting from a plasma current density of the form

\[\boldsymbol{j} = \sigma_I\hat{\boldsymbol{\varphi}}j(r) = \sigma_I\hat{\boldsymbol{\varphi}} j_0\left[ 1 - \left(\frac{r}{r_0}\right)^{q_{a2}} \right].\]

where \(j_0\) is the central current density. The parameter \(q_{a1}\) is interpreted as the on-axis safety factor, and is related to the central current density through

\[q_{a1} = \frac{2B_0}{\mu_0 j_0 R_{\rm m}}.\]

From this, we can also relate the total plasma current \(I\) the central safety factor \(q_0\equiv q(0)\equiv q_{a1}\):

\[I = \frac{2\pi B_0}{\mu_0 q_0 R_{\rm m}} r_0^2 \left( 1 - \frac{2}{n+2} \right),\]

or conversely

\[q_0 = \frac{2\pi B_0}{\mu_0 R_{\rm m} I} r_0^2 \left( 1 - \frac{2}{n+2} \right),\]

where \(\mu_0\approx 4\pi\times 10^{-7}\,\text{H/m}\) is the vacuum permeability. With \(B_0\), \(\mu_0\), \(R_{\rm m}\) and \(r_0\) in SI units, and the total plasma current in mega-ampere, the central safety factor can be computed from

\[q_0 = \frac{5n}{n+2} \frac{B_0 r_0^2}{R_{\rm m} I}.\]

More information about how to configure an analytical magnetic field can be found at @MagneticField.

Numerical magnetic field

A file defining a numerical magnetic field must contain tables of the radial, vertical and toroidal magnetic field components, as well as wall contours and some meta data. Using this information, SOFT interpolates in the magnetic field data to evaluate its components in arbitrary points of space, and uses the wall contours to block out radiation.

SOFT is able to load three different types of files: HDF5, MAT and SDT. The former two can be generated by a wide range of tools and programming languages, while the latter is a very simple text-based format developed specifically for SOFT and designed to allow exporting data to at text file on systems with limited access to more advanced libraries. Common to all these file formats is that they allow us to store data in variables. SOFT therefore expects any given input file to contain certain variables that provide the information necessary to run a simulation The below table lists the variables that SOFT look for in a magnetic field file.

Variables of numerical magnetic field

Variable	Mandatory	Type	Description
mf Bphi	Yes	nz-by-nr matrix	Toroidal field component
mf Br	Yes	nz-by-nr matrix	Radial field component
mf Bz	Yes	nz-by-nr matrix	Vertical field component
mf desc	Yes	String	(Meta) Description of data
mf maxis	Yes	2-vector	Location of magnetic axis
mf name	Yes	String	(Meta) Name of magnetic field data
mf Psi	No	nz-by-nr matrix	Poloidal magnetic flux
mf r	Yes	nr-vector	Radial grid
mf separatrix	Yes [1]	2-by-many vector	Last closed flux surface contour
mf verBphi	No	nr-vector	Verification array for Bphi
mf verBr	No	nr-vector	Verification array for Br
mf verBz	No	nr-vector	Verification array for Bz
mf wall	Yes [1]	2-by-many vector	Tokamak wall contour
mf z	Yes	nz-vector	Vertical grid

Beware that some tools handle data in column-major format, whereas SOFT uses column-major format. You may therefore have to transpose certain data to have the correct shape. The vectors ``verXXX`` can be used to ensure that the magnetic field components have the proper format.

[1]	(1, 2) At least one of the separatrix and wall variables must be present in the file.

Bphi

Br

Bz

Shape:nz-by-nr Mandatory:Yes

Magnetic field components in an \(RZ\) plane. These are given as matrices with \(n_z\)-by-\(n_r\) points (i.e. \(n_z\) rows and \(n_r\) columns), where element \(ij\) corresponds to height \(i\) on the vertical grid mf z and radius \(j\) on the radial grid mf r.

Components are given in units of Tesla.

desc

name

Type:Strings Mandatory:Yes

These variables contain meta information about the magnetic field data. While they must be present in the file, there are not requirements on their format. They should be used to keep track of what datasets magnetic field fiels come from.

maxis

Shape:2-by-1 Mandatory:Yes

Radial and vertical coordinates of the magnetic axis: \((R_{\rm axis}, Z_{\rm axis})\).

Psi

Shape:nz-by-nr Mandatory:No

The poloidal magnetic flux, \(\Psi\), given on the same \(R\)-\(Z\) grid as the magnetic field components. If the magnetic flux is provided, SOFT can evaluate the Jacobian for the guiding-center transformation using an analytical expression instead of numerically (which is forced otherwise).

Numerical evaluation of the Jacobian is associated with instabilities in certain parts of phase space and its accuracy can at most be computed to a relative error equal to the square root of the tolerance used for solving the guiding-center orbit.

r

Shape:nr-by-1 Mandatory:Yes

z

Shape:nz-by-1 Mandatory:Yes

Vectors that define the grid on which the magnetic field components are given. The whole grid must be a valid meshgrid, and this enforced by demanding it to be specified in terms of these grid vectors instead.

The magnetic field components are evaluated on this grid, so that component \(x\) of the magnetic field is

\[B_{ij} = B_x\left( r_j, z_i \right),\]

where \(i\) is the row index and \(j\) is the column index and

separatrix

wall

Shape:2-by-many Mandatory:At least one of separatrix and wall

Contains the contour line for the separatrix/wall. The first row contains the \(R\)-coordinates of the contour and the second row contains the \(Z\)-coordinates.

At least one of these must be provided in the magnetic field file. If only one of them is provided, that contour is used both as wall and separatrix. If both are present, SOFT uses the wall to filter out radiation that comes from behind a wall and verify that particles do not collide with the walls. The separatrix is used to define the normalized minor radius.

verBphi

verBr

verBz

Shape:nr-by-1 vector Mandatory:No

These vectors can be used to verify that the magnetic field components have the correct format. Since Matlab tends to be inconsistent with whether a matrix is stored in row-major or column-major form, these vectors can be used to ensure that SOFT reads the magnetic field components with a radial dependence along the column index, and vertical dependence along the row index.

These vectors may have at most as many elements as the number of columns in the magnetic field matrices. There may be fewer elements.

Dominant particles

Running with MPI

SOFT2 is parallelized in two levels. Shared memory parallelization using the OpenMP framework is a core part of SOFT that is always enabled. By default, SOFT will run computations on as many threads as reported by OpenMP, which permits full utilization of single-node computer systems, such as laptops and desktops.

To run on distributed computer systems, however, the Message Passing Interface (MPI) must be used to achieve full parallelization. MPI runs one SOFT2 process on each node of the computer system, and allows them to communicate. Each process can (and will by default) in turn run computations on several OpenMP threads. In contrast to OpenMP, MPI is not enabled in SOFT by default and requires an MPI library to be installed.

In what follows we will describe the steps necessary to run SOFT2 with MPI and the options available for simulations.

Compiling

Before compiling SOFT2, make sure that an MPI library is installed on your system. Note that MPI is a standard, meaning that there are several implementations available. SOFT2 doesn’t care which implementation you use, and so you are free to choose from any of the available. Some popular implementations are

• OpenMPI
• MPICH
• Intel MPI

If you are using the Intel C++ compiler, you should prefer to also use Intel MPI.

To compile SOFT2 with MPI, simply run cmake with the flag -DWITH_MPI=ON:

cmake .. -DWITH_MPI=ON

This is the only modification required when compiling SOFT2, and you should otherwise adhere to the steps described in Compiling.

Note

Note that MPI support need only be configured for SOFT2 and not for softlib. Thus, softlib should be built in exactly the same way, regardless whether MPI support is desired or not.

Setting up a simulation

Any SOFT2 simulation that can be run without MPI should also be possible to run when MPI support is enabled, without any modification to the simulation script. To make the most of the MPI parallelization, SOFT2 also provides a few additional options when compiled with MPI support. These options are however primarily for advanced users, and unnecessary in most cases.

Specify order of phase space

When running with MPI, SOFT2 parallelizes the computation by dividing one of the three phase space parameters (radius and momentum parameters) evenly among the MPI processes. Which parameter to divide is determined by the option @ParticleGenerator mpi_distribute_mode, and by default this is set to the parameter with the largest number of grid points.

Green’s function output

The only type of @RadiationOutput that has a special MPI option is the (green) module. This is because the Green’s function module is typically the most memory-intensive module, and a special option has therefore been added to help better utilize computer clusters.

By default, (green) stores a version of the Green’s function in each MPI process and on each thread. At the end of the run, all local Green’s functions are reduced to a single Green’s function on the main thread of the root MPI process, which is then written to disk. This mode can also be manually specified by setting mpi_mode=contiguous in the (green) module.

To better utilize the available memory however, the Green’s function can be split into several independent pieces, with each MPI process working on its own piece of the function. At the end of the run, each process then writes its piece to a separate file. To use the Green’s function, the separate pieces should in principle be loaded and put in sequential order in memory (although one can think of other techniques for using the Green’s function that doesn’t require all pieces to reside in the same memory space). This mode is designed for generating very large Green’s functions, and is selected via the option mpi_mode=chunked in the (green) module.

Important

Make sure that the phase space parameter that is divided among the MPI processes is a part of the Green’s function, i.e. is in the Green’s functions format string. Which parameter that is divided among the MPI processes is set using the @ParticleGenerator mpi_distribute_mode option.

Running

To run SOFT2 in MPI mode, simply use the mpirun command:

$ mpirun -n N /path/to/soft inputfile

where N is the number of MPI processes to launch. Note that MPI only provides any benefits when run on a computer cluster, in which case a scheduler (such as SLURM) must be used. An example configuration for a SLURM system is presented below, and a similar configuration would most likely be necessary with other scheduler systems as well.

Using SLURM

SLURM is one of the more popular so-called schedulers which is used to submit jobs to a computer cluster. The example below illustrates how a SLURM configuration for a SOFT2 run could look like. Make sure to thoroughly read the documentation for the cluster you are running on though, as desired settings may vary from system to system.

#!/bin/bash
#
# Number of nodes to run on.
#SBATCH -N 32
#
# Number of processor cores. If each node has 32 cores, and you have
# requested 32 nodes above, then you should specify 32*32 = 1024 here.
#SBATCH -n 1024
#
# Desired memory in MB per processor core.
#SBATCH --mem-per-cpu=4000
#
# Maximum execution time of the job.
#SBATCH --time=00:30:00
#
# Queue to submit job to (these names are usually specific to the cluster).
#SBATCH -p sched_express
#
# Name of file to direct stdout/stderr to.
#SBATCH -o soft-%j.out

module load gsl hdf5 mpi

# Run SOFT2
mpirun /path/to/soft inputfile

Examples

This page collects a number of example SOFT simulations that you can use to learn how to run SOFT. The example pages show you the expected result and gives the necessary parameters to reproduce the result. You can then try to write a SOFT configuration script that properly reproduces the image, spectrum or figure. Each example comes with a complete solution as well.

All examples, along with plotting scripts, can be found in the examples/ subdirectory of the SOFT2 repository.

Orbits

Orbits

Example of how to generate guiding-center and particle orbits using SOFT.

Example of orbits computed using SOFT. The left figure shows guiding-center orbits, while the right figure shows particle orbits.

Important points

SOFT always calculates orbits, but to explicitly output the orbits calculated, the @Orbits module must be used.

One-file example

The following is a one-file example of a configuration file that will generate the guiding-center orbits in the figure at the top of this page:

# Simulate orbits in a circular tokamak
# magnetic field.
##############################

magnetic_field     = mf;
tools              = orbits;
include_drifts     = yes;

# Configuration of EAST-like magnetic equilibrium
@MagneticField mf (analytical) {
    B0     = 5;     # On-axis field strength (T)
    Rm     = 0.68;  # Major radius (m)
    rminor = 0.22;  # Minor radius (m)

    # Safety-factor
    qtype  = linear;
    qa1    = 2;
    qa2    = 1;
    sigmaB = ccw;
}

# Phase space grid
@ParticleGenerator PGen {
    a      = 0.1, 0.9, 9;  # Minor radius
    gamma  = 10, 10, 1;    # Energy (mc^2) (approx. 30 MeV)
    xi     = 0.3, 0.3, 1;  # Cosine of pitch angle

    position = guiding-center;
}

# Orbit generator
@ParticlePusher PPusher {
    nt       = 100;        # Number of timesteps per orbit (resolution parameter)
}

@Orbits orbits {
    output = "data/guiding-center.mat";
}

Configuration scripts

The following configuration scripts are available in the SOFT2 GitHub repository and will generate the two figures at the top of this page.

File	Description
baseline	Baseline configuration applied to all simulations.
gc	Simulate guiding-center orbits.
particle	Simulate particle orbits.

Orbits

This example illustrates how guiding-center and particle orbits can be computed using the @Orbits module of SOFT.

Simulation time: Short

Check out: Orbits.

Radiation images

Basic synchrotron image

A basic example of a script for generating a synchrotron radiation image.

Important points

As always, we must configure the following three modules:

• @MagneticField — The magnetic field to use
• @ParticleGenerator — Generates particles in our phase space
• @ParticlePusher — Computes orbits

In addition to these, we must also set up the tool to use. To generate a synchrotron image, we need to use the @Radiation module. This module further requires us to configure the following three sub-modules:

• @Detector — Detector measuring the radiation
• @RadiationModel — Model to use for radiation (i.e. (cone) or (angdist))
• @RadiationOutput — Defining the type of output to generate (i.e. (image))

Example configuration

The following generates a Green’s function \(G(p_\parallel, p_\perp)\):

# Generate a basic single-particle
# synchrotron radiation image.
##################################

magnetic_field     = mf;
tools              = rad;

# Configuration of EAST-like magnetic equilibrium
@MagneticField mf (analytical) {
    B0     = 5;     # On-axis field strength (T)
    Rm     = 0.68;  # Major radius (m)
    rminor = 0.22;  # Minor radius (m)

    # Safety-factor
    qtype  = linear;
    qa1    = 2;
    qa2    = 1;

    sigmaB = ccw;
}

# Phase space grid
@ParticleGenerator PGen {
    a      = 0.0, 0.95, 600; # Normalized minor radius
    p      = 60, 60, 1;
    thetap = 0.2, 0.2, 1;
}

# Orbit generator
@ParticlePusher PPusher {
    nt = 2000;        # Number of timesteps per orbit (resolution parameter)
}

# Radiation tool
@Radiation rad {
    detector = det;

    ntoroidal   = 7000;    # No. of toroidal sections in tokamak (resolution parameter)
    model       = cone;    # Radiation model to use
    output      = image;   # List of configuration of output
}

# Detector properties
# Set up a tangentially viewing HXR camera.
@Detector det {
    aperture     = 0.006;
    position     = 0.68, -0.68, 0;
    direction    = 0, 1, 0;
    vision_angle = 0.78 fov;
    spectrum     = no;
}

# Radiation model
@RadiationModel cone (cone) {
    emission = synchrotron;
}

@RadiationOutput image (image) {
    pixels = 600;
    output = "data/image.mat";
}

Synchrotron images from Zhou et al. (2014)

This example reproduces Figs. 5, 6 and 7 of [Zhou et al. 2014]. Many other figures of that paper can also be reproduced with slight modifications to these examples.

Figure 6 from |zhou2014|_, simulated using SOFT2. This figure illustrates how the radiation spot varies with the magnetic field safety factor.

Tokamak and detector properties

[Zhou et al. 2014] use an analytical magnetic field with circular flux surfaces, which happens to also be implemented in SOFT. The parameters of the magnetic field are as follows:

Parameter	Value
Magnetic field on-axis	2 T
Major radius	185 cm
Minor radius	45 cm
Toroidal field direction	Counter-clockwise (CCW)

For the detector, the parameters used in [Zhou et al. 2014] are (note the directions of the axes in Fig. 4, which are different from how the corresponding axes are defined in SOFT)

Parameter	Value
Position	\(148\hat{x} + 185\hat{y}\,\text{cm}\)
Viewing direction	\(\hat{y}\)

Configuration scripts

Since all synchrotron images of [Zhou et al. 2014] are generated in very similar setups, we can benefit from splitting the configuration into multiple files. To this end, we have divided the problem into three levels.

The first level configures everything that is the same in all figures and is represented by the file baseline in the top-level directory. Next, all subfigures of each figure share some characteristics which are configured in the files named FigXbase (with X replaced by the figure number). These files include the baseline file at the top of the script. Finally, each individual subfigure has its own configuration script. These configuration scripts include the corresponding FigXbase script at the top.

Note

Tip 1 — While giving the camera a spectral range corresponding to the visible range (for example between 380-750 nm, as was used in the experiments) is more physically accurate, it makes comparison between SOFT simulations and [Zhou et al. 2014] more difficult, since [Zhou et al. 2014] only consider the spot shape and not the distribution of radiation across the spot.

File	Description
baseline	Baseline configuration applied to all simulations.
Fig5/Fig5base	Baseline for all subfigures of Fig. 5
Fig5/a	Configuration of Fig. 5(a).
Fig5/b	Configuration of Fig. 5(b).
Fig5/c	Configuration of Fig. 5(c).
Fig5/d	Configuration of Fig. 5(d).
Fig5/e	Configuration of Fig. 5(e).
Fig5/f	Configuration of Fig. 5(f).
Fig5/g	Configuration of Fig. 5(g).
Fig5/h	Configuration of Fig. 5(h).
Fig6/Fig6base	Baseline for all subfigures of Fig. 5
Fig6/e	Configuration of Fig. 6(e).
Fig6/f	Configuration of Fig. 6(f).
Fig6/e	Configuration of Fig. 6(g).
Fig6/h	Configuration of Fig. 6(h).
Fig7/Fig7base	Baseline for all subfigures of Fig. 7
Fig7/a	Configuration of Fig. 7(a).
Fig7/b	Configuration of Fig. 7(b).
Fig7/c	Configuration of Fig. 7(c).
Fig7/d	Configuration of Fig. 7(d).

Full angular distribution of radiation

Although most of the results obtained with SOFT that has been published in the literature have used the cone model, SOFT is perfectly capable of taking the full angular distribution of radiation into account, i.e. NOT assume that all radiation is emitted exactly along the velocity vector. While this model is more accurate than the cone model, it is also significantly slower.

Resulting image when taking the full angular distribution of radiation into account. The spot then usually appears somewhat thicker than when simulated with the cone model.

Important points

The only thing that differs when taking the full angular distribution of radiation into account (in contrast to just using the cone model) is that a different @RadiationModel must be used. In essence, all we need to do is set the secondary type of the radiation model to angdist:

@RadiationModel ourModel (angdist) {
    ...
}

And this is really all there is to it. The type of radiation can be specified as usual.

There are some additional things to keep in mind when simulating the full angular distribution of radiation though. These are

• Toroidal quadrature threshold: If you set @Radiation torquad to maximize (which is also the default) in the @Radiation module, make sure that @Radiation torthreshold is set to a value greater than zero. Otherwise very little is gained from using torquad = maximize.
• Detector resolution: the full angular distribution of radiation requires explicit numerical integration over the detector surface. As such, a quadrature rule to use, as well as the number of points to discretize the detector surface with, can be selected. The quadrature to use is selected using the @RadiationModel(angdist) qrule2d option, while the number of points on the detector surface is selected with @RadiationModel(angdist) nsamples.

Example

The following is a complete configuration script that can be used to generate the synchrotron image seen at the top of this page. This script can also be found in the examples/ subdirectory of the SOFT2 GitHub repository.:

# Configuration for simulation of bremsstrahlung
# radiation images.
#
# Based on Zhou et al., PoP 21 (2014)
# https://doi.org/10.1063/1.4881469
##############################

magnetic_field     = mf;
tools              = rad;

# Configuration of analytical magnetic equilibrium
@MagneticField mf (analytical) {
    B0     = 5;     # On-axis field strength (T)
    Rm     = 0.68;  # Major radius (m)
    rminor = 0.22;  # Minor radius (m)

    # Safety-factor (overriden in other scripts)
    qtype  = constant;
    qa1    = 1;
    sigmaB = ccw;
}

# Phase space grid
@ParticleGenerator PGen {
    a      = 0, 0.95, 400;  # Normalized minor radius
    gamma  = 60, 60, 1;     # Energy (mc^2) (approx. 30 MeV)
    thetap = 0.2, 0.2, 1;   # Pitch angle (rad)

    progress = 20;
}

# Orbit generator
@ParticlePusher PPusher {
    nt       = 5000;        # Number of timesteps per orbit (resolution parameter)
}

# Radiation tool
@Radiation rad {
    detector    = "det";      # Name of detector configuration (quotations optional)
    ntoroidal   = 7000;       # No. of toroidal sections in tokamak (resolution parameter)
    model       = angdist;    # Radiation model to use
    output      = image topview;  # List of configuration of output

    ####################################
    # NOTE: This is required to acquire
    # the speedup from the 'maximize'
    # quadrature rule
    ####################################
    torthreshold = 1e-3;
}

# Detector properties
# Set up a tangentially viewing HXR camera.
@Detector det {
    aperture     = 0.006;
    position     = 0.68, -0.68, 0;
    direction    = 0, 1, 0;
    vision_angle = 0.78 fov;
    spectrum     = no;
}

@RadiationModel angdist (angdist) {
    emission = synchrotron;

    # Number of points (in each direction) on detector
    # (=> in total there will be 'nsamples*nsamples'
    # points covering the detector)
    #
    # This is the resolution parameter corresponding to
    # the integral over the detector surface
    nsamples = 4;
}
@RadiationOutput image (image) {
    pixels = 600;
    output = "data/image.mat";
}
@RadiationOutput topview (topview) {
    pixels = 1000;
    output = "data/topview.mat";
}

Basic synchrotron image

An example of a basic synchrotron radiation image.

Simulation time: Short

Check out: Basic synchrotron image

Synchrotron images from Zhou et al. (2014)

An example showing how to reproduce Figs. 5, 6 and 7 of Zhou et al., PoP 21 (2014).

Simulation time: Short

Check out: Synchrotron images from Zhou et al. (2014).

Full angular distribution of radiation

Run SOFT, taking the full angular distribution of radiation into account.

Simulation time: Medium

Check out Full angular distribution of radiation.

Bremsstrahlung images

Generate bremsstrahlung images

Simulation time: Short

Special output

Spectrometer output

Green’s function output

Green’s function output

SOFT is capable of generating so-called Green’s functions.

Important points

For Green’s function output, the (green) module must be used. It requires at least two parameters to be specified: @RadiationOutput(green) format (specifying the dependences of the Green’s function), and @RadiationOutput(green) output (name of output file).

The format option is a string consisting of any sequence of the following six characters:

Format	Description
1	(Alphabetically) first momentum parameter.
2	(Alphabetically) second momentum parameter.
i	Vertical pixel dimension.
j	Horizontal pixel dimension.
r	Radial parameter.
w	Radiation spectrum.

Dominant particles

Example configuration

The following generates a Green’s function \(G(p_\parallel, p_\perp)\):

# Generate a momentum-space
# Green's function.
##############################

magnetic_field     = mf;
tools              = rad;

# Configuration of EAST-like magnetic equilibrium
@MagneticField mf (analytical) {
    B0     = 5;     # On-axis field strength (T)
    Rm     = 0.68;  # Major radius (m)
    rminor = 0.22;  # Minor radius (m)

    # Safety-factor
    qtype  = linear;
    qa1    = 2;
    qa2    = 1;
}

# Phase space grid
@ParticleGenerator PGen {
    a      = 0.0, 0.95, 20; # Normalized minor radius
    ppar   = -2, -165, 40;     # Parallel momentum (mc)
    pperp  = 1, 15, 40;    # Perpendicular momentum (mc)

    progress = 10;
}

# Orbit generator
@ParticlePusher PPusher {
    nt = 2000;        # Number of timesteps per orbit (resolution parameter)
}

# Radiation tool
@Radiation rad {
    detector = det;

    ntoroidal   = 7000;    # No. of toroidal sections in tokamak (resolution parameter)
    model       = cone;    # Radiation model to use
    output      = green;   # List of configuration of output
}

# Detector properties
# Set up a tangentially viewing HXR camera.
@Detector det {
    aperture     = 0.006;
    position     = 0.68, -0.68, 0;
    direction    = 0, 1, 0;
    vision_angle = 0.78 fov;
    spectrum     = no;
}

# Radiation model
@RadiationModel cone (cone) {
    emission = synchrotron;
}

@RadiationOutput green (green) {
    format = 12;
    output = "data/green.mat";
}

Green’s functions

Generate a SOFT Green’s function. This example illustrates how to generate a simple momentum-space Green’s function for computing total radiated power, but can be readily generalized for any set of Green’s function parameters.

Simulation time: Medium

Check out Green’s function output

Advanced input

Numeric magnetic field

Distribution function

Tips & Suggestions

SOFT gives a ‘Step-size underflow in Runge-Kutta integrator’ error.

This error is usually due to one of two reasons: (i) something is wrong with the magnetic field you are running with, or (ii) the integrator tolerance @Equation tolerance is too strict.

The former reason is possible only when using a numerical magnetic field and is indicative of a problem in one or several of the magnetic field components. You should check to make sure that the magnetic field actually looks fine before proceeding.

The latter reason occurs if @Equation tolerance is set to a too small value. Very few applications require an accuracy in the computed quantities to machine epsilon, and for both performance and stability reasons, @Equation tolerance should therefore be set to a value significantly larger than machine epsilon. On most systems, machine epsilon is around \(2\cdot 10^{-16}\).

Indices and tables

• Index
• Module Index
• Search Page