Newsletter #3, January 2024
Happy New Year and welcome to the third edition of the DREAM newsletter. Since the last edition, a number of new features and fixes have made their way into the code and are now ready to be used. The year 2024 is set to be another successful year for the DREAM community, with a growing user base, a number of publications either in preparation or already submitted, and a number of new projects starting up!
Among the features listed in this edition of the newsletter you will find the possibility of compiling DREAM as a Python library. This feature has been long in the making, but finally it was decided to merge it into master. Although the functionality is still somewhat experimental, it can enable much more fine-grained control over our simulations than ever before and could turn out to be a powerful tool in optimizing DREAM simulations.
DREAM also recently hit somewhat of a milestone: the master branch now contains more than 50,000 lines of code (counting both C++ and Python)! At the time of writing, I count 50,407 lines of useful code, and I am looking forward to seeing how much it will grow in the next year :)
// Mathias
I. New features
i. Flux from f_hot to n_re as other quantity
by Mathias Hoppe
Since early versions of DREAM, it has been possible to store certain invidual
components of the runaway rate to the output file, as for example
other/fluid/gammaDreicer, other/fluid/gammaHottail and
other/fluid/GammaAva. These are however specialized to particular fluid
generation rates, and do not account for the runaways originating from the hot
electron grid. Instead, we have historically had to deduce the generation rate
from the hot electron grid from the all-encompassing other quantity
other/fluid/runawayRate, which is simply the numerical time derivative of
n_re. This becomes problematic when radial transport is included in the
simulation, since we also do not keep track of the transport rate of runaways
(except for the rate of lost runaways, scalar/radialloss_n_re).
A new other quantity has now been added, named other/fluid/gammaFhot, which
is the density of runaways transferred from f_hot to n_re. In
simulations involving the hot electron grid, this quantity should replace
gammaDreicer and gammaHottail in plots of the runaway rate
constituents.
ii. Kinetic tritium and Compton runaway source terms
by Ida Ekmark The runaway source terms modelling tritium decay and inverse Compton scattering runaway generation have previously been fully fluid—even in simulations with the kinetic grids enabled. Now, kinetic operators for these generation mechanisms have been developed.
The kinetic operators provide energy resolution for the generated runaways and
can be particularly useful in scenarios where the critical momentum changes
rapidly during the simulation. Although the operators are “kinetic”, they can
also be added to fluid equations (i.e. to the equation for n_re) and will
then be integrated numerically.
The kinetic tritium and Compton operators must be manually enabled via the following line in your settings file:
import DREAM.Settings.Equations.RunawayElectrons as REs
ds.eqsys.n_re.setTritium(REs.TRITIUM_MODE_KINETIC)
ds.eqsys.n_re.setCompton(REs.COMPTON_RATE_ITER_DMS_KINETIC)
iii. Residual convergence monitoring
by Mathias Hoppe and Ida Ekmark
During wide parameter scans, it was found that the DREAM Newton solver may in rare cases believe itself to have converged to a solution, while the “solution” actually does not satisfy the system of equations. The fact that this is at all possible comes from the convergence condition used, which requires that
for each physical unknown quantity individually. Normally, approaches zero as the candidate solution approaches the true solution, but in rare situations, in which the jacobian becomes near singular, may approach zero without the candidate solution being close to the true solution.
While we currently have no method for avoiding situations like these, we have implemented an additional routine in DREAM which monitors the residual vector . Since, by definition of a “solution” to our system of equations, the residual must satisfy , the new routine therefore emits a warning if is “far” from zero. Specifically, this is done by comparing each element of the residual to a pre-defined “typical scale length” of each equation and requiring that each element be strictly smaller than this scale length.
This check is done at the end of Newton iteration, when the correct solution
is believed to have been found. The result of the check is stored to the output
file under the solver category—a flag indicating whether the residual was
converged or not—along with the maximum residual value (normalized to the
scale length) obtained in that time step.
Note that while the monitor has been observed to guard against false positives,
the monitoring does also occasionally report false negatives. Before discarding
your simulations after DREAM telling you they are diverged, you should therefore
manually check the residual error (e.g. by calling
solver.plotResidualMaxError() in the CLI (py/cli/cli.py)) and judge
whether the large residual seems to be an actual problem or not. In particular,
if the residual is too large over many time steps, this is an indication that
your simulation is most likely diverged. Occasional large residuals can however
be tolerated and need not be incorrect. Also, since the “typical scale length”
is fixed, solutions with unusually large values for seemingly diverged unknowns
may in fact not be diverged at all.
iv. Convergence progress to output
The Newton solver has, since long, a “verbose” setting which makes it print information about the convergence of each unknown quantity to the terminal, colour coded to indicate whether the unknown quantity is considered converged or not. With this new feature, it is possible to also store this information to the output file. To enable this functionality, the user must specify
ds.solver.setSaveConvergenceInfo(True)
in the settings. The resulting output is stored under solver/convergence,
with the 2-norms of the solution x and step dx given in each iteration
and time step. The solver/convergence group also contains the absolute and
relative tolerances used for determining convergence.
v. DREAM as a Python library
by Mathias Hoppe
Normally when we run DREAM, we use the Python interface to generate an HDF5
file with all desired settings and run the program build/iface/dreami with
the HDF5 file as its argument (most people perhaps just call the Python function
runiface(), but it just does all of this under the hood).
A nice feature of the old kinetic solver CODE was that it could plot the
evolution of the distribution function in real time, as the simulation was
progressing. Development of a graphical tool capable of visualizing the
various unknowns throughout DREAM simulations, in real time, as started a few
years ago and required a tighter integration between the Python code and the
C++ code. Since DREAM was designed as a library already from the beginning,
this only required a replacement of the code under iface, as well as some
new code for translating a Python dict into a C++ DREAM::Settings
object. Although this new Python interface was completed a few years ago, the
graphical interface was not, and so both have remained on a separate development
branch. With the recent merge, the new Python interface will therefore be
available to users while the graphical interface remains under development.
So what is so great about a tighter integration between the Python and C++ code? In principle, it allows user-defined Python functions to be executed in the C++ code, which opens a world of new possibilities. Already, the clearest example is perhaps the possibility to defer specification of the max simulation time to user-defined Python routine which is called after each time step. The Python routine has access to the values of all unknown quantities and can, based on that, make a decision about whether to stop the simulation or proceed. This, for example, allows the user to run the simulation until the maximum RE current is achieved, or until the ohmic current reaches a desired value. The new Python interface also allows monitoring of unknown quantities throughout the simulation.
To use the new Python interface, DREAM must first be compiled as a Python library. This can be done while also compiling the conventional DREAM interface, so there is no trade-off to enabling the new interface:
$ cmake .. -DDREAM_BUILD_PYFACE=ON
$ make
Once the DREAM Python library has been compiled, it can be used using the code
below. Two examples of how to use the library are also available under
examples/pyface and examples/pyface-terminate.
import sys
sys.path.append('../../') # To import 'dreampyface'
sys.path.append('../../build/dreampyface/cxx/') # To import C++ library
import dreampyface # Automatically imports C++ library
from DREAM import DREAMSettings
ds = DREAMSettings()
...
def callback(sim):
print(f'Simulation time: {sim.getCurrentTime()}')
# Register a Python function to be called during the simulation
dreampyface.register_callback_timestep_finished(callback)
# Prepare a simulation object
s = dreampyface.setup_simulation(ds)
uinfo = s.unknowns.getInfo()
for uq, info in uinfo.getInfo():
print('{uq:12s} {info["nelements"]:8d} {info["description"]}')
# Run the simulation
do = s.run()
# do is a normal DREAMOutput object
II. Bug fixes
i. Minor bug with numerical equilibrium
by Mathias Hoppe
A small bug was discovered in the flux surface averaging module. The bug manifested itself in numerical geometry in special cases where the two bounce points were found to coincide (i.e. on the magnetic axis). Since the quadrature routines used in GSL do not allow the limits of integration to be identical, they would throw an error and stop the simulation. Hence, this bug would cause simulations not to start and so does not affect simulations which run to completion.
The bug was fixed in pull request #255 and is now available on master.
ii. Make RR heat diffusion work with infinite major radius
by Oskar Vallhagen
After modifying the Rechester-Rosenbluth heat diffusion operator (as noted in the previous newsletter) to include a major radius factor, it was discovered that the special case in which an infinite major radius is specified, in analytic toroidal geometry, does not work. This has now been fixed by letting when .
This bug was fixed in pull request #251 and is now available on master.
iii. Lower limit of AMJUEL ionization losses
by Oskar Vallhagen and Ida Ekmark
The way DREAM handles ionization losses caused the losses to some times be smaller than the ionization energy, which should be the lowest possible energy lost per ionization event, when using AMJUEL ionization coefficients. A lowest possible ionization energy loss has now been implemented.
This bug was fixed in pull request #259 and is now available on master.
iv. Fix warning in CMake scripts for PETSc
by Mathias Hoppe
With newer versions of CMake, warnings are emitted when the FindPETSc.cmake
script is called (during the build process). The warnings are related to the
default settings used within CMake, and it was possible to simply increase the
default compatibility version number.
This bug was fixed in pull request #264 and is now available on master.