Newsletter #6, January 2025

Welcome to the first newsletter of 2025! In the few months that have passed since the last newsletter, several things have happened to the code, as you might be able to tell from the long table of contents. Perhaps most importantly though, the DREAM Developer Council has decided to start making regular releases of the code alongside every newsletter. As usual, the master branch will be updated as soon as new changes have been approved in a pull request, but with every newsletter we will now also tag the latest version of the master branch with a version number. This might make it easier to communicate the DREAM version used in publications.

I. Announcements

i. DREAM release versions

DREAM has so far not utilized the concept of release versions, and we have recommended users to always work on the latest version of the master branch. In recent discussions in the developer council, it has however been proposed that regular release versions of the code are published in order to more easily refer to different code versions. Each new release will be published on the DREAM GitHub page alongside each newsletter, and the version will be referred to by year number and release number. The first official release version of DREAM will therefore have version number 25.1, corresponding to this newsletter, as it is the first release of the code in 2025. Subsequent issues of the newsletter will hence be accompanied by code releases with unit increments in the minor version number.

II. New features

i. Option for Walkowiak collision fequencies

by Jedrzej Walkowiak

The Walkowiak model for the collision frequencies has been available in DREAM for some time, however a named option for using these frequencies has been missing in the Python interface. Now, it is possible to use the Walkowiak frequencies by adding

import DREAM.Settings.Collisions as Collisions

ds.collisions.collfreq_type = Collisions.COLLFREQ_TYPE_WALKOWIAK

These collision frequencies can be used for $Z>8$, and should be used for $Z>18$ (i.e. elements heavier than Ar).

This option was added in pull request #357 and is now available on master.

ii. Convert between different types of input current densities

by Mathias Hoppe

The equations involving current densities in DREAM are written such that quantities of the form $j/B$ appear. Furthermore, under the approximations made in DREAM, the poloidal variation of $j(\theta)\propto B(\theta)$, and so the quantity $j/B$ can be considered a flux function.

When prescribing an initial current density profile, DREAM therefore assumes that the profile is given in the point where $B=B_{\rm min}$. This is rather specific to DREAM, and other codes tend to use other definitions for the current density. To facilitate compabitibility with other codes, additional options have been added for importing other types of current densities directly into the C++ code.

At present, the following alternative current density definitions are supported, and it is possible to specify current density profiles definied according to these directly to DREAM, without prior conversion:

Toroidal current density

$$\left\langle \boldsymbol{j}\cdot\nabla\phi \right\rangle = \frac{j_{\parallel,{\rm min}}}{B_{\rm min}}\left[ G\left\langle\frac{1}{R^2}\right\rangle \right]$$

CORSICA current density

$$\bar{J}_\phi = \frac{j_{\parallel,{\rm min}}}{B_{\rm min}} \frac{\left\langle B^2 \right\rangle}{G\left\langle 1/R^2 \right\rangle}$$

For a detailed description of how to use this functionality, please refer to the DREAM online documentation.

This functionality was added in pull request #355 and is now available on master.

iii. Show ‘energyloss_T_cold’ in power balance plot

by Mathias Hoppe

The power balance plot (T_cold.plotEnergyBalance()) shows most of the terms involved in the power balance equation, but has previously not included the power lost through the plasma boundary. With this change, also the quantity energyloss_T_cold is included in this plot when the spatially integrated power is visualized.

This functionality was added in pull request #360 and is now available on master.

iv. Store equilibrium quantities in output file

by Mathias Hoppe

DREAM has previously not stored quantities describing the magnetic equilibrium used in a simulation in the output file. Visualizations in a poloidal plane have therefore relied on equilibrium data stored in external files, which is not always guaranteed to be available when the output file is transferred to different locations.

With this modification, DREAM now saves flux surface contours and some associated quantities to all output files under the group do.grid.eq. Along with this, some SPI plotting and post-processing routines have also been updated to utilize the new equilibrium data.

This functionality was added in pull request #329 and is now available on master.

v. Halo region heat losses

by Lorenzo Votta

A new source term describing heat losses in the halo region during a VDE has been added. The source term was derived by Kiramov et al and is implemented in DREAM as

$$\frac{\partial W_{\rm cold}}{\partial t} = \frac{2}{3}\frac{q_\parallel}{L_\parallel}\Theta\left(r > r_{\rm LCFS}\right),$$

with

$$q_\parallel = 2\gamma_{\rm sh} n_e c_{\rm s} T_e.$$

This functionality was added in pull request #384 and is now available on master.

vi. Open DREAMOutput from pathlib.Path

by Mathias Hoppe

Until now, the name of a DREAM output file had to be specified as a string. With this change, it is now possible to use a pathlib.Path object instead.

This functionality was added in pull request #366 and is now available on master.

vii. Improved post-processing of SPI shards

by Mathias Hoppe

Due to the lack of high-quality equilibrium information in the output file, it has previously not been possible to accurately visualize the SPI shard positions on a poloidal cross-section, or calculate the shard arrival time. This addition uses the newly added equilibrium information (see iv above) to calculate and visualize the shard positions in the actual magnetic geometry. To visualize the pellet shards in a poloidal cross-section, use the x_p.plotAtTime() method:

do = DREAMOutput('output.h5')
# Plot the shards at the specified time step
# (all parameters are optional; default values indicated here)
do.eqsys.x_p.plotAtTime(t=-1, shards=None)

This method also allows one to only plot a subset of the shards, via the shards parameter.

A similar routine, x_p.plotTrajectoryPoloidal(), which plot the trajectory of a set of shards in the poloidal plane, has also been added:

x_p.plotTrajectoryPoloidal(shards=None)

It is now also possible to easily calculate the arrival time of the shards using the method x_p.arrivalTime(), which returns the time index at which the first shard arrives to the outermost flux surface. The method can be applied to a subset of the shards by specifying the shards parameter.

This functionality was added in pull request #366 and is now available on master.

viii. Custom random number generator for SPI shards

by Mathias Hoppe

The SPI shard positions and velocities are specified in the Python interface, and are typically set with the help of the numpy.random random number generator. In certain context, the numpy.random generator is not suitable, such as when running several DREAM simulations in parallel within the same process (which can be the case in certain Bayesian inference frameworks). In such cases, it is preferable to specify a custom random number generator.

With this new functionality, all SPI methods which draw random numbers now take the additional optional argment random which specifies the module to use for generating random numbers. The module is expected to be interface-compatible with the numpy random number interface.

This change was introduced in pull request #366 and is now available on master.

ix. Add routines for downloading TCV equilibrium

by Mathias Hoppe

A set of routines for download TCV magnetic equilibrium data in the format required by DREAM was added to the eqget tool. The routines rely on the tcvpyget library, developed by Mathias.

This functionality was implemented in pull request #369 and is now available on master.

x. Adaptive current gradient-driven transport

by Mathias Hoppe and Lorenzo Votta

A common problem in disruption simulations is the formation of ohmic current channels which excessively heat the plasma locally and lead to a type of runaway solutions. These were shown already by Putvinski to be valid solutions of the power balance equation, and so are to be expected in certain situations. In experiments, it is believed that such current channels would drive MHD instabilities which rapidly flatten them.

In an attempt to emulate the MHD activity driven by thin current channels, and access scenarios otherwise inaccessible by codes such as DREAM, an adaptive current gradient-driven transport model was implemented. The model monitors the current gradient everywhere in the plasma and triggers transport in $n_{\rm re}$, $\psi_{\rm p}$, and $T_{\rm cold}$ if it exceeds a user-prescribed value. The transport remains active until the gradient drops below a specified value (90% of the initial value by default).

The transport can be activated for the three physical quantities separately, or all at once with consistent parameters. Details about this can be found in the online documentation.

This functionality was implemented in pull request #364 and is now available on master.

xi. Ionization time as other quantity

by Mathias Hoppe

The adaptive ionization time stepper has become a popular time stepping option in DREAM. It relies on the estimated ionization time scale

$$\tau_{\rm ioniz} = \frac{n_{\rm cold}}{|\partial n_{\rm cold}/\partial t|}.$$

With this addition, the ionization time scale is now evaluated at all radii and stored as an other quantity (if requested). It is available under fluid/tIoniz. The shortest estimated ionization time scale, which is the value actually used by the ionization-based time stepper, is also calculated and stored as scalar/tIoniz.

This functionality was implemented in pull request #378 and is now available on master.

xii. Code information in output

by Mathias Hoppe

Until now, DREAM output files have not contained any metadata which clearly identifies the version of the code used to generate it. With this update, DREAM will now store code information under a separate group /code in the output file. Specifically, the following information will be saved:

• commit: Commit hash for the version of DREAM used to generate the output file.
• datetime_commit: Date and time of DREAM commit.
• datetime_simulation: Date and time at which the simulation finished.
• refspec: DREAM refspec for the git repository (this includes e.g. current branch name).

This functionality was implemented in pull request #391 and is now available on master.

III. Interface changes

i. Add dedicated switch for custom radial grid

by Mathias Hoppe

A previous interface change (issue #5, II.i) broke the custom radial grid, which would be activated if nr = 0. In this fix, the custom radial grid is granted a dedicated flag which informs the code whether or not it should be used. It therefore no longer depends on the nr parameter.

This change was introduced in pull request #393 and is now available on master.

ii. Add parameter ‘Ninj’ to output

by Mathias Hoppe

When calculating SPI assimilation fractions, it is necessary to know the number of injected particles. Prior to this change, it was necessary to either hard-code this value in your post-processing scripts, or manually invert the shard radius. With this change, the number of injected particles is stored as a dummy parameter in the settings file and will be propagated through to the settings section of the output file (if saved). This enables easy access to the number of injected particles via do.settings.eqsys.spi.Ninj.

This change was introduced in pull request #366 and is now available on master.

iii. Specify electric field when initializing f_re

by Mathias Hoppe

The initial value for the runaway distribution function, $f_{\rm re}$, can be taken as an analytical avalanche distribution. So far, the parameters of this distribution are taken from the initial values of electric field, temperature and ion densities. In some cases, such as when $E=0$ (which can happen erroneously when initializing it from $j_{\rm tot}$, as reported in issue #367), the avalanche distribution would be flat in $\xi_0$, which is undesirable. It was therefore proposed that the electric field should be decoupled from the other parameters and that the user should specify the electric field value for which to initialize the analytical avalanche distribution. This is also motivated by the fact that in scenarios where plasma parameters vary in time, the avalanche distribution function does not necessarily correspond to one evaluated for the instantaneous value of $E$.

The analytical avalanche distribution function can now be initialized via

r_init = np.linspace(0, a)
E_init = 0.05 * (1-0.5*r_init/a)   # V/m
ds.eqsys.f_re.setInitialAvalancheDistribution(E=E_init, r=r_init)

Note that the radial grid option is optional. If a uniform electric field profile is desired, the electric field value can be given as a scalar.

This change was introduced in pull request #368 and is now available on master.

IV. Bug fixes

i. Incorrect moments in KineticQuantity

by Ida Ekmark

The calculation of the plasma current from a distribution function was missing a geometrical factor, and general moments of kinetic quantities contained an indexing bug which, both of which have now been resolved. If you have relied on the Python interface for calculating the total plasma current or other moments of the distribution function (in a simulation with more than radial grid point) your results may contain some errors.

This bug was fixed in pull request #371 and is now available on master.

ii. Use actual pressure derivative from EQDSK file

by Ida Ekmark

The script EQDSK.py was previously calculating the derivative $\mathrm{d}p/\mathrm{d}\psi$ of the pressure from the spline fitted to the input pressure profile. This was however found to cause problems when certain COCOS numbers were used, and so since data for this derivative is typically found in the EQDSK file, the Python script now loads this data directly from the file.

This bug was fixed in pull request #389 and is now available on master.

iii. Divide toroidal flux by 2π

by Lorenzo Votta

A factor of $1/2\pi$ was missing in the calculation of the toroidal magnetic flux in the C++ kernel of DREAM. This quantity mostly does not affect calculations, aside from calculations involving hyper-resistive diffusion, for which the hyper-resistive diffusion coefficient $\Lambda$ should now be multiplied by a factor $2\pi$.

This bug was fixed in pull request #383 and is now available on master.

iv. Bug in jacobian of IonFluidRunawayIonizationTerm

by Mathias Hoppe

An indexing bug was found in the jacobian for the new IonFluidRunawayIonizationTerm which caused the Newton solver to converge more slowly when the term was included.

This bug was fixed in pull request #373 and is now available on master.

v. Bug in radial grid of ion transport

by Mathias Hoppe

When applying ion transport without specifying the radial grid, the Python interface tries to generate an automatic radial grid on the interval [0, a]. If a (the plasma minor radius) has not been defined at the point when the call is made, this will result in an array of zeros.

Since the calcTransportCoefficientExpdecaySingleChargeState() method does not prescribe any radial variation of the transport coefficient, there is no point in having a radial grid with more than one grid point. Line 1166-&ellip; therefore been adjusted to just set r = [0] if no radial grid has been specified by the user.

This bug was fixed in pull request #375 and is now available on master.

vi. Make it possible to run SPI simulations through dreampyface

by Mathias Hoppe

A number of issues were found in dreampyface and the SPI module, which caused problems when running SPI simulations. These problems where such that the simulations would crash because of invalid input, and so previous simulations which finished succesfully should not have been affected by this bug.

This bug was fixed in pull request #366 and is now available on master.

vii. Generalize SPI drift radius

by Mathias Hoppe

The original implementation of the SPI drift effect assumed that the pellet would travel in, or close to, the midplane. In situtations where this is not the case, such as when the pellet is injected off-midplane, the calculated drift could be incorrect and even take unphysical values. With this change, the SPI drift radius is not calculated more directly by determining radial flux surface label for the flux surface to which the plasmoid drifts.

This bug was fixed in pull request #366 and is now available on master.

viii. Avoid dummy SPI settings

by Lorenzo Votta

In cases where todict() was called on a DREAMSettings object without SPI, the nbrShiftGridCell and TDrift parameters would be assigned dummy values. This meant that if SPI settings were later set on this object, the dummy value would apply to the very first shard and cause weird behaviour in the simulation. This fix resolves this issue so that the nbrShiftGridCell and TDrift parameters are consistent with the number of calls which set SPI parameters.

This bug was fixed in pull request #398 and is now available on master.

ix. Tritium generation not enabled

by Lorenzo Votta

Since 78fb299 (2023-04-25) the call ds.eqsys.n_re.setTritium(True) has not actually enabled tritium generation, but rather set the tritium option to TRITIUM_MODE_NEGLECT. This bug fix resolves the problem and allows one to call ds.eqsys.n_re.setTritium(True) to enable tritium generation.

(It should be noted that this call was not possible prior to the above commit, although the code seems to indicate that it should have been made possible)

This bug was fixed in pull request #396 and is now available on master.

x. NumPy array handling in dreampyface

by Mathias Hoppe

A problem with how NumPy arrays are handled in dreampyface (DREAM compiled as a Python library) was discovered. Arrays using non-default strides would be incorrectly loaded. Such arises can arise from certain NumPy operations which transform the shape and structure of the array without actually modifying the array in memory (e.g. transpose).

With this fix, dreampyface should take the array stride (number of bytes between consecutive elements of the array) into account.

This bug was fixed in pull request #395 and is now available on master.

V. Recent publications

• C. Wang et al, Nucl. Fusion 65, 016012 (2025)
• O. Vallhagen et al, submitted to J. Plasma Phys. (2024)
• M. Hoppe et al, submitted to Plasma Phys. Control. Fusion (2024)
• P. Halldestam et al, submitted to J. Plasma Phys. (2025)