Overview of basic building blocks

Particle indexing: ParticleIndexer and ParticleIndexerArray

One of the main parts of the code is the ParticleIndexer struct. It contains the starting and ending indices of the particles of a single species in a single cell. During collisions, the number of the particles of a certain species in a cell may increase, either due to inelastic processes, or due to particle splitting in variable-weight collisions. These new particles are "tacked onto" the end of the whole array of particles of the species. This two-group design allows efficient handling of newly created particles during collisions without disturbing the ordered indices of existing particles Thus, each ParticleIndexer instance potentially tracks two blocks of particles that are all in a single cell: the ones that were there before collisions were performed, and the new ones that got created during collisions and are at the end of the array.

particle_indexer_trim

The figure above shows a graphical representation of the meaning of some of the fields of the ParticleIndexer struct. Assume we have a set of 10 particles in 3 different cells (different cells are colour-coded), with particles 1, 2, 3, 7, 8, 9, 10 located in cell 1, particles 4, 5 located in cell 2, and particle 6 located in cell 3. The ParticleIndexer struct corresponding to cell 1 will then have the following values for its fields:

particle_indexer_cell_1.start1 = 1particle_indexer_cell_1.end1 = 3particle_indexer_cell_1.start2 = 7particle_indexer_cell_1.end2 = 10

The struct also has additional utility fields:

  • ParticleIndexer.n_group1 = ParticleIndexer.end1 - ParticleIndexer.start1 + 1 stores the number of particles in the first block
  • ParticleIndexer.n_group2 = ParticleIndexer.end2 - ParticleIndexer.start2 + 1 stores the number of particles in the second block
  • ParticleIndexer.n_local = ParticleIndexer.n_group1 + ParticleIndexer.n_group2 stores the total number of particles of this species in this cell.

The fields of the ParticleIndexer struct are not really intended for manual manipulation unless one is developing new sorting or merging routines, in which case one might have to set the starting and ending fields oneself. Most of the time, the existing sorting/collision/merging routines take care of keeping the particle indexers updated and correct.

Since we are usually interested in multi-species and/or multi-dimensional simulations, an array of ParticleIndexer instances is needed to keep track of the particles of different species in different cells. For this purpose Merzbild.jl provides the ParticleIndexerArray struct. It the following fields:

  • ParticleIndexerArray.n_cells: number of grid cells in the simulation
  • ParticleIndexerArray.n_species: number of species in the simulation
  • ParticleIndexerArray.indexer: a 2-dimensional array of ParticleIndexer instances with dimensions n_cells*n_species
  • ParticleIndexerArray.n_total: a 1-dimensional vector of length n_species with per-species total particle counts
  • ParticleIndexerArray.contiguous: a 1-dimensional vector of length n_species describing whether indexing is contiguous for each species (see the section on Particle buffers and contiguous indexing)
  • ParticleIndexerArray.index_last: a 1-dimensional vector of length n_species holding the index of the last particle of each species.

IMPORTANT: the indexing layout always assumes that any particles in any of the index ranges pointed to by a group2 come after all the particles in any of the index ranges pointed to by a group2. So the first index of a particle in any of the group2 ranges is always larger than the last index of any of the particles in any of the group1 ranges.

To iterate over the particles of a specific species in a specific cell, one uses an instance of the ParticleIndexerArray (called pia in the code by convention):

for i in pia.indexer[cell,species].start1:pia.indexer[cell,species].end1    # do something with particles[species][i]endif pia.indexer[cell,species].start2 > 0    for i in pia.indexer[cell,species].start2:pia.indexer[cell,species].end2        # do something with particles[species][i]    endend

Some utility functions are available for updating particle indexer arrays for developing new functionality: Merzbild.update_particle_indexer_new_particle!. It is assumed that if the number of particles in the second block pointed to by a ParticleIndexer instance is 0, then the value of start2 is set to a value smaller than 0, so to iterate over the particles in the second block one can either do a check similar to the pia.indexer[cell,species].start2 > 0 check seen above, or check if pia.indexer[cell,species].n_group2 <= 0. See also the section on Particle buffers and contiguous indexing. The function pretty_print_pia can be used to print the particle indexer array in a human-readable format.

Particles: Particle and ParticleVector

Now that we can index particles, we need to create some lists of particles to index. For that, we need to define what a particle is. For this purpose, a Particle{D} struct is available in the code. It has the following fields:

  • w: the computational weight of the particle (in a fixed-weight DSMC simulation, this is the $F_{num}$ parameter)
  • v: the 3-dimensional velocity vector of the particle
  • x: the D-dimensional position vector of the particle

For spatially homogeneous simulations, D can be set to 0, so that the location of the particles is not tracked at all. Each species has its own list of particles associated with it, so a particles variable in the simulation could have the following the type Vector{Vector{Particle{D}}}. Then particles[species_1] would correspond to the list of all particles of chemical species species_1. pia.indexer[cell,species_1] would then be used to index the particles of species_1 in a specific cell cell. If no D is specified explicitly, the dimension of the position vector defaults to 3.

The drawback of using Vector{Vector{Particle{D}}} is that for non-spatially homogeneous simulations, the particles need to be sorted after each convection step, and this would involving constantly re-writing the position and velocity vectors. To reduce the computational cost of sorting, an additional abstraction layer is added via the struct ParticleVector{D} that is intended to be used instead of a simple Vector{Particle{D}} instance. An instance of ParticleVector{D} has the following fields:

  • particles: the underlying vector of particles with a D-dimensional position.
  • index: the sorted indices of the particles
  • cell: used in particle sorting to keep track of new assigned cells
  • buffer: a LIFO queue used to track which particles from particles are not being used (i.e. allocated in memory but not present in the simulation)
  • nbuffer: the number of elements in buffer

Given a ParticleVector{D} instance pv, one can still transparently access a particle at index i as pv[i]. This access operation however uses the sorted index list to get the actual index of the particle, so pv[i] is equivalent to pv.particles[pv.index[i]]. During particle sorting, only the indices in index are shuffled around, which is computationally cheaper than sorting the particles directly. If no D is specified explicitly, the dimension of the position vector of the particles defaults to 3.

particle_vector_trim

The figure above shows the how the underlying particles (top array) are indexed via the help of the index field of the ParticleVector struct, and finally how the indices of the particles in a specific cell are indexed by the ParticleIndexer instance. In the example above, the different colors correspond to different grid cells. It can be seen that the array of Particle instances ("underlying Particle vector") is unsorted. The indices (ParticleVector.index) are however sorted. Thus, ParticleVector.index[1] points to the 1-st particle in the underlying Particle array, and ParticleVector.index[2] points to the 4-th particle. ParticleVector.index[6] points to the 2-nd particle, and so on. It may be that due to creation of new particles during collisions some particles in the same cell are at the very end of the array (as they just got created). In the example above, particles 9 and 10 are in the same cell as particles 1, 4, 5, 7, 8. Therefore the corresponding ParticleIndexer instance for this cell will have start1=1, end1=5 (pointing to the first block of indices), and start2=9, end2=9 (pointing to the second block of indices).

For all practical purposes, it should be sufficient to deal with particles in a ParticleVector instance pv by directly accessing them as pv[i]; the index and cell fields need to be changed only when one is writing new sorting routines.

Important: over the course of a simulation, as particles get re-sorted, deleted, created, or merged, indexing may become sub-optimal; and particle access becomes cache-inefficient. This can be fixed by calling restore_particle_ordering!, it is up to the user to decide how frequently this should be done to balance the cost of re-indexing and particle access in other routines. Example of such usage can be found in 1D simulations in simulations/1D (for example, in simulations/1D/couette_varweight_octree.jl). Re-ordering indexing every 10 timesteps seems to be a good starting point.

For a more in-depth overview of particle indexing, especially relevant for multi-dimensional simulations where particles might be frequently deleted due to merging or outflow, the reader is referred to the documentation on contiguous indexing, which describes additional concepts associated with the indexing approach used in Merzbild.jl and provides some tips on improving performance in the listed scenarios.

To summarize, for multi-species multi-dimensional simulations, one needs to use a Vector{ParticleVector} instance, each element of which (of type ParticleVector) corresponds to a vector of particles of a specific species. If one does not intend to sort particles (i.e. a 0-D simulation), one can in theory use a simple Vector{Vector{Particle}} instance, however, some of the particle sampling routines assume a ParticleVector instance (as they write the cell index data).

Physical data: Species

To sample particles, we need to know some properties of the chemical species that they represent. This is done via instances of the Species struct, which have the following fields:

  • name: the name of the species
  • mass: the mass of the species
  • charge: the charge of the species (in terms of elementary charge)
  • charge_div_mass: the charge of the species (in Coulombs) divided by the mass of the species

Data about chemical species can be loaded by using the load_species_data function, which reads a TOML file with the relevant species information.

Merzbild.jl exports a MERZBILD_DATA_PATH variable that points to the data directory containing data bundled with the package. For example, one can write

particles_data_path = joinpath(MERZBILD_DATA_PATH, "particles.toml")species_data = load_species_data(particles_data_path, "Ar")

to construct the path to the particles.toml file provided with Merzbild.jl and load data for Argon from the file.

Sampling particles

In order to sample particles of a certain species from a specific distribution, Merzbild.jl provides several functions. For fixed-weight DSMC simulations, one can use the sample_particles_equal_weight! to sample equal-weight particles from either a Maxwell or a BKW distribution. For variable-weight DSMC simulations, one can of course also use the same function (there is no inherent difference between fixed- and variable-weight DSMC simulations in Merzbild).

However, in case variable-weight particles are allowed, one can also sample them from a distribution by evaluating that distribution of a discrete velocity grid. This functionality is provided by the sample_on_grid! function and its more specific version sample_maxwellian_on_grid! that assumes a Maxwell distribution. Sampling on a grid can lead to very high particle counts (for example, sampling on a 30x30x30 grid will lead to creation of up to 27000 particles in each physical cell), so care must be taken, and particles might need to be merged immediately after having been created.

Computing macroscopic physical properties: PhysProps

Now that we have a vector of particles, we can compute some macroscopic properties (density, velocity, etc.). To store and use these properties where they might be needed, the PhysProps struct is provided. An instance of PhysProps has the following fields:

  • ndens_not_Np: a boolean value used to distinguish between the meanings of the n field (see below) and ensure consistency
  • n_cells: number of grid cells
  • n_species: number of species in the simulation
  • lpa: vector of length n_species storing the lengths of the particle arrays (i.e. how many elements have been allocated, actual particle counts may be less)
  • np: array with dimensions n_cells*n_species, stores the number of particles of each species in each grid cell
  • n: array with dimensions n_cells*n_species, stores either the number of physical particles of each species in each grid cell or the number density of each species in each grid cell, see below for explanation
  • v: array with dimensions 3*n_cells*n_species, stores the x, y, and z components of the macroscopic velocity of each species in each grid cell
  • T: array with dimensions n_cells*n_species, stores the temperature of each species in each grid cell

One can see that the definition of the n field is somewhat ambiguous. It can store either: (a) the count of physical particles, or (b) the number density (particles/volume). The ndens_not_Np flag is true when storing number density, false when storing particle counts. To distinguish between these two cases, the following convention is assumed, one can provide a value of ndens_not_NP during instantiation (by default it is false, i.e. the number of physical particles is computed and not the number density). We can simply instantiate a PhysProps instance like this: props = PhysProps(pia).

The compute_props! function computes the macroscopical physical properties of all species in all cells in the simulation. This computes only the number of particles in a cell, regardless of the value of the ndens_not_Np field.

There is an optimized version of this function, which assumes the particles are only indexed by the first group of a ParticleIndexer instance: compute_props_sorted!. This is the case immediately after sorting the particles on a grid. If a grid is passed as a parameter, it will compute either the number of particles in a cell or the number density in a cell depending on the value of ndens_not_NP field of the PhysProps instance passed to the function.

The avg_props! function can also be used to time-average physical properties; this requires use of two PhysProps instances, one of which holds the values of the physical properties at the current timestep, and the other one holds the values of the averaged physical properties. Similarly to the previous case, trying to average one PhysProps instance "into" another PhysProps instance with a different value of the ndens_not_NP field will raise an error.

Writing output: NCDataHolder

Finally, once the properties have been computed, we need to output them. Merzbild.jl uses the binary NetCDF format for its output, specifically the NetCDF4 version of the format. The functionality is provided by the NetCDF.jl and NCDatasets.jl libraries. The testing code relies on NCDatasets.jl as originally Merzbild.jl used the library throughout the whole code; however, the NetCDF.jl library was found to be significantly faster when used for output, so all the non-testing code has been ported to NetCDF.jl. The NCDatasets.jl dependency will be dropped in a future version once tests are ported to NetCDF.jl.

In order to write the physical properties to a NetCDF file, one needs to create an NCDataHolder instance. It keeps track of the variables being output, the file handle, and the number of times data has been written to a file. It also keeps track of whether the ndens output variable actually corresponds to the number of particles in a cell or to the number density in a cell, based on the flag in the PhysProps instance used to initialize the NCDataHolder instance. (In a future version of Merzbild.jl, it is planned to change the output variable name based on the ndens_not_Np flag value in the PhysProps instance used in the constructor.)

We can construct an NCDataHolder instance like this:

NCDataHolder("output_filename.nc", species_data, phys_props)

Additional global attributes can be written to the file if they are passed as a dictionary using the global_attributes keyword of the constructor. Additionally, it is possible to pass a list of variable names that should be skipped during output (i.e. they won't be written to the NetCDF file), a more detailed description is given in the public API reference.

In order to write the physical properties to a file at a certain timestep, we need to call

write_netcdf(ds, phys_props, timestep)

Here ds is the NCDataHolder instance, and timestep is the number of the timestep of the simulation (it also gets written to the output). Additionally, the sync_freq keyword parameter can be set to control how often the file is flushed (synchronized) to disk.

Finally, at the end of the simulation we need to close the file, so we call

close_netcdf(ds)

Computing surface properties due to particle-surface interactions: SurfProps

Details on calculation of surface properties due to particle-surface interactions can be found in the section on 1D DSMC simulations. Time-averaging works exactly the same as for PhysProps, via use of the avg_props! function. Output is performed similarly, via an NCDataHolderSurf struct and calls to write_netcdf.

Example: bringing it all together

An example of particle sampling, property computation, and output for a 0-D single-species gas is presented here.

using Merzbildusing Random# set our random seed for reproducibilityseed = 1Random.seed!(seed)rng = Xoshiro(seed)# load particle dataparticles_data_path = joinpath(MERZBILD_DATA_PATH, "particles.toml")species_data = load_species_data(particles_data_path, "Ar")# init particle vector for a 1000 particlesn_particles = 1000particles = [ParticleVector{0}(n_particles)]# set number density to 1e23ndens = 1e23# and compute the particle weight FnumFnum = ndens / n_particles# set particle temperatureT = 500.0# create the particle indexer array for 1 cell and 1 species, assume 0 particles since we haven't# sampled any particles yetpia = ParticleIndexerArray(0)# for D=3 this would sample particles in a [0.0, 1.0]x[0.0, 1.0]x[0.0, 1.0] cell; for 0-D particles the position is not trackedsample_particles_equal_weight!(rng, particles[1], pia, 1, 1, n_particles, species_data[1].mass, T, Fnum,                               0.0, 1.0, 0.0, 1.0, 0.0, 1.0; distribution=:Maxwellian)# create struct for computation of physical propertiesphys_props = PhysProps(pia)# create struct for netCDF outputds = NCDataHolder("output.nc", species_data, phys_props)# compute physical properties compute_props!(particles, pia, species_data, phys_props)# and output them (t=0)write_netcdf(ds, phys_props, 0)close_netcdf(ds)

Summary

Now we have an overview of how to

  1. Create a structure to hold particles
  2. Index the vector of particles
  3. Load species' data
  4. Sample particles from a distribution
  5. Compute macroscopic physical properties
  6. Output these properties to disk

In the next section, an overview of how to simulate collisions will be given.