Add automatic padding for the FullGridCellList

Project.toml

@@ -1,7 +1,7 @@
 name = "PointNeighbors"
 uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 authors = ["Erik Faulhaber <erik.faulhaber@uni-koeln.de>"]
-version = "0.4.6-dev"
+version = "0.4.6"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/cell_lists/full_grid.jl

@@ -42,14 +42,16 @@ end
 
 supported_update_strategies(::FullGridCellList) = (SemiParallelUpdate, SerialUpdate)
 
-function FullGridCellList(; min_corner, max_corner, search_radius = 0.0,
+function FullGridCellList(; min_corner, max_corner,
+                          search_radius = zero(eltype(min_corner)),
                           backend = DynamicVectorOfVectors{Int32},
                           max_points_per_cell = 100)
-    # Pad domain to avoid 0 in cell indices due to rounding errors.
+    # Add one layer in each direction to make sure neighbor cells exist.
+    # Also pad domain a little more to avoid 0 in cell indices due to rounding errors.
     # We can't just use `eps()`, as one might use lower precision types.
     # This padding is safe, and will give us one more layer of cells in the worst case.
-    min_corner = SVector(Tuple(min_corner .- 1.0f-3 * search_radius))
-    max_corner = SVector(Tuple(max_corner .+ 1.0f-3 * search_radius))
+    min_corner = SVector(Tuple(min_corner .- 1.001f0 * search_radius))
+    max_corner = SVector(Tuple(max_corner .+ 1.001f0 * search_radius))
 
     if search_radius < eps()
         # Create an empty "template" cell list to be used with `copy_cell_list`
@@ -120,8 +122,10 @@ end
 function push_cell!(cell_list::FullGridCellList, cell, particle)
     (; cells) = cell_list
 
+    @boundscheck check_cell_bounds(cell_list, cell)
+
     # `push!(cell_list[cell], particle)`, but for all backends
-    pushat!(cells, cell_index(cell_list, cell), particle)
+    @inbounds pushat!(cells, cell_index(cell_list, cell), particle)
 
     return cell_list
 end
@@ -134,17 +138,21 @@ end
 @inline function push_cell_atomic!(cell_list::FullGridCellList, cell, particle)
     (; cells) = cell_list
 
+    @boundscheck check_cell_bounds(cell_list, cell)
+
     # `push!(cell_list[cell], particle)`, but for all backends.
     # The atomic version of `pushat!` uses atomics to avoid race conditions when `pushat!`
     # is used in a parallel loop.
-    pushat_atomic!(cells, cell_index(cell_list, cell), particle)
+    @inbounds pushat_atomic!(cells, cell_index(cell_list, cell), particle)
 
     return cell_list
 end
 
 function deleteat_cell!(cell_list::FullGridCellList, cell, i)
     (; cells) = cell_list
 
+    @boundscheck check_cell_bounds(cell_list, cell)
+
     # `deleteat!(cell_list[cell], i)`, but for all backends
     deleteatat!(cells, cell_index(cell_list, cell), i)
 end
@@ -170,8 +178,10 @@ end
     return cells[cell_index(cell_list, cell)]
 end
 
-@inline function is_correct_cell(cell_list::FullGridCellList, cell_coords, cell_index_)
-    return cell_index(cell_list, cell_coords) == cell_index_
+@inline function is_correct_cell(cell_list::FullGridCellList, cell, cell_index_)
+    @boundscheck check_cell_bounds(cell_list, cell)
+
+    return cell_index(cell_list, cell) == cell_index_
 end
 
 @inline index_type(::FullGridCellList) = Int32
@@ -188,7 +198,21 @@ function max_points_per_cell(cells::DynamicVectorOfVectors)
     return size(cells.backend, 1)
 end
 
-# Fallback when backend is a `Vector{Vector{T}}`
+# Fallback when backend is a `Vector{Vector{T}}`. Only used for copying the cell list.
 function max_points_per_cell(cells)
     return 100
 end
+
+@inline function check_cell_bounds(cell_list::FullGridCellList, cell::Tuple)
+    (; linear_indices) = cell_list
+
+    if !all(cell[i] in axes(linear_indices)[i] for i in eachindex(cell))
+        error("particle coordinates are NaN or outside the domain bounds of the cell list")
+    end
+end
+
+@inline function check_cell_bounds(cell_list::FullGridCellList, cell::Integer)
+    (; cells) = cell_list
+
+    checkbounds(cells, cell)
+end

src/nhs_grid.jl

@@ -185,7 +185,7 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::Abstra
     return neighborhood_search
 end
 
-# WARNING! Undocumented, experimental feature:
+# WARNING! Experimental feature:
 # By default, determine the parallelization backend from the type of `x`.
 # Optionally, pass a `KernelAbstractions.Backend` to run the KernelAbstractions.jl code
 # on this backend. This can be useful to run GPU kernels on the CPU by passing

src/vector_of_vectors.jl

@@ -59,9 +59,10 @@ end
 @inline function pushat!(vov::DynamicVectorOfVectors, i, value)
     (; backend, lengths) = vov
 
+    # Outer bounds check
     @boundscheck checkbounds(vov, i)
 
-    # Activate new entry in column `i`
+    # Activate new entry in column `i`. Inner bounds check is included here.
     backend[lengths[i] += 1, i] = value
 
     return vov
@@ -76,15 +77,17 @@ end
 @inline function pushat_atomic!(vov::DynamicVectorOfVectors, i, value)
     (; backend, lengths) = vov
 
+    # Outer bounds check
     @boundscheck checkbounds(vov, i)
 
     # Increment the column length with an atomic add to avoid race conditions.
     # Store the new value since it might be changed immediately afterwards by another
     # thread.
-    new_length = Atomix.@atomic lengths[i] += 1
+    new_length = @inbounds Atomix.@atomic lengths[i] += 1
 
     # We can write here without race conditions, since the atomic add guarantees
     # that `new_length` is different for each thread.
+    # Inner bounds check is included here.
     backend[new_length, i] = value
 
     return vov
@@ -101,10 +104,10 @@ end
 
     # Replace value to delete by the last value in this column
     last_value = backend[lengths[i], i]
-    backend[j, i] = last_value
+    @inbounds backend[j, i] = last_value
 
     # Remove the last value in this column
-    lengths[i] -= 1
+    @inbounds lengths[i] -= 1
 
     return vov
 end

test/cell_lists/full_grid.jl

@@ -0,0 +1,47 @@
+@testset "`FullGridCellList`" verbose=true begin
+    # Test that `update!` throws an error when a particle is outside the bounding box
+    @testset "`update!` bounds check" begin
+        @testset "$(N)D" for N in 1:3
+            min_corner = fill(0.0, N)
+            max_corner = fill(10.0, N)
+            search_radius = 1.0
+
+            cell_list = FullGridCellList(; search_radius, min_corner, max_corner)
+
+            # Introduce the same rounding errors for this to pass
+            @test cell_list.min_corner == fill(-1.001f0, N)
+            @test cell_list.max_corner == fill(10.0 + 1.001f0, N)
+
+            nhs = GridNeighborhoodSearch{N}(; cell_list, search_radius)
+            y = rand(N, 10)
+            error_nan = ErrorException("particle coordinates contain NaNs")
+            error_string_bounds = "particle coordinates are outside the domain bounds of the cell list"
+            error_bounds = ErrorException(error_string_bounds)
+
+            y[1, 7] = NaN
+            @test_throws error_nan initialize!(nhs, y, y)
+
+            y[1, 7] = -0.01
+            @test_throws error_bounds initialize!(nhs, y, y)
+
+            y[1, 7] = 10.01
+            @test_throws error_bounds initialize!(nhs, y, y)
+
+            y[1, 7] = 0.0
+            @test_nowarn_mod initialize!(nhs, y, y)
+            @test_nowarn_mod update!(nhs, y, y)
+
+            y[1, 7] = 10.0
+            @test_nowarn_mod update!(nhs, y, y)
+
+            y[1, 7] = NaN
+            @test_throws error_nan update!(nhs, y, y)
+
+            y[1, 7] = 10.01
+            @test_throws error_bounds update!(nhs, y, y)
+
+            y[1, 7] = -0.01
+            @test_throws error_bounds update!(nhs, y, y)
+        end
+    end
+end

test/nhs_grid.jl

@@ -75,7 +75,7 @@
         @test PointNeighbors.cell_coords(coords1, nothing, cell_list, (1.0, 1.0)) ==
               (typemax(Int), typemin(Int)) .+ 1
         @test PointNeighbors.cell_coords(coords2, nothing, cell_list, (1.0, 1.0)) ==
-              (typemax(Int), 0) .+ 1
+              (typemax(Int), 1) .+ 1
         @test PointNeighbors.cell_coords(coords3, nothing, cell_list, (1.0, 1.0)) ==
               (typemax(Int), typemin(Int)) .+ 1
     end

test/unittest.jl

@@ -5,4 +5,5 @@
     include("nhs_trivial.jl")
     include("nhs_grid.jl")
     include("neighborhood_search.jl")
+    include("cell_lists/full_grid.jl")
 end;