Ions

Ions play an important role in runaway electron simulations. They are often the main source of pitch angle scattering, cause electrons to lose energy through the production of bremsstrahlung and provide a source of free electrons as they ionize. The ions also affect the energy balance of the plasma as they further ionize and recombine, which is particularly important during tokamak disruptions.

In DREAM, the user can specify an arbitrary number of ions to include in the simulation. During the simulation DREAM also keeps track of the charge states of the ions which enters into calculations of the plasma energy balance, as well as the interaction of electrons with partially ionized atoms.

Warning

Since ions and electrons are specified separately (the latter by prescribing either the electron distribution function or cold electron density), it is possible to set up a simulation in which ions and electrons are unbalanced, thus breaking quasi-neutrality. In kinetic mode, DREAM will always preserve quasi-neutrality by adding as many cold electrons as necessary, which may not be what the user desires.

To ensure that the number of electrons and ions agree, use the getFreeElectronDensity() method in the Ions class after specifying the ions for the simulation, and use this density to initialize the electrons in the simulation.

Ion models

DREAM provides three different modes for evolving ions. It is possible to evolve different ions species with different modes in the same simulation:

  1. Prescribe the evolution in time and space of the ion densities and charge states.

  2. Assume that ion charge states are distributed as to be in equilibrium in every time step of the simulation.

  3. Evolve the ion charge states using an ion rate equation.

In mode (1) it is up to the user to provide the full evolution of the ions. Mode (3) evolves the ion charge states by solving the ion rate equation

\[\begin{split}\frac{\partial n_i^{(j)}}{\partial t} = \left( I_i^{(j-1)} n_{\rm cold} + \mathcal{I}_i^{(j-1)} \right) n_i^{(j-1)} - \left( I_i^{(j)} n_{\rm cold} + \mathcal{I}_i^{(j)} \right) n_i^{(j)}\\ + R_i^{(j+1)} n_i^{(j+1)} n_{\rm cold} - R_i^{(j)} n_i^{(j)} n_{\rm cold}\end{split}\]

where \(n_i^{(j)}\) denotes the density of ion species \(i\) in charge state \(j\), \(I_i^{(j)}\) is the rate at which ion species \(i\) ionizes from charge state \(j\) to charge state \(j+1\), \(R_i^{(j)}\) is the rate at which ion species \(i\) recombines from charge state \(j\) to charge state \(j-1\), and \(\mathcal{I}_i^{(j)}\) is the rate at which ion species \(i\) in charge state \(j\) is ionized due to collisions with fast electrons. The ionization and recombination rates \(I_i^{(j)}\) and \(R_i^{(j)}\) are taken from the OPEN-ADAS database <https://open.adas.ac.uk/>.

In ion mode (2), the above ion rate equation is also solved, but with \(\partial n_i^{(j)} / \partial t = 0\) such that the equilibrium solution is sought. In this case, the ion densities can be controlled (via prescribing the values or radial transport) by the total ion density \(n_i = \sum_j n_i^{(j)}\).

The available ion modes are:

Name

Description

IONS_PRESCRIBED

The ion densities and charge states are prescribed in time and space.

IONS_EQUILIBRIUM

Ion charge states are assumed to be in equilibrium.

IONS_DYNAMIC

Ion charge states are evolved according to an ion rate equation.

The above modes are the only three models available for evolving ions in DREAM. However, to simplify initialization of ion densities, the following pseudo-modes are also available when adding ions to the simulation:

Name

Description

IONS_DYNAMIC_NEUTRAL

Ion densities are evolved dynamically, with all ions located in the neutral charge state initially.

IONS_DYNAMIC_FULLY_IONIZED

Ion densities are evolved dynamically, with all ions located in the fully ionized charge state initially.

IONS_PRESCRIBED_NEUTRAL

Ion densities are prescribed, with all ions located in the neutral charge state.

IONS_PRESCRIBED_FULL_IONIZED

Ion densities are prescribed, with all ions located in the fully ionized charge state.

IONS_EQUILIBRIUM_NEUTRAL

Ion densities are assumed in equilibrium at every time step, with all ions located in the neutral charge state initially.

IONS_EQUILIBRIUM_FULLY_IONIZED

Ion densities are evolved in equilibrium at every time step, with all ions located in the fully ionized charge state initially.

Adding ions

Ion species are added to a DREAM simulation one by one using the addIon() method. When adding a new ion species, one must specify the name of the species (which must be unique), the atomic charge \(Z\), which equation to use for evolving the ion densities, as well as any initial or prescribed ion densities. Since DREAM uses a user-specified name parameter to keep track of ion species, it is possible to add several ion populations of the same species and evolve them independently of each other.

When initializing ion densities it is also possible to specify which charge state \(Z_0\) should have a non-zero density. This is done using the adIon() input parameter Z0. If Z0 is not specified, and the iontype parameter does not specify the charge state, the user must provide the input density n as an array with each element of the first dimension representing the density in the corresponding charge state.

When adding a tritium ion species the tritium parameter should also be set to True. This allows DREAM to correctly identify the species as tritium and use its density when calculating the tritium decay runaway rate.

Warning

If you would like to add ions to a restart simulation then you will not be able to initialize the ion densities from the previous output. Instead you will have to prescribe the initial ion densities again.

The reason why ion densities cannot be initialized from a previous DREAM output if new ions are added is because DREAM is not able to determine if ions are added, removed or renamed. This makes it impossible for DREAM to keep track of which ion densities in the previous output map to which ions in the new input.

Example

Prescribed ions

The following example illustrates how to prescribe the density to be constant and uniform for an ion species:

import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
...
ds.eqsys.n_i.addIon(name='D', Z=1, iontype=Ions.IONS_PRESCRIBED, Z0=1, n=2e19)

For more advanced use cases, one can also prescribe the ion densities for each charge state, at a set of radii and at a set of time points.

import numpy as np
import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
...
# Ion name and atomic charge
name = 'Ar'
Z    = 18

# Define radial and time grids
r = np.linspace(r0, r1, nr)
t = np.linspace(t0, t1, nt)

# Construct ion densities
n = np.array([...])      # Shape (Z+1, nt, nr)

# n[0,:] = Neutral charge state
# ...
# n[Z,:] = Fully ionized charge state

ds.eqsys.n_i.addIon(name=name, Z=Z, iontype=Ions.IONS_PRESCRIBED, n=n, r=r, t=t)

Note

All ion densities must be prescribed on the same radial/time grid. Scalar ions densities can still be prescribed, must be so after prescribing at least one ion species with the non-scalar radial/time grid.

Dynamic ions

Dynamic ions can be added similarly to prescribed ions, but they are instead evolved using the ion rate equations. The following example illustrates how to add an ion species that is evolved using the ion rate equations and is initialized with a uniform radial density profile with all ions in the fully ionized charge state:

import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
...
ds.eqsys.n_i.addIon(name='D', Z=1, iontype=Ions.IONS_DYNAMIC, Z0=1, n=2e19)

Alternatively, the ion charge state densities can be directly prescribed. The following example illustrates how to prescribe the initial ion charge state density profiles and evolve them using the ion rate equations:

import numpy as np
import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
...
# Ion name and atomic charge
name = 'Ar'
Z    = 18

# Define radial grid
r = np.linspace(r0, r1, nr)

# Construct ion densities
n = np.array([...])      # Shape (Z+1, nr)

# n[0,:] = Neutral charge state
# ...
# n[Z,:] = Fully ionized charge state

ds.eqsys.n_i.addIon(name=name, Z=Z, iontype=Ions.IONS_DYNAMIC, n=n, r=r)

Note

All ion densities must be prescribed on the same radial/time grid. Scalar ions densities can still be prescribed, must be so after prescribing at least one ion species with the non-scalar radial/time grid.

Ion temperature

When the temperature in the plasma is solved for self-consistently, the default mode is to model only the electron temperature evolution due to radiation, heating and transport. It is possible to also include the evolution of the temperature of each ion species, where different charge states of the same species are assumed to have the same temperature. Ions and electrons exchange energy via elastic collisions, taking the form of a rate equation

\[\frac{\partial W_i}{\partial t} = \sum_j Q_{ij} + Q_{ie},\]

where the sum over \(j\) is taken over all ion species, and \(Q_{ij}\) is the collisional energy transfer integrated over Maxwellians of different density \(N_i = \sum_j n_i^{(j)}\) and heat \(W_i = 3 e N_i T_i / 2\),

\[Q_{ij} = \sqrt{\frac{3}{\pi}} \frac{Z_i^2 Z_j^2 e^4 \ln\Lambda_{ij} N_i N_j\sqrt{m_i N_i m_j N_j} }{4\pi \varepsilon_0^2} \frac{N_i W_j - N_j W_i}{(m_j N_j W_i + m_i N_i W_j)^{3/2}}.\]

In addition, the cold-electron temperature will pick up a similar contribution,

\[\left( \frac{\partial W_\mathrm{cold}}{\partial t} \right)_Q = \sum_i Q_{ei}\]

where the anti-symmetry of \(Q_{ij} = -Q_{ji}\) ensures that the sum of the heat equations for all ions and electrons exhibits energy conservation by these elastic collisions.

The initialization and behavior of ion temperatures are controlled via the T argument in Ions::addIon(..., T). It behaves as follows:

  • If at least one T is explicitly set, DREAM will add the additional quantities N_i and W_i to the equation system and evolve them as nontrivial unknowns.

  • Species which are not explicitly set will be initialized to T=0.

  • If the type for T_cold is set to TYPE_SELFCONSISTENT, the ion and electron heat W_i and W_cold will be evolved according to the equations above. If TYPE_PRESCRIBED it will be given by its initial value W_i=constant.

Example

In the below example, we consider a scenario where a hydrogenic population from before a disruption has the same initial temperature as the electrons, whereas an injected hydrogenic species and Neon impurity is introduced at effectively zero temperature. The ions and electrons will be allowed to exchange energy via collisions. All ion and electron initial temperatures are taken to be uniform in radius.

import numpy as np
import DREAM.Settings.Equations.IonSpecies as Ions
import DREAM.Settings.Equations.ColdElectronTemperature as T_cold

ds = DREAMSettings()

...

T_initial = 4e3 #eV

ds.eqsys.T_cold.setPrescribedData(T_initial)
ds.eqsys.n_i.addIon(name='D',     Z=1,  T=T_initial, iontype=Ions.IONS_DYNAMIC_FULLY_IONIZED, n=1e20)
ds.eqsys.n_i.addIon(name='Ne',    Z=10, iontype=Ions.IONS_DYNAMIC_NEUTRAL, n=1e19)
ds.eqsys.n_i.addIon(name='D_inj', Z=1,  iontype=Ions.IONS_DYNAMIC_NEUTRAL, n=1e21)

ds.eqsys.T_cold.setType(ttype=T_cold.TYPE_SELFCONSISTENT)

Start from coronal equilibrium

Coronal equilibrium refers to a distribution of charge states of an ion species such that ionization and recombination exactly balance for all charge states. It can be shown (see doc/notes) that the density \(n_i^{(l)}\) of species \(i\) in charge state l is given by

\[n_i^{(l)} = n_i\left( 1 + \sum_{j=0}^{l-1}\prod_{k=j+1}^l \frac{R_i^{(k)}}{I_i^{(k-1)}} + \sum_{j=l+1}^Z\prod_{k=l}^{j-1} \frac{I_i^{(k)}}{R_i^{(k+1)}} \right)^{-1},\]

where \(n_i=\sum_l n_i^{(l)}\).

DREAM simulations can be launched with ions in coronal equilibrium initially. This is specified via a flag init_equil=True to the addIon() call, along with the desired species density \(n_i\). The code will then automatically evaluate the equilibrium state at startup, iterating until \(n_{\rm free}\) (upon which the ionization/recombination rates depend).

Example

To initiate one ion species in coronal equilibrium, the following syntax is used:

import numpy as np
import DREAM.Settings.Equations.IonSpecies as Ions
import DREAM.Settings.Equations.ColdElectronTemperature as T_cold

ds = DREAMSettings()

...

ds.eqsys.n_i.addIon(name='D',  Z=1,  iontype=Ions.IONS_DYNAMIC_FULLY_IONIZED, n=1e20)
ds.eqsys.n_i.addIon(name='Ne', Z=10, iontype=Ions.IONS_DYNAMIC, n=1e19, init_equil=True)

Note

Initialization in coronal equilibrium is only possible for dynamically evolved ion species.

Atomic data

Various types of atomic data are downloaded from ADAS and NIST during the configuration phase of DREAM (prior to compilation). The data is downloaded using the scripts tools/get_adas.py and tools/get_nist.py, which should be automatically invoked by CMake.

ADAS

The Atomic Data and Analysis Structure (ADAS) is an interconnected set of computer codes and data collections for modelling the radiating properties of ions and atoms in plasmas. It can address plasmas ranging from the interstellar medium through the solar atmosphere and laboratory thermonuclear fusion devices to technological plasmas. ADAS assists in the analysis and interpretation of spectral emission and supports detailed plasma models.

—ADAS website, https://www.adas.ac.uk/about.php

DREAM utilizes data for four types of coefficients from the OpenADAS database (https://open.adas.ac.uk), name ACD (effective recombination coefficients), SCD (effective ionization coefficients), PLT (line power driven by excitation of dominant ions) and PRB (continuum and line power driven by recombination and bremsstrahlung of dominant ions). Data is downloaded for the elements specified in the file tools/elements.json. The configuration file is given in JSON format with one entry per element in a key-value format. The key should be the official name of the element and the value is the year in which the dataset to use was published.

The get_adas.py script can also be run manually, separately from the DREAM CMake script to generate a C/C++ file of ADAS data. The script is then invoked from the command line and accepts the following arguments:

Name

Description

--cachedir DIR

Use DIR to store raw data downloaded from ADAS. This can be re-used later to avoid making an HTTP request to Open-ADAS.

--elements FILE

Load elements to fetch data for from the JSON file FILE.

--no-cache

Do not store raw ADAS data locally.

--no-compile

Do not generate C++ source files with the rate coefficients.

-o, --output

Name of output C++ source file to generate.

--type-int

C++ type to use for integers.

--type-real

C++ type to use for real numbers.

Opacity and AMJUEL

The ADAS-coefficients are calculated in the low density limit, corresponding to a completely transparent plasma, an assumption which might not always be valid. At high densities, the plasma can become opaque to line radiation resulting from transistions involving the ground state. For example, this might be the case for the Lyman radiation from hydrogen isotopes after a massive injection of the type planned for ITER disruption mitigation. In that case, it might be more accurate to use coefficients where the bound-bound transitions involving the ground state are removed. For hydrogen isotopes, such coefficients are available in the AMJUEL database. Thus, there are two opacity modes available:

Name

Description

ION_OPACITY_MODE_TRANSPARENT

Use ADAS coefficients, corresponding to a completely transparent plasma

ION_OPACITY_MODE_GROUND_STATE_OPAQUE

Use coefficients disregarding bound-bound transitions involving the ground state, assuming a complete opacity to those lines

Note

The ION_OPACITY_MODE_GROUND_STATE_OPAQUE is currently only supported for hydrogen species.

Example

The following example illustrates how to create a deuterium ion species whith radiation, ionization and recombination calculated assuming opacity to Lyman radiation:

ds.eqsys.n_i.addIon(name='D', Z=1, iontype=Ions.IONS_DYNAMIC_NEUTRAL, n=40e20, opacity_mode=Ions.ION_OPACITY_MODE_GROUND_STATE_OPAQUE)

NIST

DREAM uses ionization and binding energy values tabulated in the database of the american National Institute of Standards (NIST) database. The values are downloaded by making a request via the form https://physics.nist.gov/PhysRefData/ASD/ionEnergy.html and are then embedded into a C++ source file. The elements to download are specifed by name in an array in the main() function of the script tools/get_nist.py.

While the tools/get_nist.py script should be automatically invoked when running CMake for DREAM, it can also be run manually. The script accepts the following command-line arguments:

Name

Description

--cachedir DIR

Use DIR to store raw data downloaded from ADAS. This can be re-used later to avoid making an HTTP request to Open-ADAS.

--ionization

Downloads ionization energy data instead of the default binding energy data.

--no-cache

Do not store raw ADAS data locally.

--no-compile

Do not generate C++ source files with the rate coefficients.

-o, --output

Name of output C++ source file to generate.

--type-int

C++ type to use for integers.

--type-real

C++ type to use for real numbers.

Free electron density

When running kinetic simulations, ẗhe density of all ion species, as well as of electrons, are locally conserved. The net charge of the plasma is also strictly conserved. Therefore, in order to satisfy quasi-charge neutrality, \(n_e = \sum_{ij} Z_{0j}n_i^{(j)}\), the initial plasma must do so.

When setting the initial electron distribution, the ions have been equipped with the function getFreeElectronDensity() to help the user set the correct density value of the electrons.

Example

An example of how to initialize the electron distribution to a Maxwellian with the correct value for the electron density is provided by the following:

T_initial = 100 # eV
n0_initial, rn0 = ds.eqsys.n_i.getFreeElectronDensity()

ds.eqsys.f_hot.setInitialProfiles(n0=n_initial, rn0=rn0, T0=T_initial)

Warning

The method getFreeElectronDensity() must be called after all initial ion densities have been prescribed.

Ionization models

When prescribing initial ion densities as DYNAMIC, their charge states will be evolved via rate equations including the effects of ionization and recombination. In DREAM, cold electrons will contribute ionization via ADAS rate coefficients which have been averaged over a Maxwellian distribution. There is also the option to include a kinetic model of ionization due to the fast (non-thermal) electrons, valid up to arbitrarily high energies. The kinetic ionization model is described in N A Garland et al, Phys Plasmas 27, 040702 (2020), but where DREAM has fitted parameters in the model in order to reproduce the ADAS data as accurately as possible.

The ionization model to be used is controlled with the help of three settings:

Name

Description

IONIZATION_MODE_FLUID

Only cold electrons contribute to ionization, via ADAS rate coefficients.

IONIZATION_MODE_KINETIC

All electrons in the f_hot and f_re distributions contribute to ionization, employing the full jacobian.

IONIZATION_MODE_KINETIC_APPROX_JAC

All electrons in the f_hot and f_re distributions contribute to ionization, but only the FLUID jacobian is used

In a simulation with many ion species (for example including one hydrogen species and one argon species, yielding N=21), the jacobian matrix in mode KINETIC will pick up a large number of non-zero elements, which will generally make simulations extremely heavy. Therefore, it is recommended to include kinetic ionization effects via KINETIC_APPROX_JAC unless this clearly leads to convergence issues.

Example

Kinetic ionization can be activated in a DREAM simulation as follows:

import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
ds.eqsys.n_i.setIonization(Ions.IONIZATION_MODE_KINETIC_APPROX_JAC)

Particle injection

Particles can be injected through the plasma edge at a specified rate. The user specifies the number of particles per unit time which are to be injected, which can either be specified as constant in time or as following a specified time evolution. The user can choose to either inject particles of a specific charge state, or to inject particles of multiple charge states simultaneously.

Example

The ion source is enabled via a single line:

from DREAM import DREAMSettings
import DREAM.Settings.Equations.IonSpecies as Ions
import numpy as np

ds = DREAMSettings()
...
ds.eqsys.n_i.addIon(name='D', Z=1, iontype=Ions.IONS_DYNAMIC, Z0=1, n=5e19)
# Enable the ion source for species 'D'
ds.eqsys.n_i.addIonSource('D', dNdt=1e19)    # Add 1e19 deuterium atoms per second

The charge state of the injected particles is \(Z_0=0\), but can be modified with the Z0 argument to addIonSource():

...
ds.eqsys.n_i.addIonSource('D', dNdt=1e19, Z0=1)  # Add 1e19 deuterium ions per second

The full time evolution of the ion source can also be provided:

...
t = np.linspace(0, 1)
dNdt = 5e18 * ((1-t)**2 + 1)
ds.eqsys.n_i.addIonSource('D', dNdt=dNdt, t=t, Z0=1)

Finally, particles in several charge states may be added simultaneously:

...
Z = 1
t = np.linspace(0, 1)
dNdt = np.zeros((Z+1, t.size))
# D(Z0=0)
dNdt[0,:] = 5e18 * ((1-t)**2 + 1)
dNdt[0,:] = 5e18 * (t**2 + 1)

Note

The injection rate specified is the number of particles per unit time which is to be injected. This means that the number dNdt specified equals

\[dNdt = \frac{1}{\Delta t}\int \left[ n_i^{(j)}(r,t+\Delta t) - n_i^{(j)}(r,t) \right]V'\,\mathrm{d}r\]

i.e. the change in the volume integral of the ion density per unit time. Thus, the unit of dNdt is simply \(\mathrm{s}^{-1}\) (or particles per second).

Ion/neutral transport

Todo

Document the ion/neutral transport functionality

Tritium

Runaway electrons can be generated when tritium decays into helium-3 via the beta decay process

\[\mathrm{T} \to \ ^3_2\mathrm{He} + \mathrm{e}^- + \bar{\nu}_{\rm e}.\]

The corresponding runaway rate is given by

\[\left( \frac{\mathrm{d} n_{\rm RE}}{\mathrm{d} t} \right)_{\rm T} \approx \ln 2 \frac{n_{\rm T}}{\tau_{\rm T}} F_\beta\left( \gamma_{\rm c} \right),\]

where \(n_{\rm T}\) is the tritium density, \(\tau_{\rm T} = 4800\pm 8\) days is the tritium half-life, and \(F_\beta(\gamma_{\rm c})\) denotes the fraction of beta electrons generated with an energy above the critical energy \(\gamma_{\rm c}\) for runaway to occur.

Tritium runaway generation is enabled with the n_re 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, as described above under Adding ions.

Note

It is possible to include multiple tritium populations in the simulation, simply by providing tritium=True when adding each of them.

Example

The following example illustrates how to enable the tritium decay runaway mechanism in a DREAM simulation:

import DREAM.Settings.Equations.IonSpecies as Ions

ds = DREAMSettings()
...
# Include source term in equation for n_re
ds.eqsys.n_re.setTritium(True)

# Add tritium ion species to list of ions
ds.eqsys.n_i.addIon('T', Z=1, iontype=Ions.IONS_DYNAMIC, n=2e19, tritium=True)

Connection to SPI

An ion species can be connected to a Shattered Pellet Injection (SPI), so that a specified fraction of the ablated material is added to the density of this species. The newly added material is added with the equilibrium charge state distribution, reflecting the fact that the ionization takes place very rapidly in the usually rather hot plasma into which the pellet is injected. The most general way to connect an ion species to an SPI is to set the argument SPIMolarFraction to a vector specifying the fraction of the particles in each shard added to the ds.eqsys.spi object which consists of the ion species being created. By default, SPIMolarFraction=-1, telling DREAM that this species is not part of an SPI. Additionally, if one wants the pellet to consist of a certain isotope of a certain species, this can also be specified by the argument isotope to the addIon()-method. The default value of isotope is 0, meaning the naturally occuring mix of isotopes.

Example

The following example illustrates how to add ions in a way that makes the injected pellet consist of 95% deuterium and 5% neon.

...
nShardD=1000 # Number of shards
...

# Add ions connected to the shards using the SPIMolarFraction argument
ds.eqsys.n_i.addIon(name='D_inj', Z=1, isotope=2, iontype=Ions.IONS_DYNAMIC_NEUTRAL, n=1e0, SPIMolarFraction=0.95*np.ones(nShard))
ds.eqsys.n_i.addIon(name='Ne', Z=10, iontype=Ions.IONS_DYNAMIC_NEUTRAL, n=1e0, SPIMolarFraction=0.05*np.ones(nShard))

Note

To avoid numerical errors arrising when the density of an ion species is zero, the initial ion density must be set to a small but finite value also for the species related to an SPI, in this case n=1e0.

In some cases, however, it might be more convenient to create the ion species connected to an SPI by passing a pointer to the ds.eqsys.n_i-object to a helper function in the ds.eqsys.spi-class, se the documentation of the SPI settings.

Class documentation

class DREAM.Settings.Equations.Ions.Ions(settings, ionization=1)

Bases: DREAM.Settings.Equations.UnknownQuantity.UnknownQuantity

__init__(settings, ionization=1)

Constructor.

addIon(name, Z, iontype=1, Z0=None, isotope=0, SPIMolarFraction=- 1, opacity_mode=1, charged_diffusion_mode=1, charged_prescribed_diffusion=None, rChargedPrescribedDiffusion=None, tChargedPrescribedDiffusion=None, neutral_diffusion_mode=1, neutral_prescribed_diffusion=None, rNeutralPrescribedDiffusion=None, tNeutralPrescribedDiffusion=None, charged_advection_mode=1, charged_prescribed_advection=None, rChargedPrescribedAdvection=None, tChargedPrescribedAdvection=None, neutral_advection_mode=1, neutral_prescribed_advection=None, rNeutralPrescribedAdvection=None, tNeutralPrescribedAdvection=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, diffusion_initial_all_cs=None, diffusion_final_all_cs=0, diffusion_offset_all_cs=0, advection_initial_all_cs=None, advection_final_all_cs=0, advection_offset_all_cs=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, init_equil=False, T=None, n=None, r=None, t=None, tritium=False, hydrogen=False)

Adds a new ion species to the plasma.

Parameters
  • name (str) – Name by which the ion species will be referred to.

  • Z (int) – Ion charge number.

  • isotope (int) – Ion mass number.

  • iontype (int) – Method to use for evolving ions in time.

  • Z0 (int) – Charge state to populate (used for populating exactly one charge state for the ion).

  • n – Ion density (can be either a scalar, 1D array or 2D array, depending on the other input parameters)

  • SPIMolarFraction (float) – Molar fraction of the SPI injection (if any). A negative value means that this species is not part of the SPI injection

  • r (numpy.ndarray) – Radial grid on which the input density is defined.

  • T – Ion initial temperature (can be scalar for uniform temperature, otherwise 1D array matching r in size)

  • r – Radial grid on which the input density and temperature is defined.

  • t (numpy.ndarray) – Time grid on which the input density is defined.

  • tritium (bool) – If True, the ion species is treated as Tritium.

  • hydrogen (bool) – If True, the ion species is treated as Hydrogen (single proton).

addIonSource(species, dNdt=None, t=None, Z0=0)

Add a source term for the specified ion species.

Parameters
  • species – Name of the ion species to add this source term to.

  • dNdt – Number of particles to add per unit time.

  • t – Time grid associated with dNdt (if any).

  • Z0 – For scalar or 1D dNdt, the charge state for which to add the source term.

changeRadialGrid(r)

Change the radial grid used for the ion species.

fromdict(data)

Load settings from the specified dictionary.

Parameters

data (dict) – Dictionary containing all settings to load.

getChargedAdvectionModes()

Returns a list of ion charged advection modes for the various ion species contained by this object.

getChargedDiffusionModes()

Returns a list of ion charged diffusion modes for the various ion species contained by this object.

getCharges()

Returns a list of the charges of the various ion species contained by this object.

getFreeElectronDensity(t=0)

Returns the plasma free electron density at the given time index, based on the prescribed/initialized ion densities.

Parameters

t (int) – Index of time for which to retrieve the free electron density.

getHydrogenSpecies()

Returns a list of names of the ion species which are treated as Hydrogen.

getIndex(species)

Return the index of the ion species with the specified name.

getIon(i=None)

Returns the ion species with the specified index or name.

Parameters

i – Index or name of ion species to retrieve.

getIsotopes()

Returns a list of the isotopes of the various ion species contained by this object.

getNeutralAdvectionModes()

Returns a list of ion neutral advection modes for the various ion species contained by this object.

getNeutralDiffusionModes()

Returns a list of ion neutral diffusion modes for the various ion species contained by this object.

getOpacityModes()

Returns a list of ion opacity modes for the various ion species contained by this object.

getSPIMolarFraction()

Returns a list of the SPI molar fractions of the various ion species contained by this object.

getTritiumSpecies()

Returns a list of names of the ion species which are treated as Tritium.

getTypes()

Returns a list of ion types for the various ion species contained by this object.

setAdvectionInterpolationMethodCharged(ad_int=1, ad_jac=2, fluxlimiterdamping=1.0)

Sets the interpolation method that is used in the charged advection terms of the transport equation.

Parameters
  • 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).

setAdvectionInterpolationMethodNeutral(ad_int=1, ad_jac=2, fluxlimiterdamping=1.0)

Sets the interpolation method that is used in the neutral advection terms of the transport equation.

Parameters
  • 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).

setChargedAdvection(species, mode, Ar=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Ar_0=None, Ar_f=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial advection for charged particles.

Parameters
  • species – Species to apply transport to.

  • mode – Type of transport to prescribe.

  • Ar – Advection coefficient to prescribe.

  • r – Radial grid on which Ar is given.

  • t – Time grid on which Ar is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Ar_0 – Initial value of advection coefficient when decaying exponentially.

  • Ar_f – Final value of advection coefficient when decaying exponentially.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setChargedDiffusion(species, mode, Drr=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Drr_0=None, Drr_f=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial diffusion for charged particles.

Parameters
  • species – Species to apply transport to.

  • mode – Type of transport to prescribe.

  • Drr – Diffusion coefficient to prescribe.

  • r – Radial grid on which Drr is given.

  • t – Time grid on which Drr is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Drr_0 – Initial value of diffusion coefficient when decaying exponentially.

  • Drr_f – Final value of diffusion coefficient when decaying exponentially.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setIonType(index, ttype)

Modifies the type of equation used for the specified ion species.

Parameters
  • index – Index or name of ion species to set type for.

  • ttype (int) – Type of equation to use for evolving the ion species.

setIonization(ionization=1)

Sets which model to use for ionization.

Parameters

ionization (int) – Flag indicating which model to use for ionization.

setNeutralAdvection(species, mode, Ar=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Ar_0=None, Ar_f=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial advection for neutral particles.

Parameters
  • species – Species to apply transport to.

  • mode – Type of transport to prescribe.

  • Ar – Advection coefficient to prescribe.

  • r – Radial grid on which Ar is given.

  • t – Time grid on which Ar is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Ar_0 – Initial value of advection coefficient when decaying exponentially.

  • Ar_f – Final value of advection coefficient when decaying exponentially.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setNeutralDiffusion(species, mode, Drr=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Drr_0=None, Drr_f=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial diffusion for neutral particles.

Parameters
  • species – Species to apply transport to.

  • mode – Type of transport to prescribe.

  • Drr – Diffusion coefficient to prescribe.

  • r – Radial grid on which Drr is given.

  • t – Time grid on which Drr is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Drr_0 – Initial value of diffusion coefficient when decaying exponentially.

  • Drr_f – Final value of diffusion coefficient when decaying exponentially.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

shiftTimeTranspCoeffs(tShift)

Shift the time grids for the ion transport coefficients by an amount tShift. This is needed between restarts.

Parameters

tShift (float) – Amount of time the time grids for the ion transport coefficient should be shifted

todict()

Returns a Python dictionary containing all settings of this Ions object.

verifySettings()

Verify that all settings are consistent.

class DREAM.Settings.Equations.IonSpecies.IonSpecies(settings, name, Z, ttype=0, Z0=None, isotope=0, SPIMolarFraction=- 1.0, opacity_mode=1, charged_diffusion_mode=1, charged_prescribed_diffusion=None, rChargedPrescribedDiffusion=None, tChargedPrescribedDiffusion=None, neutral_diffusion_mode=1, neutral_prescribed_diffusion=None, rNeutralPrescribedDiffusion=None, tNeutralPrescribedDiffusion=None, charged_advection_mode=1, charged_prescribed_advection=None, rChargedPrescribedAdvection=None, tChargedPrescribedAdvection=None, neutral_advection_mode=1, neutral_prescribed_advection=None, rNeutralPrescribedAdvection=None, tNeutralPrescribedAdvection=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, diffusion_initial_all_cs=None, diffusion_final_all_cs=0, diffusion_offset_all_cs=0, advection_initial_all_cs=None, advection_final_all_cs=0, advection_offset_all_cs=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, init_equil=False, T=None, n=None, r=None, t=None, interpr=None, interpt=None, tritium=False, hydrogen=False)

Bases: object

__init__(settings, name, Z, ttype=0, Z0=None, isotope=0, SPIMolarFraction=- 1.0, opacity_mode=1, charged_diffusion_mode=1, charged_prescribed_diffusion=None, rChargedPrescribedDiffusion=None, tChargedPrescribedDiffusion=None, neutral_diffusion_mode=1, neutral_prescribed_diffusion=None, rNeutralPrescribedDiffusion=None, tNeutralPrescribedDiffusion=None, charged_advection_mode=1, charged_prescribed_advection=None, rChargedPrescribedAdvection=None, tChargedPrescribedAdvection=None, neutral_advection_mode=1, neutral_prescribed_advection=None, rNeutralPrescribedAdvection=None, tNeutralPrescribedAdvection=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, diffusion_initial_all_cs=None, diffusion_final_all_cs=0, diffusion_offset_all_cs=0, advection_initial_all_cs=None, advection_final_all_cs=0, advection_offset_all_cs=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, init_equil=False, T=None, n=None, r=None, t=None, interpr=None, interpt=None, tritium=False, hydrogen=False)

Constructor.

Parameters
  • settings (DREAMSettings) – Parent DREAMSettings object.

  • name (str) – Name by which the ion species will be referred to.

  • Z (int) – Ion charge number.

  • isotope (int) – Ion mass number.

  • ttype (int) – Method to use for evolving ions in time.

  • Z0 (int) – Charge state to populate with given density.

  • n (float) – Ion density (can be either a scalar, 1D array or 2D array, depending on the other input parameters)

  • SPIMolarFraction (float) – Molar fraction of the SPI injection (if any). A negative value means that this species is not part of the SPI injection

  • init_equil (bool) – Initialize ion species in coronal equilibrium.

  • T (float) – Ion initial temperature (can be scalar for uniform temperature, otherwise 1D array matching r in size)

  • r (numpy.ndarray) – Radial grid on which the input density is defined.

  • t (numpy.ndarray) – Time grid on which the input density is defined.

  • interpr (numpy.ndarray) – Radial grid onto which ion densities should be interpolated.

  • interpt (numpy.ndarray) – Time grid onto which ion densities should be interpolated.

  • tritium (bool) – If True, this ion species is treated as Tritium.

  • hydrogen (bool) – If True, this ion species is treated as Hydrogen.

calcTransportCoefficientExpdecayAllChargedStates(t_exp, c0, cf=0, co=0, t_start=0, r=None, t=None)

Construct transport coefficients which decay exponentially in time, for all charge states of this species.

calcTransportCoefficientExpdecaySingleChargeState(t_exp, c0, cf=0, co=0, t_start=0, r=None, t=None)

Contsruct a transport coefficient which decays exponentially in time according to

c = co + cf + (c0-cf)*exp(-(t-t_start)/t_exp).

The coefficient is set to ‘co’ for t <= t_start.

getChargedAdvectionMode()

Returns the charged advection mode to use for evolving the ion densities for this species.

getChargedDiffusionMode()

Returns the charged diffusion mode to use for evolving the ion densities for this species.

getChargedPrescribedAdvection()

Returns the prescribed charged advection coefficient array for this ion species.

getChargedPrescribedDiffusion()

Returns the prescribed charged diffusion coefficient array for this ion species.

getDensity()

Returns the prescribed density array for this ion species.

getInitialSpeciesDensity()

Returns the initial density of all charge states combined (used when initializing the species in coronal equilibrium).

getIsotope()
getName()

Returns the name of this ion species.

getNeutralAdvectionMode()

Returns the neutral advection mode to use for evolving the ion densities for this species.

getNeutralDiffusionMode()

Returns the neutral diffusion mode to use for evolving the ion densities for this species.

getNeutralPrescribedAdvection()

Returns the prescribed neutral advection coefficient array for this ion species.

getNeutralPrescribedDiffusion()

Returns the prescribed neutral diffusion coefficient array for this ion species.

getOpacityMode()

Returns the opacity mode to use for evolving the ion densities for this species.

getR()

Returns the radial grid on which the ion densities are defined.

getRChargedPrescribedAdvection()

Returns the radial grid for the prescribed charged advection coefficient array for this ion species.

getRChargedPrescribedDiffusion()

Returns the radial grid for the prescribed charged diffusion coefficient array for this ion species.

getRNeutralPrescribedAdvection()

Returns the radial grid for the prescribed neutral advection coefficient array for this ion species.

getRNeutralPrescribedDiffusion()

Returns the radial grid for the prescribed neutral diffusion coefficient array for this ion species.

getSPIMolarFraction()
getSourceDensity()

Returns the (time evolution of the) source term used to seed this ion species with new particles.

getSourceTime()

Returns the time grid on which the ion source is defined.

getSourceType()

Returns the type of the source term used.

getTChargedPrescribedAdvection()

Returns the time grid for the prescribed charged advection coefficient array for this ion species.

getTChargedPrescribedDiffusion()

Returns the time grid for the prescribed charged diffusion coefficient array for this ion species.

getTNeutralPrescribedAdvection()

Returns the time grid for the prescribed neutral advection coefficient array for this ion species.

getTNeutralPrescribedDiffusion()

Returns the time grid for the prescribed neutral diffusion coefficient array for this ion species.

getTemperature()

Returns the initial temperature array to use for evolving the ion heat of this species

getTime()

Returns the time grid on which the ion densities are defined.

getType()

Returns the type of equation to use for evolving the ion densities for this species.

getZ()

Returns the atomic charge for this ion species.

initializeToEquilibrium()

Returns a flag indicating whether or not this ion species should be initialized according to coronal equilibrium.

initialize_charged_prescribed_advection(charged_prescribed_advection=None, rChargedPrescribedAdvection=None, tChargedPrescribedAdvection=None, interpr=None, interpt=None)

Prescribes the evolution of the charged advection coefficients for this ion species.

initialize_charged_prescribed_diffusion(charged_prescribed_diffusion=None, rChargedPrescribedDiffusion=None, tChargedPrescribedDiffusion=None, interpr=None, interpt=None)

Prescribes the evolution of the charged diffusion coefficients for this ion species.

initialize_dynamic(n=None, r=None, init_equil=False, interpr=None)

Evolve ions according to the ion rate equation in DREAM.

initialize_dynamic_charge_state(Z0, n=None, r=None, interpr=None, init_equil=False)

Evolve the ions dynamically, initializing them all to reside in the specified charge state Z0.

initialize_dynamic_fully_ionized(n=None, r=None, interpr=None, init_equil=False)

Evolve the ions dynamically, initializing them all as fully ionized.

initialize_dynamic_neutral(n=None, r=None, interpr=None, init_equil=False)

Evolve the ions dynamically, initializing them all as neutrals.

initialize_equilibrium(n=None, r=None, interpr=None)

Evolve ions according to the equilibrium equation in DREAM.

initialize_neutral_prescribed_advection(neutral_prescribed_advection=None, rNeutralPrescribedAdvection=None, tNeutralPrescribedAdvection=None, interpr=None, interpt=None)

Prescribes the evolution of the neutral advection coefficients for this ion species.

initialize_neutral_prescribed_diffusion(neutral_prescribed_diffusion=None, rNeutralPrescribedDiffusion=None, tNeutralPrescribedDiffusion=None, interpr=None, interpt=None)

Prescribes the evolution of the neutral diffusion coefficients for this ion species.

initialize_prescribed(n=None, r=None, t=None)

Prescribes the evolution for this ion species.

initialize_prescribed_charge_state(Z0, n=None, r=None, t=None, interpr=None, interpt=None)

Prescribe the ions to all be situated in the specified charge state Z0.

initialize_prescribed_fully_ionized(n=None, r=None, t=None, interpr=None, interpt=None)

Prescribe the ions to be fully ionized.

initialize_prescribed_neutral(n=None, r=None, t=None, interpr=None, interpt=None)

Prescribe the ions to be neutral.

initialize_source(n, t=None, Z0=0)

Initialize the ion source term associated with this species.

isHydrogen()

Returns True if this ion species is a hydrogen species.

isTritium()

Returns True if this ion species is a tritium species.

setChargedAdvection(mode, Ar=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Ar_0=None, Ar_f=0, Ar_o=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial advection for charged particles.

Parameters
  • mode – Type of transport to prescribe.

  • Ar – Advection coefficient to prescribe.

  • r – Radial grid on which Ar is given.

  • t – Time grid on which Ar is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Ar_0 – Initial value of advection coefficient when decaying exponentially.

  • Ar_f – Final value of advection coefficient when decaying exponentially.

  • Ar_o – Value of advection coefficient before exponential decay.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setChargedDiffusion(mode, Drr=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Drr_0=None, Drr_f=0, Drr_o=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial diffusion for charged particles.

Parameters
  • mode – Type of transport to prescribe.

  • Drr – Diffusion coefficient to prescribe.

  • r – Radial grid on which Drr is given.

  • t – Time grid on which Drr is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Drr_0 – Initial value of diffusion coefficient when decaying exponentially.

  • Drr_f – Final value of diffusion coefficient when decaying exponentially.

  • Drr_o – Value of diffusion coefficient before exponential decay.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setNeutralAdvection(mode, Ar=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Ar_0=None, Ar_f=0, Ar_o=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial advection for neutral particles.

Parameters
  • mode – Type of transport to prescribe.

  • Ar – Advection coefficient to prescribe.

  • r – Radial grid on which Ar is given.

  • t – Time grid on which Ar is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Ar_0 – Initial value of advection coefficient when decaying exponentially.

  • Ar_f – Final value of advection coefficient when decaying exponentially.

  • Ar_o – Value of advection coefficient before exponential decay.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setNeutralDiffusion(mode, Drr=None, r=None, t=None, t_transp_expdecay_all_cs=None, t_transp_start_expdecay_all_cs=0, Drr_0=None, Drr_f=0, Drr_o=0, r_expdecay_all_cs=None, t_expdecay_all_cs=None, interpr=None, interpt=None)

Set ion radial diffusion for neutral particles.

Parameters
  • mode – Type of transport to prescribe.

  • Drr – Diffusion coefficient to prescribe.

  • r – Radial grid on which Drr is given.

  • t – Time grid on which Drr is given.

  • t_transp_expdecay_all_cs – e-folding time of transport coefficient decay (exponential decay).

  • t_transp_start_expdecay_all_cs – Start time of exponential decay.

  • Drr_0 – Initial value of diffusion coefficient when decaying exponentially.

  • Drr_f – Final value of diffusion coefficient when decaying exponentially.

  • Drr_o – Value of diffusion coefficient before exponential decay.

  • r_expdecay_all_cs – Radial grid on which the coefficient should be defined.

  • t_expdecay_all_cs – Time grid on which the coefficient should be defined.

  • interpr – Radial grid onto which ion transport coefficients should be interpolated.

  • interpt – Time grid onto which ion transport coefficients should be interpolated.

setSPIMolarFraction(SPIMolarFraction)
setTemperature(T)

Sets the ion temperature from an input value T. For scalar T, sets a uniform radial profile, otherwise requires the T profile to be given on the r grid which is provided to the IonSpecies constructor.

verifySettings()

Verify that the settings of this ion species are correctly set.