merzbild · knstmrd · Apr 9, 2026 · Feb 6, 2026 · Feb 6, 2026 · Feb 9, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,17 @@
 # Changelog
 
+## v0.7.9
+* `restore_particle_ordering!` added, this restores optimal indexing of particles and can lead to simulation speed-ups due to improved cache usage
+* Minor optimizations in octree merging
+* Simplified elastic VHS collision code via a unified `collide_2particles_vhs!` function, minor speed-ups
+* New keyword parameter `dw_tol` in collision routines that sets tolerance in weight difference under which particle
+collisions are treated as equal weight collisions and no additional particle splitting is performed
+* Minor improvements in readability and speed of electron-neutral collision routines
+* `@inbounds` added to acceleration routine
+* `ntc_equal_weight` added for slightly faster collisions in equal-weight simulations
+* Documentation improvements
+* Improved test coverage
+
 ## v0.7.8
 * Fixed ionization simulations without event splitting
 * Tests for simulations with elastic and ionizing electron-neutral collisions (without event splitting)

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Merzbild"
 uuid = "de01fcb4-c117-45a4-a951-7e39b0f12516"
 authors = ["Georgii Oblapenko <kunstmord@kunstmord.com>", "Leo Basov <basov.leo@gmail.com>"]
-version = "0.7.8"
+version = "0.7.9"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

diff --git a/docs/src/contiguous_indexing.md b/docs/src/contiguous_indexing.md
@@ -106,6 +106,11 @@ the value of `pia.contiguous` and does nothing if that value is set to `true`.
 **Summary**: one needs to restore continuity of particle indexing if particles are deleted and created in a simulation, otherwise this
 might lead to erroneous results.
 
+## Optimal particle indexing
+Over time, particle indexing becomes quite fragmented, leading to poor cache performance. Especially in non-0D simulations this affects the performance
+of sorting and convection routines. This can be fixed by calling [`restore_particle_ordering!`](@ref) **after** the particles have been sorted.
+Restoring particle indexing every 10 timesteps seems to be a good balance between the cost of the re-indexing operation and the performance gain.
+
 ## Debugging
 Several utility functions are available to verify/help debug simulations.
 

diff --git a/docs/src/overview_1d.md b/docs/src/overview_1d.md
@@ -67,7 +67,17 @@ sort_particles!(gridsorter, grid, particles[species_id], pia, species_id)
 Since sorting indices only can lead to increase fragmentation of the particle layout in memory,
 a function [`count_disordered_particles`](@ref) is available that counts the number of 
 non-continuously laid out particles; this can serve as a metric as to whether the underlying
-particles (and not just their indices) need to be re-sorted (currently not implemented).
+particles (and not just their indices) need to be re-sorted.
+The fragmentation of particles indices can be fixed by calling [`restore_particle_ordering!`](@ref) after the sorting routine;
+calling it every 10 timesteps or so gives a good balance between cost of re-indexing and speed-up due to improved memory access.
+
+The re-indexing can be called as follows:
+```julia
+restore_particle_ordering!(particles[species_id], index_inv_map)
+```
+
+Here `index_inv_map` is a pre-allocated array of integers of length equal to the number of particles in the simulation (if it is smaller,
+it will be resized in by the `restore_particle_ordering!` function).
 
 ## Creating boundary conditions
 Next, we need to create boundary conditions for the left and right walls.
@@ -169,6 +179,7 @@ n_particles = ppc * nx
 particles = [ParticleVector(n_particles)]
 pia = ParticleIndexerArray(grid.n_cells, 1)
 gridsorter = GridSortInPlace(grid, n_particles)
+index_inv_map = zeros(Int64, n_particles)
 
 # sample particles
 # Fnum * ppc = Np in cell = ndens * V_cell
@@ -227,6 +238,11 @@ for t in 1:n_timesteps
     # sort particles
     sort_particles!(gridsorter, grid, particles[1], pia, 1)
 
+    # restore indexing/ordering every 10 timesteps
+    if t%10 == 0
+        restore_particle_ordering!(particles[1], index_inv_map)
+    end
+
     # compute props and do averaging
     if (t >= avg_start)
         compute_props_sorted!(particles, pia, species_data, phys_props)

diff --git a/docs/src/overview_blocks.md b/docs/src/overview_blocks.md
@@ -112,6 +112,12 @@ For all practical purposes, it should be sufficient to deal with particles in a
 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!`](@ref),
+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](@ref "Contiguous indexing"), which describes additional

diff --git a/docs/src/overview_fixedweight.md b/docs/src/overview_fixedweight.md
@@ -33,11 +33,11 @@ passed to the collision routine, which will use it as required.
 
 ## NTC-specific data: CollisionFactors
 The No-Time Counter (NTC) collision algorithm requires an estimate of ``(\sigma g w)_{max}`` for each species pair
-for each cell in the flow. The `CollisionFactors` data structure stores the values of this factor, as well as 
+for each cell in the flow. The [`CollisionFactors`](@ref) data structure stores the values of this factor, as well as 
 other NTC-related parameters (number of collisions performed, number of collision partners).
 One needs to initialize a 3-dimensional array of `CollisionFactors` of shape 
 `n_species x n_species x n_cells` in order to store the required factors for each species' pair for all cells in the flow;
-this is done by calling `create_collision_factors_array(n_species, n_cells)`.
+this is done by calling [`create_collision_factors_array(n_species, n_cells)`](@ref).
 Some initial estimate for the values of ``(\sigma g w)_{max}``; the simplest way to do so for a fixed-weight DSMC simulation
 is to call the `estimate_sigma_g_w_max!(collision_factors, interactions, species_data, T_list, Fnum; mult_factor=1.0)`
 function.
@@ -61,8 +61,11 @@ Now that the interaction data has been loaded, the `CollisionData` instance to s
 has been instantiated, the 3-dimensional array of `CollisionFactors` has been created, and the ``(\sigma g w)_{max}``
 values precomputed, one can perform collisions.
 
-Single-species elastic VHS collisions can be performed by calling
-`ntc!(rng, collision_factors, collision_data, interaction, particles, pia, cell, species, Δt, V)`.
+Single-species elastic VHS collisions for fixed-weight particles can be performed by calling
+[`ntc_equal_weight!`](@ref).
+A more generic function [`ntc!`](@ref)
+is available, which checks whether particles have equal or non-equal weights (see [Variable-weight DSMC simulations](@ref)); for
+purely equal-weight simulations it is however somewhat slower, thus it is recommended to use [`ntc_equal_weight!`](@ref).
 Here `collision_factors` is the specific instance of `CollisionFactors`, i.e. a specific element of the
 3-dimensional array of `CollisionFactors`. `Δt` is the timestep, and `V` is the volume of the physical cell
 (for spatially homogeneous simulations this can be set to 1.0).
@@ -145,11 +148,11 @@ for ts in 1:n_t  # loop over time
         for s1 in s2:n_species  # loop over second species
             if (s1 == s2)
                 # collisions between particles of same species
-                ntc!(rng, collision_factors[s1,s1,1], collision_data,
+                ntc_equal_weight!(rng, collision_factors[s1,s1,1], collision_data,
                      interaction_data, particles[s1], pia, 1, s1, Δt, V)
             else
                 # collisions between particles of different species
-                ntc!(rng, collision_factors[s1,s2,1], collision_data,
+                ntc_equal_weight!(rng, collision_factors[s1,s2,1], collision_data,
                      interaction_data, particles[s1], particles[s2],
                      pia, 1, s1, s2, Δt, V)
             end

diff --git a/docs/src/overview_varweight.md b/docs/src/overview_varweight.md
@@ -71,10 +71,14 @@ pia = ParticleIndexerArray(n_sampled)
 ```
 
 ## Colliding variable-weight particles
-As mentioned above, the existing DSMC routines take care of the required particle splitting, so no specific adaptation is required
-for the variable-weight case. However, one needs to estimate ``(\sigma g w)_{max}``. The simplest approach is to
+As mentioned above, the existing DSMC routines, unless they have `equal_weight` explicitly in their name,
+take care of the required particle splitting. For variable-weight NTC simulations, one needs to estimate ``(\sigma g w)_{max}``.
+The simplest approach is to
 use `estimate_sigma_g_w_max` function, but it requires a (fixed) value of `Fnum` - so one can simply compute it as
-`ndens/n_sampled` (i.e. the number density divided by the total number of sampled particles).
+`ndens/n_sampled` (i.e. the number density divided by the total number of sampled particles). If we merge after sampling,
+an effective `Fnum` is then computed as `ndens/n_post_merge` (i.e. number density divided by the total number of post-merge particles).
+The [`ntc!`](@ref) function has a keyword parameter `dw_tol` that governs when particle collisions are treated as equal-weight so that small
+discrepancies in particle weights do not cause unnecessary splitting.
 
 ## Merging variable-weight particles
 In order to merge variable-weight particles, one needs to set up a merging algorithm and the associated

diff --git a/docs/src/reference_internal.md b/docs/src/reference_internal.md
@@ -39,6 +39,8 @@ Merzbild.sample_maxwellian!
 ```@docs
 Merzbild.compute_n_coll_single_species
 Merzbild.compute_n_coll_two_species
+Merzbild.collide_2particles_vhs!
+Merzbild.collide_2particles_vhs_equal_weight!
 Merzbild.compute_vhs_factor
 Merzbild.compute_com!
 Merzbild.compute_g!
@@ -149,6 +151,7 @@ Merzbild.AbstractNCDataHolder
 
 ## Parallel computations
 ```@docs
+Merzbild.swap_particles_true_index!
 Merzbild.swap_particles!
 Merzbild.update_swap_indexing!
 Merzbild.push_particles!

diff --git a/docs/src/reference_public.md b/docs/src/reference_public.md
@@ -23,6 +23,7 @@ ParticleIndexerArray(n_particles::Integer)
 ParticleIndexerArray(n_particles::T) where T<:AbstractVector
 ParticleIndexerArray(grid, species_data::Array{Species}) 
 squash_pia!
+restore_particle_ordering!
 count_disordered_particles
 check_unique_index
 check_pia_is_correct
@@ -97,6 +98,7 @@ estimate_sigma_g_w_max!
 estimate_sigma_g_w_max_ntc_n_e!
 estimate_sigma_g_max!
 ntc!
+ntc_equal_weight!
 ntc_n_e!
 ntc_n_e_es!
 swpm!

diff --git a/scripts/README.md b/scripts/README.md
@@ -32,4 +32,4 @@ python scripts/plot_1D.py --file scratch/data/couette_0.0005_50_500.0_300.0_1000
 ```
 
 ## Paper reproducibility scripts
-Scripts to post-process simulation results for various papers can be found in the `reproducibility` subdirectory.
+Scripts to post-process simulation results for various papers can be found in the `reproducibility` subdirectory.
diff --git a/simulations/0D/BKW/bkw.jl b/simulations/0D/BKW/bkw.jl
@@ -3,7 +3,6 @@
 # 20k particles, julia --project=. simulations/basic/bkw.jl  4.73s user 1.75s system 133% cpu 4.874 total
 # 50k particles, no print julia --project=. simulations/basic/bkw.jl  5.54s user 1.65s system 134% cpu 5.357 total
 # 200k particles, 10 moments, no print julia --project=. simulations/basic/bkw.jl  13.29s user 1.78s system 114% cpu 13.211 total
-
 using Merzbild
 using Random
 
@@ -77,7 +76,7 @@ function run(seed)
     V::Float64 = 1.0
 
     @time for ts in 1:n_t
-        ntc!(rng, collision_factors, collision_data, interaction_data, particles[1],
+        ntc_equal_weight!(rng, collision_factors, collision_data, interaction_data, particles[1],
              pia, 1, 1, Δt, V)
 
         compute_props!(particles, pia, species_data, phys_props)

diff --git a/simulations/0D/basic/basic_collisions_multispecies.jl b/simulations/0D/basic/basic_collisions_multispecies.jl
@@ -60,9 +60,9 @@ function run(seed)
         for s2 in 1:n_species
             for s1 in s2:n_species
                 if (s1 == s2)
-                    ntc!(rng, collision_factors[s1,s1,1], collision_data, interaction_data, particles[s1], pia, 1, s1, Δt, V)
+                    ntc_equal_weight!(rng, collision_factors[s1,s1,1], collision_data, interaction_data, particles[s1], pia, 1, s1, Δt, V)
                 else
-                    ntc!(rng, collision_factors[s1,s2,1], collision_data, interaction_data, particles[s1], particles[s2],
+                    ntc_equal_weight!(rng, collision_factors[s1,s2,1], collision_data, interaction_data, particles[s1], particles[s2],
                          pia, 1, s1, s2, Δt, V)
                 end
             end

diff --git a/simulations/1D/couette_benchmarking.jl b/simulations/1D/couette_benchmarking.jl
@@ -53,14 +53,16 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, n_timesteps, avg_star
 
     n_avg = n_timesteps - avg_start + 1
 
+    index_inv_map = zeros(Int64, n_particles)
+
     for t in 1:n_timesteps
         if t % 1000 == 0
             println(t)
         end
 
         # collide particles
         for cell in 1:grid.n_cells
-            @timeit "collide" @inbounds ntc!(rng, collision_factors[1, 1, cell],
+            @timeit "collide" @inbounds ntc_equal_weight!(rng, collision_factors[1, 1, cell],
                                    collision_data, interaction_data, particles[1], pia, cell, 1, Δt, grid.cells[cell].V)
         end
 
@@ -72,6 +74,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, n_timesteps, avg_star
 
         # compute props and do I/O
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
+
         if (t >= avg_start)
             @timeit "props compute" compute_props_sorted!(particles, pia, species_data, phys_props)
             @timeit "avg physprops" avg_props!(phys_props_avg, phys_props, n_avg)

diff --git a/simulations/1D/couette_fp.jl b/simulations/1D/couette_fp.jl
@@ -54,6 +54,8 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, output_freq, n_timest
     write_grid("scratch/data/couette_$(L)_$(nx)_grid.nc", grid)
 
     n_avg = n_timesteps - avg_start + 1
+
+    index_inv_map = zeros(Int64, n_particles)
 
     for t in 1:n_timesteps
         if t % 1000 == 0
@@ -71,6 +73,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, output_freq, n_timest
         # sort particles
         @timeit "sort" sort_particles!(gridsorter, grid, particles[1], pia, 1)
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
+
         # compute props and do I/O
 
         if (t < avg_start)

diff --git a/simulations/1D/couette_fp_vw.jl b/simulations/1D/couette_fp_vw.jl
@@ -68,6 +68,8 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, merge_threshold, merge_tar
 
     n_avg = n_timesteps - avg_start + 1
 
+    index_inv_map = zeros(Int64, n_particles)
+
     for t in 1:n_timesteps
         if t % 1000 == 0
             println(t)
@@ -92,6 +94,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, merge_threshold, merge_tar
         # sort particles
         @timeit "sort" sort_particles!(gridsorter, grid, particles[1], pia, 1)
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
+
         # compute props and do I/O
 
         if (t < avg_start)

diff --git a/simulations/1D/couette_varweight_nnls.jl b/simulations/1D/couette_varweight_nnls.jl
@@ -118,6 +118,8 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
 
     n_avg = n_timesteps - avg_start + 1
 
+    index_inv_map = zeros(Int64, n_particles)
+
     for t in 1:n_timesteps
         if t % 1000 == 0
             println("$t, # of particles=$(pia.n_total[1])")
@@ -152,6 +154,9 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
         # sort particles
         @timeit "sort" sort_particles!(gridsorter, grid, particles[1], pia, 1)
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
 
         # compute props and do I/O
         if (t < avg_start)

diff --git a/simulations/1D/couette_varweight_octree.jl b/simulations/1D/couette_varweight_octree.jl
@@ -79,6 +79,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
 
     n_avg = n_timesteps - avg_start + 1
 
+    index_inv_map = zeros(Int64, n_particles)
+
+    index_inv_map = zeros(Int64, n_particles)
+
     for t in 1:n_timesteps
         if t % 1000 == 0
             println(t)
@@ -107,7 +111,11 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
 
         # count % of particles where indexing is disordered
         if t % 1000 == 0
-            println(count_disordered_particles(particles[1], pia, 1) / pia.n_total[1] * 100.0)
+            @timeit "disordered count" println(count_disordered_particles(particles[1], pia, 1) / pia.n_total[1] * 100.0)
+        end
+
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
         end
 
         # compute props and do I/O

diff --git a/simulations/1D/couette_varweight_octree_swpm.jl b/simulations/1D/couette_varweight_octree_swpm.jl
@@ -77,6 +77,8 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
 
     n_avg = n_timesteps - avg_start + 1
 
+    index_inv_map = zeros(Int64, n_particles)
+
     for t in 1:n_timesteps
         if t % 1000 == 0
             println(t)
@@ -108,6 +110,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc_sampled, merge_threshold, m
             println(count_disordered_particles(particles[1], pia, 1) / pia.n_total[1] * 100.0)
         end
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
+
         # compute props and do I/O
         if (t < avg_start)
             if (t % output_freq == 0)

diff --git a/simulations/1D/couette_with_surface_quantities.jl b/simulations/1D/couette_with_surface_quantities.jl
@@ -66,6 +66,8 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, output_freq, n_timest
     write_netcdf(ds, phys_props, 0)
     write_grid("scratch/data/couette_$(L)_$(nx)_grid.nc", grid)
 
+    index_inv_map = zeros(Int64, n_particles)
+
     n_avg = n_timesteps - avg_start + 1
 
     for t in 1:n_timesteps
@@ -98,6 +100,10 @@ function run(seed, T_wall, v_wall, L, ndens, nx, ppc, Δt, output_freq, n_timest
             end
         end
 
+        if t%10 == 0
+            @timeit "restore ordering" restore_particle_ordering!(particles[1], index_inv_map)
+        end
+
         # compute props and do I/O
 
         if (t < avg_start)