fix nbatches update bug, and tiny-system issue of batches

CHANGELOG.md

@@ -32,10 +32,12 @@ Version 0.10.0-DEV
 - ![INFO][badge-info] Remove deprecated code, improve code coverage.
 - ![INFO][badge-info] Remove internal coordinate swapping by size.
 - ![INFO][badge-info] Internal implementation of neighborlist moved to ParticleSystem interface.
-- ![INFO][badge-info] The internal representation of coordinantes is done with `ParticleSytemPositions` array type - which is not a subtype of abstract array. 
+- ![INFO][badge-info] The internal representation of coordinantes is done with `ParticleSystemPositions` array type - which is not a subtype of abstract array. 
 - ![INFO][badge-info] Add performance test.
 - ![INFO][badge-info] Default name for output in ParticleSystem is `default_output_name`, but it can be accessed though `sys.output`, as before.
 - ![BUGFIX][badge-bugfix] Fix cross-computation (`pairwise!(f, x, sys)`) with `NonPeriodicCell` failing to find pairs when particle coordinates span negative values or large coordinate ranges / do not modify input coordinates in `NonPeriodicCell`. This bug was never released, it was created and fixed in the development version.
+- ![BUGFIX][badge-bugfix] If the number of particles is reduced by updating and becomes smaller than nbatches, there was a crash. Fixed. Also fixed automatic updating of nbatches when first set of positions change (bug never released).
+- ![BUGFIX][badge-bugfix] Fix `nbatches` not being updated for the `:map` phase of `CellListPair` systems when the particle count changes (bug never released).
 
 Version 0.9.17
 --------------

src/internals/CellLists.jl

@@ -216,13 +216,17 @@ function update_number_of_batches!(
         _nbatches::Union{NumberOfBatches, Nothing} = nothing;
         parallel = true
     ) where {N, T}
-    # For CellListPair, compute nbatches independently for each set
-    # Build: based on each set's particle count
-    # Map: based on product of particle counts for cross-interactions
-    # If _nbatches is nothing (default), use auto=(true, true) for proper recalculation
-    # Otherwise, respect user-specified values
     if isnothing(_nbatches)
-        _nbatches = NumberOfBatches((true, true), (0, 0))
+        # Default call (e.g. when particle count changed): reconstruct auto flags from stored
+        # state. build_cell_lists.auto is correctly preserved across updates. map_computation.auto
+        # is always stored as false in individual lists (to prevent per-list recalculation based
+        # on individual particle counts), so we use build_auto as the pair-level map auto flag.
+        # This is correct for all documented use cases: (0,0) all-auto and (n,m) all-manual.
+        build_auto = first(cl.ref_list.nbatches.build_cell_lists)
+        _nbatches = NumberOfBatches(
+            (build_auto, build_auto),
+            (last(cl.ref_list.nbatches.build_cell_lists), last(cl.ref_list.nbatches.map_computation))
+        )
     end
     auto = (first(_nbatches.build_cell_lists), first(_nbatches.map_computation))
     n_build_ref = last(_nbatches.build_cell_lists)
@@ -756,6 +760,7 @@ function UpdateCellList!(
         @sync for ibatch in eachindex(aux.idxs, aux.lists)
             @spawn begin
                 prange = aux.idxs[ibatch]
+                isempty(prange) && return
                 aux.lists[ibatch] = reset!(aux.lists[ibatch], box, length(prange))
                 xt = @view(x[prange])
                 aux.lists[ibatch] = add_particles!(xt, box, prange[begin] - 1, aux.lists[ibatch])

test/API/ParticleSystem.jl

@@ -293,6 +293,15 @@ end
             CellListMap._nbatches_map_computation(5000 * 10000)
     )
     @test CellListMap.get_nbatches(sys) == expected
+
+    # If the setting is manual, no changes are expected
+    sys = ParticleSystem(xpositions = x, ypositions = y, cutoff = 0.1, unitcell = [1, 1, 1], output = 0.0, nbatches=(2,2))
+    @test CellListMap.get_nbatches(sys) == (2,2)
+    x = rand(SVector{3, Float64}, 5000)
+    y = rand(SVector{3, Float64}, 10000)
+    update!(sys; xpositions=x, ypositions=y)
+    @test CellListMap.get_nbatches(sys) == (2,2)
+
 end
 
 @testitem "ParticleSystem - validate coordinates" begin

test/internals/CellLists.jl

@@ -4,25 +4,25 @@
     x = rand(SVector{3,Float64}, 100)
     y = rand(SVector{3,Float64}, 10)
     box = CellListMap.Box([1, 1, 1], 0.1)
-    cl = CellListMap.CellList(PSP(x), box; nbatches = (2, 4))
+    cl = CellListMap.CellList(PSP(x), box; nbatches=(2, 4))
     @test CellListMap.nbatches(cl) == (2, 4)
     @test CellListMap.nbatches(cl, :build) == 2
     @test CellListMap.nbatches(cl, :map) == 4
-    cl = CellListMap.CellList(PSP(x), PSP(y), box; nbatches = (2, 4))
+    cl = CellListMap.CellList(PSP(x), PSP(y), box; nbatches=(2, 4))
     @test CellListMap.nbatches(cl) == (2, 4)
     @test CellListMap.nbatches(cl, :build) == 2
     @test CellListMap.nbatches(cl, :map) == 4
     cl = CellListMap.CellList(PSP(x), PSP(y), box)
     # For CellListPair, build batches are independent, but map batches should match
     @test CellListMap.nbatches(cl.ref_list, :map) == CellListMap.nbatches(cl.target_list, :map)
-    cl = CellListMap.CellList(PSP(x), box; nbatches = (2, 4), parallel = false)
+    cl = CellListMap.CellList(PSP(x), box; nbatches=(2, 4), parallel=false)
     @test CellListMap.nbatches(cl) == (1, 1)
     # The automatic set of number of batches for this small system:
     if Threads.nthreads() == 10
-        x = rand(SVector{3, Float64}, 10)
-        sys = ParticleSystem(positions = x, cutoff = 0.1, unitcell = [1, 1, 1], output = 0.0)
+        x = rand(SVector{3,Float64}, 10)
+        sys = ParticleSystem(positions=x, cutoff=0.1, unitcell=[1, 1, 1], output=0.0)
         @test CellListMap.get_nbatches(sys) == (4, 4)
-        x = rand(SVector{3, Float64}, 10000)
+        x = rand(SVector{3,Float64}, 10000)
         resize!(sys.xpositions, length(x))
         sys.xpositions .= x
         pairwise!((pair, out) -> out += pair.d2, sys)
@@ -32,33 +32,54 @@
     end
 end
 
-@testitem "automatic nbatches update on UpdateCellList!" setup=[Testing] begin
+@testitem "automatic nbatches update on UpdateCellList!" setup = [Testing] begin
     using CellListMap, StaticArrays
     # 3-arg UpdateCellList! (without preallocated aux) should update nbatches
     # when the number of particles changes
-    x = rand(SVector{3, Float64}, 2)
+    x = rand(SVector{3,Float64}, 2)
     box = CellListMap.Box([1, 1, 1], 0.1)
     cl = CellListMap.CellList(PSP(x), box)
     nb_initial = CellListMap.nbatches(cl)
-    x = rand(SVector{3, Float64}, 10000)
+    x = rand(SVector{3,Float64}, 10000)
     cl = CellListMap.UpdateCellList!(PSP(x), box, cl)
     expected = (
         CellListMap._nbatches_build_cell_lists(10000),
         CellListMap._nbatches_map_computation(10000),
     )
     @test CellListMap.nbatches(cl) == expected
     # nbatches should not change when particle count is unchanged
-    x .= rand(SVector{3, Float64}, 10000)
+    x .= rand(SVector{3,Float64}, 10000)
     cl = CellListMap.UpdateCellList!(PSP(x), box, cl)
     @test CellListMap.nbatches(cl) == expected
+
+    # Test updating to very few particles, in parallel runs (caused crash in a bug)
+    sys = ParticleSystem(
+        xpositions=rand(3, 10),
+        ypositions=rand(3, 100),
+        cutoff=0.2,
+        output=0.0,
+        unitcell=[1, 1, 1],
+        parallel=true,
+        nbatches=(4, 4),
+    )
+    @test CellListMap.get_nbatches(sys) == (4,4)
+    pairwise!((pair, sumd) -> sumd += pair.d, sys)
+    update!(sys; xpositions=rand(3,1))
+    @test pairwise!((pair, sumd) -> sumd += pair.d, sys) > -1 
+    update!(sys; xpositions=rand(3,100), ypositions=rand(3,1))
+    @test pairwise!((pair, sumd) -> sumd += pair.d, sys) > -1 
+    update!(sys; xpositions=rand(3,1), ypositions=rand(3,1))
+    @test pairwise!((pair, sumd) -> sumd += pair.d, sys) > -1
+    @test CellListMap.get_nbatches(sys) == (4,4)
+
 end
 
-@testitem "nbatches symbol variants, show, and CellListPair" setup=[Testing] begin
+@testitem "nbatches symbol variants, show, and CellListPair" setup = [Testing] begin
     using CellListMap, StaticArrays
-    x = rand(SVector{3, Float64}, 100)
-    y = rand(SVector{3, Float64}, 50)
+    x = rand(SVector{3,Float64}, 100)
+    y = rand(SVector{3,Float64}, 50)
     box = CellListMap.Box([1, 1, 1], 0.1)
-    cl = CellListMap.CellList(PSP(x), box; nbatches = (2, 4))
+    cl = CellListMap.CellList(PSP(x), box; nbatches=(2, 4))
     # Full symbol names
     @test CellListMap.nbatches(cl, :map_computation) == 4
     @test CellListMap.nbatches(cl, :build_cell_lists) == 2
@@ -82,57 +103,57 @@ end
     @test_throws ArgumentError CellListMap.set_idxs!(idxs, 100, 3)
 end
 
-@testitem "UpdateCellList! for CellListPair" setup=[Testing] begin
+@testitem "UpdateCellList! for CellListPair" setup = [Testing] begin
     using CellListMap, StaticArrays
-    x = rand(SVector{3, Float64}, 100)
-    y = rand(SVector{3, Float64}, 50)
+    x = rand(SVector{3,Float64}, 100)
+    y = rand(SVector{3,Float64}, 50)
     box = CellListMap.Box([1, 1, 1], 0.1)
     cl = CellListMap.CellList(PSP(x), PSP(y), box)
     # 3-arg update, serial path
-    x2 = rand(SVector{3, Float64}, 100)
-    y2 = rand(SVector{3, Float64}, 50)
-    cl = CellListMap.UpdateCellList!(PSP(x2), PSP(y2), box, cl; parallel = false)
+    x2 = rand(SVector{3,Float64}, 100)
+    y2 = rand(SVector{3,Float64}, 50)
+    cl = CellListMap.UpdateCellList!(PSP(x2), PSP(y2), box, cl; parallel=false)
     @test cl.ref_list.n_real_particles == 100  # x (first arg) has 100
     @test cl.target_list.n_real_particles == 50  # y (second arg) has 50
     # 3-arg update, parallel path
-    x3 = rand(SVector{3, Float64}, 100)
-    y3 = rand(SVector{3, Float64}, 50)
-    cl = CellListMap.UpdateCellList!(PSP(x3), PSP(y3), box, cl; parallel = true)
+    x3 = rand(SVector{3,Float64}, 100)
+    y3 = rand(SVector{3,Float64}, 50)
+    cl = CellListMap.UpdateCellList!(PSP(x3), PSP(y3), box, cl; parallel=true)
     @test cl.ref_list.n_real_particles == 100  # x (first arg) has 100
     @test cl.target_list.n_real_particles == 50  # y (second arg) has 50
 end
 
-@testitem "celllists - validate coordinates" setup=[Testing] begin
+@testitem "celllists - validate coordinates" setup = [Testing] begin
     using StaticArrays
-    x = rand(SVector{3, Float64}, 100)
+    x = rand(SVector{3,Float64}, 100)
     x[50] = SVector(1.0, NaN, 1.0)
-    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, unitcell=[1,1,1], cutoff=0.1, output=0.0)
+    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
     x[50] = SVector(1.0, 0.5, 1.0)
-    sys = ParticleSystem(xpositions=x, unitcell=[1,1,1], cutoff=0.1, output=0.0)
+    sys = ParticleSystem(xpositions=x, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
     sys.xpositions[1] = SVector(1.0, NaN, 1.0)
     @test_throws "Invalid coordinates" pairwise!(f1, sys)
     x[50] = SVector(1.0, 0.5, 1.0)
-    y = rand(SVector{3, Float64}, 100)
+    y = rand(SVector{3,Float64}, 100)
     y[50] = SVector(1.0, NaN, 1.0)
-    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, ypositions=y, unitcell=[1,1,1], cutoff=0.1, output=0.0)
+    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, ypositions=y, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
     y[50] = SVector(1.0, 0.5, 1.0)
-    sys = ParticleSystem(xpositions=x, ypositions=y, unitcell=[1,1,1], cutoff=0.1, output=0.0)
+    sys = ParticleSystem(xpositions=x, ypositions=y, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
     sys.ypositions[1] = SVector(1.0, NaN, 1.0)
     @test_throws "Invalid coordinates" pairwise!(f1, sys)
     x = rand(3, 100)
     x[2, 50] = NaN
-    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, unitcell=[1,1,1], cutoff=0.1, output=0.0)
+    @test_throws "Invalid coordinates" ParticleSystem(xpositions=x, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
 end
 
 @testitem "CellList dimension mismatch error" begin
     using CellListMap, StaticArrays
     PSP = CellListMap.ParticleSystemPositions
     # 3D positions with 2D box should throw DimensionMismatch
-    x = rand(SVector{3, Float64}, 100)
+    x = rand(SVector{3,Float64}, 100)
     box = CellListMap.Box([1.0, 1.0], 0.1)
     @test_throws DimensionMismatch CellListMap.CellList(PSP(x), box)
     # 2D positions with 3D box should throw DimensionMismatch
-    x = rand(SVector{2, Float64}, 100)
+    x = rand(SVector{2,Float64}, 100)
     box = CellListMap.Box([1.0, 1.0, 1.0], 0.1)
     @test_throws DimensionMismatch CellListMap.CellList(PSP(x), box)
 end
@@ -143,7 +164,7 @@ end
     # The function adjusts cell indices at the boundary cells (lcell and nc - lcell + 1)
 
     # Create a box where we can control the boundary conditions
-    box = CellListMap.Box([1.0, 1.0, 1.0], 0.1, lcell = 1)
+    box = CellListMap.Box([1.0, 1.0, 1.0], 0.1, lcell=1)
 
     # Test lower boundary case: particle at cell index == lcell should be moved to lcell + 1
     cartesian_index = CartesianIndex(box.lcell, box.lcell + 1, box.lcell + 2)