Merzbild.jl internal API reference

map_cont_index(particle_indexer, i)

Maps a continuous index in the range [0, n_local-1] to an index in the particle array given a ParticleIndexer instance describing how particle indices are split across 2 groups.

Positional arguments

• particle_indexer: the ParticleIndexer instance
• i: the index to map

map_cont_index(pia, cell, species, i)

Maps a continuous index in the range [0, n_local-1] to an index in the particle array for particles of a specific species in a specific cell.

Positional arguments

• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particles are located
• species: the index of the particles' species
• i: the index to map

update_particle_indexer_new_particle!(pia, cell, species)

Update a ParticleIndexerArray instance when a particle of a given species in a given cell is created. This places the particle index in the 2-nd group of particle indices in the ParticleIndexer instance.

Positional arguments

• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is created
• species: the index of the species of which the particle is created

Base.getindex(pv::ParticleVector{D}, i) where D

Returns the underlying particle in a ParticleVector instance with index i.

Is usually called as ParticleVector[i].

Positional arguments

• pv: ParticleVector instance
• i: the index of the particle to be selected

Base.setindex!(pv::ParticleVector{D}, p::Particle{D}, i::Integer) where D

Set the underlying particle in a ParticleVector instance with index i to a new particle.

Is usually called as ParticleVector[i] = p.

Positional arguments

• pv: ParticleVector instance
• p: the Particle instance to write
• i: the index of the particle to be written to

Base.length(pv::ParticleVector{D}) where D

Returns the length of a ParticleVector instance.

Is usually called as length(ParticleVector).

Positional arguments

• pv: ParticleVector instance

Base.resize!(pv::ParticleVector{D}, n::Integer) where D

Resize a ParticleVector instance, taking care of the indices, buffer, and creating placeholder new particles with weight 0, velocity 0, and position 0.

Positional arguments

• pv: ParticleVector instance
• n: the new length of the ParticleVector instance (i.e. the length of all the vector fields of the instance)

add_particle!(pv::ParticleVector{D1}, position, w, v::SVector{3,Float64}, x::SVector{D2,Float64}) where {D1,D2}

Create a new particle in a ParticleVector instance at position position. The ParticleIndexer/ParticleIndexerArray instances should be updated independently. See update_particle_buffer_new_particle! for more information regarding how the buffer of the ParticleVector is updated. This should not be used to update an existing particle. Particles at positions before position should exist in the ParticleVector array. Any components of x that exceed the x dimension of pv are discarded.

Positional arguments

• pv: ParticleVector instance
• position: the position in the index vector to which to write the index of the new particle
• w: the computational weight of the particle to create
• v: the velocity of the particle to create

update_particle_buffer_new_particle!(pv::ParticleVector{D}, position) where D

Update the buffer in a ParticleVector instance when a new particle is created. This writes the index of the new particle (the last index stored in the active part of the buffer) to the index vector at position position, and reduces the length of the active part of the buffer by 1.

Positional arguments

• pv: ParticleVector instance
• position: the position in the index vector to which to write the index of the new particle

update_particle_buffer_new_particle!(pv::ParticleVector{D}, pia, species) where D

Update the buffer in a ParticleVector instance when a new particle is created at the end of the particle array, and reduces the length of the active part of the buffer by 1. This assumes that the pia structure has already an updated particle count (that accounts for the newly created particle), as the index of the new particle taken from the buffer is written to pv.index.[pia.index_last[species]].

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• species: the index of the species of which a new particle is created

update_buffer_index_new_particle!(pv::ParticleVector{D}, pia, cell, species) where D

Update a ParticleIndexerArray and the buffer in a ParticleVector instance when a particle of a given species in a given cell is created. The particle index is added to the second group of particles pointed to by the ParticleIndexer. See the documentation of update_particle_indexer_new_particle! and update_particle_buffer_new_particle!

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is created
• species: the index of the species of which the particle is created

delete_particle!(pv::ParticleVector{D}, pia, cell, species, i) where D

Delete particle with index i of species species in cell cell and update the particle indexers and buffers accordingly. This changes the ordering of the non-deleted particles in the cell.

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is deleted
• species: the index of the species of which the particle is deleted
• i: the index of the particle to delete

delete_particle_end!(pv::ParticleVector{D}, pia, cell, species) where D

Delete the last particle of species species in cell cell: if particles are present in the 2nd group of the indices stored in the ParticleIndexer instance, it will delete the last particle in that group; otherwise it will delete the last particle in the 1st group of particles pointed to by the ParticleIndexer instance. If no particles are present in the cell, the function does nothing. This does not set the value of the contiguous field of pia to false even if pia becomes discontinuous; this has to be done outside of this function.

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is deleted
• species: the index of the species of which the particle is deleted

delete_particle_end_group1!(pv::ParticleVector{D}, pia, cell, species) where D

Delete particle with index pia.indexer[cell, species].end1of speciesspeciesin cellcelland update the particle indexers and buffers accordingly (i.e. delete the last particle in the 1st group of particles of a given species in a given cell). This also sets the weight of the deleted particle to 0. If no particles are present in the 1st group of particles, the function does nothing. This does not set the value of thecontiguousfield ofpiatofalseeven ifpia` becomes discontinuous; this has to be done outside of this function.

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is deleted
• species: the index of the species of which the particle is deleted

delete_particle_end_group2!(pv::ParticleVector{D}, pia, cell, species) where D

Delete particle with index pia.indexer[cell, species].end2of speciesspeciesin cellcelland update the particle indexers and buffers accordingly (i.e. delete the last particle in the 2nd group of the particles of a given species in a given cell). This also sets the weight of the deleted particle to 0. If no particles are present in the 2nd group of particles, the function does nothing. This does not set the value of thecontiguousfield ofpiatofalseeven ifpia` becomes discontinuous; this has to be done outside of this function.

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is deleted
• species: the index of the species of which the particle is deleted

delete_batch_end!(pv::ParticleVector{D}, pia, cell, species, n) where D

Delete the last n particles of species species in cell cell: first any particles in the 2nd group of the indices stored in the ParticleIndexer instance will be deleted; followed by particles in the 1st group of particles pointed to by the ParticleIndexer instance (starting at the end of each group) If no particles are present in the cell, the function does nothing. This does not set the value of the contiguous field of pia to false even if pia becomes discontinuous; this has to be done outside of this function. This updates the indexing, the buffer, n_total and index_last. It is equivalent to calling delete_particle_end! n times, but does the bookkeeping once for the whole batch and resolves index_last with at most one scan over the cells. Note: this assumes that n <= n_group1 + n_group2.

Positional arguments

• pv: ParticleVector instance
• pia: the ParticleIndexerArray instance
• cell: the index of the cell in which the particle is deleted
• species: the index of the species of which the particle is deleted

delete_batch_end_group1!(pv::ParticleVector{D}, indexer, n) where D

Delete n particles from the end of the 1st group of particles of a given species in a given cell indexed by indexer and update the particle indexer and buffer accordingly. This also sets the weight of the deleted particles to 0. If n is 0, the function does nothing. This does not set the value of the contiguous field of pia to false even if pia becomes discontinuous; this has to be done outside of this function. Note: 1) This does not check that n is not larger than the number of particles in the 1st group of particles.

1. This also does not update index_last or n_total; use delete_batch_end! to have those updated.

Positional arguments

• pv: ParticleVector instance
• indexer: the ParticleIndexer instance
• n: the number of particles to delete

delete_batch_end_group2!(pv::ParticleVector{D}, indexer, n) where D

Delete n particles from the end of the 2nd group of particles of a given species in a given cell indexed by indexer and update the particle indexer and buffer accordingly. This also sets the weight of the deleted particles to 0. If n is 0, the function does nothing. This does not set the value of the contiguous field of pia to false even if pia becomes discontinuous; this has to be done outside of this function. Note: 1) This does not check that n is not larger than the number of particles in the 2nd group of particles of species in the cell. 2) This also does not update index_last or n_total; use delete_batch_end! to have those updated.

Positional arguments

• pv: ParticleVector instance
• indexer: the ParticleIndexer instance
• n: the number of particles to delete

find_index_last_after_group1_delete!(pia, cell, species)

Find and set the new value of pia.index_last[species] after the 1st group of particles of species species in cell cell has been emptied whilst holding the last index. Only the 1st groups of the cells before cell are searched, since a 1st group holding the last index implies that no 2nd group is populated; failing that, index_last is set to 0.

Positional arguments

• pia: the ParticleIndexerArray instance
• cell: the index of the cell the particles were deleted from
• species: the index of the species the particles were deleted from

find_index_last_after_group2_delete!(pia, cell, species)

Find and set the new value of pia.index_last[species] after the 2nd group of particles of species species in cell cell has been emptied whilst holding the last index. Cells before cell are searched for a non-empty 2nd group; failing that, all cells are searched for a non-empty 1st group; failing that, index_last is set to 0.

Positional arguments

• pia: the ParticleIndexerArray instance
• cell: the index of the cell the particles were deleted from
• species: the index of the species the particles were deleted from

compute_mu_ref(m_AB, vhs_o, vhs_Tref, vhs_d)

Compute reference viscosity for the VHS model.

Positional arguments

• m_AB: average mass of colliding particles
• vhs_o: exponent for the VHS potential
• vhs_Tref: reference temperature for the VHS potential
• vhs_d: diameter for the VHS potential

Returns

• Reference viscosity

compute_vss_mu_ref_factor(vss_alpha)

Compute the factor $(\alpha + 1)(\alpha + 2) / (6\alpha)$ by which the VHS reference viscosity has to be multiplied to obtain the VSS reference viscosity for the same total cross-section. It is the inverse of the ratio of the VSS and VHS viscosity cross-sections $\sigma_\mu = \sigma_T (1 - \langle\cos^2\chi\rangle)$, and is equal to 1 for $\alpha = 1$ (isotropic scattering).

Positional arguments

• vss_alpha: VSS exponent (alpha)

Returns

• the reference viscosity correction factor

References

• G.A. Bird, Eq. (4.63), Molecular gas dynamics and the direct simulation of gas flows, Clarendon Press, Oxford, 1994.

interaction_from_toml(m1, m2, interaction_toml)

Construct an Interaction instance from the masses of the two species and the dictionary of the interaction parameters of the species pair, as read from an interaction data TOML file.

The elastic scattering model is given by the optional model key (see Merzbild.parse_scattering_model); if it is not present, the VHS model is used. Both models require the vhs_d, vhs_o, and vhs_Tref keys, with the VSS model additionally requiring the vss_alpha key.

Positional arguments

• m1: molecular mass of the first species
• m2: molecular mass of the second species
• interaction_toml: dictionary of the interaction parameters of the species pair

Returns

• the constructed Interaction instance

Throws

KeyError if a parameter required by the scattering model is not found in the dictionary, ArgumentError if the scattering model is not recognized.

UnitDVGrid

Stores information about a uniform discrete 3-dimensional velocity grid with extent $[-1.0,1.0]\times[-1.0,1.0]\times[-1.0,1.0]$.

Fields

• nx: number of cells in x direction
• ny: number of cells in y direction
• nz: number of cells in z direction
• dx: grid spacing in x direction
• dy: grid spacing in y direction
• dz: grid spacing in z direction
• vx_grid: Vector of grid nodes in x direction
• vy_grid: Vector of grid nodes in y direction
• vz_grid: Vector of grid nodes in z direction

DVGrid

Stores information about a uniform discrete 3-dimensional symmetric velocity grid with extent $[-v_{x,max},v_{x,max}]\times[-v_{y,max},v_{y,max}]\times[-v_{z,max},v_{z,max}]$.

Fields

• base_grid: the underlying unit (non-scaled) UnitDVGrid uniform grid
• vx_max: extent of the grid in the x direction
• vy_max: extent of the grid in the y direction
• vz_max: extent of the grid in the z direction
• dx: grid spacing in x direction
• dy: grid spacing in y direction
• dz: grid spacing in z direction
• vx_grid: Vector of grid nodes in x direction
• vy_grid: Vector of grid nodes in y direction
• vz_grid: Vector of grid nodes in z direction

VDF

Stores the values of a function evaluated on 3-dimensional velocity grid.

Fields

• nx: the number of grid elements in the x direction
• ny: the number of grid elements in the y direction
• nz: the number of grid elements in the z direction
• w: the 3-dimensional array of values

sample_bkw!(rng, particles::ParticleVector{D}, nparticles, offset, m, T, v0) where D

Sample particle velocities from the BKW distribution with temperature T for a species with mass m at t=0 and add a velocity offset. Note: This does not update the particle weights, positions, or any indexing structures.

Positional arguments

• rng: the random number generator
• particles: the Vector-like structure holding the particles
• nparticles: the number of particles to sample
• offset: offset the starting position for writing the sampled particles to the particles array
• m: species' mass
• T: temperature
• v0: the 3-dimensional velocity to add to the sampled velocities

sample_bkw!(rng, particles::ParticleVector{D}, nparticles, m, T, v0) where D

Sample particles' velocities from the BKW distribution with temperature T for a species with mass m at t=0 and add a velocity offset. This does not update the particle weights, positions, or any indexing structures. The sampled particles are written to the start of the particles array.

Positional arguments

• rng: the random number generator
• particles: the Vector-like structure holding the particles
• nparticles: the number of particles to sample
• m: species' mass
• T: temperature
• v0: the 3-dimensional velocity to add to the sampled velocities

evaluate_distribution_on_grid!(vdf, distribution_function, grid, w_total, cutoff_v; normalize=true)

Evaluate a distribution on a discrete velocity grid, considering only points inside a sphere of a given radius (the value of the VDF at points outside of the sphere will be 0.0). The values can be re-normalized so that the distribution has the prescribed computational weight/density.

Positional arguments

• vdf: The VDF instance where the evaluated values will be stored
• distribution_function: the distribution function to be evaluated which takes the x, y, and z velocities as parameters
• grid: the DVGrid on which the distribution function is evaluated
• w_total: the total weight used in the re-normalization
• cutoff_v: the radius of the sphere used to cut off high velocities: on grid points outside the sphere the VDF will be 0

Keyword arguments

• normalize: if true, the resulting values of the VDF will be renormalized so that their sum is equal to w_total

sample_maxwellian!(rng, particles::ParticleVector{D}, nparticles, offset, m, T, v0)

Sample nparticles particles from a Maxwellian with temperature T for a species with mass m and add a velocity offset. Note: This does not update the particle weights, positions, or any indexing structures.

Positional arguments

• rng: the random number generator
• particles: the Vector-like structure holding the particles
• nparticles: the number of particles to sample
• offset: offset the starting position for writing the sampled particles to the particles array
• m: species' mass
• T: temperature
• v0: the 3-dimensional velocity to add to the sampled velocities

norm3(v)

Compute the Euclidean norm of a 3-component vector.

Bit-for-bit identical to LinearAlgebra.norm(v), but cheaper in hot loops. It must stay outside any @muladd block: contracting the sum of squares into muladd changes the result by up to 0.3 ulp relative to norm, which is enough to break the bit-exact reference-solution tests.

Positional arguments

• v: the 3-component vector

Returns

The Euclidean norm of v.

compute_n_coll_single_species(rng, collision_factors, np, Δt, V)

Compute the non-integer number of collisions between particles of same species.

Positional arguments

• rng: the random number generator
• collision_factors: the CollisionFactors holding the estimate of $(\sigma g w)_{max}$ for the species in question in the cell
• np: number of particles in the cell
• Δt: timestep
• V: cell volume

Returns

The non-integer number of collisions.

compute_n_coll_single_species(rng, collision_factors_swpm, np, w_max, G, Δt, V)

Compute the non-integer number of collisions between particles of same species for use in the SWPM method.

Positional arguments

• rng: the random number generator
• collision_factors: the CollisionFactorsSWPM instance holding the estimate of $(\sigma g)_{max}$ for the species in question in the cell
• np: number of particles in the cell
• w_max: the maximum computational weight of the particles in the cell
• G: non-negative value defining the weight transfer function
• Δt: timestep
• V: cell volume

Returns

The non-integer number of collisions.

compute_n_coll_two_species(rng, collision_factors, np1, np2, Δt, V)

Compute number of collisions between particles of different species

Positional arguments

• rng: the random number generator
• collision_factors: the CollisionFactors holding the estimate of $(\sigma g w)_{max}$ for the species in question in the cell
• np1: number of particles of the first species in the cell
• np2: number of particles of the second species in the cell
• Δt: timestep
• V: cell volume

Returns

The non-integer number of collisions.

collide_2particles!(rng, model, collision_data, collision_factors, interaction, pa_i::Particle{D}, pa_k::Particle{D},
                    particles_1::ParticleVector{D}, particles_2::ParticleVector{D}, pia, cell, species1, species2; dw_tol=1e-16) where D

Collide two particles elastically using the elastic scattering model model. Particles can be of same or different species. If particles' weights differ by less than dw_tol, an equal-weight collision is performed and no particles are split.

Positional arguments

• rng: the random number generator
• model: the AbstractScatteringModel singleton for the scattering model
• collision_data: CollisionData instance used for storing collisional quantities
• collision_factors: the CollisionFactors holding the estimate of $(\sigma g w)_{max}$ for the species in question in the cell
• interaction: the Interaction instance for the colliding species
• particles_1: ParticleVector of the particles of the first species being collided
• particles_2: ParticleVector of the particles of the second species being collided
• pia: the ParticleIndexerArray
• cell: the index of the cell in which collisions are performed
• species1: the index of the first species for which collisions are performed
• species2: the index of the second species for which collisions are performed

Keyword arguments

• dw_tol: if weights of particles differ by less than this amount, an equal-weight collision is assumed

and no particle splitting is performed

collide_2particles_equal_weight!(rng, model, collision_data, collision_factors, interaction, pa_i::Particle{D}, pa_k::Particle{D}) where D

Collide two particles elastically using the elastic scattering model model, assuming equal weights - no particle splitting is performed even if weights are unequal. Particles can be of same or different species.

Positional arguments

• rng: the random number generator
• model: the AbstractScatteringModel singleton for the scattering model
• collision_data: CollisionData instance used for storing collisional quantities
• collision_factors: the CollisionFactors holding the estimate of $(\sigma g w)_{max}$ for the species in question in the cell
• interaction: the Interaction instance for the colliding species
• pa_i: the first particle being collided
• pa_k: the second particle being collided

compute_vhs_factor(vhs_Tref, vhs_d, vhs_o, m_r)

Compute interaction-specific factor $\pi D_{VHS}^2 (2 T_{ref,VHS}/m_r)^{(\omega_{VHS} - 0.5)} \frac{1}{\Gamma(2.5 - \omega_{VHS})}$.

Positional arguments

• vhs_Tref: reference temperature for the VHS potential
• vhs_d: diameter for the VHS potential
• vhs_o: exponent for the VHS potential
• m_r: collision-reduced mass

Returns

• computed VHS cross-section factor

compute_com!(collision_data::CollisionData, interaction::Interaction, p1, p2)

Compute center of mass velocity of two particles.

Positional arguments

• collision_data: the CollisionData instance where the computed velocity will be stored
• interaction: the Interaction instance for the species of the colliding particles
• p1: the first particle
• p2: the second particle

compute_g!(collision_data::CollisionData, p1, p2)

Compute relative velocity (vector and its magnitude) of two particles.

Positional arguments

• collision_data: the CollisionData instance where the computed velocity will be stored
• p1: the first particle
• p2: the second particle

compute_g_new_ionization!(collision_data, interaction, E_i, energy_splitting)

Compute post-ionization magnitudes of the relative velocities of the impacting and newly created electrons.

Positional arguments

• collision_data: the CollisionData instance where the computed velocity will be stored
• interaction: the Interaction instance for the electron-neutral interaction which produced the ion
• E_i: ionization energy of the neutral species in electron-volt
• energy_splitting: how the energy is divided among the electrons: if ElectronEnergySplitEqual, the energy is split equally, if ElectronEnergySplitZeroE, one electron takes all of the energy

scatter!(rng, model, collision_data, interaction, p1, p2)

Scatter two particles using the scattering law of the elastic scattering model: isotropic scattering (Merzbild.scatter_vhs!) for the VHS model, VSS scattering (Merzbild.scatter_vss!) for the VSS model.

Positional arguments

• rng: the random number generator
• model: the AbstractScatteringModel singleton for the scattering model
• collision_data: the CollisionData instance which stores the center-of-mass velocity and the pre-collisional relative velocity of the particles, and to which the new post-collisional velocity will be written
• interaction: the Interaction instance for the colliding particles
• p1: the first colliding particle
• p2: the second colliding particle

scatter_vhs!(rng, collision_data, interaction, p1, p2)

Scatter two particles using VHS (isotropic) scattering.

Positional arguments

• rng: the random number generator
• collision_data: the CollisionData instance to which stores the center-of-mass velocity, the magnitude of the pre-collisional relative velocity of the particles and to which the new post-collisional velocity will be written
• interaction: the Interaction instance for the colliding particles
• p1: the first colliding particle
• p2: the second colliding particle

scatter_vss!(rng, collision_data, interaction, p1, p2)

Scatter two particles using VSS (anisotropic) scattering. The cosine of the deflection angle $\chi$ (the angle between the pre- and post-collisional relative velocities) is sampled as $\cos\chi = 2 R^{1/\alpha} - 1$, where $R$ is a uniformly distributed random number and $\alpha$ is the VSS exponent of the interaction; the azimuthal angle $\varepsilon$ is sampled uniformly in $[0, 2\pi)$. The post-collisional relative velocity is then computed by rotating the pre-collisional one.

The power $R^{1/\alpha}$ is evaluated as exp2(vss_inv_alpha * log2(R)) for the same reason as in Merzbild.sigma_vhs.

If the relative velocity is nearly aligned with the x-axis, the rotation is performed in the $(y,z)$ plane instead, as the general expressions become ill-conditioned.

Positional arguments

• rng: the random number generator
• collision_data: the CollisionData instance which stores the center-of-mass velocity, the pre-collisional relative velocity of the particles (vector and magnitude), and to which the new post-collisional velocity will be written
• interaction: the Interaction instance for the colliding particles
• p1: the first colliding particle
• p2: the second colliding particle

References

• G.A. Bird, Eq. (2.22), Molecular gas dynamics and the direct simulation of gas flows, Clarendon Press, Oxford, 1994.
• K. Koura, H. Matsumoto, Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential. Phys. Fluids A, 1991.

scatter_electron_vhs!(rng, particle_electron, g_new)

Scatter an electron using VHS (isotropic) scattering and re-scale its relative velocity to g_new. This DOES NOT add the velocity of the center of mass to the electron.

Positional arguments

• rng: the random number generator
• particle_electron: the electron particle to scatter off of the neutral particle
• g_new: the magnitude of the post-collisional relative velocity

scatter_ionization_electrons_and_ion!(rng, collision_data, p_e1, p_e2, p_ion, mass_ratio)

Scatter electrons and ion after an ionization reaction using VHS (isotropic) scattering.

Positional arguments

• rng: the random number generator
• collision_data: the CollisionData instance to which stores the center-of-mass velocity the magnitude of the pre-collisional relative velocity of the electron and the neutral, and the post-collisional magnitudes of the velocities of the electrons
• p_e1: the first electron particle to scatter off of the neutral
• p_e2: the second electron particle to scatter off of the neutral
• p_ion: the ion produced in the ionization reaction
• mass_ratio: ratio of the electron mass to the ion mass

References

• K. Nanbu, Eqns. (47)-(53b), IEEE Trans. Plasma. Sci., 2000

sigma(model, interaction, g)

Compute the total elastic collision cross-section for a given scattering model. The VSS model uses the same power law as the VHS model.

Positional arguments

• model: the AbstractScatteringModel singleton tag of the scattering model
• interaction: the Interaction instance
• g: the relative velocity of the collision

Returns

The value of the computed cross-section.

sigma_vhs(interaction, g)

Computes the VHS cross-section.

The power law is evaluated as exp2(vhs_exp * log2(g)) rather than g^vhs_exp: Float64^Float64 carries the logarithm in double-double precision so that the result stays below 1 ulp for any exponent, which costs roughly twice the work. Since vhs_exp is a small fixed exponent, the plain decomposition is accurate to ~2 ulp.

Positional arguments

• interaction: the Interaction instance
• g: the relative velocity of the collision

Returns

The value of the computed cross-section.

@scattering_barrier model call

Turn the Merzbild.ScatteringModel enum value model into the corresponding scattering model singleton and insert it into the function call call as the second positional argument (that is, directly after rng).

The macro expands into an if/elseif chain with a literal singleton in each arm, so that each arm is a separate call with a concrete singleton type. This acts as a function barrier: the branch is taken once per call of a collision routine, and everything below it is compiled for one fixed model, with the cross-section and scattering functions inlined and no dispatch left in the collision loop.

Example

@scattering_barrier interaction_l.model ntc!(rng, collision_factors, collision_data, interaction,                                             particles, pia, cell, species, Δt, V; dw_tol=dw_tol)

Tuple of all (enum value, singleton tag) pairs of the implemented elastic scattering models, used by Merzbild.@scattering_barrier to generate the enum-to-tag conversion. Any newly added scattering model has to be listed here.

parse_scattering_model(name)

Convert the name of an elastic scattering model, as written in an interaction data TOML file, to the corresponding Merzbild.ScatteringModel enum value. The comparison is case-insensitive; the recognized names are "VHS" and "VSS".

Positional arguments

• name: the name of the scattering model

Returns

The ScatteringModel enum value.

Throws

ArgumentError if the model name is not recognized.

compute_tabulated_cs_constant_continuation(tabulated_cs_data, E_coll)

Compute an energy-dependent cross-section from tabulated data using linear interpolation. If the energy is smaller than the minimum energy used for the tabulation, return first element of the table; if energy is larger than the maximum energy used for the tabulation, return last element of the table.

Positional arguments

• tabulated_cs_data: a TabulatedCSData instance with the tabulated energies and cross-section values
• E_coll: the collision energy for which to compute the cross-section

Returns

The value of the cross-section.

compute_tabulated_cs_zero_continuation(tabulated_cs_data, E_coll)

Compute an energy-dependent cross-section from tabulated data using linear interpolation. If the energy is outside of the range of energies used for the tabulation, returns 0.0.

Positional arguments

• tabulated_cs_data: a TabulatedCSData instance with the tabulated energies and cross-section values
• E_coll: the collision energy for which to compute the cross-section

Returns

The value of the cross-section.

compute_cross_sections_only!(computed_cs, interaction, g, electron_neutral_interactions, neutral_species_index, extend::CSExtend)

Compute electron-impact ionization and excitation cross-sections, and electron-neutral elastic scattering crosss-sections, and return collision energy in eV. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• computed_cs: the vector of ComputedCrossSection instances in which the computed values will be stored
• interaction: the Interaction instance describing the electron-neutral interaction being considered
• g: the magnitude of the relative collision velocity
• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data
• neutral_species_index: the index of the neutral species being considered

Returns

The electron-neutral collision energy in eV.

compute_cross_sections!(computed_cs, interaction, g, electron_neutral_interactions, neutral_species_index; extend::CSExtend=CSExtendConstant)

Compute electron-impact ionization and excitation cross-sections, electron-neutral elastic scattering cross-sections, total collision cross-section, probabilities of the different processes, and return collision energy in eV. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance. Out-of-tabulation-range values are treated according to what extend method is used.

Positional arguments

• computed_cs: the vector of ComputedCrossSection instances in which the computed values will be stored
• interaction: the Interaction instance describing the electron-neutral interaction being considered
• g: the magnitude of the relative collision velocity
• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data
• neutral_species_index: the index of the neutral species being considered

Keyword arguments

• extend: enum of CSExtend type that sets how out-of-range energy values are treated when computing cross-sections

Returns

The electron-neutral collision energy in eV.

get_cs_total(electron_neutral_interactions, computed_cs, neutral_species_index)

Get the value of a computed total collision cross-section for electron-neutral interactions. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data
• computed_cs: the vector of ComputedCrossSection instances holding the computed cross-section values
• neutral_species_index: the index of the neutral species being considered

Returns

The value of the total collision cross-section.

get_cs_elastic(electron_neutral_interactions, computed_cs, neutral_species_index)

Get the value of a computed electron-neutral elastic scattering cross-section. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data
• computed_cs: the vector of ComputedCrossSection instances holding the computed cross-section values
• neutral_species_index: the index of the neutral species being considered

Returns

The value of the elastic scattering cross-section.

get_cs_ionization(electron_neutral_interactions, computed_cs, neutral_species_index)

Get the value of a computed electron-impact ionization cross-section. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data
• computed_cs: the vector of ComputedCrossSection instances holding the computed cross-section values
• neutral_species_index: the index of the neutral species being considered

Returns

The value of the electron-impact ionization scattering cross-section.

get_ionization_threshold(electron_neutral_interactions, neutral_species_index)

Get the ionization threshold energy. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• electron_neutral_interactions the ElectronNeutralInteractions instance storing the tabulated cross-section data
• neutral_species_index: the index of the neutral species being considered

Returns

The ionization threshold energy of a specific neutral species.

get_electron_energy_split(electron_neutral_interactions, neutral_species_index)

Get the way electron energy is split during ionization for a specific electron-neutral interaction. Here neutral_species_index is the index of the neutral species being considered in the overall array of the Species instances that includes all species in the simulation and was used to construct the ElectronNeutralInteractions instance.

Positional arguments

• electron_neutral_interactions the ElectronNeutralInteractions instance storing the tabulated cross-section data
• neutral_species_index: the index of the neutral species being considered

Returns

The ElectronEnergySplit value for the specific electron-neutral interaction.

sample_normal_rands!(rng, collision_data_fp, n_local)

Sample 3*n_local normally distributed random numbers (with mean 0 and variance 1), one for each velocity component, and write them to a CollisionDataFP instance.

Positional arguments

• rng: the random number generator
• collision_data_fp: CollisionDataFP instance used for storing collisional quantities
• n_local: the number of particles to sample the numbers for

scale_norm_rands!(collision_data_fp, particles, indexer, local_w)

Scale sampled normally distributed random numbers velocity component-wise, so that their weight-averaged means for each component are exactly 0, and their weight-averaged variances are exactly 1. Weighting by the particles' computational weights is required so that the Langevin step conserves momentum (and hence energy) exactly for variable-weight particles; for equal-weight particles this reduces to the plain (unweighted) normalization.

Positional arguments

• collision_data_fp: CollisionDataFP instance used for storing collisional quantities
• particles: ParticleVector of the particles being collided
• indexer: the ParticleIndexer for the cell/species being collided (pia.indexer[cell, species])
• local_w: the total computational weight of the particles in the cell

compute_relaxation_time(interaction, species_data, species, V, es_old, local_w)

Compute the relaxation time for a single-species flow, as given by $\tau = 2\mu / p$.

Positional arguments

• interaction: the Interaction instance for the self-interaction of the species in question
• species_data: the vector of SpeciesData
• species: the index of the species for which collisions are performed
• V: cell volume
• es_old: the total kinetic energy of the particles in the cell, divided by the species' molecular mass
• local_w: the total computational weight of the particles in the cell

Returns

The elastic collision relaxation time.

TabulatedCSData

Structure for storing tabulated cross-section data as a function of relative collision energy

Fields

• n_vals: number of values stored
• E: array of relative collision energies
• sigma: array of cross-section values at these energies
• ΔE: how much energy is lost (gained) in the collision in case it is inelastic

ElasticScattering

Structure to hold data on an elastic scattering cross-section.

Fields

• data: a TabulatedCSData instance holding the cross-section values
• scattering: the ScatteringLaw model for the scattering of the particles

ExcitationSink

Structure to hold data on electron-impact electronic excitation cross-sections for a specific species.

Fields

• n_reactions: number of electron-impact electronic excitation reactions
• data: an array of TabulatedCSData instances (of length n_reactions) holding the cross-section values for each reaction
• scattering: the ScatteringLaw model for the scattering of the particles

Ionization

Structure to hold data on an electron-impact ionization cross-section.

Fields

• data: a TabulatedCSData instance holding the cross-section values
• scattering: the ScatteringLaw model for the scattering of the electrons
• split: the ElectronEnergySplit model for energy splitting across the primary and secondary electrons

Find a chemical species in an LXCat-format XML; returns the first instance of the species found and the id of the "Groups" element in which the species was found. The following conditions need to be satisfied for a tuple (i,j) for it to be returned: tag(xml_data[i]) == "Groups" and attributes(xml_data[i][j])["id"] == species_name.

Positional arguments

• xml_data: the LXCat-format XML data to be searched
• species_name: the name of the species to search for

Returns

Tuple (i,j) denoting the index of the group and group element in which the species is located

load_elastic_data(xml_data)

Load electron-neutral elastic scattering data from the part of an LXCat-format XML file for a specific species. The following should hold for the data to be loaded: tag(xml_data[i]) == "Processes", attributes(xml_data[i][j])["type"] == "Ionization". This will load the first set of data found for the ionization process.

Positional arguments

• xml_data: the LXCat-format XML data to be searched

Returns

TabulatedCSData structure containing the electron-neutral elastic scattering cross-section data.

Throws

DataMissingException if data not found or not all required data present.

load_ionization_data(xml_data)

Load electron-impact ionization data from the part of an LXCat-format XML file for a specific species. The following should hold for the data to be loaded: tag(xml_data[i]) == "Processes", attributes(xml_data[i][j])["type"] == "Ionization". This will load the first set of data found for the ionization process.

Positional arguments

• xml_data: the LXCat-format XML data to be searched

Returns

TabulatedCSData structure containing the electron-impact ionization cross-section data.

Throws

DataMissingException if data not found or not all required data present.

GridCell{D}

Struct for keeping track of merging-related quantities in a velocity grid cell with a D-dimensional position vector.

Fields

• np: number of particles in cell
• w: total computational weight of particles in cell
• v_mean: mean velocity of particles in cell
• v_std_sq: variance of velocity of particles in cell
• x_mean: mean position of particles in cell
• x_std_sq: variance of position of particles in cell
• particle_index1: index of first particle in the cell (first as in the first particle that is found to belong to the cell)
• particle_index2: index of second particle in the cell (second as in the second particle that is found to belong to the cell)
• w1: the post-merge weight to assign to the first particle in the cell
• w2: the post-merge weight to assign to the second particle in the cell
• v1: the post-merge velocity to assign to the first particle in the cell
• v2: the post-merge velocity to assign to the second particle in the cell
• x1: the post-merge position to assign to the first particle in the cell
• x2: the post-merge position to assign to the second particle in the cell

compute_velocity_extent!(merging_grid::GridN2Merge{D}, cell, species, species_data, phys_props::PhysProps) where D

Compute extent of velocity grid based on temperature in the cell.

Positional arguments:

• merging_grid: the grid merging (GridN2Merge) instance for which to compute the extent
• cell: the index of the cell in which the temperature was computed
• species: the index of the species for which the temperature was computed
• species_data: the array of Species data
• phys_props: the PhysProps instance containing the computed temperature

compute_velocity_extent!(merging_grid::GridN2Merge{D}, vx_extent, vy_extent, vz_extent) where D

Compute extent of velocity grid based on explicitly set extents.

Positional arguments:

• merging_grid: the grid merging (GridN2Merge) instance for which to compute the extent
• vx_extent: lower and upper bounds of the grid extent in the x velocity direction
• vy_extent: lower and upper bounds of the grid extent in the y velocity direction
• vz_extent: lower and upper bounds of the grid extent in the z velocity direction

compute_grid_index(merging_grid::GridN2Merge{D}, v) where D

Compute index of cell on the merging grid in which a velocity is located (the last 8 indices correspond to the velocity octants outside the grid).

Positional arguments:

• merging_grid: the grid merging (GridN2Merge) instance defining the velocity space grid
• v: the velocity vector

Returns:

Index of cell on the merging grid.

clear_merging_grid!(merging_grid::GridN2Merge{D}) where D

Resets all data for a merging grid instance.

Positional arguments:

• merging_grid: the grid merging (GridN2Merge) instance defining the velocity space grid

compute_grid!(merging_grid::GridN2Merge{D}, particles::ParticleVector{D}, pia, cell, species) where D

Compute all the required cell properties for a grid-based merge.

Positional arguments:

• merging_grid: the grid merging (GridN2Merge) instance defining the velocity space grid
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the cell index
• species: the species index

compute_new_particles!(rng, merging_grid::GridN2Merge{D}, particles::ParticleVector{D}, pia, cell, species) where D

Compute new particles based on the grid cell properties without checking particle locations. So particles may end up outside of the domain.

Positional arguments:

• rng: the random number generator instance
• merging_grid: the grid merging (GridN2Merge) instance defining the velocity space grid
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the cell index
• species: the species index

compute_new_particles!(rng, merging_grid::GridN2Merge{D}, particles::ParticleVector{D}, pia, cell, species, grid::Grid1DUniform) where D

Compute new particles based on the grid cell properties; placing out-of-domain particles back into the domain.

Positional arguments:

• rng: the random number generator instance
• merging_grid: the grid merging (GridN2Merge) instance defining the velocity space grid
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the cell index
• species: the species index
• grid: the Grid1DUniform grid

compute_new_particles!(rng, octree::OctreeMerge{D,2}, particles::ParticleVector{D}, pia, cell, species) where D

Compute post-merge particles with particles based on octree bin properties without checking or setting particle locations (for spatially homogeneous merging). N:2 merging in each bin.

Positional arguments:

• rng: the random number generator instance
• octree: the OctreeMerge instance
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the cell index
• species: the species index

compute_new_particles!(rng, octree::OctreeMerge{D,2}, particles::ParticleVector{D}, pia, cell, species, grid::Grid1DUniform) where D

Compute post-merge particles particles based on octree bin properties; placing out-of-domain particles back into the domain. N:2 merging in each bin.

Positional arguments:

• rng: the random number generator instance
• octree: the OctreeMerge instance
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the cell index
• species: the species index
• grid: the Grid1DUniform grid

compute_new_particles!(rng, octree::OctreeMerge{D,1}, particles::ParticleVector{D}, pia, cell, species) where D

Compute post-merge particles with particles based on octree bin properties without checking or setting particle locations (for spatially homogeneous merging). N:1 merging in each bin. It computes properties of all particles. And after the first merge, it scales the particles' velocities based on these properties to achieve conservation.

Positional arguments:

• rng: the random number generator instance
• octree: the OctreeMerge instance
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

compute_new_particles!(rng, octree::OctreeMerge{D,1}, particles::ParticleVector{D}, pia, cell, species, grid::Grid1DUniform) where D

Compute post-merge particles with particles based on octree bin properties, N:1 merging in each bin; placing out-of-domain particles back into the domain.. It computes properties of all particles. And after the first merge, it scales the particles based on these properties to achieve conservation.

Positional arguments:

• rng: the random number generator instance
• octree: the OctreeMerge instance
• particles: the ParticleVector instance of the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

write_back_to_particles!(octree::OctreeMerge{D,M}, particles::ParticleVector{D}, pia, cell, species) where {D,M}

Write back merged particles to the particle array (0D case - position of particles not set). It will only be used in the function computenewparticles! without a grid argument, which means no boundary handling is needed.

Positional arguments:

• octree: the OctreeMerge instance
• particles: the ParticleVector instance containing the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

Returns

The current particle index after writing.

write_back_to_particles!(octree::OctreeMerge{D,M}, particles::ParticleVector{D}, pia, cell, species, grid::Grid1DUniform) where {D,M}

Write back merged particles to the particle array with boundary clamping (1D case). Particles positioned outside the grid boundaries (grid.min_x, grid.max_x) are clamped to the boundary. It will only be used in the function computenewparticles! with grid::Grid1DUniform argument.

Positional arguments:

• octree: the OctreeMerge instance
• particles: the ParticleVector instance containing the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• grid: the Grid1DUniform instance

Returns

The current particle index after writing.

Relative tolerance below which a post-merge variance is treated as having collapsed to zero, see variance_scaling.

variance_scaling(var_before::SVector{N,Float64}, var_post::SVector{N,Float64}) where N

Scale the the components of a variance vector var_before by the square root of the ratio of the components of it and a vector var_post, checking for zero/negative values in var_post and non-finite values in the scaling factor.

A component of var_post is treated as zero if it is not larger than variance_scaling_rel_tol times the corresponding component of var_before, and a scaling factor of 1.0 is returned for it. A plain var_post > 0.0 check is not enough: when all post-merge particles collapse onto the mean (e.g. an N:1 merge of a single octree bin), var_post is zero only in exact arithmetic, and is in practice a round-off residue of order eps()^2 * var_before. Scaling by sqrt(var_before / var_post) would then amplify that residue by ~1/eps() and destroy the conservation of the mean, instead of restoring a variance that cannot be restored by scaling in the first place.

Positional arguments

• var_before: variance vector components of which to scale
• var_post: variance vector components by which to scale

vx_sign(octant)

Return sign of velocity vx of an octant in velocity space.

Positional arguments

• octant: index of the octant in velocity space

Returns

Sign of the x-velocity corresponding to the octant.

vy_sign(octant)

Return sign of velocity vy of an octant in velocity space.

Positional arguments

• octant: index of the octant in velocity space

Returns

Sign of the y-velocity corresponding to the octant.

vz_sign(octant)

Return sign of velocity vz of an octant in velocity space.

Positional arguments

• octant: index of the octant in velocity space

Returns

Sign of the z-velocity corresponding to the octant

base_multi_index_moments()

Base multi indices corresponding to conservation of mass, momentum, and energy components.

Positional arguments

None

Returns

A vector of multi-indices (3-tuples) which corresponding to mass, momentum, and directional energy.

compute_w_total_v0!(nnls_merging, particles::ParticleVector{D}, pia, cell, species) where D

Compute total computational weight of particles and mean velocity, as well as velocity bounds of the set of particles in each velocity direction.

Positional arguments

• nnls_merging: the NNLSMerge instance where the computed values will be stored
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

fill_powers!(powers, v, v0)

Fill a table of the powers of the components of the centered velocity / position vector v - v0, so that powers[j, e+1] holds (v[j] - v0[j])^e. Only the first size(powers, 1) components are considered, and powers up to size(powers, 2) - 1 are computed. Evaluating the moments off such a table avoids the runtime-exponent ^ calls that dominate the cost of building the LHS matrix.

Positional arguments

• powers: the table of powers to fill
• v: the velocity / position vector
• v0: the mean velocity / position vector

ccm_vel(vel_powers, mim)

Compute an unweighted central velocity moment from a table of powers filled by fill_powers!.

Positional arguments

• vel_powers: the table of powers of the centered velocity components
• mim: the 3-dimensional multi-index

Returns

Computed unweighted central moment.

ccm_pos(pos_powers, mim, D)

Compute an unweighted central spatial moment from a table of powers filled by fill_powers!. Only the first D components of the multi-index are used, matching the dimensionality of the particle position vectors.

Positional arguments

• pos_powers: the table of powers of the centered position components
• mim: the 3-dimensional multi-index
• D: the dimensionality of the position vectors

Returns

Computed unweighted central moment.

compute_lhs_and_rhs!(nnls_merging::NNLSMerge{D}, lhs_matrix, vel_pos_matrix, particles::ParticleVector{D}, pia, cell, species) where D

Compute LHS matrix and RHS vector for NNLS merging. Returns the pre-merge number of particles.

Positional arguments

• nnls_merging: the NNLSMerge instance where the RHS vector will be stored
• lhs_matrix: the matrix of size n_total_conserved x n_particles, where n_total_conserved is the number of conserved moments and n_particles is the pre-merge number of particles
  • any fictitious particles
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

Returns

The pre-merge number of particles.

compute_lhs_and_rhs_rate_preserving!(nnls_merging::NNLSMerge{D}, lhs_matrix, vel_pos_matrix,
                                     interaction, electron_neutral_interactions, computed_cs,
                                     particles::ParticleVector{D}, pia, cell, species, neutral_species_index, extend) where D

Compute LHS matrix and RHS vector for the rate-preserving NNLS merging (for electrons). Approximate elastic scattering and electron-impact ionization rates are conserved.

Positional arguments

• nnls_merging: the NNLSMerge instance where the RHS vector will be stored
• lhs_matrix: the matrix of size n_total_conserved x n_particles, where n_total_conserved is the number of conserved moments and n_particles is the pre-merge number of particles
  • any fictitious particles
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• interaction: the Interaction instance describing the electron-neutral interaction being considered
• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data used to compute the rates
• computed_cs: the vector of ComputedCrossSection instances in which the computed values will be stored
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• neutral_species_index: the index of the neutral species which is the collision partner in the electron-neutral collisions for which approximate rates are being preserved.
• extend: enum of CSExtend type that sets how out-of-range energy values are treated when computing cross-sections

compute_lhs_and_rhs_rate_preserving!(nnls_merging::NNLSMerge{D}, lhs_matrix, vel_pos_matrix,
                                     interaction, electron_neutral_interactions, computed_cs,
                                     particles::ParticleVector{D}, particles_neutral::ParticleVector{D}, pia, cell, species, neutral_species_index, extend) where D

Compute LHS matrix and RHS vector for the rate-preserving NNLS merging (for electrons). Exact elastic scattering and electron-impact ionization rates are conserved.

Positional arguments

• nnls_merging: the NNLSMerge instance where the RHS vector will be stored
• lhs_matrix: the matrix of size n_total_conserved x n_particles, where n_total_conserved is the number of conserved moments and n_particles is the pre-merge number of particles
  • any fictitious particles
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• interaction: the Interaction instance describing the electron-neutral interaction being considered
• electron_neutral_interactions: the ElectronNeutralInteractions instance storing the tabulated cross-section data used to compute the rates
• computed_cs: the vector of ComputedCrossSection instances in which the computed values will be stored
• particles: the ParticleVector instance containing the particles that are being merged
• particles_neutrals: the ParticleVector instance containing the neutral collision partner particles
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• neutral_species_index: the index of the neutral species which is the collision partner in the electron-neutral collisions for which approximate rates are being preserved.
• extend: enum of CSExtend type that sets how out-of-range energy values are treated when computing cross-sections

scale_lhs_rhs_variance!(nnls_merging::NNLSMerge{D}, lhs_matrix, lhs_ncols) where D

Scale the LHS and RHS of the NNLS system using the computed variances of the particles velocities in the $x$, $y$, $z$ directions. If any of the variances is smaller than 1e-6, then $v_{ref}$ is used. Each moment with multi-index (i,j,k) is scaled by $(1/Ev[1])^{i}(1/Ev[2])^{j}(1/Ev[3])^{k}$.

Positional arguments

• nnls_merging: the NNLSMerge instance
• lhs_matrix: the matrix of the LHS
• lhs_ncols: number of columns in the LHS matrix

scale_lhs_rhs_vref!(nnls_merging::NNLSMerge{D}, lhs_matrix, lhs_ncols) where D

Scale the LHS and RHS of the NNLS system using the reference velocity $v_{ref}$. Each moment is scaled by $(1/v_{ref})^{n_{tot}}$, where $n_{tot}$ is the total order of the moment (i.e. for a moment with multi-index (i,j,k) the total order is i+j+k).

Positional arguments

• nnls_merging: the NNLSMerge instance
• lhs_matrix: the matrix of the LHS
• lhs_ncols: number of columns in the LHS matrix

scale_lhs_rhs_spatial_variance!(nnls_merging::NNLSMerge{D}, lhs_matrix, lhs_ncols) where D

Scale the spatial moments in the LHS and RHS of the NNLS system using the computed variances of the particles positions in the $x$, $y$, $z$ directions. If any of the variances is smaller than 1e-6, then no scaling is done. Each spatial moment with multi-index (i,j,k) is scaled by $(1/Ex[1])^{i}(1/Ex[2])^{j}(1/Ex[3])^{k}$, for non-3-dimensional particle vectors only the first D components are considered and scaled.

Positional arguments

• nnls_merging: the NNLSMerge instance
• lhs_matrix: the matrix of the LHS
• lhs_ncols: number of columns in the LHS matrix

scale_lhs_rhs!(nnls_merging::NNLSMerge{D}, lhs_matrix, scaling, lhs_ncols) where D

Scale the LHS and RHS of the NNLS system using either the reference velocity or the computed variances of the velocity in each direction.

Positional arguments

• nnls_merging: the NNLSMerge instance
• lhs_matrix: the matrix of the LHS
• scaling: how to scale entries (:vref or :variance)
• lhs_ncols: number of columns in the LHS matrix

scale_lhs_rhs_rate_preserving!(nnls_merging, lhs_matrix, ref_k_elastic, ref_k_ion, scaling, lhs_ncols) where D

Scale the LHS and RHS of the NNLS system for the rate-preserving electron merging using the reference velocity $v_{ref}$ and reference elastic scattering and ionization cross-sections. Each moment is scaled by $(1/v_{ref})^{n_{tot}}$, where $n_{tot}$ is the total order of the moment (i.e. for a moment with multi-index (i,j,k) the total order is i+j+k). The entries in the LHS and RHS corresponding to the rates are scaled by ref_k_elastic, ref_k_ion.

Positional arguments

• nnls_merging: the NNLSMerge instance
• lhs_matrix: the matrix of the LHS
• ref_k_elastic: the reference elastic collision rate coefficient
• ref_k_ion: the reference ionization rate coefficient
• scaling: how to scale entries (:vref or :variance)
• lhs_ncols: number of columns in the LHS matrix

compute_post_merge_particles_nnls!(nnls_merging::NNLSMerge{3}, x::Vector{Float64}, particles::ParticleVector{3},
                                   pia, cell, species, lhs_ncols,
                                   vel_pos_matrix,
                                   max_err, w_threshold, work_index, column_norms)

Compute post-merge particles based on the solution of the NNLS problem. This will replace the particles with the post-merge ones and delete any extraneous particles.

Positional arguments

• nnls_merging: the NNLSMerge instance
• x: the solution vector of the NNLS system
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• lhs_ncols: the number of columns in the LHS matrix
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• max_err: maximum allowed value of the residual of the NNLS system
• w_threshold: the relative (w.r.t the total computational weight of the particles being merge) value of the computational weight below which particles are discarded
• work_index: the index of the NNLSWorkspace used to solve the NNLS system
• column_norms: the vector of the column-wise norms of the LHS matrix

compute_post_merge_particles_nnls!(nnls_merging::NNLSMerge{2}, x::Vector{Float64}, particles::ParticleVector{2},
                                   pia, cell, species, lhs_ncols,
                                   vel_pos_matrix,
                                   max_err, w_threshold, work_index, column_norms)

Compute post-merge particles based on the solution of the NNLS problem. This will replace the particles with the post-merge ones and delete any extraneous particles.

Positional arguments

• nnls_merging: the NNLSMerge instance
• x: the solution vector of the NNLS system
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• lhs_ncols: the number of columns in the LHS matrix
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• max_err: maximum allowed value of the residual of the NNLS system
• w_threshold: the relative (w.r.t the total computational weight of the particles being merge) value of the computational weight below which particles are discarded
• work_index: the index of the NNLSWorkspace used to solve the NNLS system
• column_norms: the vector of the column-wise norms of the LHS matrix

compute_post_merge_particles_nnls!(nnls_merging::NNLSMerge{1}, x::Vector{Float64}, particles::ParticleVector{1},
                                   pia, cell, species, lhs_ncols,
                                   vel_pos_matrix,
                                   max_err, w_threshold, work_index, column_norms)

Compute post-merge particles based on the solution of the NNLS problem. This will replace the particles with the post-merge ones and delete any extraneous particles.

Positional arguments

• nnls_merging: the NNLSMerge instance
• x: the solution vector of the NNLS system
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• lhs_ncols: the number of columns in the LHS matrix
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• max_err: maximum allowed value of the residual of the NNLS system
• w_threshold: the relative (w.r.t the total computational weight of the particles being merge) value of the computational weight below which particles are discarded
• work_index: the index of the NNLSWorkspace used to solve the NNLS system
• column_norms: the vector of the column-wise norms of the LHS matrix

compute_post_merge_particles_nnls!(nnls_merging::NNLSMerge{1}, x::Vector{Float64}, particles::ParticleVector{1},
                                   pia, cell, species, lhs_ncols,
                                   vel_pos_matrix,
                                   max_err, w_threshold, work_index, column_norms)

Compute post-merge particles based on the solution of the NNLS problem. This will replace the particles with the post-merge ones and delete any extraneous particles.

Positional arguments

• nnls_merging: the NNLSMerge instance
• x: the solution vector of the NNLS system
• particles: the ParticleVector instance containing the particles that are being merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged
• lhs_ncols: the number of columns in the LHS matrix
• vel_pos_matrix: the matrix of size 6 x n_particles, where n_particles is the pre-merge number of particles
  • any fictitious particles, where their velocities and positions will be stored
• max_err: maximum allowed value of the residual of the NNLS system
• w_threshold: the relative (w.r.t the total computational weight of the particles being merge) value of the computational weight below which particles are discarded
• work_index: the index of the NNLSWorkspace used to solve the NNLS system
• column_norms: the vector of the column-wise norms of the LHS matrix

OctreeCell

Struct holding computed bin properties needed for refinement of an octree bin.

Fields

• np: number of particles in cell
• w: total computational weight of particles in cell
• v_min: vector of the per-component lower bounds of the velocities in the cell
• v_max: vector of the per-component upper bounds of the velocities in the cell
• depth: level of refinement the cell is at
• can_be_refined: whether the cell can be refined further

OctreeFullCell{D,M}

Struct holding computed bin properties required to merge the particles in a bin down to M particles.

Fields

• v_mean: mean velocity of particles in cell
• v_std_sq: variance of velocity of particles in cell
• x_mean: mean position of particles in cell
• x_std_sq: variance of position of particles in cell
• particle_indices: indices of particles in the cell
• weights: the post-merge weights to assign to the particles in the cell
• velocities: the post-merge velocities to assign to the particles in the cell
• positions: the post-merge positions to assign to the particles in the cell

fill_bins(Nbins)

Fill the octree bins structs with zero data, used as a utility function for initialization.

Positional arguments:

• Nbins: number of OctreeCell bins to create

Returns: An array of Nbins OctreeCell instances filled with zeros.

fill_full_bins(::Val{D}, ::Val{M}, Nbins) where {D,M}

Fill the octree bins full structs with zero data, used as a utility function for initialization.

Positional arguments:

• D: dimension of position vectors of particles to be merged
• M: number of post-merge particles in each bin
• Nbins: number of OctreeFullCell bins to create

Returns: An array of Nbins OctreeFullCell instances filled with zeros.

clear_octree!(octree)

Reset octree before doing a new merge.

Positional arguments:

• octree: the OctreeMerge instance

resize_octree_buffers!(octree, n_particles)

Check and resize octree buffers if needed to accommodate a larger number of particles. The size of the buffers is set to n_particles + DELTA_PARTICLES if the their sizes are smaller than n_particles.

Positional arguments:

• octree: the OctreeMerge instance
• n_particles: the current number of particles being merged

compute_octant(particle_v, v_middle)

Compute octant of particle velocity relative to a v_middle` The order of the octants is:

1. - - -2. + - -3. - + -4. + + -5. - - +6. + - +7. - + +8. + + +

Positional arguments:

• particle_v: the velocity of the particle
• v_middle: the velocity relative to which the octant is computed

Returns:

The octant number

bin_bounds_inherit!(octree, bin_id, v_min_parent, v_max_parent, v_middle, octant)

Compute new bin bounds of one of the 8 octant sub-bins inheriting bounds of parent bin. The order of the octants is

1. - - -2. + - -3. - + -4. + + -5. - - +6. + - +7. - + +8. + + +

Positional arguments:

• octree: the OctreeMerge instance
• bin_id: the index of the bin for which the velocity bounds are being recomputed
• v_min_parent: the velocity vector of the lower bound of the velocities of the parent bin
• v_max_parent: the velocity vector of the upper bound of the velocities of the parent bin
• v_middle: the velocity across which the split is being performed
• octant: the octant of the parent bin to which this bin corresponds

bin_bounds_recompute!(octree, bin_id, bs, be, particles::ParticleVector{D}) where D

Recompute bin bounds based on particle velocities by setting them to the smallest and largest velocities of the particles in each velocity direction.

Positional arguments:

• octree: the OctreeMerge instance
• bin_id: the index of the bin for which the velocity bounds are being recomputed
• bs: index of the first particle in bin
• be: index of the last particle in bin
• particles: the ParticleVector instance of the particles to be merged

compute_v_mean!(octree, bs, be, particles::ParticleVector{D}) where D

Compute mean velocity of particles in a bin.

Positional arguments:

• octree: the OctreeMerge instance
• bs: index of the first particle in bin
• be: index of the last particle in bin
• particles: the ParticleVector instance of the particles to be merged

compute_v_median!(octree, bs, be, particles::ParticleVector{D}) where D

Compute median velocity of particles in a bin. NOTE: allocates memory and is probably not fully correct!

Positional arguments:

• octree: the OctreeMerge instance
• bs: index of the first particle in bin
• be: index of the last particle in bin
• particles: the ParticleVector instance of the particles to be merged

get_new_bin_id(i, bin_id, Nbins)

Get index of a newly created bin once a bin with index bin_id is split into 8 sub-bins. The bin_id index is re-used for the 1-st sub-bin, and the other 7 sub-bins are tacked onto the end of the list of bins (so they have indices Nbins + 1, Nbins + 2, ..., where Nbins was the total number of octree bins before the split).

Positional arguments

• i: index of sub-octant (ranging from 1 to 8) of the split bin for which to return a new bin index
• bin_id: index of bin being split
• Nbins: total number of octree bins before the split

Returns

Index of a newly created bin corresponding to a bin created from sub-octant i of bin bin_id.

split_bin!(octree::OctreeMerge{D,M}, bin_id, particles::ParticleVector{D}) where {D,M}

Sort particles into sub-bins of a bin with index bin_id (by splitting it into octants), keeping track of which sub-bins particles end up in. Also sets the velocity bounds of the new bins.

Positional arguments:

• octree: the OctreeMerge instance
• bin_id: octree bin index
• particles: the ParticleVector instance of the particles to be merged

compute_bin_props!(octree::OctreeMerge{D,M}, bin_id, particles::ParticleVector{D}) where {D,M}

Compute properties in a bin required for merging: total computational weight, mean velocity and position, standard deviation of particle velocities and positions.

Positional arguments:

• octree: the OctreeMerge instance
• bin_id: octree bin index
• particles: the ParticleVector instance of the particles to be merged

get_bin_post_merge_np(octree::OctreeMerge{D,M}, bin_id) where {D,M}

Get number of post-merge particles in a bin: M if the number of particles in the bin is >= M, otherwise the number of particles in the bin is returned.

Positional arguments:

• octree: the OctreeMerge instance
• bin_id: octree bin index

Returns

The number of post-merge particles in a single octree bin (0, 1, or 2).

init_octree!(octree::OctreeMerge{D,M}, particles::ParticleVector{D}, pia, cell, species) where {D,M}

Initialize the top bin in an octree by copying particle indices and setting bin bounds.

Positional arguments

• octree: the OctreeMerge instance
• particles: the ParticleVector instance containing the particles to be merged
• pia: the ParticleIndexerArray instance
• cell: the index of the grid cell in which particles are being merged
• species: the index of the species being merged

compute_octree!(octree::OctreeMerge{D,M}, particles::ParticleVector{D}, target_np) where {D,M}

Perform refinement of an octree for N:2 merging until target number of particles reached or nothing left to refine.

Positional arguments

• octree: the OctreeMerge instance
• particles: the ParticleVector instance containing the particles to be merged
• target_np: the target post-merge number of particles; the post-merge number of particles will not exceed this value but may be not exactly equal to it

Cell1D

Cell element of a 1-D grid

Fields

• xlo: coordinate of the left end of the element
• xhi: coordinate of the right end of the element
• V: cell volume
• inv_V: inverse of cell volume

Cell1D(xlo, xhi, V)

Create cell of 1-D uniform grid

Positional arguments

• xlo: coordinate of the left end of the element
• xhi: coordinate of the right end of the element
• V: cell volume

get_cell(grid1duniform::Grid1DUniform, x_pos::SVector{D,Float64}) where D

Find in which cell of a 1-D uniform grid the coordinate is located

Positional arguments

• grid1duniform: the 1-D uniform grid
• x_pos: the D-dimensional coordinate vector for which the cell index is to be determined (only the first component is used)

Returns

Index of cell the particle is located in

convect_single_particle!(rng, grid::Grid1DUniform, bc_list, particle::Particle{D}, species, Δt) where D

Convect a singe particle on a 1-D uniform grid.

Positional arguments

• rng: the random number generator
• grid: the grid on which the convection is performed
• bc_list: the Tuple of boundary conditions (left and right wall)
• particle: the particle to be convected
• species: the index of the species being convected
• Δt: the convection timestep

convect_single_particle!(rng, grid::Grid1DUniform, bc_list, particle::Particle{D}, species, surf_props::SurfProps, mass, Δt) where D

Convect a singe particle on a 1-D uniform grid, updating surface properties if it collides with a wall.

Positional arguments

• rng: the random number generator
• grid: the grid on which the convection is performed
• bc_list: the Tuple of boundary conditions (left and right wall)
• particle: the particle to be convected
• species: the index of the species being convected
• surf_props: the SurfProps struct where the computed surface properties will be stored
• mass: the molecular mass of the species
• Δt: the convection timestep

convect_single_particle_periodic!(grid::Grid1DUniform, particle::Particle{D}, Δt) where D

Convect a singe particle on a 1-D uniform grid assuming a periodic grid.

Positional arguments

• grid: the grid on which the convection is performed
• particle: the particle to be convected
• Δt: the convection timestep

set_x(x::SVector{D,Float64}, val::Float64) where D

Set the x component of a Particle position vector.

Positional arguments

• x: position vector
• val: value to write in the 1-st component of the position vector

Returns

New position vector.

AbstractFieldBC1D

Abstract type for boundary conditions imposed on the electrostatic field in 1-D simulations. Not to be confused with AbstractBC, which describes particle-surface interaction.

poisson_unknowns(n_cells, bc_left, bc_right)

Compute the number of unknowns of the discrete 1-D Poisson system and the offset between the index of an unknown and the index of the corresponding node, for a given pair of boundary conditions. Throws an error for the boundary condition combinations which are not admissible (Neumann on both sides, which is singular and leaves no way of fixing the gauge, and a periodic boundary condition paired with a non-periodic one).

Positional arguments

• n_cells: number of cells of the grid
• bc_left: the AbstractFieldBC1D boundary condition on the left boundary
• bc_right: the AbstractFieldBC1D boundary condition on the right boundary

Returns

Tuple of the number of unknowns and the node offset.

set_neumann_diagonal_left!(poisson_solver, bc_left)

Replace the first row of the discrete Poisson matrix by the finite-volume integration over the left half-cell if the left boundary condition is a Neumann one, otherwise do nothing. The diagonal entry is halved with respect to an interior row, as the control volume of the boundary node is $\Delta x / 2$.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• bc_left: the AbstractFieldBC1D boundary condition on the left boundary

set_neumann_diagonal_right!(poisson_solver, bc_right)

Replace the last row of the discrete Poisson matrix by the finite-volume integration over the right half-cell if the right boundary condition is a Neumann one, otherwise do nothing. The diagonal entry is halved with respect to an interior row, as the control volume of the boundary node is $\Delta x / 2$.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• bc_right: the AbstractFieldBC1D boundary condition on the right boundary

factorize_poisson!(poisson_solver)

Pre-compute the forward elimination of the tridiagonal matrix of the Thomas algorithm and store the inv_m, cp, and a_inv_m coefficients. These depend only on the matrix and not on the right-hand side, so with a fixed matrix only the right-hand side recurrence and the back substitution have to be performed per solve.

This assumes that the matrix never changes during a simulations, i.e. PIC grid and BC types are fixed (the prescribed values of the potential/electric fields at Dirichlet/Neumann BCs can vary).

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance to be factorized

poisson_rhs_left!(poisson_solver, field_props, bc_left)

Apply the correction of the first row of the right-hand side of the discrete Poisson system stemming from the left boundary condition: for a Dirichlet boundary condition, the eliminated known value of the potential in the boundary node is added to the right-hand side; for a Neumann boundary condition, the row is that of the finite-volume integration over the left half-cell, with a halved charge density and the prescribed value of the electric field.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the charge density
• bc_left: the AbstractFieldBC1D boundary condition on the left boundary

poisson_rhs_right!(poisson_solver, field_props, bc_right)

Apply the correction of the last row of the right-hand side of the discrete Poisson system stemming from the right boundary condition, mirroring poisson_rhs_left!.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the charge density
• bc_right: the AbstractFieldBC1D boundary condition on the right boundary

assemble_poisson_rhs!(poisson_solver, field_props)

Assemble the right-hand side of the discrete Poisson system from the charge density stored in field_props, applying the boundary corrections. In the periodic case, the mean charge density over the unique nodes is subtracted (otherwise the singular periodic system has no solution at all) and stored in field_props.net_charge_density; a non-zero value of the latter means that the simulation carries a net charge, which is unphysical in a periodic domain.

The charge density is not modified.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the charge density

solve_tridiagonal!(poisson_solver)

Solve the tridiagonal system using the pre-computed factorization of the matrix (see factorize_poisson!) and the right-hand side stored in poisson_solver.rhs, storing the result in poisson_solver.x.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance

apply_gauge!(poisson_solver, field_props)

Shift the potential in the unique nodes 1:n_cells of a periodic domain so that its mean is zero. On a periodic uniform grid this gauge is identical to requiring a zero mean of the cell-averaged potential, as each unique node contributes to exactly two cells.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the potential

finalize_potential_left!(field_props, bc_left)

Write the known value of the potential into the left boundary node for a Dirichlet boundary condition; for a Neumann boundary condition the boundary node is an unknown of the system and nothing is done.

Positional arguments

• field_props: the ElectrostaticFieldProps instance holding the potential
• bc_left: the AbstractFieldBC1D boundary condition on the left boundary

finalize_potential_right!(field_props, bc_right)

Write the known value of the potential into the right boundary node for a Dirichlet boundary condition, mirroring finalize_potential_left!.

Positional arguments

• field_props: the ElectrostaticFieldProps instance holding the potential
• bc_right: the AbstractFieldBC1D boundary condition on the right boundary

finalize_potential!(poisson_solver, field_props)

Write the values of the potential which are not part of the tridiagonal solution into the nodes: the prescribed values for Dirichlet boundary conditions, and, in the periodic case, the pinned value $\phi = 0$ in node n_cells (which has to be written every timestep, as a stale value left over from the previous timestep would corrupt both the gauge and the electric field in the neighbourhood of that node), followed by the gauge shift and the mirroring of node 1 into the last node.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the potential

electric_field_left!(poisson_solver, field_props, bc_left)

Compute the x-component of the electric field in the left boundary node: a one-sided difference of the potential for a Dirichlet boundary condition, and the prescribed value for a Neumann boundary condition.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the potential and the electric field
• bc_left: the AbstractFieldBC1D boundary condition on the left boundary

electric_field_right!(poisson_solver, field_props, bc_right)

Compute the x-component of the electric field in the right boundary node: a one-sided difference of the potential for a Dirichlet boundary condition, and the prescribed value for a Neumann boundary condition.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the potential and the electric field
• bc_right: the AbstractFieldBC1D boundary condition on the right boundary

compute_electric_field!(poisson_solver, field_props)

Compute the x-component of the electric field $E_x = -d\phi/dx$ in the nodes from the potential, using central differences in the interior nodes and the boundary condition-specific treatment of the boundary nodes.

Positional arguments

• poisson_solver: the PoissonSolver1DUniform instance
• field_props: the ElectrostaticFieldProps instance holding the potential and the electric field

AbstractBC

Abstract type for boundary conditions

specular_reflection_x!(particle, Δt)

Perform specular reflection of a particle in the x direction.

Positional arguments

• particle: the Particle instance for which the velocity is reflected
• Δt: the time left for the particle to move

Returns

Unchanged value of Δt.

diffuse_reflection_x!(rng, particle, wall_reflection_v_sq, wall_normal_sign, wall_v, Δt)

Perform diffuse reflection of a particle, assuming the wall is orthogonal to the x axis.

Positional arguments

• rng: the random number generator
• particle: the Particle instance for which the velocity is reflected
• wall_reflection_v_sq: the squared thermal velocity of the species reflected at the wall temperature
• wall_normal_sign: sign of the wall normal
• wall_v: wall velocity vector
• Δt: the time left for the particle to move

Returns

Unchanged value of Δt.

update_surface_incident!(particle, species, surf_props, surface_element_id)

Update surface properties for surface element surface_element_id for an incident particle.

Positional arguments

• particle: the particle hitting the surface before its velocity is updated
• species: the species of the particle
• surf_props: the SurfProps instance
• surface_element_id: id of the surface element which the particle has impacted

update_surface_reflected!(particle, species, surf_props, surface_element_id)

Update surface properties for surface element surface_element_id for a reflected particle.

Positional arguments

• particle: the particle hitting the surface after its velocity is updated
• species: the species of the particle
• surf_props: the SurfProps instance
• surface_element_id: id of the surface element which the particle has impacted

surface_props_scale!(species, species_data, surf_props, Δt)

Scale computed surface properties using the molecular mass of species, the inverse of the timestep, and the inverse surface area.

Positional arguments

• species: the species of the particle
• species_data: the vector of SpeciesData of the chemical species in the flow
• surf_props: the SurfProps instance
• Δt: the timestep over which the surface properties were computed

apply_bc!(rng, particle::Particle{D}, wallbc::MaxwellWallBC1D, surface_normal::SVector{3,Float64}, Δt) where D

Apply a Maxwell 1D boundary condition.

Positional arguments

• rng: the random number generator
• particle: the Particle instance for which the velocity is reflected
• wallbc: the MaxwellWallBC1D instance
• surface_normal: the vector of the surface normal
• Δt: the time left for the particle to move

Returns

Unchanged value of Δt.

apply_bc!(rng, particle::Particle{D}, wallbc::FullyDiffuseBC1D, surface_normal::SVector{3,Float64}, Δt) where D

Apply a fully diffuse 1D boundary condition.

Positional arguments

• rng: the random number generator
• particle: the Particle instance to which the boundary condition is applied
• wallbc: the FullyDiffuseBC1D instance
• surface_normal: the vector of the surface normal
• Δt: the time left for the particle to move

Returns

Unchanged value of Δt.

apply_bc!(rng, particle::Particle{D}, wallbc::FullySpecularBC1D, surface_normal::SVector{3,Float64}, Δt) where D

Apply a fully specular 1D boundary condition.

Positional arguments

• rng: the random number generator
• particle: the Particle instance to which the boundary condition is applied
• wallbc: the FullySpecularBC1D instance
• surface_normal: the vector of the surface normal
• Δt: the time left for the particle to move

Returns

Unchanged value of Δt.

apply_bc_dispatched!(rng, particle::Particle{D}, bcs::Tuple, i::Int, surface_normal::SVector{3,Float64}, Δt::Float64) where D

Generates a function that dispatches to the correct boundary condition function based on the index i. For use in a particle convection routine.

Positional arguments:

• rng: the random number generator
• particle: the Particle instance for which the boundary condition is applied
• bcs: Tuple of boundary conditions
• i: index of boundary condition to apply
• surface_normal: the vector of the surface normal
• Δt: the time left for the particle to move

Returns

Potentially changed value of Δt.

AbstractNCDataHolder

Abstract type that holds NetCDF-output related data for I/O

swap_particles_true_index!(pv1::ParticleVector{D}, pv2::ParticleVector{D}, i, j) where D

Swap particles pv1.particles[i] and pv2.particles[j] in two ParticleVector instances. This does not update any associated indices or buffers. This uses the underlying ("true") indices of the particles.

Positional arguments

• pv1: the first ParticleVector
• pv2: the second ParticleVector
• i: the underlying index of the particle in pv1
• j: the underlying index of the particle in pv2

swap_particles!(pv1::ParticleVector{D}, pv2::ParticleVector{D}, i, j) where D

Swap particles pv1[i] and pv2[j] in two ParticleVector instances. This does not update any associated indices or buffers.

Positional arguments

• pv1: the first ParticleVector
• pv2: the second ParticleVector
• i: index of the particle in pv1
• j: index of the particle in pv2

update_swap_indexing!(chunk_exchanger, pia_chunks, species, i, j, s_ci_ij, e_ci_ij, s_ji, n_swap)

Updates index bookkeeping in chunk_exchanger for particles that were swapped from chunk i to chunk j.

After particles are swapped, this function records where in chunk j those particles from chunk i were placed, cell by cell. It does not perform the actual swapping. For correct bookkeeping it thus needs to be called twice, with all arguments dependent on i and j symmetrically swapped.

Positional arguments

• chunk_exchanger: the ChunkExchanger instance to track post-swap and post-push indices
• pia_chunks: Vector of ParticleIndexerArray instances for each chunk
• species: the particle species being redistributed
• i: source chunk index
• j: destination chunk index
• s_ci_ij: start cell index in chunk i for the swap
• e_ci_ij: end cell index in chunk i for the swap
• s_ji: starting index in chunk j's particle array where swapped particles from chunk i were placed
• n_swap: total number of particles that were successfully swapped

Returns

• s_ci_ij2: the last cell index in the iteration (i.e., where the last particle was swapped)
• offset_ij: number of particles swapped from s_ci_ij2, useful if the swap was interrupted mid-cell

push_particles!(chunk_exchanger, particles_chunks::Vector{Vector{ParticleVector{D}}}, pia_chunks, species, i, j, offset_ij, s_ci_ij2, e_ci_ij) where D

Pushes particles of the specified species from chunk i to the end of chunk j, handling partial (unfinished) swaps.

This function finalizes a particle redistribution process between two spatial chunks in a parallel or chunked particle simulation. It ensures all remaining particles from chunk i that need to be moved to chunk j are properly transferred and indexed, even if the previous transfer process via swapping particles was interrupted mid-cell. This also updates the buffer of the source ParticleVector as particles are removed from it.

Positional arguments

• chunk_exchanger: the ChunkExchanger instance to track post-swap and post-push indices
• particles_chunks: Vector of Vector of ParticleVector (per chunk and per species, i.e. particles_chunks[chunk_id][species] is the correct order of access)
• pia_chunks: Vector of ParticleIndexerArray instances for each chunk
• species: the particle species being redistributed
• i: Source chunk index
• j: Destination chunk index
• offset_ij: number of particles already transferred from a cell where swapping was performed but incomplete (0 if no incomplete swapping)
• s_ci_ij2: cell index in chunk i where the swapping was interrupted, or where the push should start in case no swapping was performed
• e_ci_ij: final cell index to transfer from chunk i to chunk j

push_particles_to_end!(pv_i::ParticleVector{D}, pv_j::ParticleVector{D}, s2, e2, offset) where D

Copy the particles pv_i[s2+offset:e2+offset] to the positions s2:e2 of pv_j, taking new particles in pv_j from its buffer, and adding the freed particles of pv_i to its buffer. This does not update any ParticleIndexer/ParticleIndexerArray indexing.

Positional arguments

• pv_i: the source ParticleVector
• pv_j: the destination ParticleVector
• s2: first index in pv_j to write to
• e2: last index in pv_j to write to
• offset: offset between the indices in pv_j and the indices of the particles in pv_i

Speed of light, m/s

Electron-Volt, K

Electron-Volt, J

1.0/eV[J], 1/J (or equivalent to how much 1 J is equal to expressed in eV)

$2 \pi$

Electron mass divided by 1 eV, kg/J

Positive and negative direction signs

Elementary charge, C

create_position_vector_D(rng, ::Val{D}, xlo, xhi, ylo, yhi, zlo, zhi) where D

Helper function to create D-dimensional position vectors for particle sampling in a 0/1/2/3-dimensional box with extent xlo:xhi, ylo:yhi, zlo:zhi; positions are uniformly distributed.

Positional arguments

• rng: the random number generator
• ::Val{D}: the dimension of the position vector
• xlo - lower x-extent of the box
• xhi - upper x-extent of the box
• ylo - lower y-extent of the box
• yhi - upper y-extent of the box
• zlo - lower z-extent of the box
• zhi - upper z-extent of the box

compute_thermal_velocity(m, T)

Compute the thermal velocity $\sqrt(2kT/m)$.

Positional arguments

• m: the molecular mass of the species
• T: the temperature

Returns

The thermal velocity

binary_search(x, val)

Binary search for value val in a sorted array x. Finds the position mid such that x[mid] < val < x[mid+1].

Positional arguments

• x: the sorted array to be searched
• val: the value to search for

Returns

• -1 if val < x[1]
• 0 if val > x[end]
• Otherwise, returns the index mid satisfying x[mid] < val < x[mid+1]

linear_interpolation(x, y, val, pos, lower_limit, upper_limit)

Perform linear interpolation of a tabulated function y(x), given a parameter value val, a sorted array of values of x and corresponding values of y, as well as the index of the closest values of x to val. In case val is outside of the limits of x, return placeholder values. That is, given that x[pos] < val < x[pos+1], we want to interpolate y(val) given that y[x[pos]] = y[pos], y[x[pos+1]] = y[pos+1].

Positional arguments

• x: the vector of values of the parameter of the function to be interpolated
• y: the vector of the corresponding function values
• val: the value of the parameter at which the function is to be interpolated
• pos: the index of the element in x such that x[pos] < val < x[pos+1], or -1 if val < x[1], 0 if val > x[end]
• lower_limit: the value to return if val < x[1]
• upper_limit: the value to return if val > x[end]

Returns

• The linearly interpolated value of y(val) if x[1] <= val <= x[end]
• lower_limit, if val < x[1]
• upper_limit, if val > x[end]

compute_mixed_moment(particles, pia, cell, species, powers; sum_scaler=1.0, res_scaler=1.0)

Compute mixed velocity moment of particles in a cell: $\sum_i w_i v_{x,i}^{p_x} v_{y,i}^{p_y} v_{z,i}^{p_z}$.

Positional arguments

• particles: the Vector of ParticleVectors containing all the particles in a simulation
• pia: the ParticleIndexerArray instance
• cell: the cell in which the moment is being computed
• species: the species for which the mixed moment is being computed
• powers: a Vector of the powers to which the x-, y-, and z-components of the velocity are to be raised

Keyword arguments

• sum_scaler: during the summation over the particles, scale each summand by this factor potentially reduce round-off issues
• res_scaler: scaling factor by which to multiply the result at the end

scale_columns!(matrix, column_norms)

Scale the columns of the matrix to have unit L2 norm and store the inverse of the original norm in a vector.

Positional arguments

• matrix: the matrix of the LHS
• column_norms: the vector in which to store the computed inverses of the original column-wise norms

Algorithm NNLS: NONNEGATIVE LEAST SQUARES

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 15, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM A * X = B SUBJECT TO X .GE. 0

CONSTRUCTION AND/OR APPLICATION OF A SINGLE HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

UnsafeVectorView only works for isbitstype types. For other types, we're already allocating lots of memory elsewhere, so creating a new View is fine.

This function looks type-unstable, but the isbitstype(T) test can be evaluated by the compiler, so the result is actually type-stable.

Fallback for non-contiguous arrays, for which UnsafeVectorView does not make sense.

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 15, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

Views in Julia still allocate some memory (since they need to keep a reference to the original array). This type allocates no memory and does no bounds checking. Use it with caution.

COMPUTE ORTHOGONAL ROTATION MATRIX.. The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A2+B2)) (-S,C) (-S,C)(B) ( 0 ) COMPUTE SIG = SQRT(A2+B2) SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT SIG MAY BE IN THE SAME LOCATION AS A OR B .

CONSTRUCTION AND/OR APPLICATION OF A SINGLE HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

apply_householder_sweep!(A::Matrix{T}, up::T, ju, nsetp, m, idx, iz1, iz2) where T

Apply the Householder transformation held in rows nsetp:m of column ju of A to rows nsetp:m of every column in set Z, i.e. the columns idx[iz1], ..., idx[iz2].

This is a batched replacement for calling apply_householder! once per column. It is equivalent to that loop, but the reflector scalars (the up * u[1] sign test and the reciprocal 1 / (up * u[1])) are the same for every column, so they are computed once outside the column loop instead of once per column, and the columns are processed four at a time so that the reflector stays in registers and the four dot products overlap.

Positional arguments

• A: the matrix being triangularized
• up: the first component of the Householder vector, as returned by construct_householder!
• ju: index of the column holding the Householder vector
• nsetp: number of columns already moved to set P (the transformation acts on rows nsetp:m)
• m: number of rows of A
• idx: the index vector of the NNLS workspace
• iz1, iz2: the range of positions in idx holding the set-Z columns