merzbild · knstmrd · May 19, 2026 · Apr 9, 2026 · Apr 9, 2026 · May 9, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,12 @@
 # Changelog
 
+## v0.7.10
+* Documentation improvements
+* Improved test coverage
+* `convect_particles_and_compute_cell!` function added that computes new cell indices during the convection, a new version of `sort_particles!`
+that takes advantage of these pre-computed cell indices has also been added
+* `squash_pia!` now also squashes particle cell indices
+
 ## 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

diff --git a/PAPER_REPRODUCIBILITY.md b/PAPER_REPRODUCIBILITY.md
@@ -32,7 +32,7 @@ For the BKW test case (see the commented-out part on "multiple runs with ensembl
 * `simulations/0D/BKW/bkw_varweight_nnls.jl` - for the variable-weight simulations using NNLS merging
 
 For the ionization test case (**when run over all ensembles, these produce VERY LARGE amounts of data, 100s of GBs**);
-external data from [LXCat](https://us.lxcat.net/home/) is required for the cross-sections (IST-Lisbon database) (see notes above on exact XML format required):
+external data from [LXCat](https://us.lxcat.net/home/) is required for the cross-sections (IST-Lisbon database) (see notes below on exact XML format required):
 * `simulations/0D/ionization/0D_ionization_1neutralspecies_es.jl` - for the variable-weight simulations using octree N:2 merging
 (uncomment lines below comment "Uncomment set-up below ..." to get the full set-up running over all parameters and ensembles; setting `paramset` to [12000, 6000] will produce results used as reference values)
 * `simulations/0D/ionization/0D_ionization_1neutralspecies_nnls_es.jl` - for the variable-weight simulations using different versions of NNLS merging

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.9"
+version = "0.7.10"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

diff --git a/data/particles.toml b/data/particles.toml
@@ -16,4 +16,12 @@ charge = 1.0
 
 ["e-"]
 mass = 9.109383632e-31
-charge = -1.0
+charge = -1.0
+
+["Xe"]
+mass = 2.181e-25
+charge = 0.0
+
+["Xe+"]
+mass = 2.181e-25
+charge = 1.0
diff --git a/data/test_neutral_electron_data.xml b/data/test_neutral_electron_data.xml
@@ -103,7 +103,7 @@
     </Database>
 
 
-    <Database name="Constant elastic and ionizaton cross-sections, Ar" id="ConstantDB">
+    <Database name="Constant elastic and ionization cross-sections, Ar" id="ConstantDB">
         <Permlink>
         knstmrd.github.io
         </Permlink>
@@ -194,4 +194,95 @@
         </Groups>
     </Database>
 
+    <Database name="constant value elastic and zero-value ionization cross-sections, Xe" id="ZeroDB">
+        <Permlink>
+        knstmrd.github.io
+        </Permlink>
+        <Description>
+        These data are constant (sigma=1e-19 for elastic cross-section, sigma=0.0 for ionization cross-section).
+        </Description>
+
+        <Contact>
+        email georgii.oblapenko@fastmail.com
+        </Contact>
+
+        <HowToReference>
+        Merzbild.jl reference
+        </HowToReference>
+        <MerzbildUrl>
+        https://github.com/merzbild/Merzbild.jl
+        </MerzbildUrl>
+
+        <Groups>
+            <Group id="Xe">
+                <Description>
+                Const E-sigma relationship</Description>
+                <Processes>
+                    <Process class="Scattering Cross Sections" type="Elastic">
+                        <Species>
+                            <Reactant>e</Reactant><Reactant>He</Reactant>
+                            <Product>E</Product><Product>He</Product>
+                        </Species>
+
+                        <Reaction>
+                        E + Xe -&gt; E + Xe
+                        </Reaction>
+
+                        <Parameters>
+                            <mM>
+                            1.37396e-5
+                            </mM>
+                            <Parameter>
+                            complete set
+                            </Parameter>
+                        </Parameters>
+                        <Comment>
+                        Constant, sigma = 1e-19
+                        </Comment>
+                        <DataX type="Energy" units="eV" size="4">
+                        0.000000e+0 1.0e2 5.0e2 1.0e3
+                        </DataX>
+                        <DataY type="Cross section" units="m2" size="4">
+                        1.0e-19 1.0e-19 1.0e-19 1.0e-19
+                        </DataY>
+                    </Process>
+
+                    <Process class="Scattering Cross Sections" type="Ionization">
+                    <Species>
+                        <Reactant>e</Reactant><Reactant>Ar</Reactant>
+                        <Product>E</Product><Product>E</Product><Product>Ar^+</Product>
+                    </Species>
+
+                    <Reaction>
+                    E + Xe -&gt; E + E + Xe^+
+                    </Reaction>
+
+                    <Parameters>
+                        <E units="eV">
+                        1.212984e+1
+                        </E>
+                        <Parameter>
+                        complete set
+                        </Parameter>
+                    </Parameters>
+                    <Comment>
+                        Constant
+                    </Comment>
+                    <Updated>
+                    2024-10-17
+                    </Updated>
+
+                    <DataX type="Energy" units="eV" size="3">
+                    1.212984e+1 3.0e+1 1.0e+2
+                    </DataX>
+
+                    <DataY type="Cross section" units="m2" size="3">
+                    0.000000e+0 0.0e-20 0.0e-20
+                    </DataY>
+                    </Process>
+                </Processes>
+            </Group>
+        </Groups>
+    </Database>
+
 </lxcat>
diff --git a/docs/src/modelling_ionization.md b/docs/src/modelling_ionization.md
@@ -16,7 +16,12 @@ Examples of synthetic data used for testing can be found in `data/test_neutral_e
 
 In general, the XML file has to have the following structure
 ```xml
+<?xml version="1.0" encoding="UTF-8"?>
 <lxcat version="" created="" message="example structure">
+    <References>
+        <Source>source. Generated on ---. All rights reserved.</Source>
+        <Reference>- ??? database, www.lxcat.net, retrieved on ---</Reference>
+    </References>
     <Database name="Linear dependence of cross-section on energy, He" id="LinearDB">
         <Groups>
             <Group id="He">
@@ -96,8 +101,8 @@ pair; each one is species-specific, therefore a vector for all the neutral speci
 [`create_computed_crosssections`](@ref) and passing an [`ElectronNeutralInteractions`](@ref) instance.
 
 ## Performing ionizing collisions
-Electron-neutral collisions with elastic scattering and ionization reactions are modelled by the [`ntc_n_e!`] and
-[`ntc_n_e_es!`] functions. The latter implements the Event Splitting
+Electron-neutral collisions with elastic scattering and ionization reactions are modelled by the [`ntc_n_e!`](@ref) and
+[`ntc_n_e_es!`](@ref) functions. The latter implements the Event Splitting
 of [Oblapenko et al. (2022)](https://doi.org/10.1016/j.jcp.2022.111390) and reduces the level of stochastic noise,
 but is suitable only for variable-weight simulations with particle merging.
 The function implements the collision mechanics as described in [Nanbu (2000)](https://doi.org/10.1109/27.887765);

diff --git a/docs/src/overview_1d.md b/docs/src/overview_1d.md
@@ -122,7 +122,7 @@ end
 
 ## Performing convection
 Having set up the grid and boundary conditions, we can convect particles.
-This is done by calling the `convect_particles!` function.
+This is done by calling the [`convect_particles!`](@ref) function.
 The convection should be followed by particle sorting before any computations of physical properties are done.
 
 ```julia
@@ -135,6 +135,20 @@ a `SurfProps` instance needs to be passed:
 convect_particles!(rng, grid, boundaries, particles[species_id], pia, species_id, species_data, surf_props, Δt)
 ```
 
+## Convection and sorting with precomputed particle/cell indices
+The approach described above assumes that during the convection process, the indices of the cells the particles find themselves in
+are not computed; therefore the call to `sort_particles!` requires passing in the `grid` instance, and the `sort_particles!`
+calls a `get_cell` function internally. For 1-D uniform grids, this is an efficient operation, but for other grid types,
+the computation of the cell index given only the particle position can be more expensive than keeping track of the cell index inside the convection routine.
+
+To this purpose, one can call [`convect_particles_and_compute_cell!`](@ref), which will also set the values of the `cell` array of the `ParticleVector`
+instance. One can then call another version of `sort_particles!` that does **not** take the grid as a parameter and instead uses the pre-computed `cell` values
+to sort the particles:
+
+```
+sort_particles!(gridsorter, particles[species_id], pia, species_id)
+```
+
 ## Bringing it all together
 Now we can combine all the pieces to set up a simulation of a single-species Couette flow in a channel
 with a width of 0.5 mm, discretized with 50 cells. The y-velocity of the left wall is assumed to be -500 m/s,

diff --git a/docs/src/reference_public.md b/docs/src/reference_public.md
@@ -172,6 +172,7 @@ sort_particles!
 ## Particle movement
 ```@docs
 convect_particles!
+convect_particles_and_compute_cell!
 ```
 
 ## Particle-surface interactions

diff --git a/src/Merzbild.jl b/src/Merzbild.jl
@@ -61,7 +61,7 @@ export estimate_sigma_g_w_max_ntc_n_e!, ntc_n_e!, ntc_n_e_es!
 export ParticleVector
 export AbstractGrid, Grid1DUniform, write_grid
 export GridSortInPlace, sort_particles!
-export MaxwellWallBC, MaxwellWalls1D, convect_particles!
+export MaxwellWallBC, MaxwellWalls1D, convect_particles!, convect_particles_and_compute_cell!
 export pretty_print_pia
 export ChunkExchanger, exchange_particles!, reset!, sort_particles_after_exchange!
 export count_disordered_particles, check_pia_is_correct, check_unique_index, check_unique_buffer

diff --git a/src/convection/convection_1D.jl b/src/convection/convection_1D.jl
@@ -205,4 +205,105 @@ function convect_particles!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D
     surface_props_scale!(species, surf_props, species_data, Δt)
 end
 
+
+"""
+    convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, Δt)
+
+Convect particles on a 1-D uniform grid and write post-convection cell index to `particles.cell`.
+
+# Positional arguments
+* `rng`: the random number generator
+* `grid`: the grid on which the convection is performed
+* `boundaries`: the `MaxwellWalls1D` struct describing the boundaries (it is assumed that the wall with index 1 is the left wall and
+    the wall with index 2 is the right wall)
+* `particles`: the `ParticleVector` of particles to be convected
+* `pia`: the `ParticleIndexerArray` instance
+* `species`: the index of the species being convected
+* `species_data`: the vector of `Species` data
+* `Δt`: the convection timestep
+"""
+function convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, Δt)
+    # @inbounds @simd for i in 1:pia.n_total[species]
+
+    @inbounds if pia.contiguous[species]
+        @inbounds n_tot = pia.n_total[species]
+        @inbounds @simd for i in 1:n_tot
+            convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt)
+            particles.cell[i] = get_cell(grid, particles[i].x)
+        end
+    else
+        for cell in 1:grid.n_cells
+            @inbounds s = pia.indexer[cell, species].start1
+            @inbounds e = pia.indexer[cell, species].end1
+
+            @inbounds @simd for i in s:e
+                convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt)
+                particles.cell[i] = get_cell(grid, particles[i].x)
+            end
+
+            @inbounds if pia.indexer[cell, species].n_group2 > 0
+                @inbounds s = pia.indexer[cell, species].start2
+                @inbounds e = pia.indexer[cell, species].end2
+
+                @inbounds @simd for i in s:e
+                    convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt)
+                    particles.cell[i] = get_cell(grid, particles[i].x)
+                end
+            end
+        end
+    end
+end
+
+"""
+    convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, surf_props::SurfProps, particles, pia, species, species_data, Δt)
+
+Convect particles on a 1-D uniform grid and write post-convection cell index to `particles.cell`, computing surface properties if particles hit a surface.
+
+# Positional arguments
+* `rng`: the random number generator
+* `grid`: the grid on which the convection is performed
+* `boundaries`: the `MaxwellWalls1D` struct describing the boundaries (it is assumed that the wall with index 1 is the left wall and
+    the wall with index 2 is the right wall)
+* `particles`: the `ParticleVector` of particles to be convected
+* `pia`: the `ParticleIndexerArray` instance
+* `species`: the index of the species being convected
+* `species_data`: the vector of `Species` data
+* `surf_props`: the `SurfProps` struct where the computed surface properties will be stored
+* `Δt`: the convection timestep
+"""
+function convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, surf_props::SurfProps, Δt)
+    # @inbounds @simd for i in 1:pia.n_total[species]
+
+    clear_props!(surf_props)
+    @inbounds if pia.contiguous[species]
+        @inbounds n_tot = pia.n_total[species]
+        @inbounds @simd for i in 1:n_tot
+            convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt)
+            particles.cell[i] = get_cell(grid, particles[i].x)
+        end
+    else
+        for cell in 1:grid.n_cells
+            @inbounds s = pia.indexer[cell, species].start1
+            @inbounds e = pia.indexer[cell, species].end1
+
+            @inbounds @simd for i in s:e
+                convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt)
+                particles.cell[i] = get_cell(grid, particles[i].x)
+            end
+
+            @inbounds if pia.indexer[cell, species].n_group2 > 0
+                @inbounds s = pia.indexer[cell, species].start2
+                @inbounds e = pia.indexer[cell, species].end2
+
+                @inbounds @simd for i in s:e
+                    convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt) 
+                    particles.cell[i] = get_cell(grid, particles[i].x)
+                end
+            end
+        end
+    end
+
+    surface_props_scale!(species, surf_props, species_data, Δt)
+end
+
 end