DREAM includes support for introducing ions to the plasma through Shattered Pellet Injections (SPI). The user can initialise an arbitrary number of shards with given initial sizes, positions, velocities and compositions. This class deals with the settings for the initialization and evolution of the shard sizes, positions and velocities as well as the distribution of the ablated material around the shards. The compositions are set by adding an ion species for every species the pellet is composed of and connecting this ion species to the relevant shards, see the documentation for the Ions class.
Page overview
There are six types of modes with various options for the SPI modelling:
the ablation mode, which determines the model used for the evolution of the shard radii
the velocity mode, which determines the motion of the shards
the deposition mode, which determines the spread of the ablated material around the shards
the shift mode, which determines the mode used for the shift of the ablated material due to the plasmoid drift effect
the heat absorbtion mode, which determines the absorbtion of heat flowing into the neutral cloud surrounding the shards
the cloud radius mode, which determines the calculation of the radius of the neutral cloud. This radius is used to determine the cross section for the heat absorbtion (if included), and also determines the length scale of the deposition kernel (if the deposition kernel used has a finit width).
This section describes the various modes available in all these five categories.
The available ablation modes are
Name |
Description |
|---|---|
|
No ablation (default) |
|
Ablation according to the Neutral Gas Shielding (NGS) formula
presented by P. Parks at TSDW 2017, calculated using the fluid
density and temperature ( |
|
Similar to |
The available velocity modes are
Name |
Description |
|---|---|
|
The shards do not move (default) |
|
Prescribed (and constant) shard velocities |
The available deposition modes are
Name |
Description |
|---|---|
|
No deposition (default) |
|
Delta function source (averaged over the current time step and over the grid cells) |
|
Same as |
|
Gaussian deposition kernel in the radial coordinate, with 1/e-length scale equal to the radius of the neutral cloud. NOTE: This kernel is not averaged over the time step, so it should be used with caution if the shards travel a longer distance than the radius of the neutral cloud during a single time step! |
The available shift modes are
Name |
Description |
|---|---|
|
No shift (default) |
|
Calculate the shift using the analytical model derived by Vallhagen et al, JPP 2023 |
The available heat absorbtion modes are
Name |
Description |
|---|---|
|
No heat absorbtion, only the energy required to ionize the ablated material to the equilibrium distribution of ionisation states is subtracted from the background plasma |
|
All heat flowing through a disc of radius \(r_\mathrm{cld}\) is subtracted from the background plasma as a delta function source around each shard, where \(r_\mathrm{cld}\) is the radius of the neutral cloud (no additional sink for the ionization of the ablated material) |
|
Same as |
Note
If the ablated material is immediately deposited locally at the shard positions, they immediately carry the energy absorbed in the cloud back to the plasma, except for the energy needed for the ionization. In such a situation, HEAT_ABSORBTION_MODE_NEGLECT would be the most physically appropriate heat absorbtion mode.
The available cloud radius modes are
Name |
Description |
|---|---|
|
Cloud radius not used (default) |
|
Constant prescribed cloud radius |
|
Currently gives simply \(r_\mathrm{cld}=10r_\mathrm{p}\), which is a very rough approximation |
The most general way to add shards to the SPI is to directly provide a vector containing the shard sizes, initial positions and velocities to the setInitialData() method. The initial positions and velocities are given in cartesian coordinates, with the origin at the magnetic axis and the xy-plane coinciding with the poloidal cross section. The vectors specifying the initial positions and velocities should have the format \(\boldsymbol{x}_\mathrm{p}=(x_\mathrm{p,1},y_\mathrm{p,1},z_\mathrm{p,1},x_\mathrm{p2},y_\mathrm{p2},z_\mathrm{p2},...)\). The example below shows how to initialise one shard with radius 1 cm on the horisontal mid-plane 2.15 m from the magnetic axis, traveling at a speed of 200 m/s directly towards the magnetic axis.
ds = DREAMSettings()
...
ds.eqsys.spi.setInitialData(rp=0.01, xp=np.array([2.15,0,0]), vp=np.array([-200,0,0]))
There are, however, a number of helper-functions implemented to more easily set up injections with some standard distributions for the shard parameters. These functions are covered in the rest of this section.
To set nShard new shard positions to a singel scattering point (x0,y0,z0), the setShardPositionSinglePoint() can be used:
ds.eqsys.spi.setShardPositionSinglePoint(nShard=nShard, shatterPoint=np.array([x0,y0,z0]), add=True)
The last argument add determines wether a new set of shards should be added to the existing ones (add=True, default) or if the shard position vector should be cleared (add=False).
To select nShard new shard velocities uniformly within a magnitude range v0-DeltaV,v0+DeltaV and directions chosen uniformly over an nDim dimensional cone of opening angle alpha/2 and axis anti-parallell with the x-axis, the setShardVelocitiesUniform() can be used:
ds.eqsys.setShardVelocitiesUniform(nShard=nShard, abs_vp_mean=v0, abs_vp_diff=DeltaV, alpha_max=alpha, nDim=nDim)
If nDim=1, all shards simply move along the x-axis. If nDim=2, the direction is chosen along an arc spanning an angle alpha and if nDim=3 the direction is chosen over an ordinary 3-dimensional cone with opening angle alpha/2. Similarly to setShardPositionSinglePoint(), there is an additional argument add which is set to True by default.
In some cases, such as when making a staggered injection, it is practical to create some of the shards and add the corresponding ion species already from the beginning, but not set the shards into motion until a later restart. In this way, one does not have to re-set the ion species when making a restart, but can simply initialise the ion species from the previous output. This can be done by using the additional parameter shards in setShardVelocitiesUniform(), which is a slice specifying the indices for the shards whose velocities should be updated. If one, for example, wants to change the velocities of the nShard2 last shards added, this could be done as
ds.eqsys.setShardVelocitiesUniform(abs_vp_mean=v0, abs_vp_diff=DeltaV, alpha_max=alpha, nDim=nDim, shards=slice(-nShard2,None))
When shards is not None, the nShard parameter is automatically set to the number of shards specified by the shards-parameter, and add is set to False.
There are a number of helper methods implemented to select shard sizes from the Bessel-like statistical distribution found in P. Parks GA report:
where \(K_0\) is the zeroth modified Bessel function of the second kind. The inverse characteristic shard size \(k_\mathrm{p}\) is related to the particle content, composition and degree of shattering of the pellet according to
where \(n_\mathrm{p}\) is the solid particle density of the pellet, \(N_\mathrm{s}\) is the number of shards into which the pellet is shattered and \(N_\mathrm{inj}\) is the total number of particles contained in the pellet.
If the desired \(k_\mathrm{p}\) is already known, one can sample nShard shards from the above distribution as
rp=ds.eqsys.spi.sampleRpDistrParksStatistical(N=nShard, kp=kp)
If one instead wants to specify \(N_\mathrm{s}\), \(N_\mathrm{inj}\) and the pellet composition, one can instead use the setRpParksStatistical()-method. The following lines of code illustrate how to add a pellet containing \(10^{24}\) particles shattered into 1000 shards, consisting of 5% neon and 95% deuterium, with radii selected from the above distribution. The Ion species connected to this pellet are also added by the setRpParksStatistical()-method by passing a pointer to the ds.eqsys.n_i-object.
Ninj=1e24 # Total number of injected particles
nShard=1000 # Number of shards into which the pellet is shattered
Zs=[1,10] # List of charge numbers of the species the pellet is composed of
isotopes=[2,0] # List of isotopes, 0 meaning naturally occuring mix
molarFractions=[0.95,0.05] # List of molar fractions specifying the pellet composition
ionNames=['D_inj','Ne_inj'] # List of names of the ion species connected to this pellet
ds.eqsys.spi.setRpParksStatistical(nShard=nShard, Ninj=Ninj, Zs=Zs, isotopes=isotopes, molarFractions=molarFractions, ionNames=ionNames, n_i=ds.eqsys.n_i)
As earlier, an extra argument add=False resets the shard radii instead of adding new ones.
The code block below illustrates how to set up staggered Deuterium-Neon injections similar to those investigated in O. Vallhagens MSc thesis, using the wrapper function setParamsVallhagenMSc() to set the initial positions, velocities and radii of the shards on the same line.
radius=[0,2] # Span of the radial grid
radius_wall=2.15 # Location of the wall
...
# Settings for the first deuterium SPI
nShardD=1742 # Number of shards
NinjD=2e24 # Number of atoms
alpha_maxD=0.17 # Divergence angle
abs_vp_meanD=800 # Mean shard speed
abs_vp_diffD=0.2*abs_vp_meanD # Width of the uniform shard speed distribution
ds.eqsys.spi.setParamsVallhagenMSc(nShard=nShardD, Ninj=NinjD, Zs=[1], isotopes=[2], molarFractions=[1], ionNames=['D_inj'], n_i=ds.eqsys.n_i, shatterPoint=np.array([0,0,radius_wall]), abs_vp_mean=abs_vp_meanD, abs_vp_diff=abs_vp_diffD, alpha_max=alpha_maxD)
# Settings for the second neon SPI
nShardNe=50 # Number of shards
NinjNe=1e23 # Number of atoms
alpha_maxNe=0.17 # Divergence angle
abs_vp_meanNe=200 # Mean shard speed
abs_vp_diffNe=0.2*abs_vp_meanNe # Width of the uniform shard speed distribution
# Initialise neon shards with zero velocity, and set the velocity before the restart when the neon injection should actually start
ds.eqsys.spi.setParamsVallhagenMSc(nShard=nShardNe, Ninj=NinjNe, Zs=[10], isotopes=[0], molarFractions=[1], ionNames=['Ne_inj'], n_i=ds.eqsys.n_i, shatterPoint=np.array([0,0,radius_wall]), abs_vp_mean=0, abs_vp_diff=0, alpha_max=alpha_maxNe)
...
runiface(ds, 'output_D_inj.h5')
...
ds2=DREAMSettings(ds)
ds2.fromOutput('output_D_inj.h5', ignore=['v_p','x_p'])
...
# Set the shards of the second injection into motion and advance them until
# the fastest shards reach the plasma edge
do=DREAMOutput('output_D_inj.h5')
ds2.eqsys.spi.vp=do.eqsys.v_p.data[-1,:].flatten()
ds2.eqsys.spi.xp=do.eqsys.x_p.data[-1,:].flatten()
ds2.eqsys.spi.setShardVelocitiesUniform(abs_vp_mean=abs_vp_meanNe,abs_vp_diff=abs_vp_diffNe,alpha_max=alpha_maxNe,shards=slice(-nShardNe,None))
t_edge=(radius_wall-radius[1])/np.max(-ds2.eqsys.spi.vp[-3*nShardNe::3])
ds2.eqsys.spi.xp[-3*nShardNe:]=ds2.eqsys.spi.xp[-3*nShardNe:]+ds2.eqsys.spi.vp[-3*nShardNe:]*t_edge
runiface(ds2, 'output_Ne_inj.h5')
The effect of plasmoid drifts can be accounted for using the analytical model given by equation (A4) in Vallhagen et al, JPP 2023. To use this model, the following model parameters have to be speified: representative temperature \(T\) during the drift duration; the initial temperature \(T_0\) close to the pellet (before the cloud has drifted away from the pellet); the characteristic half-thickness of the drifting cloud \(\Delta y\) (should be similar to the radius of the neutral cloud around the pellet); the major radius \(R_\mathrm{m}\) at the magnetic axis (only used if the major radius is otherwise infinite in the simulation); the average charge states \(Z_\mathrm{avg}^{(i)}\) inside the drifting cloud for all ion species include in the pellets. Note that the average charge states can not be calculated using the ADAS rates because the conditions in the drifting cloud, especially the density and optical thickness, are very different from the validity range and assumptions in ADAS, and we therefore take user-given estimates for them. These estimates should be given as a look-up table in the form of the lists ZavgArray, Zs and isotopes, where ZavgArray provides the average charge states to be used for ion species with atomic numbers Zs and isotopes isotopes.
Below is an example showing how to include the analytical drift model for a mixed D+Ne SPI with typical settings. The values of the model parameters are motivated by e.g. numerical studies by Matsuyama and experimental studies by Müller et al. In particular, these studies indicate that the representative temperature \(T\) is significantly different for purely hydrogenic pellets and impurity doped pellets, as the doped pellets will have much stronger radiation losses. Therefore, different temperatures can be set
for different shards as long as they are in the correct order. The following example includes nShardD
deuterium shards (atomic number 1, isotope 2) with \(T=30\,\rm eV\) and nShardNe number of neon doped shards (atomic number 10, isotope 0 i.e. naturally occuring mix) with \(T=5\,\rm eV\).
ds.eqsys.spi.setShiftParamsAnalytical(shift = SHIFT_MODE_ANALYTICAL, T = np.concatenate([30*np.ones(nShardD), 5*np.ones(nShardNe)]), T0 = 2, delta_y = 0.0125, Rm = R0 , ZavgArray =[1 , 2], Zs, [1,10], isotopes = [2,0]):
DREAM.Settings.Equations.SPI.SPI(settings, rp=None, vp=None, xp=None, t_delay=None, VpVolNormFactor=1, rclPrescribedConstant=0.01, velocity=1, ablation=1, deposition=1, heatAbsorbtion=1, cloudRadiusMode=1, magneticFieldDependenceMode=1, abl_ioniz=1, shiftMode=1, TDrift=None, T0Drift=None, DeltaYDrift=None, RmDrift=None, ZavgDriftArray=None, ZsDrift=None, isotopesDrift=None)¶Bases: DREAM.Settings.Equations.UnknownQuantity.UnknownQuantity
__init__(settings, rp=None, vp=None, xp=None, t_delay=None, VpVolNormFactor=1, rclPrescribedConstant=0.01, velocity=1, ablation=1, deposition=1, heatAbsorbtion=1, cloudRadiusMode=1, magneticFieldDependenceMode=1, abl_ioniz=1, shiftMode=1, TDrift=None, T0Drift=None, DeltaYDrift=None, RmDrift=None, ZavgDriftArray=None, ZsDrift=None, isotopesDrift=None)¶Constructor.
settings (DREAMSettings) – Parent DREAMSettings object.
rp (numpy.ndarray) – Initial shard radii.
vp (numpy.ndarray) – Initial shard velocities (cartesian coordinates)
xp (numpy.ndarray) – Initial shard positions (cartesian coordinates)
VpVolNormFactor (float) – Factor used to renormalize the value of VpVol used for calculating the voluma of the flux tubes (eg major radius in case of cylindrical geometry)
rclPrescribedConstant (float) – Constant, prescribed radius of the neutral cloud surrounding each pellet shard (only applicable if cloudRadiusMode=CLOUD_RADIUS_MODE_PRESCRIBED_CONSTANT)
velocity (int) – Model used for the shard velocities
ablation (int) – Model used for shard ablation
deposition (int) – Model used for the deposition of the ablated material
heatAbsobtion (int) – Model used for absorbtion of heat flowing into the neutral clouds
cloudRadiusMode (int) – Mode used for calculating the radius of the neutral clouds
magneticFieldDependenceMode (int) – Mode used for calculating the magnetic field dependence of the albation
shift (int) – Model used for determining the cloud drift
fromdict(data)¶Set all options from a dictionary.
resetTimeDelay()¶rpDistrParksStatistical(rp, kp)¶Evaluates the shard size distribution function referred to as the ‘statistical model’ in P. Parks 2016 GA report (DOI:10.2172/1344852)
sampleRpDistrParksStatistical(N, kp)¶Samples N shard radii according to the distribution function given by rpDistrParksStatistical()
setAblIoniz(abl_ioniz)¶Specifies which model to use for calculating the charge state distribution with which the recently ablated material is deposited
setAblation(ablation)¶Specifies which model to use for calculating the ablation rate.
setCloudRadiusMode(cloudRadiusMode)¶Specifies which model to use for calculating the radius of the the neutral pellet cloud
setDeposition(deposition)¶Specifies which model to use for calculating the deposition of ablated material.
setHeatAbsorbtion(heatAbsorbtion)¶Specifies which model to use for calculating the heat absorbtion in the neutral pellet cloud
setInitialData(rp=None, vp=None, xp=None, t_delay=None, nbrShiftGridCell=None, TDrift=None)¶setMagneticFieldDependenceMode(magneticFieldDependenceMode)¶Specifies which model to use for calculating the magnetic field dependence of the ablation
setParamsVallhagenMSc(nShard, Ninj, Zs, isotopes, molarFractions, ionNames, shatterPoint, abs_vp_mean, abs_vp_diff, alpha_max, t_delay=0, nDim=2, add=True, opacity_modes=None, nbrShiftGridCell=0, TDrift=None, **kwargs)¶Wrapper for setRpParksStatistical(), setShardPositionSinglePoint() and setShardVelocitiesUniform(), which combined are used to set up an SPI-scenario similar to those in Oskar Vallhagens MSc thesis (available at https://hdl.handle.net/20.500.12380/302296).
setRclPrescribedConstant(rclPrescribedConstant)¶setRpParksStatistical(nShard, Ninj, Zs, isotopes, molarFractions, ionNames, opacity_modes=None, add=True, n=1.0, charged_advection_modes=None, charged_prescribed_advections=None, rChargedPrescribedAdvections=None, tChargedPrescribedAdvections=None, neutral_advection_modes=None, neutral_prescribed_advections=None, rNeutralPrescribedAdvections=None, tNeutralPrescribedAdvections=None, charged_diffusion_modes=None, charged_prescribed_diffusions=None, rChargedPrescribedDiffusions=None, tChargedPrescribedDiffusions=None, neutral_diffusion_modes=None, neutral_prescribed_diffusions=None, rNeutralPrescribedDiffusions=None, tNeutralPrescribedDiffusions=None, **kwargs)¶sets (or adds) nShard shards with radii distributed accordin to rpDistrParksStatistical(), with the characteristic inverse shard size kp calculated from the given pellet and shattering parameters. Also updates the ion settings with the appropriate molar fractions contributing to each ion species
nShard (int) – Number of shards into which the pellet is shattered
Ninj (float) – Numbr of particles contained in the pellet
Zs (list) – List of charge numbers for every ion species the pellet consists of
isotopes (list) – List of isotopes for every ion species the pellet consists of
molarFractions (numpy.ndarray) – Molar fraction with which each ion species contribute
ionNames (list) – List of names for the ion species to be added and connected to the ablation of this pellet
opacity_modes (list) – List of opacity modes for every ion species the pellet consists of. If ‘None’, this argument is omitted when adding the ion species, so that the default settings (transparent) is used
add (bool) – If ‘True’, add the new pellet shards to the existing ones, otherwise existing shards are cleared
the inverse characteristic shard size kp
setShardPositionSinglePoint(nShard, shatterPoint, add=True)¶Sets self.xp to a vector of the (x,y,z)-coordinates of nShard initial pellet shard positions starting from the single point shatterPoint
nShard (int) – Number of shards
shatterPoint (numpy.ndarray) – (x,y,z)-coordinates for the starting point of the shards to be set
add (bool) – If ‘True’, add the new pellet shard positions to the existing ones, otherwise existing shards are cleared
setShardVelocitiesUniform(nShard, abs_vp_mean, abs_vp_diff, alpha_max, t_delay=0, nDim=2, add=True, shards=None)¶Sets self.vp to a vector storing the (x,y,z)-components of nShard shard velosities, assuming a uniform velocity distribution over a nDim-dimensional cone whose axis is anti-parallell to the x-axis. TODO: implement support for an arbitrary axis?
nShard (int) – Number of shards
abs_vp_mean (float) – Mean of the magnitude of the shard velocities
abs_vp_diff (float) – width of the uniform distribution of the magnitude of the shard velocities
alpha_max (float) – Span of divergence angle (ie twice the opening angle of the cone)
nDim (int) – number of dimensions into which the shards should be spread
add (bool) – If ‘True’, add the new pellet shard velocities to the existing ones, otherwise existing shards are cleared
shards (slice) – indices of existing shards whose velocities should be updated. If not ‘None’, add is set to ‘False’ and nShard is set to the number of indices to be updated
setShift(shift)¶Specifies which model to use for calculating the ablation cloud drift
setShiftParamsAnalytical(shift=3, TDrift=None, T0Drift=0, DeltaYDrift=0, RmDrift=- 1, ZavgDriftArray=[0.0], ZsDrift=[0], isotopesDrift=[0], add=True)¶Specifies model parameters to be used for calculating the shift. Apart from the shift mode-argument, the parameters below apply to SHIFT_MODE_ANALYTICAL
shift (int) – Model used for determining the cloud drift
T0Drift (float) – cloud temperature close to the pellet (before the cloud has drifted away from the pellet)
TDrift (numpy.ndarray) – representative cloud temperature during the drift motion for each shard
DeltaYDrift (float) – characteristic half-thickness of the drifting cloud (should be similar to the radius of the neutral cloud around the pellet)
RmDrift (float) – major radius of the magnetic axis, only used if the major radius is otherwise infinite in the simulation
ZavgDriftArray (list) – average charge states inside the drifting cloud of all drifting ion species. These can not be calculated using the ADAS rates because the conditions in the drifting cloud, especially the density and optical thickness, are very different from the validity range and assumptions in ADAS, and we therefore take user-given estimates for them. Note that his list does NOT neccessarily have the same shape as the list of atomic numbers and isotopes included in the simulation, but instead the ZavgDriftArray-list and the ZsDrift and isotopesDrift-lists below will instead be used to look up the average charge state inside the drifting cloud for all the simulated ion species included in the pellet.
ZsDrift (list) – atomic numbers of all the drifting ion species, corresponding to the average charge states listed in the ZavgDriftArray-list above
isotopesDrift (list) – isotopes of all the drifting ion species, corresponding to the average charge states listed in the ZavgDriftArray-list above
setShiftParamsPrescribed(shift=2, nbrShiftGridCell=None, add=True)¶setVelocity(velocity)¶Specifies mode to calculate shard velocities.
setVpVolNormFactor(VpVolNormFactor)¶shiftTimeDelay(tShift)¶todict()¶Returns a Python dictionary containing all settings of this SPI object.
verifySettings()¶Verify that the settings of this unknown are correctly set.
verifySettingsPrescribedInitialData()¶