DistributionFunction

The DistributionFunction class contains all shared settings for the hot and runaway electron distribution functions. The distribution functions are automatically enabled whenever their corresponding grids (hottailgrid and runawaygrid) are enabled.

The distribution functions are solved for using the 3D Fokker-Planck equation

\[\frac{\partial f}{\partial t} + eE\left[ \frac{1}{p^2}\frac{\partial}{\partial p}\left( p^2\xi f \right) + \frac{1}{p}\frac{\partial}{\partial\xi}\left[ \left( 1-\xi^2 \right) f \right] \right] = C\left\{ f \right\} + \frac{1}{\mathcal{V}'}\frac{\partial}{\partial r} \left[ \mathcal{V}'\left( -Af + D\frac{\partial f}{\partial r} \right) \right] + S,\]

where the second term on the LHS represents electric field acceleration along magnetic field lines, \(C\{f\}\) denotes the Fokker-Planck collision operator, the second term on the RHS represents advective-diffusive radial transport with arbitrary advection and diffusion coefficients \(A\) and \(D\), and \(S\) is a collection of source terms (most notably the production of secondary runaway electrons through the avalanche mechanism and the particle source term which ensures quasi-charge neutrality). Additional terms not shown here can also be added, as described further down on this page.

For the hot electron distribution function, it is also possible to solve the pitch angle-averaged 2D Fokker-Planck equation by setting nxi=1 on the hottailgrid. In this case, the equation is reduced with the ordering \(\nu_D \gg eE/p \gg \mathrm{remainder}\), in which case the pitch-angle distribution takes an analytic form with which the equation is analytically averaged. The approximation becomes increasingly better for higher plasma effective charge and lower energy. It is inadequate for describing Dreicer seed generation or the shape of the runaway tail, but can provide a fair approximation for the runaway rate of other mechanisms.

Advection interpolation

In the finite volume method, which is used for discretizing most equations in DREAM (particularly the Fokker-Planck equation), advection terms are discretized as

\[\left[ \frac{1}{\mathcal{V}'}\frac{\partial}{\partial z^k}\left( \mathcal{V}' A_k f \right) \right]_i = \frac{\mathcal{V}_{i+1/2}' A_{k,i+1/2} f_{i+1/2} - \mathcal{V}_{i-1/2}' A_{k,i-1/2} f_{i-1/2}} {\mathcal{V}_i' \Delta z^k}\]

where half-indices are located on cell faces, while the distribution function \(f\) is given at cell centres, corresponding to integer indices. Thus, to obtain a practical discretization for the advection term it becomes necessary to interpolate in \(f\). This is generally achieved by introducing a set of interpolation coefficients \(\left\{ \delta^{(i)}_j \right\}\) and letting

\[f_{i-1/2} = \sum_{j=-2}^1 \delta^{(i)}_j f_{i+j}.\]

The interpolation coefficients can be chosen in a number of ways, and one of the simplest choices imaginable is \(\delta^{(i)}_{-2} = \delta^{(i)}_{1} = 0\), \(\delta^{(i)}_{-1} = \delta^{(i)}_{0} = 1/2\). This is the so-called CENTRED interpolation scheme, which is used by default in DREAM. For strongly advective problems, however, it is usually desirable to choose these coefficients more carefully.

Advective problems can famously yield unstable solutions containing spurious, unphysical oscillations. These oscillations result from the unphysical nature of the interpolation done with schemes such as the CENTRED described above and can in fact be avoided with other interpolation schemes. DREAM provides a number of such positivity-preserving interpolation schemes, which avoid oscillations and ensure that the resulting distribution function is everywhere positive. The simplest schemes are the linear first and second-order upwind schemes, which are computationally cheap but less accurate. For most cases, the non-linear flux limiter schemes are to be preferred.

Note

All flux limiter schemes require the non-linear solver to be used.

List of schemes

The following interpolation schemes are provided by DREAM:

Name	Non-linear?	Reference
AD_INTERP_CENTRED	No
AD_INTERP_UPWIND	No
AD_INTERP_UPWIND_2ND_ORDER	No
AD_INTERP_DOWNWIND	No
AD_INTERP_QUICK	No	Leonard (1980)
AD_INTERP_MUSCL	Yes	van Leer (1979)
AD_INTERP_SMART	Yes	Gaskell and Lau (1988)
AD_INTERP_OSPRE	Yes	Waterson and Deconinck (1995)
AD_INTERP_TCDF	Yes	Zhang et al (2015)

Two useful systematic reviews of advection interpolation schemes and flux limiter methods are given by N P Waterson & H Deconinck (2007) and Zhang et al (2015).

Note

For most applications where flux limiters are desired, we recommend the TCDF flux limiter scheme.

For flux limiter schemes, it is also possible to specify a flux limiter damping parameter. This parameter, here denoted \(\eta\), is used when updating the interpolation coefficients in the non-linear iteration process to under-relax the algorithm as this may accelerate convergence. In each iteration the new interpolation coefficient is set to

\[\delta = \delta_{\rm prev} + \eta\left( \delta_{\rm calc} - \delta_{\rm prev} \right),\]

where \(\delta_{\rm prev}\) is the old value of the interpolation coefficient and \(\delta_{\rm calc}\) is the calculated new value.

Advection jacobian mode

The flux limiter interpolation schemes are non-linear as the coefficients depend on the gradients of the unknown in which they interpolate. There are three different methods of treating the jacobian of the interpolation coefficients:

Name	Description
AD_INTERP_JACOBIAN_LINEAR	Default. Neglects any contribution from the unknown to the jacobian
AD_INTERP_JACOBIAN_FULL	Fully accounts for the unknown contribution to the jacobian
AD_INTERP_JACOBIAN_UPWIND	Uses the jacobian corresponding to UPWIND interpolation coefficients

Note

FULL is only active for the smooth TCDF and OSPRE flux limiter schemes, for all others it defaults to LINEAR. UPWIND requires a non-linear solver as the jacobian does not match the function vector, and it will sometimes take many iterations to converge, but each iteration is faster due to fewer non-zero elements in the matrix.

Example

The following example illustrates how to use the TCDF flux limiter with full jacobian in all dimensions on the hot electron grid:

import DREAM.Settings.Equations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
ds.eqsys.f_hot.setAdvectionInterpolationMethod(
   ad_int=DistFunc.AD_INTERP_TCDF,
   ad_jac=DistFunc.AD_INTERP_JACOBIAN_FULL)

It is also possible to specify different interpolation methods in different dimensions. Using a flux limiter scheme may be more important in the \(p\) dimension due to the strong advective flows due to the electric field, whereas the pitch dimension may be dominated by (diffusive) pitch-angle scattering, and the radial dimension has no fluxes at all. In this case, it may be sufficient to apply the flux limiter to the \(p\) dimension only:

import DREAM.Settings.Equations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
ds.eqsys.f_hot.setAdvectionInterpolationMethod(ad_int_p1=DistFunc.AD_INTERP_TCDF)

The other dimensions default to the AD_INTERP_CENTRED scheme.

Warning

Although flux limiters will generally produce smooth and stable solutions, they are NOT at all guaranteed to be accurate. You should always verify that the results you obtain are correct by varying the various resolution parameters of your simulation (usually nr, np, nxi and nt). Convergence of the non-linear solver is not guaranteed with flux limiters, as sometimes there may be two valid solutions that the system oscillates between. It can help to increase resolution, try fluxlimiterdamping, another jacobian mode or another flux limiter.

Boundary conditions at pMax

Numerical solution of the Fokker-Planck equation on a finite momentum grid requires boundary condition to specified at \(p=p_{\rm max}\), the upper boundary in the momentum dimension. DREAM supports three different forms for this boundary condition, corresponding to

1. setting \(f(p>p_{\rm max}) = 0\) (default)

2. ensuring that the flux out across \(p=p_{\rm max}\) is the same as the flux into the last cell on the momentum grid.

3. extrapolating the flux out across \(p=p_{\rm max}\) based on the flux into the two last cells on the momentum grid.

Which boundary condition should I use? You can choose whichever of the boundary conditions you prefer. The conditions (1) and (2) may be marginally cheaper to use, but all are equally arbitrary. If the boundary condition appears to affect the result of your simulations, you should rather increase \(p_{\rm max}\) so that it is located sufficiently far from where the interesting physics are happening.

In some cases, the boundary condition may give rise to numerical instabilities on the boundary. In these cases it is well motivated to try an alternative boundary condition. If none of the boundary conditions are able to stabilize the solution near the boundary, you should proceed to applying a flux limiter scheme.

List of boundary conditions

Name	Description
BC_F_0	Assume \(f(p>p_{\rm max}) = 0\) (default)
BC_PHI_CONST	Assume \(\Phi^{(p)}_{N_p+1/2} = \Phi^{(p)}_{N_p-1/2}\)
BC_DPHI_CONST	Assume \(\Phi^{(p)}_{N_p+1/2} = \left( 1 - \delta \right)\Phi^{(p)}_{N_p-1/2} + \delta\Phi^{(p)}_{N_p-3/2}\)

Note

When using both hot and runaway electron grids, the boundary condition at \(p_{\rm hot} = p_{\rm hot,max}\) is automatically set to allow particles to flow freely across the boundary between the grids. It is then only necessary to specify the boundary condition at the \(p_{\rm RE} = p_{\rm RE,max}\) boundary on the runaway grid.

Note

All particles leaving the hot electron grid through \(p=p_{\rm max}\) are automatically moved to the runaway density n_re. When particles leave the runaway electron grid through \(p=p_{\rm max}\) they are still counted as part of n_re, and thus this quantity may deviate from the density moment of f_re.

Example

The following example shows how to use the BC_PHI_CONST:

import DREAM.Settings.Equations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
ds.eqsys.f_hot.setBoundaryCondition(DistFunc.BC_PHI_CONST)

Initialization

An initial condition is required when solving the Fokker-Planck equation and in DREAM this can be either in the form of a Maxwellian with prescribed density and temperature profiles, or as an arbitrary prescribed function.

As Maxwellian

DREAM can initialize the distribution function to a Maxwell-Jüttner distribution function:

\[f(r,p,\xi) = \frac{n_0(r)}{4\pi\Theta_0(r) K_2\left(1/\Theta_0(r)\right)} e^{-\gamma/\Theta_0(r)},\]

where \(n_0(r)\) is the initial electron density profile, \(\Theta_0(r) = T_0(r)/m_ec^2\) is the initial electron temperature profile normalized to the electron rest mass, \(\gamma\) is the Lorentz factor and \(K_2(x)\) is the second-order modified Bessel function of the second kind.

To prescribe the initial distribution function as a Maxwellian, the user should call the method setInitialProfiles(), providing the density and temperature profiles determining the shape of the initial distribution function at every radius.

Example

The following example illustrates how to prescribe the initial hot electron distribution function as a Maxwellian:

ds = DREAMSettings()
...
# Radial grid (may be different for n0 and T0)
r = np.linspace(r0, r1, nr)

# Density profile
n0 = np.array([...])     # shape (nr,)  [m^-3]
# Temperature profile
T0 = np.array([...])     # shape (nr,)  [eV]

ds.eqsys.f_hot.setInitialProfiles(n0=n0, T0=T0, rn0=r, rT0=r)

If uniform density and/or temperature profiles are desired, n0 and/or T0 can be given as scalars and their corresponding radial grid parameters omitted from the call to setInitialProfiles():

ds = DREAMSettings()
...
# Uniform density of n=2x10^19 m^-3 and temperature T=1 keV
ds.eqsys.f_hot.setInitialProfiles(n0=2e19, T0=1000)

Arbitrary distribution function

Both f_hot and f_re can also be initialized with arbitrary distribution functions. The functions should be represented as 3D arrays with the shape (nr, np2, np1) (for \(p/\xi\) coordinates, np1=np and np2=nxi). Initialization is done via a call to the method setInitialValue().

Example

The following example illustrates how to prescribe an arbitrary initial hot electron distribution function:

ds = DREAMSettings()
...
r  = np.linspace(r0, r1, nr)
p  = np.linspace(p0, p1, np)
xi = np.linspace(xi0, xi1, nxi)

f0 = np.array([...])     # shape (nr, nxi, np)

ds.eqsys.f_hot.setInitialValue(f=f0, r=r, p=p, xi=xi)

Magnetic ripple

Due to the finite number of toroidal field coils used for tokamaks, the toroidal symmetry of the devices is somewhat broken. The perturbed magnetic field can lead to both increased particle transport, as well as increased pitch-angle scattering.

In DREAM, we can account for the increased pitch-angle scattering using the diffusion operator derived by Kurzan et al (1995), which for DREAM takes the form

\[D^{\xi_0\xi_0} = \frac{\pi}{32}\frac{eB}{m_e} \frac{B_{\rm min}}{B}\frac{\xi^2}{\xi_0^2}\frac{1-\xi_0^2}{p^2} \frac{v_\parallel}{c}\left( \frac{\delta B_{nm}}{B} \right)^2 H_{nm}(p\xi),\]

with

\[\begin{split}H_{nm}(p\xi) = \begin{cases} \frac{1}{\Delta p_{nm}}, \qquad & \left|p_\parallel - p_{nm} \right| < \Delta p_{nm},\\ 0, \qquad & \text{otherwise} \end{cases}\end{split}\]

and resonant momentum and resonance width

\[\begin{split}p_{nm} &= \frac{e G(r)}{m_e cn N_{\rm c}},\\ \Delta p_{nm} &= \sqrt{\frac{\delta B_{nm}}{B} p_\parallel p_\perp}.\end{split}\]

where \(G(r)\) gives the toroidal magnetic field variation and \(N_{\rm c}\) denotes the number of toroidal field coils. Since the derivation assumes small pitch angles and slab geometry, the geometric factor \(\xi^2B_{\rm min}/\xi_0^2B = 1\).

We also allow the user to employ an alternative model for the resonance width, where instead of a sharp boundary we model the resonance with a smooth Gaussian region

\[H_{nm} = \frac{2}{\sqrt{\pi}\Delta p_{nm}} \exp\left[ \frac{(p\xi - p_{nm})^2}{\Delta p_{nm}} \right]\]

where the characteristic width is evaluated at \(p_\parallel=p_{nm}\),

\[\Delta p_{nm} = p_{nm}\sqrt{\frac{\delta B_{nm}}{B}\sqrt{1-\xi^2}}.\]

The model used for the resonance can be controlled using the following options:

Name	Model for the resonance width
RIPPLE_MODE_NEGLECT	The ripple resonance is ignored altogether.
RIPPLE_MODE_BOX	Constant in \(\left|p_\parallel - p_{nm} \right| < \Delta p_{nm}\).
RIPPLE_MODE_GAUSSIAN	Gaussian envelope function

Input data

To include the magnetic ripple pitch scattering in simulation, it is necessary to provide the number of magnetic field coils \(N_{\rm c}\), as well as the perturbation magnitude \(\delta B_{nm}/B\) at one or more mode numbers \((n,m)\) (toroidal, poloidal). The magnetic perturbations can be specified as functions of time and radius, and will be interpolated as necessary.

In cylindrical (radial grid) mode, the function \(G(r)\) becomes infinite. However, this is compensated by the fact that a large-aspect ratio tokamak would have an infinite number of toroidal field coils, \(N_{\rm c}\). To still simulate the ripple in a cylindrical field, we therefore introduce the distance between magnetic field coils

\[\Delta_{\rm coils} = \frac{2\pi R_0}{N_{\rm c}}.\]

The resonant momentum then takes the form

\[p_{nm} = \frac{e G(r) \Delta_{\rm coils}}{2\pi m_ec n R_0},\]

with \(G(r)/R_0\) forming the finite toroidal magnetic field stregth.

Resonant momentum

When including the ripple pitch scattering in a simulation, it is possible to save the resonance momentum \(p_{nm}\) calculated at each time step, and for each pair of mode numbers, as an other quantity. The available ripple quantities are

Name of OtherQuantity	Description
fluid/ripple_m	Poloidal mode numbers corresponding to each resonance.
fluid/ripple_n	Toroidal mode numbers corresponding to each resonance.
fluid/f_hot_ripple_pmn	Resonant momentum for f_hot.
fluid/f_re_ripple_pmn	Resonant momentum for f_re.

Example

This example shows how to include the magnetic ripple pitch scattering mechanism in a simulation of the hot electron distribution function:

import DREAM.Settings.Equations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
# Number of toroidal field coils
nCoils = 16

# Include three resonances
dB_B   = [1e-4,5e-5,2e-5]
m      = [1,1,1]
n      = [1,2,3]

# Apply settings
ds.radialgrid.setRipple(m=m, n=n, dB_B=dB_B, ncoils=nCoils)
ds.eqsys.f_hot.setRippleMode(DistFunc.RIPPLE_MODE_BOX)

In cylindrical grid mode, the parameter deltaCoil should be specified instead of ncoils.

The arrays m and n should always have the same number of elements, which should also be the same number of elements as in the first dimension of dB_B.

To prescribe a set of spatiotemporally varying perturbations, replace the perturbations in the example above with 2D arrays and specify time and radial grid vectors:

import numpy as np

ds = DREAMSettings()

t = np.linspace(0, 1, 10)
r = np.linspace(0, 0.5, 20)

dB_B_11 = np.array([...])
dB_B_21 = np.array([...])
dB_B_31 = np.array([...])

dB_B = [dB_B_11, dB_B_21, dB_B_31]
m    = [1,1,1]
n    = [1,2,3]

ds.runawaygrid.setRipple(m=m, n=n, dB_B=dB_B, r=r, t=t, ncoils=nCoils)
ds.eqsys.f_hot.setRippleMode(DistFunc.RIPPLE_MODE_GAUSSIAN)

Approximate ion jacobian

When running DREAM with a self-consistent ion charge state evolution, the Newton-solver jacobian of the collision operator with respect to these densities will produce a large number of non-zero elements due to the sometimes large number of ion states. For example, with argon impurities in a deuterium plasma there are 21 ion densities, each giving a full distribution-grid contribution to the jacobian matrix.

Since the partially screened collision operator depends only logarithmically on ion charge state in the relativistic limit, a good approximation may often be obtained by neglecting the ion contribution to the jacobian. DREAM supports such an approximation using the fullIonJacobian setting, which can be True or False.

Example

In the below example we disable the ion jacobian only in the runaway distribution.

ds = DREAMSettings()

...

ds.eqsys.f_re.enableIonJacobian(False)

Note

The approximate ion jacobian can reduce the number of non-zero elements in the jacobian by more than a factor of two, and produce performance gains of orders-of-magnitude in certain situations. However, as with all approximations to the jacobian, there is a risk of deteriorated stability.

Prescribed distribution function

It is possible to run DREAM with a prescribed distribution function, which can be useful for studying the effect of different distribution function moments on the fluid equations (e.g. the fast electron impact ionization on the ion densities). The distribution function can be prescribed in time, space, momentum and pitch. Only one function call is needed to prescribe the distribution function:

ds.eqsys.f_re.prescribe(f, t, r, xi, p)

It is also possible to prescribe the distribution function in cylindrical momentum coordinates, i.e.:

ds.eqsys.f_re.prescribe(f, t, r, pperp=pperp, ppar=ppar)

Note

Prescribing the distribution function is currently only possible for the runaway electron distribution function. The hot electron distribution gives rise to additional moments which would require additional care to define properly in prescribed mode.

Example

from DREAM.Formulas.Distribution import getAvalancheDistribution

ds = DREAMSettings()
...
# Grid definition
np = 100
nxi = 40
pMin, pMax = 1, 100

ds.runawaygrid.setEnabled(True)
ds.runawaygrid.setNp(np)
ds.runawaygrid.setNxi(nxi)
ds.runawaygrid.setPmin(pMin)
ds.runawaygrid.setPmax(pMax)

# Distribution function
f   = np.zeros((nt, nr, nxi, np))
pp  = np.linspace(pMin, pMax, np+1)
xip = np.linspace(-1, 1, nxi+1)
p   = 0.5 * (pp[1:] + pp[:-1])
xi  = 0.5 * (xip[1:] + xip[:-1])

# Avalanche distribution parameters
E = 2    # E/Ec
Ztot = 8 # Total plasma charge (c.f. [Embreus JPP 84 (2018)])
nre = 1e16

f[0,0,:] = getAvalancheDistribution(p=p, xi=xi, E=E, Z=Ztot, nre=nre)

ds.eqsys.f_re.prescribe(f=f, t=[0], r=[0], xi=xi, p=p)

Radial Transport

A radial transport term can be added to the Fokker-Planck equation. The general transport term takes the form

\[\nabla\cdot\boldsymbol{\Phi}_{\rm transport} = \frac{1}{\mathcal{V}'}\frac{\partial}{\partial r}\left[ \mathcal{V}'\left( A f + D\frac{\partial f}{\partial r} \right) \right],\]

where \(A\) and \(D\) are the advection and diffusion coefficients respectively. DREAM provides two different ways of prescribing these coefficients:

1. Let the user prescribe \(A\) and/or \(D\) directly.

2. Use the model derived by Rechester and Rosenbluth and let the user prescribe the magnetic perturbation level \(\delta B/B\).

Prescribed coefficients

The user can directly prescribe the coefficients \(A=A(t,r,p,\xi)\) and \(D=D(t,r,p,\xi)\) by calling the methods transport.prescribeAdvection() and transport.prescribeDiffusion(). The coefficients can be either scalars (in which case they are assumed to be constant and uniform in all parameters) or 4D arrays, representing the coefficient’s dependence on time, radius, momentum and pitch, in that order.

Example

The following example shows how to prescribe advective-diffusive transport for the hot electron distribution function:

ds = DREAMSettings()
...
t  = np.linspace(t0, t1, nt)
r  = np.linspace(r0, r1, nr)
p  = np.linspace(p0, p1, np)
xi = np.linspace(xi0, xi1, nxi)

A = np.array([...])  # shape (nt, nr, np, nxi)
D = np.array([...])  # shape (nt, nr, np, nxi)

ds.eqsys.f_hot.transport.prescribeAdvection(ar=A, t=t, r=r, p=p, xi=xi)
ds.eqsys.f_hot.transport.prescribeDiffusion(drr=D, t=t, r=r, p=p, xi=xi)

Contrary to what the example may seem to suggest, the coefficients do not need to share the same time, radius, momentum and pitch grids. The grids are also not required to be the same as those used internally by DREAM for the simulation, as the provided coefficients will be interpolated as needed.

If a constant coefficient is desired, the following more compact syntax may be used:

ds.eqsys.f_hot.transport.prescribeDiffusion(10)

which will prescribe a constant diffusion coefficient with the value \(10\,\mathrm{m}^2/\mathrm{s}\).

Rechester-Rosenbluth

Rechester and Rosenbluth derived a diffusion operator for describing the radial transport due to turbulent perturbations of the magnetic field. Given a magnetic field perturbation strength \(\delta B/B\), the diffusion coefficient for the Rechester-Rosenbluth operator takes the form

\[D = \pi q R_0\left(\frac{\delta B}{B}\right)^2 \left|v_\parallel\right|,\]

where \(q(r)\) is the tokamak safety factor, \(R_0\) is the tokamak major radius and \(v_\parallel\) is the particle speed along the magnetic field line. In DREAM, the user provides the magnetic perturbation \(\delta B/B\) as a function of time and radius.

Note

For simplicity in this equation term, DREAM does not calculate the normalized safety factor \(q(r)R_0\) dynamically from the current profile. Instead, in the implementation of this operator it is set to \(qR_0=1\,\mathrm{m}\), effectively requiring the user to provide this dependence in the input parameter \(\delta B/B\).

Example

The following example illustrates how to prescribe a Rechester-Rosenbluth diffusion coefficient in DREAM:

ds = DREAMSettings()
...
# Grids on which the perturbation is defined
t = np.linspace(t0, t1, nt)
r = np.linspace(r0, r1, nr)

# Magnetic perturbation
dB_B = np.array([...])   # shape (nt, nr)

ds.eqsys.f_hot.transport.setMagneticPerturbation(dBB=dB_B, t=t, r=r)

If a constant and uniform perturbation level is desired, the more compact syntax

ds.eqsys.f_hot.transport.setMagneticPerturbation(dBB=1e-3)

can also be used. In this example, the magnetic perturbation \(\delta B/B = 10^{-3}\) is used in every time step and at all radial positions.

Frozen current mode

In the frozen current mode, a target plasma current \(I_{\rm p}\) is prescribed by the user and an associated diffusive radial transport with coefficient \(D_0\) is introduced. In each time step, DREAM will then adjust \(D_0\) such that \(I_{\rm p}\) exactly matches its prescribed value. This is particularly useful for simulations of experiments, where the plasma current is often accurately known, while the radial transport is generally poorly constrained by experimental measurements. This feature is inspired by a similar feature in the kinetic solver LUKE.

Frozen current mode can be enabled on any transportable (particle) quantity by calling setFrozenCurrentMode(). The call takes two mandatory parameters: (i) the type of transport to model, and (ii) the plasma current target. Additionally, two optional parameters may be specified: (iii) the time vector corresponding to the prescribed plasma current, if the latter is to vary in time, and (iv) the maximum permitted value of \(D_0\) (default: 1000 m/s^2).

If the prescribed plasma current is higher than what can be achieved with a diffusive transport (i.e. if \(I_{\rm p}\) is below its prescribed value in the absence of radial transport), the diffusion coefficient is set to zero and the target plasma current is ignored.

List of options

The radial diffusion coefficient used in the frozen current mode is generally written on the form

\[D_{rr}(r,p,\xi_0) = D_0h(r,p,\xi_0)\]

where the function \(h(r,p,\xi_0)\) can take different forms. The following forms for \(h(r,p,\xi_0)\) are currently supported in DREAM:

Name	\(h(r,p,\xi_0)\)
FROZEN_CURRENT_MODE_DISABLED	\(0\)
FROZEN_CURRENT_MODE_CONSTANT	\(1\)
FROZEN_CURRENT_MODE_BETAPAR	\(\beta_\parallel = v_\parallel/c\)

Example

The following example illustrates how to use the frozen current mode:

from DREAM import DREAMSettings
import DREAM.Settings.Transport as Transport

ds = DREAMSettings()
...
Ip = 800e3   # Plasma current target (A)
ds.eqsys.f_hot.transport.setFrozenCurrentMode(
    Transport.FROZEN_CURRENT_MODE_BETAPAR,
    Ip_presc=Ip
)

If one wishes to prescribe a time-evolving plasma current, the following call can be made instead:

from DREAM import DREAMSettings
import DREAM.Settings.Transport as Transport

ds = DREAMSettings()
...
# Time vector
t  = np.linspace(0, 2, 40)
# Plasma current target (A)
Ip = 800e3*np.linspace(0.2, 1, t.size)
ds.eqsys.f_hot.transport.setFrozenCurrentMode(
    Transport.FROZEN_CURRENT_MODE_BETAPAR,
    Ip_presc=Ip, Ip_presc_t=t
)

In poorly confined plasmas, or situations where strong transport leads to unstable simulations, one may want to adjust the maximum diffusion coefficient:

from DREAM import DREAMSettings
import DREAM.Settings.Transport as Transport

ds = DREAMSettings()
...
Ip = 800e3   # Plasma current target (A)
ds.eqsys.f_hot.transport.setFrozenCurrentMode(
    Transport.FROZEN_CURRENT_MODE_BETAPAR,
    Ip_presc=Ip, D_I_max=200
)

Boundary conditions

The general transport term requires a boundary condition to specified at the plasma boundary \(r=a\). DREAM currently supports two different boundary conditions at this boundary:

1. Assume that no particles can cross the plasma boundary (default).

2. Assume that \(f(r>a)=0\) and that fluxes are set accordingly.

Condition (1) is a reflective boundary condition which naturally conserves the total number of particles in \(f\). It prevents particles from leaving the plasma.

Condition (2) generally results in a net outflux of particles, and allows particles to be lost from the system. However, to conserve quasi-neutrality, DREAM compensates for the lost electrons by adding additional cold electrons.

List of options

The following boundary conditions can be applied:

Name	Description
BC_CONSERVATIVE	Assume that no particles can leave the plasma through \(r=a\). (default)
BC_F_0	Assume that \(f(r>a) = 0\).
BC_DF_CONST	Assume that \(\partial f/\partial r = \mathrm{const.}\) on plasma edge.

Example

The following example illustrates how to apply a lossy boundary condition to the hot electron distribution function:

import DREAM.Settings.TransportSettings as Transport

ds = DREAMSettings()
...
ds.eqsys.f_hot.transport.setBoundaryCondition(Transport.BC_F_0)

Output

When applying transport to the kinetic equation, two diagnostic quantities are saved as an other quantity for each distribution function:

Name of OtherQuantity	Description
scalar/radialloss_f_re	Density moment: Rate at which runaways leave the plasma due to transport
scalar/energyloss_f_re	Energy moment: Rate at which runaway energy leaves the plasma due to transport
scalar/radialloss_f_hot	Density moment: Rate at which hot electrons leave the plasma due to transport
scalar/energyloss_f_hot	Energy moment: Rate at which hot electrons leave the plasma due to transport

Runaway avalanching

Secondary runaway electrons (produced through the so-called avalanche mechanism) can be included in the simulation of the distribution function through a Rosenbluth-Putvinski source term (Rosenbluth and Putvinski 1997). This source term is however enabled on the n_re settings object using the setAvalanche() method with the option RunawayElectrons.AVALANCHE_MODE_KINETIC.

Since the source term diverges in the origin, creating an infinite number of knock-ons, a cutoff parameter pCutAvalanche is introduced, defined as the minimum momentum to which the source is applied.

Warning

The momentum cutoff must be chosen to be smaller than the critical momentum for runaway generation, otherwise the cutoff will act as the effective critical momentum, thus artificially limiting the avalanche growth rate.

Example

Kinetic avalanche runaways can be included in the simulation thusly:

import DREAM.Settings.Equations.RunawayElectrons as Runaways

ds = DREAMSettings()
...
ds.eqsys.n_re.setAvalanche(Runaways.AVALANCHE_MODE_KINETIC, pCutAvalanche=0.1)

This will add the Rosenbluth-Putvinski source term to all kinetic equations in the simulation.

Note

When activating kinetic avalanche, the source integrated over momenta outside the hottail grid, i.e. for all \(p>p_\mathrm{max}\), will be added directly to n_re. This way, the total number of knock-ons created is independent of grid parameters.

Synchrotron losses

Synchrotron radiation provides one of the most important energy loss mechanisms for runaway electrons. DREAM implements the model of Decker et al (2016), which uses the radiation time scale

\[\tau_{\rm rad}^{-1} = \frac{e^4B^2}{6\pi\epsilon_0 m_e^3c^3}.\]

Synchrotron losses are disabled by default.

List of options

Name	Description
SYNCHROTRON_MODE_INCLUDE	Include synchrotron losses in simulation
SYNCHROTRON_MODE_NEGLECT	Do not include synchrotron losses in simulation

Example

Synchrotron losses can be enabled via a call to setSynchrotronMode():

import DREAM.Settings.Equations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
ds.eqsys.f_re.setSynchrotronMode(DistFunc.SYNCHROTRON_MODE_INCLUDE)

# To explicitly disable, call
#ds.eqsys.f_re.setSynchrotronMode(DistFunc.SYNCHROTRON_MODE_NEGLECT)

Alternatively, a bool may be used to enable/disable synchrotron losses:

ds = DREAMSettings()
...
ds.eqsys.f_re.setSynchrotronMode(True)

# To explicitly disable, call
#ds.eqsys.f_re.setSynchrotronMode(False)

Time-varying magnetic field

Due to the conservation of magnetic moment, a time-varying magnetic field will introduce an effective force acting on the electrons which transfers momentum from the parallel to the perpendicular direction (or oppositie direction, depending on whether the magnetic field increases or decreases). In terms of the magnetic moment \(\mu=mc^2p_\perp^2/2B\), the normalized momenta are

\[\begin{split}p_\parallel &= \sqrt{p^2-\left(\frac{2\mu B}{mc^2}\right)^2},\\ p_\perp &= \sqrt{\frac{2\mu B}{mc^2}},\end{split}\]

from which the resulting advection coefficient can be derived:

(1)\[A^\xi = -\frac{mc}{2B}\frac{1-\xi^2}{\xi}\frac{\mathrm{d}B}{\mathrm{d}t}.\]

Since, at the time of writing (and implementing this term), DREAM does not support actual time-varying magnetic fields, this term only emulates the effect of having a time-varying magnetic field. As such, we replace \(\mathrm{d}B/\mathrm{d}t\) with the simpler parameter \(\mathrm{d}B_0/\mathrm{d}t\), which corresponds to the time variation of the on-axis magnetic field strength and is exact in cylindrical geometry where \(B(R)=B_0R_0/R\) holds. The bounce-averaged form of (1) can then be written

\[\left\{ A^\xi \right\} = -\frac{mc}{2B}\frac{\mathrm{d}B_0}{\mathrm{d}t} \left\{\frac{1-\xi^2}{\xi}\right\},\]

where curly braces denote a bounce-average.

Warning

This operator will not actually change the magnetic field strength used by DREAM — it only emulates the effective force felt by electrons due to the time-varying magnetic field.

Example

from DREAM import DREAMSettings
import DREAM.Settings.Eqyations.DistributionFunction as DistFunc

ds = DREAMSettings()
...
# On-axis magnetic field strength (T)
B0 = 5.2
# Time-derivative of magnetic field strength (T/s)
dB0dt = 3.0

ds.eqsys.f_hot.setTimeVaryingB(True)
ds.radialgrid.setTimeVaryingB(dB0dt_B0=dB0dt/B0)

Alternatively, the time-derivative can be made to vary in time itself:

from DREAM import DREAMSettings
import DREAM.Settings.Eqyations.DistributionFunction as DistFunc
import numpy as np

ds = DREAMSettings()
...
# Maximum time to give dB/dt for (s)
tMax = 2.0
# On-axis magnetic field strength (T)
B0 = 5.2
# Time-derivative of magnetic field strength (T/s)
dB0dt = lambda t : 3.0  + 0.5*tMax
# Time points to give dB/dt in
t = np.linspace(0, tMax)

ds.eqsys.f_hot.setTimeVaryingB(True)
ds.radialgrid.setTimeVaryingB(dB0dt_B0=dB0dt(t)/B0)

Analytical distribution

DREAM supports the option of replacing the numerical kinetic description of the hot and runaway electrons by an approximate analytical model for their distribution functions. This is controlled by the setting

Name	Description
DISTRIBUTION_MODE_NUMERICAL	The distribution is resolved on a kinetic (finite-volume) grid
DISTRIBUTION_MODE_ANALYTICAL	The distribution is modelled with an analytical function

Typically, in the analytical mode, f_hot is described by a hot-tail slowing-down distribution and f_re by a quasi-steady state avalanche distribution.

Example

ds = DREAMSettings()

# disable kinetic grid
ds.hottailgrid.setEnabled(False)

# enable the analytical distribution function
ds.f_hot.enableAnalyticalDistribution()

...

Tips

Note

Are you experiencing unstable solutions? Try to apply one of the alternative advection term interpolation schemes available in DREAM. These can efficiently suppress numerical oscillations and stabilize solutions. If the non-linear solver fails to converge, increasing resolution sometimes helps, in particular decreasing the time step. If a floating-point exception is thrown, the matrix may have become singular to working precision, which is often resolved by using a different linear solver, or increasing grid and/or time resolution.

Class documentation

class DREAM.Settings.Equations.DistributionFunction.DistributionFunction(settings, name, grid, f=[0], initr=[0], initp=[0], initxi=[0], initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=2, ad_jac_p1=2, ad_jac_p2=2, mode=1, fluxlimiterdamping=1.0)

Bases: DREAM.Settings.Equations.UnknownQuantity.UnknownQuantity

__init__(settings, name, grid, f=[0], initr=[0], initp=[0], initxi=[0], initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=2, ad_jac_p1=2, ad_jac_p2=2, mode=1, fluxlimiterdamping=1.0)

Constructor.

enableAnalyticalDistribution(mode=True)

Enables/disables the use of an analytical distribution function to represent the electron population

enableIonJacobian(includeJacobian)

Enables/disables the ion jacobian in the kinetic equation.

Parameters

includeJacobian (bool) – Flag indicating whether the ion jacobian will be added. True by default, False to disable.

fromdict(data)

Load data for this object from the given dictionary.

Parameters

data (dict) – Dictionary to load distribution function from.

prescribe(f, t, r, xi=None, p=None, pperp=None, ppar=None)

Prescribe the time evolution of this distribution function instead of solving for it using a kinetic equation.

setAdvectionInterpolationMethod(ad_int=None, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac=None, ad_jac_r=2, ad_jac_p1=2, ad_jac_p2=2, fluxlimiterdamping=1.0)

Sets the interpolation method that is used in the advection terms of the kinetic equation. To set all three components, provide ad_int and/or ad_jac. Otherwise the three components can use separate interpolation methods.

Parameters

• ad_int (int) – Interpolation method to use for all coordinates.

• ad_int_r (int) – Interpolation method to use for the radial coordinate.

• ad_int_p1 (int) – Interpolation method to use for the first momentum coordinate.

• ad_int_p2 (int) – Interpolation method to use for the second momentum coordinate.

• ad_jac (int) – Jacobian interpolation mode to use for all coordinates.

• ad_jac_r (int) – Jacobian interpolation mode to use for the radial coordinate.

• ad_jac_p1 (int) – Jacobian interpolation mode to use for the first momentum coordinate.

• ad_jac_p2 (int) – Jacobian interpolation mode to use for the second momentum coordinate.

• fluxlimiterdamping (float) – Damping parameter used to under-relax the interpolation coefficients during non-linear iterations (should be between 0 and 1).

setBoundaryCondition(bc)

Sets the boundary condition at p=pmax. For ‘f_hot’, this boundary condition is only used when ‘f_re’ is disabled.

Parameters

bc (int) – Flag specifying which boundary condition to use.

setInitialProfiles(n0, T0, rn0=None, rT0=None)

Sets the initial density and temperature profiles of the electron population.

Parameters

• rn0 – Radial grid on which the density is given.

• n0 – Electron density profile.

• rT0 – Radial grid on which the temperature is given.

• T0 – Electron temperature profile.

setInitialValue(f, r, p=None, xi=None, ppar=None, pperp=None)

Set the initial value of this electron distribution function. Only one of the pairs (p, xi) and (ppar, pperp) of momentum grids need to be given.

Parameters

• f – Array representing the distribution function value on the grid (must have size (nr, nxi, np) or (nr, npperp, nppar))

• r – Radial grid on which the initial distribution is given.

• p – Momentum grid.

• xi – Pitch grid.

• ppar – Parallel momentum grid.

• pperp – Perpendicular momentum grid.

setRippleMode(mode)

Enables/disables inclusion of pitch scattering due to the magnetic ripple.

Parameters

mode (int) – Flag indicating whether or not to include magnetic ripple effects.

setSynchrotronMode(mode)

Sets the type of synchrotron losses to have (either enabled or disabled).

Parameters

mode (int) – Flag indicating whether or not to enable synchrotron losses (may be bool).

setTimeVaryingB(mode)

Enables/disable the time-varying magnetic field strength operator.

Parameters

mode (int) – Flag indicating whether or not to include the time-varying magnetic field operator.

todict()

Returns a Python dictionary containing all settings of this DistributionFunction object.

Returns

a dictionary, containing all settings of this object, which can be directly given to DREAM.

verifyInitialDistribution()

Verifies that the initial distribution function has been set correctly and consistently.

verifyInitialProfiles()

Verifies that the initial density and temperature profiles are set correctly.

verifySettings()

Verify that the settings of this unknown are correctly set.

HotElectronDistribution

The HotElectronDistribution class contains settings specific to the hot electron distribution, f_hot. It contains the following methods:

Momentum threshold

When resolving the thermal population on the grid, the question arises how to split the electron population into cold (which enter into the temperature evolution) and hot electrons (corresponding to the remainder). This is relevant to the split of current into j_ohm and j_hot, density into n_cold and n_hot, and also in kinetic equation terms where the contribution from hot electrons are added as an integral moment of the distribution.

In DREAM, this is handled using the pThreshold parameter, which can be set using three different modes:

Name	Definition of hot electrons
HOT_REGION_P_MODE_MC	\(\int_\mathrm{pCutoff}^\mathrm{pMax} \mathrm{d}p\)
HOT_REGION_P_MODE_THERMAL	\(\int_{\mathrm{pCutoff}\sqrt{2T_\mathrm{cold}/m_e c^2}}^\mathrm{pMax} \mathrm{d}p\)
HOT_REGION_P_MODE_THERMAL_SMOOTH	\(\int_0^\mathrm{pMax} \mathrm{d}p ~h(p,\,T_\mathrm{cold})\)

Note

The default settings are HOT_REGION_P_MODE_THERMAL with pCutoff = 10, corresponding to a sharp cutoff at ten thermal momenta.

All modes assume an isotropic definition of the hot electrons, where the hot region is independent of angle. Mode MC corresponds to a sharp cutoff with pThreshold the momentum cutoff in units of \(m_e c\), and THERMAL gives it in units of the thermal momentum. In SMOOTH mode, the envelope function is defined as:

\[h(p,\,T_\mathrm{cold}) = \tanh\left[\left(p-\mathrm{pCutoff}\sqrt{2T_\mathrm{cold}/m_e c^2}\right)/\Delta p\right]\]

This smooth envelope function allows the jacobian of the momentum region with respect to the cold electron temperature to be evaluated, which may improve the rate of convergence in the non-linear solver.

Warning

The step width \(\Delta p\) in SMOOTH mode is currently defined as the momentum width of the cell containing the cutoff. It is due for a rework, as a failure mode presently occurs when using non-uniform momentum grids with a large jump in step width near to the thermal population, in which case a significant fraction of thermal electrons will be counted as hot.

Example

The following example shows how to change hot-electron definition to a smooth envelope function located at an elevated cutoff momentum:

import DREAM.Settings.Equations.HotElectronDistribution as FHot

ds = DREAMSettings()
ds.eqsys.f_hot.setHotRegionThreshold(pThreshold=20, pMode=FHot.HOT_REGION_P_MODE_THERMAL_SMOOTH)

This could for example be employed in an attempt to stabilize a hypothetical unstable self-consistent simulation with a rapid temperature drop when the thermal Maxwellian is resolved on-grid.

Particle source

When running with COLLISION_FREQUENCY_MODE_FULL, the entire electron population is followed on the kinetic grid. A particle source term \(S = \mathrm{S\_particle} * S_0(p)\) is added to the hottail kinetic equation which ensures that the correct density evolution is followed; although the kinetic equation is locally density conserving, such a source term is needed in order to compensate for changes in the ion charge states by ionization and recombination, and also for the radial transport term which is typically constructed such that it causes a net flow of charge which must be compensated for.

Currently, the source function shape \(S_0(p)\) is chosen as the Maxwell-Jüttner distribution evaluated at the instantaneous cold-electron temperature. It is normalized such that its density moment equals unity, so that the unknown quantity S_particle represents the net rate at which density is added to (or removed from) the distribution.

The particle source is controlled by the following settings:

Name	Model for S_particle
PARTICLE_SOURCE_ZERO	\(\mathrm{S\_particle} = 0\) (the particle source is neglected)
PARTICLE_SOURCE_IMPLICIT	\(\int f_\mathrm{hot} \mathrm{d}p = n_\mathrm{cold} + n_\mathrm{hot}\) (S_particle is such that f_hot has the correct density moment)
PARTICLE_SOURCE_EXPLICIT	\(\mathrm{S\_particle} = \Big(\mathrm{d}n_\mathrm{free}/\mathrm{d}t\Big)_\mathrm{ions} + \text{[other sources of density]}\)

Here, ZERO deactivates the particle source and allows the density of the hot distribution function to deviate from n_cold + n_hot. IMPLICIT and EXPLICIT yield exactly the same behaviour, but due to the different formulations yield different condition numbers for the matrix.

Example

The model used for the particle source can be set by:

import DREAM.Settings.Equations.HotElectronDistribution as FHot

ds = DREAMSettings()

ds.eqsys.f_hot.setParticleSource(particleSource=FHot.PARTICLE_SOURCE_IMPLICIT):

Note

If the equation system becomes ill-conditioned when running a self-consistent simulation with COLLISION_FREQUENCY_MODE_FULL, try changing particle source mode, which may make the equation system more well-conditioned.

Analytical distribution type

Different models may be implemented to describe the hot-electron distribution. These are controlled via the following options:

Name	Model for f_hot distribution
F_HOT_DIST_MODE_NONREL	The isotropic non-relativistic slowing-down distribution of Eq (9) in Smith (2008)

Presently, only one model is implemented which is used by default, so this setting does not need to be manually set. Note, that F_HOT_DIST_MODE_NONREL depends on a normalized-time parameter tau_coll, which is added to the equation system as a new UnknownQuantity whenever the ANALYTIC mode is used for f_hot combined with this distribution type.

Class documentation

class DREAM.Settings.Equations.HotElectronDistribution.HotElectronDistribution(settings, fhot=None, initr=None, initp=None, initxi=None, initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=2, ad_jac_p1=2, ad_jac_p2=2, fluxlimiterdamping=1.0, dist_mode=1, pThreshold=7, pThresholdMode=2, particleSource=3, particleSourceShape=1)

Bases: DREAM.Settings.Equations.DistributionFunction.DistributionFunction

__init__(settings, fhot=None, initr=None, initp=None, initxi=None, initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=2, ad_jac_p1=2, ad_jac_p2=2, fluxlimiterdamping=1.0, dist_mode=1, pThreshold=7, pThresholdMode=2, particleSource=3, particleSourceShape=1)

Constructor.

fromdict(data)

Load data for this object from the given dictionary.

setHotRegionThreshold(pThreshold=7, pMode=2)

Sets the boundary ‘pThreshold’ which defines the cutoff separating ‘cold’ from ‘hot’ electrons when using collfreq_mode FULL.

setParticleSource(particleSource=3, shape=1)

Sets which model to use for S_particle if using collfreq_mode FULL, which is designed to force the density moment of f_hot to n_cold+n_hot.

ZERO: The particle source is disabled and set to zero EXPLICIT/IMPLICIT: Two in principle equivalent models, but can be more or less stable in different situations.

todict()

Returns a Python dictionary containing all settings of this HotElectronDistribution object.

class DREAM.Settings.Equations.RunawayElectronDistribution.RunawayElectronDistribution(settings, fre=[0.0], initr=[0.0], initp=[0.0], initxi=[0.0], initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=1, ad_jac_p1=1, ad_jac_p2=1, fluxlimiterdamping=1.0)

Bases: DREAM.Settings.Equations.DistributionFunction.DistributionFunction

__init__(settings, fre=[0.0], initr=[0.0], initp=[0.0], initxi=[0.0], initppar=None, initpperp=None, rn0=None, n0=None, rT0=None, T0=None, bc=2, ad_int_r=1, ad_int_p1=1, ad_int_p2=1, ad_jac_r=1, ad_jac_p1=1, ad_jac_p2=1, fluxlimiterdamping=1.0)

Constructor.

fromdict(data)

Load data for this object from the given dictionary.

setInitType(inittype)

Specifies how the runaway electron distribution function f_re should be initialized from the runaway density n_re.

Parameters

inittype (int) – Flag indicating how to initialize f_re.

setInitialValue(f, r, p=None, xi=None, ppar=None, pperp=None)

Set the initial value of this electron distribution function. Only one of the pairs (p, xi) and (ppar, pperp) of momentum grids need to be given.

Parameters

• f – Array representing the distribution function value on the grid (must have size (nr, nxi, np) or (nr, npperp, nppar))

• r – Radial grid on which the initial distribution is given.

• p – Momentum grid.

• xi – Pitch grid.

• ppar – Parallel momentum grid.

• pperp – Perpendicular momentum grid.

todict()

Returns a Python dictionary containing all settings of this RunawayElectronDistribution object.