Restructure tests (#13)

* Restructure tests

* Also combine tests in `grid_nhs.jl`

* Improve comments

* Vary initalize and update functions

* Reformat code

* Remove ugly inline functions

* Extract variables

* Update test/nhs_grid.jl

Co-authored-by: Niklas Neher <[email protected]>

* Update test/nhs_grid.jl

Co-authored-by: Niklas Neher <[email protected]>

---------

Co-authored-by: Niklas Neher <[email protected]>

Author: efaulhaber

test/neighborhood_search.jl

@@ -0,0 +1,79 @@
+# This file contains tests for the generic functions in `src/neighborhood_search.jl` and
+# tests comparing all NHS implementations against the `TrivialNeighborhoodSearch`.
+@testset verbose=true "Neighborhood Searches" begin
+    @testset verbose=true "Periodicity" begin
+        # These examples are constructed by hand and are therefore a good test for the
+        # trivial neighborhood search as well.
+        # (As opposed to the tests below that are just comparing against the trivial NHS.)
+
+        # Names, coordinates and corresponding periodic boxes for each test
+        names = [
+            "Simple Example 2D",
+            "Box Not Multiple of Search Radius 2D",
+            "Simple Example 3D",
+        ]
+
+        coordinates = [
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.39],
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.42],
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.39
+             0.14 0.34 0.12 0.06 0.13],
+        ]
+
+        periodic_boxes = [
+            ([-0.1, -0.2], [0.2, 0.4]),
+            # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test
+            # to avoid split cells.
+            ([-0.1, -0.2], [0.205, 0.43]),
+            ([-0.1, -0.2, 0.05], [0.2, 0.4, 0.35]),
+        ]
+
+        @testset verbose=true "$(names[i])" for i in eachindex(names)
+            coords = coordinates[i]
+
+            NDIMS = size(coords, 1)
+            n_particles = size(coords, 2)
+            search_radius = 0.1
+
+            neighborhood_searches = [
+                TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_particles,
+                                                 periodic_box_min_corner = periodic_boxes[i][1],
+                                                 periodic_box_max_corner = periodic_boxes[i][2]),
+                GridNeighborhoodSearch{NDIMS}(search_radius, n_particles,
+                                              periodic_box_min_corner = periodic_boxes[i][1],
+                                              periodic_box_max_corner = periodic_boxes[i][2]),
+            ]
+            neighborhood_searches_names = [
+                "`TrivialNeighborhoodSearch`",
+                "`GridNeighborhoodSearch`",
+            ]
+
+            # Run this for every neighborhood search
+            @testset "$(neighborhood_searches_names[j])" for j in eachindex(neighborhood_searches_names)
+                nhs = neighborhood_searches[j]
+
+                initialize!(nhs, coords, coords)
+
+                neighbors = [Int[] for _ in axes(coords, 2)]
+
+                for_particle_neighbor(coords, coords, nhs,
+                                      particles = axes(coords, 2)) do particle, neighbor,
+                                                                      pos_diff, distance
+                    append!(neighbors[particle], neighbor)
+                end
+
+                # All of these tests are designed to yield the same neighbor lists.
+                # Note that we have to sort the neighbor lists because neighborhood searches
+                # might produce different orders.
+                @test sort(neighbors[1]) == [1, 3, 5]
+                @test sort(neighbors[2]) == [2]
+                @test sort(neighbors[3]) == [1, 3]
+                @test sort(neighbors[4]) == [4]
+                @test sort(neighbors[5]) == [1, 5]
+            end
+        end
+    end
+end;

test/nhs_grid.jl

@@ -1,6 +1,6 @@
 @testset verbose=true "GridNeighborhoodSearch" begin
-    @testset "Coordinate Limits" begin
-        # Test the threshold for very large and very small coordinates.
+    @testset "Cells at Coordinate Limits" begin
+        # Test the threshold for very large and very small coordinates
         coords1 = [Inf, -Inf]
         coords2 = [NaN, 0]
         coords3 = [typemax(Int) + 1.0, -typemax(Int) - 1.0]
@@ -15,8 +15,8 @@
 
     @testset "Rectangular Point Cloud 2D" begin
         #### Setup
-        # Rectangular filled with equidistant spaced particles
-        # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35)
+        # Rectangle of equidistantly spaced particles
+        # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35).
         range = -0.25:0.1:0.35
         coordinates1 = hcat(collect.(Iterators.product(range, range))...)
         nparticles = size(coordinates1, 2)
@@ -28,8 +28,7 @@
         # Create neighborhood search
         nhs1 = GridNeighborhoodSearch{2}(radius, nparticles)
 
-        coords_fun(i) = coordinates1[:, i]
-        initialize_grid!(nhs1, coords_fun)
+        initialize_grid!(nhs1, coordinates1)
 
         # Get each neighbor for `particle_position1`
         neighbors1 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1)))
@@ -38,57 +37,51 @@
         coordinates2 = coordinates1 .+ [1.4, -3.5]
 
         # Update neighborhood search
-        coords_fun2(i) = coordinates2[:, i]
-        update_grid!(nhs1, coords_fun2)
+        update_grid!(nhs1, coordinates2)
 
         # Get each neighbor for updated NHS
-        neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1,
-                                                              nhs1)))
+        neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1)))
 
         # Change position
         particle_position2 = particle_position1 .+ [1.4, -3.5]
 
         # Get each neighbor for `particle_position2`
-        neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2,
-                                                              nhs1)))
+        neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2, nhs1)))
 
         # Double search radius
         nhs2 = GridNeighborhoodSearch{2}(2 * radius, size(coordinates1, 2))
-        initialize_grid!(nhs2, coords_fun)
+        initialize!(nhs2, coordinates1, coordinates1)
 
         # Get each neighbor in double search radius
-        neighbors4 = sort(collect(PointNeighbors.eachneighbor(particle_position1,
-                                                              nhs2)))
+        neighbors4 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs2)))
 
         # Move particles
         coordinates2 = coordinates1 .+ [0.4, -0.4]
 
         # Update neighborhood search
-        update_grid!(nhs2, coords_fun2)
+        update!(nhs2, coordinates2, coordinates2)
 
         # Get each neighbor in double search radius
-        neighbors5 = sort(collect(PointNeighbors.eachneighbor(particle_position1,
-                                                              nhs2)))
+        neighbors5 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs2)))
 
-        #### Verification
+        #### Verification against lists of potential neighbors built by hand
         @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33]
 
         @test neighbors2 == Int[]
 
         @test neighbors3 == [17, 18, 19, 24, 25, 26, 31, 32, 33]
 
         @test neighbors4 == [
-            9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25,
-            26, 27, 28, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47,
-            48, 49]
+            9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31,
+            32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49]
 
         @test neighbors5 == [36, 37, 38, 43, 44, 45]
     end
 
     @testset "Rectangular Point Cloud 3D" begin
         #### Setup
-        # Rectangular filled with equidistant spaced particles
-        # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y) = (0.35, 0.35, 0.35)
+        # Rectangle of equidistantly spaced particles
+        # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y, z) = (0.35, 0.35, 0.35).
         range = -0.25:0.1:0.35
         coordinates1 = hcat(collect.(Iterators.product(range, range, range))...)
         nparticles = size(coordinates1, 2)
@@ -104,8 +97,7 @@
         initialize_grid!(nhs1, coords_fun)
 
         # Get each neighbor for `particle_position1`
-        neighbors1 = sort(collect(PointNeighbors.eachneighbor(particle_position1,
-                                                              nhs1)))
+        neighbors1 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1)))
 
         # Move particles
         coordinates2 = coordinates1 .+ [1.4, -3.5, 0.8]
@@ -115,17 +107,15 @@
         update_grid!(nhs1, coords_fun2)
 
         # Get each neighbor for updated NHS
-        neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1,
-                                                              nhs1)))
+        neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1)))
 
         # Change position
         particle_position2 = particle_position1 .+ [1.4, -3.5, 0.8]
 
         # Get each neighbor for `particle_position2`
-        neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2,
-                                                              nhs1)))
+        neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2, nhs1)))
 
-        #### Verification
+        #### Verification against lists of potential neighbors built by hand
         @test neighbors1 ==
               [115, 116, 117, 122, 123, 124, 129, 130, 131, 164, 165, 166, 171, 172,
             173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229]
@@ -137,59 +127,46 @@
             173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229]
     end
 
-    @testset verbose=true "Periodicity 2D" begin
-        @testset "Simple Example" begin
-            coords = [-0.08 0.0 0.18 0.1 -0.08
-                      -0.12 -0.05 -0.09 0.15 0.39]
-
-            # 3 x 6 cells
-            nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2),
-                                            periodic_box_min_corner = [-0.1, -0.2],
-                                            periodic_box_max_corner = [0.2, 0.4])
-
-            initialize_grid!(nhs, coords)
-
-            neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i],
-                                                                  nhs)))
-                         for i in 1:5]
-
-            # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells
-            @test neighbors[1] == [1, 2, 3, 5]
-            @test neighbors[2] == [1, 2, 3]
-            @test neighbors[3] == [1, 2, 3]
-            @test neighbors[4] == [4]
-            @test neighbors[5] == [1, 5]
-
-            neighbors_loop = [Int[] for _ in axes(coords, 2)]
-
-            for_particle_neighbor(coords, coords, nhs,
-                                  particles = axes(coords, 2)) do particle,
-                                                                  neighbor,
-                                                                  pos_diff,
-                                                                  distance
-                append!(neighbors_loop[particle], neighbor)
-            end
-
-            @test sort(neighbors_loop[1]) == [1, 3, 5]
-            @test sort(neighbors_loop[2]) == [2]
-            @test sort(neighbors_loop[3]) == [1, 3]
-            @test sort(neighbors_loop[4]) == [4]
-            @test sort(neighbors_loop[5]) == [1, 5]
-        end
-
-        @testset "Rounding Up Cell Sizes" begin
-            coords = [-0.08 0.0 0.18 0.1 -0.08
-                      -0.12 -0.05 -0.09 0.15 0.42]
-
-            # 3 x 6 cells
-            nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2),
-                                            periodic_box_min_corner = [-0.1, -0.2],
-                                            periodic_box_max_corner = [0.205, 0.43])
+    @testset verbose=true "Periodicity" begin
+        # These setups are the same as in `test/neighborhood_search.jl`,
+        # but instead of testing the actual neighbors with `for_particle_neighbor`,
+        # we only test the potential neighbors (particles in neighboring cells) here.
+
+        # Names, coordinates and corresponding periodic boxes for each test
+        names = [
+            "Simple Example 2D",
+            "Box Not Multiple of Search Radius 2D",
+            "Simple Example 3D",
+        ]
+
+        coordinates = [
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.39],
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.42],
+            [-0.08 0.0 0.18 0.1 -0.08
+             -0.12 -0.05 -0.09 0.15 0.39
+             0.14 0.34 0.12 0.06 0.13],
+        ]
+
+        periodic_boxes = [
+            ([-0.1, -0.2], [0.2, 0.4]),
+            # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test
+            # to avoid split cells.
+            ([-0.1, -0.2], [0.205, 0.43]),
+            ([-0.1, -0.2, 0.05], [0.2, 0.4, 0.35]),
+        ]
+
+        @testset verbose=true "$(names[i])" for i in eachindex(names)
+            coords = coordinates[i]
+
+            nhs = GridNeighborhoodSearch{size(coords, 1)}(0.1, size(coords, 2),
+                                                          periodic_box_min_corner = periodic_boxes[i][1],
+                                                          periodic_box_max_corner = periodic_boxes[i][2])
 
             initialize_grid!(nhs, coords)
 
-            neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i],
-                                                                  nhs)))
+            neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs)))
                          for i in 1:5]
 
             # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells
@@ -198,32 +175,16 @@
             @test neighbors[3] == [1, 2, 3]
             @test neighbors[4] == [4]
             @test neighbors[5] == [1, 5]
-
-            neighbors_loop = [Int[] for _ in axes(coords, 2)]
-
-            for_particle_neighbor(coords, coords, nhs,
-                                  particles = axes(coords, 2)) do particle,
-                                                                  neighbor,
-                                                                  pos_diff,
-                                                                  distance
-                append!(neighbors_loop[particle], neighbor)
-            end
-
-            @test sort(neighbors_loop[1]) == [1, 3, 5]
-            @test sort(neighbors_loop[2]) == [2]
-            @test sort(neighbors_loop[3]) == [1, 3]
-            @test sort(neighbors_loop[4]) == [4]
-            @test sort(neighbors_loop[5]) == [1, 5]
         end
 
         @testset "Offset Domain Triggering Split Cells" begin
-            # This used to trigger a "split cell bug", where the left and right boundary
-            # cells were only partially contained in the domain.
-            # The left particle was placed inside a ghost cells, which causes it to not
-            # see the right particle, even though it is within the search distance.
+            # This test used to trigger a "split cell bug", where the left and right
+            # boundary cells were only partially contained in the domain.
+            # The left particle was placed inside a ghost cells, which caused it to not
+            # see the right particle, even though it was within the search distance.
             # The domain size is an integer multiple of the cell size, but the NHS did not
             # offset the grid based on the domain position.
-            # See https://github.com/trixi-framework/PointNeighbors.jl/pull/211
+            # See https://github.com/trixi-framework/TrixiParticles.jl/pull/211
             # for a more detailed explanation.
             coords = [-1.4 1.9
                       0.0 0.0]
@@ -244,40 +205,4 @@
             @test neighbors[2] == [1, 2]
         end
     end
-
-    @testset verbose=true "Periodicity 3D" begin
-        coords = [-0.08 0.0 0.18 0.1 -0.08
-                  -0.12 -0.05 -0.09 0.15 0.39
-                  0.14 0.34 0.12 0.06 0.13]
-
-        # 3 x 6 x 3 cells
-        nhs = GridNeighborhoodSearch{3}(0.1, size(coords, 2),
-                                        periodic_box_min_corner = [-0.1, -0.2, 0.05],
-                                        periodic_box_max_corner = [0.2, 0.4, 0.35])
-
-        initialize_grid!(nhs, coords)
-
-        neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs)))
-                     for i in 1:5]
-
-        # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells
-        @test neighbors[1] == [1, 2, 3, 5]
-        @test neighbors[2] == [1, 2, 3]
-        @test neighbors[3] == [1, 2, 3]
-        @test neighbors[4] == [4]
-        @test neighbors[5] == [1, 5]
-
-        neighbors_loop = [Int[] for _ in axes(coords, 2)]
-
-        for_particle_neighbor(coords, coords,
-                              nhs) do particle, neighbor, pos_diff, distance
-            append!(neighbors_loop[particle], neighbor)
-        end
-
-        @test sort(neighbors_loop[1]) == [1, 3, 5]
-        @test sort(neighbors_loop[2]) == [2]
-        @test sort(neighbors_loop[3]) == [1, 3]
-        @test sort(neighbors_loop[4]) == [4]
-        @test sort(neighbors_loop[5]) == [1, 5]
-    end
-end
+end;

test/runtests.jl

@@ -3,4 +3,5 @@ include("test_util.jl")
 @testset verbose=true "PointNeighbors.jl Tests" begin
     include("nhs_trivial.jl")
     include("nhs_grid.jl")
+    include("neighborhood_search.jl")
 end