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 FokkerPlanck equation
where the second term on the LHS represents electric field acceleration along magnetic field lines, \(C\{f\}\) denotes the FokkerPlanck collision operator, the second term on the RHS represents advectivediffusive 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 quasicharge 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 angleaveraged 2D FokkerPlanck 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
pitchangle 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.
Page overview
In the finite volume method, which is used for discretizing most equations in DREAM (particularly the FokkerPlanck equation), advection terms are discretized as
where halfindices 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
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 socalled
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 positivitypreserving interpolation schemes, which avoid oscillations
and ensure that the resulting distribution function is everywhere positive. The
simplest schemes are the linear first and secondorder upwind schemes, which
are computationally cheap but less accurate. For most cases, the nonlinear
flux limiter schemes are to be preferred.
Note
All flux limiter schemes require the nonlinear solver to be used.
The following interpolation schemes are provided by DREAM:
Name 
Nonlinear? 
Reference 


No 


No 


No 


No 


No 


Yes 


Yes 


Yes 
Waterson and Deconinck (1995) 

Yes 
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 nonlinear iteration process to underrelax the algorithm as this may accelerate convergence. In each iteration the new interpolation coefficient is set to
where \(\delta_{\rm prev}\) is the old value of the interpolation coefficient and \(\delta_{\rm calc}\) is the calculated new value.
The flux limiter interpolation schemes are nonlinear 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 


Default. Neglects any contribution from the unknown to the jacobian 

Fully accounts for the unknown contribution to the jacobian 

Uses the jacobian corresponding to 
Note
FULL
is only active for the smooth TCDF
and OSPRE
flux limiter
schemes, for all others it defaults to LINEAR
.
UPWIND
requires a nonlinear 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 nonzero elements in the matrix.
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) pitchangle 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 nonlinear 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.
Numerical solution of the FokkerPlanck 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
setting \(f(p>p_{\rm max}) = 0\) (default)
ensuring that the flux out across \(p=p_{\rm max}\) is the same as the flux into the last cell on the momentum grid.
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.
Name 
Description 


Assume \(f(p>p_{\rm max}) = 0\) (default) 

Assume \(\Phi^{(p)}_{N_p+1/2} = \Phi^{(p)}_{N_p1/2}\) 

Assume \(\Phi^{(p)}_{N_p+1/2} = \left( 1  \delta \right)\Phi^{(p)}_{N_p1/2} + \delta\Phi^{(p)}_{N_p3/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
.
An initial condition is required when solving the FokkerPlanck 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.
DREAM can initialize the distribution function to a MaxwellJüttner distribution function:
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 secondorder 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.
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)
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()
.
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)
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 pitchangle scattering.
In DREAM, we can account for the increased pitchangle scattering using the diffusion operator derived by Kurzan et al (1995), which for DREAM takes the form
with
and resonant momentum and resonance width
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
where the characteristic width is evaluated at \(p_\parallel=p_{nm}\),
The model used for the resonance can be controlled using the following options:
Name 
Model for the resonance width 


The ripple resonance is ignored altogether. 

Constant in \(\leftp_\parallel  p_{nm} \right < \Delta p_{nm}\). 

Gaussian envelope function 
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 largeaspect 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
The resonant momentum then takes the form
with \(G(r)/R_0\) forming the finite toroidal magnetic field stregth.
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 


Poloidal mode numbers corresponding to each resonance. 

Toroidal mode numbers corresponding to each resonance. 

Resonant momentum for 

Resonant momentum for 
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 = [1e4,5e5,2e5]
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)
When running DREAM with a selfconsistent ion charge state evolution, the Newtonsolver jacobian of the collision operator with respect to these densities will produce a large number of nonzero 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 distributiongrid 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
.
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 nonzero elements in the jacobian by more than a factor of two, and produce performance gains of ordersofmagnitude in certain situations. However, as with all approximations to the jacobian, there is a risk of deteriorated stability.
A radial transport term can be added to the FokkerPlanck equation. The general transport term takes the form
where \(A\) and \(D\) are the advection and diffusion coefficients respectively. DREAM provides two different ways of prescribing these coefficients:
Let the user prescribe \(A\) and/or \(D\) directly.
Use the model derived by Rechester and Rosenbluth and let the user prescribe the magnetic perturbation level \(\delta B/B\).
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.
The following example shows how to prescribe advectivediffusive 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 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 RechesterRosenbluth operator takes the form
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\).
The following example illustrates how to prescribe a RechesterRosenbluth 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=1e3)
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.
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:
Assume that no particles can cross the plasma boundary (default).
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 quasineutrality, DREAM compensates for the lost electrons by adding additional cold electrons.
The following boundary conditions can be applied:
Name 
Description 


Assume that no particles can leave the plasma through \(r=a\). (default) 

Assume that \(f(r>a) = 0\). 

Assume that \(\partial f/\partial r = \mathrm{const.}\) on plasma edge. 
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)
When applying transport to the kinetic equation, two diagnostic quantities are saved as an other quantity for each distribution function:
Name of OtherQuantity 
Description 


Density moment: Rate at which runaways leave the plasma due to transport 

Energy moment: Rate at which runaway energy leaves the plasma due to transport 

Density moment: Rate at which hot electrons leave the plasma due to transport 

Energy moment: Rate at which hot electrons leave the plasma due to transport 
Secondary runaway electrons (produced through the socalled avalanche mechanism)
can be included in the simulation of the distribution function through a
RosenbluthPutvinski 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 knockons,
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.
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 RosenbluthPutvinski 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 knockons
created is independent of grid parameters.
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
Synchrotron losses are disabled by default.
Name 
Description 


Include synchrotron losses in simulation 

Do not include synchrotron losses in simulation 
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)
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 


The distribution is resolved on a kinetic (finitevolume) grid 

The distribution is modelled with an analytical function 
Typically, in the analytical mode, f_hot
is described by a hottail
slowingdown distribution and f_re
by a quasisteady state avalanche
distribution.
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 nonlinear solver fails to converge, increasing resolution sometimes helps, in particular decreasing the time step. If a floatingpoint 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.
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.
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.
data (dict) – Dictionary to load distribution function from.
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.
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 underrelax the interpolation coefficients during nonlinear 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.
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.
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.
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.
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).
mode (int) – Flag indicating whether or not to enable synchrotron losses (may be bool).
todict
()¶Returns a Python dictionary containing all settings of this DistributionFunction object.
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.
The HotElectronDistribution
class contains settings specific to the hot electron
distribution, f_hot
. It contains the following methods:
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 


\(\int_\mathrm{pCutoff}^\mathrm{pMax} \mathrm{d}p\) 

\(\int_{\mathrm{pCutoff}\sqrt{2T_\mathrm{cold}/m_e c^2}}^\mathrm{pMax} \mathrm{d}p\) 

\(\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:
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 nonlinear 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 nonuniform 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.
The following example shows how to change hotelectron 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 selfconsistent simulation with a rapid temperature drop when the thermal Maxwellian is resolved ongrid.
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 MaxwellJüttner distribution evaluated
at the instantaneous coldelectron 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 


\(\mathrm{S\_particle} = 0\) (the particle source is neglected) 

\(\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) 

\(\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.
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 illconditioned when running a selfconsistent simulation
with COLLISION_FREQUENCY_MODE_FULL
, try changing particle source mode, which may make the equation
system more wellconditioned.
Different models may be implemented to describe the hotelectron distribution. These are controlled via the following options:
Name 
Model for 


The isotropic nonrelativistic slowingdown 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 normalizedtime 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.
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.
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.
inittype (int) – Flag indicating how to initialize f_re.
todict
()¶Returns a Python dictionary containing all settings of this RunawayElectronDistribution object.