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`
@@ -188,7 +190,22 @@ 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
+
+function check_domain_bounds(cell_list::FullGridCellList, y, search_radius)
+    if any(isnan, y)
+        error("particle coordinates contain NaNs")
+    end
+
+    # Require one extra layer in each direction to make sure neighbor cells exist.
+    # Convert to CPU in case this lives on the GPU.
+    min_corner = Array(minimum(y, dims = 2) .- search_radius)
+    max_corner = Array(maximum(y, dims = 2) .+ search_radius)
+
+    if any(min_corner .< cell_list.min_corner) || any(max_corner .> cell_list.max_corner)
+        error("particle coordinates are outside the domain bounds of the cell list")
+    end
+end

src/nhs_grid.jl

@@ -165,6 +165,9 @@ end
 function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::AbstractMatrix)
     (; cell_list) = neighborhood_search
 
+    check_domain_bounds(neighborhood_search.cell_list, y,
+                        search_radius(neighborhood_search))
+
     empty!(cell_list)
 
     if neighborhood_search.search_radius < eps()
@@ -185,7 +188,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
@@ -203,10 +206,16 @@ end
 # Update only with neighbor coordinates
 function update_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS},
                       y::AbstractMatrix; parallelization_backend = y) where {NDIMS}
+    check_domain_bounds(neighborhood_search.cell_list, y,
+                        search_radius(neighborhood_search))
+
     update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i);
                  parallelization_backend)
 end
 
+# This is dispatched in `cell_lists/full_grid.jl` for the `FullGridCellList`
+check_domain_bounds(_, _, _) = nothing
+
 # Serial and semi-parallel update.
 # See the warning above. `parallelization_backend = nothing` will use `Polyester.@batch`.
 function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any,

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;