# SPI¶

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.

## SPI Model Settings¶

There are five 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 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.

### Ablation Modes¶

The available ablation modes are

Name

Description

ABLATION_MODE_NEGLECT

No ablation (default)

ABLATION_MODE_FLUID_NGS

Ablation according to the Neutral Gas Shielding (NGS) formula presented by P. Parks at TSDW 2017, calculated using the fluid density and temperature (T_cold and n_cold)

ABLATION_MODE_KINETIC_NGS

Similar to ABLATION_MODE_FLUID_NGS, but expressed in terms of the mean energy and isotropic heat flux of the distribution function. NOTE: the factor converting from the total mean energy and isotropic heat flux to the corresponding values in the direction towards the shards are taken to be the same as for a Maxwellian, so this mode is somewhat approximate.

### Velocity Modes¶

The available velocity modes are

Name

Description

VELOCITY_MODE NONE

The shards do not move (default)

VELOCITY_MODE_PRESCRIBED

Prescribed (and constant) shard velocities

### Deposition Modes¶

The available deposition modes are

Name

Description

DEPOSITION_MODE_NEGLECT

No deposition (default)

DEPOSITION_MODE_LOCAL

Delta function source (averaged over the current time step and over the grid cells)

DEPOSITION_MODE_LOCAL_LAST_FLUX_TUBE

Same as DEPOSITION_MODE_LOCAL but shifted one flux tube behind the shards. This can be used to mimic the effect of a finite equilibration time in the sense that the ablated material does not immediately affect the background plasma at its own location, which weakens the feedback effect between the ablation and the resulting cooling

DEPOSITION_MODE_LOCAL_GAUSSIAN

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!

### Heat Absorbtion Modes¶

The available heat absorbtion modes are

Name

Description

HEAT_ABSORBTION_MODE_NEGLECT

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

HEAT_ABSORBTION_MODE_LOCAL_FLUID_NGS

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)

HEAT_ABSORBTION_MODE_LOCAL_FLUID_NGS_GAUSSIAN

Same as HEAT_ABSORBTION_MODE_LOCAL_FLUID_NGS but with a Gaussian-shaped sink, with 1/e-length scale equal to $$r_\mathrm{cld}$$. 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!

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_MODE_NEGLECT

CLOUD_RADIUS_MODE_PRESCRIBED_CONSTANT

CLOUD_RADIUS_MODE_SELFCONSISTENT

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.

### Shard positions¶

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

### Shard velocities¶

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:

$P(r_\mathrm{p}) = k_\mathrm{p}^2 r_\mathrm{p} K_0(k_\mathrm{p} r_\mathrm{p}),$

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

$k_\mathrm{p}=\left(\frac{6\pi^2n_\mathrm{p}N_\mathrm{s}}{N_\mathrm{inj}}\right)^{1/3},$

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.

### Example: staggered injection¶

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))

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')


## Class documentation¶

class DREAM.Settings.Equations.SPI.SPI(settings, rp=None, vp=None, xp=None, VpVolNormFactor=1, rclPrescribedConstant=0.01, velocity=1, ablation=1, deposition=1, heatAbsorbtion=1, cloudRadiusMode=1)

Bases: DREAM.Settings.Equations.UnknownQuantity.UnknownQuantity

__init__(settings, rp=None, vp=None, xp=None, VpVolNormFactor=1, rclPrescribedConstant=0.01, velocity=1, ablation=1, deposition=1, heatAbsorbtion=1, cloudRadiusMode=1)

Constructor.

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

fromdict(data)

Set all options from a dictionary.

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()

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)
setParamsVallhagenMSc(nShard, Ninj, Zs, isotopes, molarFractions, ionNames, shatterPoint, abs_vp_mean, abs_vp_diff, alpha_max, nDim=2, add=True, opacity_modes=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, **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

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

Returns

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

Parameters
• 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, 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?

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

setVelocity(velocity)

Specifies mode to calculate shard velocities.

setVpVolNormFactor(VpVolNormFactor)
todict()

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

verifySettings()

Verify that the settings of this unknown are correctly set.

verifySettingsPrescribedInitialData()