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) hot-tail 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 advective-diffusive 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 hot-tail 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 non-zero
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 non-zero 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 fluid-kinetic 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 hot-tail 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 so-called 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 helium-3 via the beta decay process
The corresponding fluid runaway rate is given by
and the kinetic source rate by
where \(n_{\rm T}\) is the tritium density, \(\tau_{\rm T} = 4800\pm 8\) days is the tritium half-life, \(f_\beta(\gamma)\) is the beta energy spectrum 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.
The following settings are used to control the tritium source mode:
Option |
Description |
|---|---|
|
Do not include generation from tritium beta decay in the simulation. |
|
Use the fluid runaway rate, model described in Martin-Solis NF (2017). |
|
Use the kinetic source rate, model described in Ekmark JPP (2024). |
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(tritium=TRITIUM_MODE_FLUID)
# 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 Klein-Nishina differential cross section \(\frac{{\rm d}\sigma}{{\rm d}\Omega}(p,E_\gamma)\) depends on the momentum and the incident-photon energy \(E_\gamma\). Using the fluid Compton generation, this source is integrated over the runaway region \(p>p_c\), analogous to the tritium and avalanche calculations. The integrated cross section is given by Equation (29) in Martin-Solis NF (2017) , and this total cross section \(\sigma(p_c,\,E_\gamma)\) depends on plasma parameters through \(p_c\). For both the fluid and kinetic Compton generation models, the net generation rates are consequently obtained as the integrals of the Klein-Nishina or 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).
Accordingly, the fluid runaway rate is given by
and the kinetic source rate by
Following Martin-Solis Eq (24), we model the photon energy spectrum with
where \(\Phi_0 = \int \Phi \,\mathrm{d}E_\gamma\) is the total photon gamma flux. For ITER, according to Martin-Solis, this value is \(\Phi_0 \approx 10^{18}\,\mathrm{m}^{-2}\mathrm{s}^{-1}\) during the nuclear phase. We use \(z = [C_1 + \ln(E_\gamma[MeV])]/C_2 + C_3(E_\gamma[MeV])^2\), where \(C_1,\ C_2\) and \(C_3\) are free positive parameters used to determine the shape of the photon flux spectrum. The Compton fitting tool can be used to fit these parameters to data of this spectrum, otherwise the default values are the ones used by Martin-Solis, i.e. \(C_1 = 1.2\), \(C_2 = 0.8\) and \(C_3 = 0\).
Note that with the kinetic source rate, the electrons will be added to the kinetic grids, when possible. If the runaway grid is not activated, the generation of runaways above the upper limit of the hot-tail grid will be integrated, as with the fluid runaway rate, from the upper limit of the hot-tail grid.
The following settings are used to control the Compton source mode:
Option |
Description |
|---|---|
|
Do not include Compton scattering in the simulation. |
|
Use the fluid runaway rate, model described in Martin-Solis NF (2017), with tuneable total photon flux. |
|
Use the kinetic source rate, model described in Ekmark JPP (2024), with tuneable total photon flux. |
|
Short-hand for |
|
Short-hand for |
The fluid Compton source can be activated following this example:
import DREAM.Settings.Equations.RunawayElectons as Runaways
ds = DREAMSettings()
Phi0 = 1e18 # total photon flux in units of m^-2 s^-1, typical of ITER
c_1 = 1.2
c_2 = 0.8
c_3 = 0.
ds.eqsys.n_re.setCompton(compton=Runaways.COMPTON_MODE_FLUID,
photonFlux=Phi0, C1=c_1, C2=c_2, C3=c_3)
If one is simulating ITER specifically, it is possible to instead just do
import DREAM.Settings.Equations.RunawayElectrons as Runaways
ds = DREAMSettings()
# Equivalent to the above example
ds.eqsys.n_re.setCompton(compton=Runaways.COMPTON_MODE_ITER_DMS_FLUID)
The photon flux can also be specified as a function of time. In this case, the photon flux should be given as a 1D array along with a equally shaped time vector:
import DREAM.Settings.Equations.RunawayElectrons as Runaways
ds = DREAMSettings()
Phi0 = 1e18
Phi = np.ones((4,)) * Phi0
t = np.zeros(Phi.shape)
t[1] = 5e-3
t[2] = 6e-3
t[3] = 1e3
Phi[2:] /= 1e3
# Set photon flux density (in units of m^-2 s^-1)
ds.eqsys.n_re.setCompton(compton=Runaways.COMPTON_MODE_FLUID, photonFlux=Phi, photonFlux_t=t)
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 angle-averaged electron distribution evaluated at the critical momentum. In this hottail model, the critical momentum is found from the angle-averaged kinetic equation in the Lorentz limit of strong pitch-angle scattering:
which for a given distribution \(f_0\) is solved numerically for the critical momentum \(p=p_c\). Here we assume that only pitch-angle scattering, electric field acceleration and collisional slowing down contributes to the electron dynamics. We have taken the non-relativistic 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 high-Z 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 low-Z limit, Section 4.1.3 in Svenningsson’s thesis) |
|
The high-Z 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 pitch-angle 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 distribution-bounce 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 mean-force 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 flux-surface 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 bounce-average of the pitch-dependent 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 anti-symmetric.
Finally, these steps allow rapid calculation of a spline of the distribution-bounce 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 unknown-quantity) dependent prefactors requiring further computation, most of which is spent on the slowing-down frequency \(\nu_s\).
We solve for the critical field as a nested optimization problem, where the outer layer is the one-dimensional 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 (23-24) 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 anti-parallel 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 anti-parallel 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 anti-parallel 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 anti-parallel (-) 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) * (1-R)
# 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 momentum-dependence of the runaway transport coefficients during avalanche-dominated 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 pitch-averaged 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 user-specified 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 zeroth-order 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 cut-off \(p_\star\), typically taking to be near the critical momentum for runaway acceleration \(p_c\sim 1/\sqrt{E/E_c-1}\).
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 four-dimensional 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)
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.
The radial diffusion coefficient used in the frozen current mode is generally written on the form
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)\) |
|---|---|
|
\(0\) |
|
\(1\) |
|
\(\beta_\parallel = v_\parallel/c\) |
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.n_re.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.n_re.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.n_re.transport.setFrozenCurrentMode(
Transport.FROZEN_CURRENT_MODE_BETAPAR,
Ip_presc=Ip, D_I_max=200
)
DREAM.Settings.Equations.RunawayElectrons.RunawayElectrons(settings, density=0, radius=0, avalanche=1, dreicer=1, compton=1, Eceff=4, pCutAvalanche=0, comptonPhotonFlux=0, C1_Compton=1.2, C2_compton=0.8, C3_Compton=0.0, tritium=1, hottail=1, lcfs_loss=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, C1_Compton=1.2, C2_compton=0.8, C3_Compton=0.0, tritium=1, hottail=1, lcfs_loss=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 under-relax the interpolation coefficients during non-linear iterations (should be between 0 and 1).
setAvalanche(avalanche, pCutAvalanche=0)¶Enables/disables avalanche generation.
setCompton(compton, photonFlux=None, photonFlux_t=array([0]), machine=None, C1=None, C2=None, C3=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).
setExtrapolateDreicer(extrapolateDreicer=False)¶Extrapolates the result from the neural network for small electric fields such that the Dreicer generation rate is continuous and has continuous derivative.
setHottail(hottail)¶Specify which model to use for hottail runaway generation
setInitialProfile(density, radius=0)¶setLCFSLoss(lcfs_loss)¶Specifies which model to use for calculating the LCFS loss term.
setLCFSLossPsiEdget0(psi_edge_t0, user_input_active=True)¶Sets the value of psi_p at the plasma edge at t = 0, used to determine the LCFS radial point. Use in restarts to keep the value from the first simulation. user_input_active should be 1 (default) or set to 0 to manually switch off the user input.
setLCFSLossTime(t_loss, radius=0)¶Sets the timescale constant t_loss for the LCFS loss term.
setNegativeRunaways(negative_re=True)¶Introduce a density of runaway electrons with negative pitch, allowing the kinetic avalanche source term to properly account for large-angle 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()¶