The RunawayElectons
object, located under .eqsys.n_re
in the
DREAMSettings object is used to specify settings for the fluid quantity
n_re
, representing the density of runaway electrons in simulations. The
quantity is always present in DREAM simulations and is evolved using the
equation
where
\(\Phi^{(p)}_{\rm hot}\) is the flux of particles from the (kinetic) hottail grid.
\(\Gamma_{\rm Ava}\) is the avalanche growth rate.
\(\gamma_{\rm Dreicer}\) is the runaway rate due to the Dreicer mechanism.
\(\gamma_{\rm hottail}\) is the runaway rate due to the hottail mechanism.
\(\gamma_{\rm tritium}\) is the runaway rate due to tritium decay.
\(\gamma_{\rm Compton}\) is the runaway rate due to Compton scattering.
and the last term represents advectivediffusive radial transport. Each term can be individually enabled/disabled using the methods described further down on this page.
Note
The Dreicer generation rate \(\gamma_{\rm Dreicer}\) should typically be disabled in simulations including the hottail grid where the thermal population is resolved kinetically. This is because the Dreicer generation will then be modelled using a kinetic equation and the Dreicer generation rate become a part of the hot electron flux \(\Phi^{(p)}_{\rm hot}\).
Note
Note that the runaway density \(n_{\rm RE}\) is never defined in terms of the runaway electron distribution function \(f_{\rm RE}\) in DREAM. The quantities are instead evolved separately and set up in such a way that their local densities should agree in the end if the simulation was set up in a consistent way.
It is possible—and perfectly okay—for \(n_{\rm RE}\) to differ from the density moment of \(f_{\rm RE}\) if the latter loses a significant number of particles at \(p=p_{\rm max}\).
Page overview
By default, the runaway density is set to zero at the beginning of a DREAM
simulation. It is however possible to start the simulation with a nonzero
runaway density profile, and this is done using the setInitialProfile()
method.
Warning
If the runaway grid is enabled, one should also make sure to appropriately initialize the runaway electron distribution function when setting a nonzero initial runaway density. The density moment of the runaway electron distribution function \(f_{\rm RE}\) should always agree with the runaway electron density \(n_{\rm RE}\).
The following example illustrates how to initialize the runaway electron density:
import numpy as np
ds = DREAMSettings()
...
# Define radial grid for the initial profile
r = np.linspace(r0, r1, nr)
# Generate the runaway electron density profile
nre = np.array([...]) # shape (nr,)
ds.eqsys.n_re.setInitialProfile(density=nre, r=r)
If a uniform initial density profile is desired one can specify the density
parameter as a scalar value instead:
import numpy as np
ds = DREAMSettings()
...
ds.eqsys.n_re.setInitialProfile(density=2e16)
DREAM includes a number of common runaway electron source terms for use in both pure fluid as well as combined fluidkinetic simulations.
Runaway electrons can produce new runaway electrons by colliding with thermal electrons and transferring a sufficent amount of energy to these for them to become runaway electrons, while the original runaway electron remains in the runaway region. This runaway mechanism is commonly known as the avalanche mechanism and will lead to an exponential growth in the number of runaway electrons.
DREAM contains three different models for avalanche generation. The first two are fluid models, while the third is the kinetic source term originally derived by Rosenbluth and Putvinski (1997).
Note
Note that also fluid growth rates can be used in simulations involving the runaway grid. In these simulations, the runaways produced with the fluid source term are added to the \(\xi=1\) cell on the runaway grid (if \(E > 0\), otherwise to the \(\xi=1\) cell).
Similarly, the kinetic source term can be used in fluid simulations and is then integrated to yield the number of fluid runaways produced.
Option 
Description 


Do not include avalanche generation in the simulation. 

Modified 

Implements Eq (14) in Hesslow NF (2019) 

Kinetic source term derived by Rosenbluth and Putvinski (1997). 
The FLUID
model modifies the FLUID_HESSLOW
model by replacing the denominator
\(\sqrt{4+\bar\nu_s\bar\nu_D}\) appearing in (14) of the Hesslow paper by
\(\sqrt{4\bar\nu_s+\bar\nu_s\bar\nu_D}\), which increases the accuracy for
nearly neutral plasmas dominated by hydrogen collisions.
Warning
The kinetic avalanche source, AVALANCHE_MODE_KINETIC
, makes explicit use
of the runaway density \(n_{\rm re}\). Hence, it is not possible to model
avalanche generation in simulations where the hottail grid is intended to
contain all, or a significant fraction of, the runaway electrons. If the
runaways should be resolved kinetically, the runaway grid should also be
enabled.
Note
Running the dream_avalanche
physics test with the –plot flag
compares the FLUID
and FLUID_HESSLOW
models with the
kinetic calculation.
Note
If you are using kinetic avalanche generation and the electric field changes sign during the simulation, read negative electric fields below and call ds.eqsys.n_re.setNegativeRunaways() on the settings object to properly account for the direction of motion of the runaways.
Todo
Describe the use of the pCutAvalanche
parameter in kinetic mode.
The first runaway mechanism to be discovered is the socalled Dreicer runaway mechanism which occurs in a plasma as soon as a sufficiently strong electric is applied. The electric field immediately accelerates all electrons with momentum \(p > p_{\rm c}\) to even higher momentum, leaving a “gap” in the electron distribution function. Since collisions seek to maintain the Maxwellian form for the distribution, the electrons soon reorganize to replace the accelerated electrons. However, the new electrons which diffused to fill in the gap now find themselves with a momentum \(p > p_{\rm c}\) and are thus also immediately accelerated to higher energy. This leads to gradual diffusive leak of particles into the runaway region.
Since the Dreicer mechanism is naturally obtained as the balance between collisions and electric field acceleration, the mechanism is automatically and always present in all kinetic simulations. To enable Dreicer generation in pure fluid simulations, the options described in the table below are available.
Option 
Description 


Do not include Dreicer generation in the simulation. 

Use the formula derived by Connor and Hastie (1975), excluding the relativistic corrections (valid in the limit \(E\gg E_{\rm c}\)). 

Use the full formula derived by Connor and Hastie (1975). 

Use the neural network constructed by Hesslow et al (2019). 
Runaway electrons can be generated when tritium decays into helium3 via the beta decay process
The corresponding runaway rate is given by
where \(n_{\rm T}\) is the tritium density, \(\tau_{\rm T} = 4800\pm 8\) days is the tritium halflife, and \(F_\beta(\gamma_{\rm c})\) denotes the fraction of beta electrons generated with an energy above the critical energy \(\gamma_{\rm c}\) for runaway to occur.
Tritium runaway generation is enabled with the setTritium()
method of the
RunawayElectrons
object, but it is necessary to also provide a tritium ion
species in the ions interface. The user must make
sure to specify tritium=True
when adding the tritium ion to the simulation
using the addIon()
method.
Note
It is possible to include multiple tritium populations in the simulation,
simply by providing tritium=True
when adding each of them.
The following example illustrates how to enable the tritium decay runaway mechanism in a DREAM simulation:
import DREAM.Settings.Equations.IonSpecies as Ions
ds = DREAMSettings()
...
# Include source term in equation for n_re
ds.eqsys.n_re.setTritium(True)
# Add tritium ion species to list of ions
ds.eqsys.ions.addIon('T', Z=1, iontype=Ions.IONS_DYNAMIC, n=2e19, tritium=True)
Runaway generation by Compton scattering is modelled via the KleinNishina differential cross section integrated over the runaway region \(p>p_c\), analogous to the tritium and avalanche calculations. This integrated cross section is given by Equation (29) in MartinSolis NF (2017). This total cross section \(\sigma(p_c,\,E_\gamma)\) depends on plasma parameters through \(p_c\). The net runaway rate is consequently obtained as the integral of the total cross section over the photon energy spectrum \(\Phi(E_\gamma)\) (where we assume that bound and free electrons are equally susceptible to Compton scattering):
Following MartinSolis Eq (24), we model the photon energy spectrum with
where \(z = [1.2 + \ln(E_\gamma[MeV])]/0.8\), and \(\Phi_0 = \int \Phi \,\mathrm{d}E_\gamma\) is the total photon gamma flux. For ITER, according to MartinSolis, this value is \(\Phi_0 \approx 10^{18}\,\mathrm{m}^{2}\mathrm{s}^{1}\) during the nuclear phase.
The following settings are used to control the Compton source mode:
Option 
Description 


Do not include Compton scattering in the simulation. 

Use the model described in MartinSolis NF (2017), with tuneable total photon flux. 
A fluid model for hottail generation is implemented based on the Master’s thesis of Ida Svenningsson (2020). DREAM mainly utilizes the theory of Section 4.2 in the thesis, using the same approximations as used in ‘ISOTROPIC’ (Nxi=1) kinetic mode.
The hottail runaway generation rate is given by
where \(p_c\) is the critical runaway momentum, \(\dot{p}_c\) its time rate of change and \(f_0\) is the angleaveraged electron distribution evaluated at the critical momentum. In this hottail model, the critical momentum is found from the angleaveraged kinetic equation in the Lorentz limit of strong pitchangle scattering:
which for a given distribution \(f_0\) is solved numerically for the critical momentum \(p=p_c\). Here we assume that only pitchangle scattering, electric field acceleration and collisional slowing down contributes to the electron dynamics. We have taken the nonrelativistic limit of the equation, and neglected screening corrections which requires a high effective charge in order to be accurate (typically satisfied during the thermal quench with highZ material injection).
The distribution function \(f_0\) is taken as the analytical hot electron distribution described in HotElectronDistribution.
The following hot tail settings are supported:
Option 
Description 


Do not include hottail generation in the simulation. 

TODO (the lowZ limit, Section 4.1.3 in Svenningsson’s thesis) 

The highZ limit from Section 4.2 in Svenningsson’s thesis 
The hottail generation can be activated if and only if f_hot
is in analytical
mode.
import DREAM.Settings.Equations.RunawayElectons as Runaways
ds = DREAMSettings()
# Set f_hot mode to analytical and provide initial profiles
ds.hottailgrid.setEnabled(False)
# rn, n0, rT, T0 = ... get profiles of density and temperature
ds.eqsys.f_hot.setInitialProfiles(rn0=rn, n0=n0, rT0=rT, T0=T0)
ds.eqsys.n_re.setHottail(Runaways.HOTTAIL_MODE_ANALYTIC_ALT_PC)
DREAM allows the effective critical electric field to be calculated following the method outlined in Hesslow PPCF (2018), which has here been generalized to the case of tokamak geometry (i.e. inhomogeneous magnetic fields). The algorithm DREAM implements is relatively involved, and therefore we briefly describe the method here.
Following Hesslow PPCF (2018) an analytical expression for the runaway pitch distribution can be obtained under the assumption that energy fluxes are negligible compared with pitch fluxes, which is expected to be valid near the maximum runaway momentum where energy losses balance the electric field acceleration.
This gives
where it is assumed that the dominant bounce averaged terms are the pitch component of the electric field term \(A_{\xi_0}\) and the pitchangle scattering coefficient \(D_{\xi_0\xi_0}\). The exponent can be written in the familar form
where in the cylindrical theory \(A = 2eE/(p\nu_D)\) and \(g = 1\xi_0\), but in the general case we have
and \(1\xi_0\) is replaced by the function
where we introduce the auxiliary function
By integrating the kinetic equation over \(\mathrm{d}\xi_0 \mathcal{V}'/V'\), it takes the form
where the distributionbounce averaged momentum flow (i.e. net acceleration) \(U(p)\) is
In the DREAM calculation we account for the contributions from electric field acceleration, synchrotron losses
and slowing down (collisional with screening effect plus meanforce bremsstrahlung loss) using the exact
advection terms that would be included in a regular SUPERTHERMAL
simulation.
The critical effective field \(E_c^\mathrm{eff}\) is then defined as the minimum value of the electric field for which there exists a real solution to \(U(p) = 0\). That is,
The critical field calculation has two main steps: initialisation of splines (typically once per simulation), and solution of the optimization problem for \(E_c^\mathrm{eff}\) (at each time step or iteration of the solver).
Since flux surface averages are significantly more expensive to evaluate than splines, and it is typically sufficient for us to resolve the critical field with a relative tolerance of approximately \(10^{3}\), it has proven effective to create 1D splines (one at each radius) over quantities involving fluxsurface or bounce averages (where relatively sparse splines are sufficient, sampling only tens of points).
First, we spline the exponent of the pitch distribution \(\xi_0 / \langle \xi \rangle\) on a uniform \(\xi_0\) grid. This allows the \(g\) function to be evaluated efficiently using routines for exact integration of a spline.
Secondly, we identify that all advection terms of interest can be factorised on the form \(A_p = a_p(p) \hat{A}_p(\xi_0,\,\theta)\) where the prefactor depends only on momentum and the remainder only on pitch and poloidal angle. Then, the bounceaverage of the pitchdependent part of each advection term \(\{\hat{A}_p\}\) is splined on a uniform pitch grid in the interval \(\xi_0 \in [0,1]\) since all advection terms of interest are either symmetric or antisymmetric.
Finally, these steps allow rapid calculation of a spline of the distributionbounce average of \(\hat{A}_p\) for each advection term, sampled in the mapped variable \(X = A^2/(1+A)^2 \in [0,\,1]\) (in which the functions are smoothly varying all the way up to the limit of \(A=\infty\) corresponding to all runaways having \(\xi=1\)), where \(A\) is the pitch distribution width parameter appearing in the exponent of \(f_\mathrm{RE}\).
This procedure allows the geometric part of the acceleration function \(U(p)\) to be evaluated essentially for free, with only the momentum (and unknownquantity) dependent prefactors requiring further computation, most of which is spent on the slowingdown frequency \(\nu_s\).
We solve for the critical field as a nested optimization problem, where the outer layer is the onedimensional root finding problem
where \(U_e\) is the extremum of \(U(p)\) with respect to \(p\) at a given \(\left\langle \boldsymbol{E}\cdot\boldsymbol{B} \right\rangle\). This is solved with an unbounded secant method (Newton’s method with the numerical derivative taken between successive iterations), where the initial guess in each solve is that \(E_c^\mathrm{eff}/E_c^\mathrm{tot}\) is constant in time.
The inner layer is the minimization problem
where the electric field is provided from the outer layer. This is solved using the bounded Brent’s method
implemented in the gsl_root
library, requiring an interval in \(p\) to be given which contains the minimum.
This is done by choosing a narrow interval \(p_\mathrm{opt} \times (1 \pm 0.02)\), where \(p_\mathrm{opt}\)
is the optimum from the previous solve. If the interval does not contain the minimum, it is expanded until a minimum
is enclosed – this is typically needed less than once in a thousand solves.
The model used for \(E_c^\mathrm{eff}\) is controlled with the following settings:
Option 
Description 


Use the approximation \(E_c^\mathrm{eff} = E_c^\mathrm{tot}\) 

Uses the analytical model Eq (2324) in Hesslow PPCF (2018) 



Full model outlined above 
Note
In typical scenarios with toroidal geometry, we observe discrepencies
of up to 2% between modes FULL
and SIMPLE
, less than 10%
with CYLINDRICAL
and over 50% with EC_TOT
.
Note
SIMPLE
takes essentially the same computation time as FULL
,
and is therefore not recommended except for benchmarking. Compared
with the simpler models EC_TOT
and CYLINDRICAL
, FULL
can
make especially fluid simulations substantially slower.
Runaway generation is generally treated correctly regardless of the sign of the sign of the electric field in DREAM. One exception occurs for the avalanche source when the electric field changes sign during the simulation. In this case, in default mode, the sign change may cause excessive avalanching in the new direction of the electric field. To resolve this issue, the density of runaways travelling antiparallel to the magnetic field (called n_re_neg in the code) can be introduced in the system. When this is done, the kinetic avalanche source is automatically modified to account for whether runaways are travelling parallel or antiparallel to the magnetic field.
Note
Tracking the direction of motion of the runaways requires the runaway grid to enabled.
The density of runaway electrons travelling antiparallel to the magnetic field is defined as
The avalanche source term in the kinetic equation is then modified to become
where the superscript sign indicates whether the source creates electrons in the parallel (+) or antiparallel () direction (corresponding to positive and negative values of \(\xi_0\) respectively). In the default source, the avalanche electrons are created in the direction instantaneously parallel to the electric field.
DREAM allows the user to prescribe transport in a few different ways to the runaway density.
Arbitrary advection and diffusion coefficients \(A_{r}=A_{r}(t,r)\) and \(D_{rr}=D_{rr}(t,r)\) can be provided to DREAM through the regular radial transport settings interface. Coefficients must obey the same rules as all other prescribed spatiotemporal quantities in DREAM, i.e. they can either be given as scalars (assumed constant in time/uniform in radius), as functions of only one of the time and radius coordinates, or as functions of both coordinates.
import numpy as np
import DREAM.Settings.Transport as Transport
ds = DREAMSettings()
...
D0 = 10
tDrr = np.linspace(0, 1)
rDrr = np.linspace(0, a)
# Generate matrix representing diffusion coefficient Drr(t,r)
T, R = np.meshgrid(tDrr, rDrr)
Drr = D0 * np.sin(np.pi*T) * (1R)
# Prescribe diffusion coefficient
ds.eqsys.n_re.transport.prescribeDiffusion(drr=Drr, t=tDrr, r=rDrr)
# Specify boundary condition (assume n_re=0 outside plasma)
ds.eqsys.n_re.transport.setBoundaryCondition(Transport.BC_F_0)
Note
A list of available boundary conditions can be found on the page DistributionFunction.
A model accounting for the momentumdependence of the runaway transport coefficients during avalanchedominated runaway generation was developed by Svensson et al. (2021). The model assumes that the distribution function takes the form
with \(\hat{A}(p)=2E/(p\nu_D\tau)\) and \(F(r,p,t)\) independent of the pitch coordinate \(\xi\). By integrating the kinetic equation over \(\xi\), an equation for the energy distribution \(F\) is obtained, taking the form
where \(U(p)\) is the pitchaveraged force felt by the electrons, \(\tau=m_ec/(eE_c)\) is the relativistic electron collision time and the pitch average is given by
The runaway density is then recovered by integrating over the momentum distribution from the userspecified parameter \(p_\star\) to \(\infty\)
Integrating the kinetic equation above in the same way yields a rate equation for the runaway density:
where the radial flux of runaways \(\Gamma_0\) is given by
Here, \(F_0\) is the zerothorder solution to \(F\), given by neglecting the radial transport, and is
The Svensson transport requires either advection or diffusion coefficients to be specified (or both) as functions of radius, momentum, pitch and (optionally) time. In addition to this, the user must also specify the value to use for the lower momentum cutoff \(p_\star\), typically taking to be near the critical momentum for runaway acceleration \(p_c\sim 1/\sqrt{E/E_c1}\).
Warning
Simulations may be sensitive to the choice of \(p_\star\) and so the user should take care and choose \(p_\star\) wisely.
Coefficients should be provided as fourdimensional arrays with shape
(nt, nr, nxi, np)
.
The transport coefficients will be interpolated in time during the simulation,
and two different methods for doing this exist: nearest interpolation (taking
the value of the coefficient closest in time to the current time) and linear
interpolation. The interpolation method is set with an argument interp1d
to
the methods DREAM.Settings.Transport.Transport.setSvenssonAdvection()
and DREAM.Settings.Transport.Transport.setSvenssonDiffusion()
.
The argument should be either Transport.INTERP1D_NEAREST
or
Transport.INTERP1D_LINEAR
.
It is also possible to interpolate in the transport coefficients using the total plasma current as the interpolating variable rather than time. This can be useful in situations where the transport coefficients have been generated by a separate code (such as a 3D MHD code) and are expected to depend rather on the plasma current than the actual time.
The following example illustrates how to use the Svensson transport module:
import DREAM.Settings.Transport as Transport
ds = DREAMSettings()
...
# Load transport coefficients from 3D MHD simulation output
# (note: 'loadTransportCoefficients()' is NOT provided by the
# DREAM Python interface)
t, r, p, xi, A, D = loadTransportCoefficients('3D_MHD_output.h5')
# Prescribe transport coefficients
ds.eqsys.n_re.transport.setSvenssonAdvection(ar=A, t=t, r=r, p=p, xi=xi, interp1d=Transport.INTERP1D_LINEAR)
ds.eqsys.n_re.transport.setSvenssonDiffusion(drr=D, t=t, r=r, p=p, xi=xi, interp1d=Transport.INTERP1D_LINEAR)
# Specify lower momentum boundary (in units of mc)
ds.eqsys.n_re.transport.setSvenssonPstar(pstar=0.3)
Optionally, the transport coefficients can be evolved as functions of total plasma current, rather than time, in the following way:
import DREAM.Settings.Transport as Transport
ds = DREAMSettings()
...
# Load transport coefficients from 3D MHD simulation output
# (note: 'loadTransportCoefficients()' is NOT provided by the
# DREAM Python interface)
t, r, p, xi, Ip, A, D = loadTransportCoefficients('3D_MHD_output.h5')
# Consider A & D as functions of Ip rather than t
# (must be called before 'setSvenssonAdvection()' and 'setSvenssonDiffusion()')
ds.eqsys.n_re.transport.setSvenssonInterp1dParam(Transport.SVENSSON_INTERP1D_PARAM_IP)
# Prescribe transport coefficients
ds.eqsys.n_re.transport.setSvenssonAdvection(ar=A, Ip=Ip, r=r, p=p, xi=xi, interp1d=Transport.INTERP1D_LINEAR)
ds.eqsys.n_re.transport.setSvenssonDiffusion(drr=D, Ip=Ip, r=r, p=p, xi=xi, interp1d=Transport.INTERP1D_LINEAR)
# Specify lower momentum boundary (in units of mc)
ds.eqsys.n_re.transport.setSvenssonPstar(pstar=0.3)
DREAM.Settings.Equations.RunawayElectrons.
RunawayElectrons
(settings, density=0, radius=0, avalanche=1, dreicer=1, compton=1, Eceff=4, pCutAvalanche=0, comptonPhotonFlux=0, tritium=False, hottail=1)¶Bases: DREAM.Settings.Equations.UnknownQuantity.UnknownQuantity
, DREAM.Settings.Equations.PrescribedInitialParameter.PrescribedInitialParameter
__init__
(settings, density=0, radius=0, avalanche=1, dreicer=1, compton=1, Eceff=4, pCutAvalanche=0, comptonPhotonFlux=0, tritium=False, hottail=1)¶Constructor.
fromdict
(data)¶Set all options from a dictionary.
setAdvectionInterpolationMethod
(ad_int=1, ad_jac=2, fluxlimiterdamping=1.0)¶Sets the interpolation method that is used in the advection terms of the transport equation.
ad_int (int) – Interpolation method to use for the radial coordinate.
ad_jac (int) – Jacobian interpolation mode to use for the radial coordinate.
fluxlimiterdamping (float) – Damping parameter used to underrelax the interpolation coefficients during nonlinear iterations (should be between 0 and 1).
setAvalanche
(avalanche, pCutAvalanche=0)¶Enables/disables avalanche generation.
setCompton
(compton, photonFlux=None)¶Specifies which model to use for calculating the compton runaway rate.
setDreicer
(dreicer)¶Specifies which model to use for calculating the Dreicer runaway rate.
setEceff
(Eceff)¶Specifies which model to use for calculating the effective critical field (used in the avalanche formula).
setHottail
(hottail)¶Specify which model to use for hottail runaway generation
setInitialProfile
(density, radius=0)¶setNegativeRunaways
(negative_re=True)¶Introduce a density of runaway electrons with negative pitch, allowing the kinetic avalanche source term to properly account for largeangle collisions with runaways moving in different directions.
setTritium
(tritium)¶Specifices whether or not to include runaway generation through tritium decay as a source term.
todict
()¶Returns a Python dictionary containing all settings of this RunawayElectrons object.
verifySettings
()¶Verify that the settings of this unknown are correctly set.
verifySettingsPrescribedInitialData
()¶