Fix neighbor rebuild storage and add guard/test

src/SPHCellList.jl

@@ -908,7 +908,7 @@ using TimerOutputs: @timeit
 
         # Produce sorting related variables
         ParticleRanges         = zeros(Int, NumberOfPoints + 1 + 1) # +1 for the last particle, +1 for dummy entry
-        UniqueCells            = zeros(CartesianIndex{Dimensions}, NumberOfPoints)
+        UniqueCells            = zeros(CartesianIndex{Dimensions}, NumberOfPoints + 1) # +1 for the dummy sentinel cell
         CellListIndices        = zeros(Int, NumberOfPoints)
         FullStencil            = ConstructStencil(Val(Dimensions))
         NeighborCellLists      = [Int[] for _ in 1:length(UniqueCells)]

src/SPHNeighborList.jl

@@ -91,6 +91,11 @@ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace,
                           ParticleRanges, UniqueCells, CellListIndices)
     ExtractCells!(Particles, InverseCutOff)
 
+    RequiredCellSlots = length(Particles.Cells) + 1
+    if length(UniqueCells) < RequiredCellSlots
+        throw(ArgumentError("UniqueCells must have at least one sentinel slot plus one slot per particle; got $(length(UniqueCells)) slots for $(length(Particles.Cells)) particles"))
+    end
+
     sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace)
     Cells = @views Particles.Cells
     UniqueCells[1] = MinCell(eltype(Cells))

test/runtests.jl

@@ -2,16 +2,18 @@ using Test
 using SPHExample
 using StaticArrays
 using StructArrays
+using LinearAlgebra: norm
 
 @testset "time stepping" begin
     pos = [SVector{2,Float64}(0.0, 0.0), SVector{2,Float64}(1.0, 0.0)]
     vel = [SVector{2,Float64}(0.0, 0.0), SVector{2,Float64}(0.0, 0.0)]
     acc = [SVector{2,Float64}(0.0, 0.0), SVector{2,Float64}(0.0, -9.81)]
     sc = SimulationConstants{Float64}()
     ker = SPHKernelInstance{2, Float64}(WendlandC2(); dx=sc.dx)
-    dt  = Δt(pos, vel, acc, sc, ker)
+    max_acc = maximum(norm, acc)
+    dt  = Δt(max_acc, sc, ker)
     @test dt > 0
-    alloc = @allocated Δt(pos, vel, acc, sc, ker)
+    alloc = @allocated Δt(max_acc, sc, ker)
     @test alloc == 0
 end
 
@@ -51,15 +53,14 @@ end
 
     for _ in 1:1000
         ResetArrays!(dρdtI, particles.Acceleration)
-        dt = Δt(particles.Position, particles.Velocity, particles.Acceleration,
-                 sc, ker)
+        dt = Δt(maximum(norm, particles.Acceleration), sc, ker)
         dt2 = dt / 2
 
         SPHExample.SPHCellList.HalfTimeStep(meta, sc, particles, pos_n, vel_n,
                                            ρ_n, dρdtI, dt2)
         LimitDensityAtBoundary!(ρ_n, sc.ρ₀, particles.MotionLimiter)
         Pressure!(press, ρ_n, sc)
-        SPHExample.SPHCellList.FullTimeStep(meta, ker, sc, particles, ∇C, ∇r, dt)
+        SPHExample.SPHCellList.FullTimeStep(meta, ker, sc, particles, vel_n, ∇C, ∇r, dt)
         DensityEpsi!(dens, dρdtI, ρ_n, dt)
         LimitDensityAtBoundary!(dens, sc.ρ₀, particles.MotionLimiter)
         SPHExample.SPHCellList.UpdateMetaData!(meta, dt)
@@ -72,3 +73,31 @@ end
     @test particles.Velocity[1][1] == 0
     @test particles.Velocity[1][2] < 0
 end
+
+@testset "neighbor rebuild with all unique cells" begin
+    D = 2
+    T = Float64
+    N = 3
+    particles = StructArray((
+        Cells = [CartesianIndex(0, 0) for _ in 1:N],
+        Position = [SVector{D,T}(0, 0), SVector{D,T}(10, 0), SVector{D,T}(20, 0)],
+    ))
+    ParticleRanges = zeros(Int, N + 2)
+    UniqueCells = zeros(CartesianIndex{D}, N + 1)
+    CellListIndices = zeros(Int, N)
+    _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(particles), N)
+
+    for _ in 1:10
+        IndexCounter = UpdateNeighbors!(particles, one(T), SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices)
+        UniqueCellsView = view(UniqueCells, 1:IndexCounter)
+        NeighborCellLists = [Int[] for _ in 1:IndexCounter]
+
+        @test IndexCounter == N + 1
+        @test ParticleRanges[IndexCounter + 1] == N + 1
+        @test CellListIndices == collect(2:(N + 1))
+        @test allunique(UniqueCellsView[2:end])
+
+        BuildNeighborCellLists!(NeighborCellLists, ConstructStencil(Val(D)), UniqueCellsView, ParticleRanges)
+        @test length(NeighborCellLists) == IndexCounter
+    end
+end