Revise periodic box API (trixi-framework#30)

* Revise periodic box API

* Allow `Vector`s in constructor of `PeriodicBox`

* Fix tests

* Fix docs for `TrivialNeighborhoodSearch`

Author: efaulhaber

src/PointNeighbors.jl

@@ -13,7 +13,7 @@ include("cell_lists/cell_lists.jl")
 include("nhs_grid.jl")
 include("nhs_neighbor_lists.jl")
 
-export foreach_point_neighbor, foreach_neighbor
+export foreach_point_neighbor, foreach_neighbor, PeriodicBox
 export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch
 export initialize!, update!, initialize_grid!, update_grid!
 

src/neighborhood_search.jl

@@ -41,14 +41,26 @@ See also [`initialize!`](@ref).
     return search
 end
 
+"""
+    PeriodicBox(; min_corner, max_corner)
+
+Define a rectangular periodic domain.
+
+# Keywords
+- `min_corner`: Coordinates of the domain corner in negative coordinate directions.
+- `max_corner`: Coordinates of the domain corner in positive coordinate directions.
+"""
 struct PeriodicBox{NDIMS, ELTYPE}
     min_corner :: SVector{NDIMS, ELTYPE}
     max_corner :: SVector{NDIMS, ELTYPE}
     size       :: SVector{NDIMS, ELTYPE}
 
-    function PeriodicBox(min_corner, max_corner)
-        new{length(min_corner), eltype(min_corner)}(min_corner, max_corner,
-                                                    max_corner - min_corner)
+    function PeriodicBox(; min_corner, max_corner)
+        min_corner_ = SVector(Tuple(min_corner))
+        max_corner_ = SVector(Tuple(max_corner))
+
+        new{length(min_corner), eltype(min_corner)}(min_corner_, max_corner_,
+                                                    max_corner_ - min_corner_)
     end
 end
 

src/nhs_grid.jl

@@ -1,6 +1,6 @@
 @doc raw"""
-    GridNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner=nothing,
-                                  periodic_box_max_corner=nothing, threaded_nhs_update=true)
+    GridNeighborhoodSearch{NDIMS}(search_radius, n_points;
+                                  periodic_box = nothing, threaded_nhs_update = true)
 
 Simple grid-based neighborhood search with uniform search radius.
 The domain is divided into a regular grid.
@@ -28,15 +28,12 @@ since not sorting makes our implementation a lot faster (although less paralleli
 - `n_points`:    Total number of points.
 
 # Keywords
-- `periodic_box_min_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in negative coordinate
-                                directions.
-- `periodic_box_max_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in positive coordinate
-                                directions.
-- `threaded_nhs_update=true`:              Can be used to deactivate thread parallelization in the neighborhood search update.
-                                This can be one of the largest sources of variations between simulations
-                                with different thread numbers due to neighbor ordering changes.
+- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a
+                            [`PeriodicBox`](@ref).
+- `threaded_nhs_update = true`: Can be used to deactivate thread parallelization in the
+                                neighborhood search update. This can be one of the largest
+                                sources of variations between simulations with different
+                                thread numbers due to neighbor ordering changes.
 
 ## References
 - M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán.
@@ -59,25 +56,19 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear
     threaded_nhs_update :: Bool
 
     function GridNeighborhoodSearch{NDIMS}(search_radius, n_points;
-                                           periodic_box_min_corner = nothing,
-                                           periodic_box_max_corner = nothing,
+                                           periodic_box = nothing,
                                            threaded_nhs_update = true) where {NDIMS}
         ELTYPE = typeof(search_radius)
         cell_list = DictionaryCellList{NDIMS}()
 
         cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_points, Threads.nthreads())
         cell_buffer_indices = zeros(Int, Threads.nthreads())
 
-        if search_radius < eps() ||
-           (periodic_box_min_corner === nothing && periodic_box_max_corner === nothing)
-
+        if search_radius < eps() || isnothing(periodic_box)
             # No periodicity
-            periodic_box = nothing
             n_cells = ntuple(_ -> -1, Val(NDIMS))
             cell_size = ntuple(_ -> search_radius, Val(NDIMS))
-        elseif periodic_box_min_corner !== nothing && periodic_box_max_corner !== nothing
-            periodic_box = PeriodicBox(periodic_box_min_corner, periodic_box_max_corner)
-
+        else
             # Round up search radius so that the grid fits exactly into the domain without
             # splitting any cells. This might impact performance slightly, since larger
             # cells mean that more potential neighbors are considered than necessary.
@@ -91,9 +82,6 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear
                                     "in each dimension when used with periodicity. " *
                                     "Please use no NHS for very small problems."))
             end
-        else
-            throw(ArgumentError("`periodic_box_min_corner` and `periodic_box_max_corner` " *
-                                "must either be both `nothing` or both an array or tuple"))
         end
 
         new{NDIMS, ELTYPE, typeof(cell_list),
@@ -103,7 +91,7 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear
     end
 end
 
-@inline Base.ndims(neighborhood_search::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS
+@inline Base.ndims(::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS
 
 @inline function npoints(neighborhood_search::GridNeighborhoodSearch)
     return size(neighborhood_search.cell_buffer, 1)
@@ -321,13 +309,8 @@ end
 
 # Create a copy of a neighborhood search but with a different search radius
 function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, x, y)
-    if nhs.periodic_box === nothing
-        search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs))
-    else
-        search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs),
-                                                    periodic_box_min_corner = nhs.periodic_box.min_corner,
-                                                    periodic_box_max_corner = nhs.periodic_box.max_corner)
-    end
+    search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs),
+                                                periodic_box = nhs.periodic_box)
 
     # Initialize neighborhood search
     initialize!(search, x, y)

src/nhs_neighbor_lists.jl

@@ -1,7 +1,6 @@
 @doc raw"""
     PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points;
-                                         periodic_box_min_corner = nothing,
-                                         periodic_box_max_corner = nothing)
+                                         periodic_box = nothing)
 
 Neighborhood search with precomputed neighbor lists. A list of all neighbors is computed
 for each point during initialization and update.
@@ -17,32 +16,24 @@ initialization and update.
 - `n_points`:    Total number of points.
 
 # Keywords
-- `periodic_box_min_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in negative coordinate
-                                directions.
-- `periodic_box_max_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in positive coordinate
-                                directions.
+- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a
+                            [`PeriodicBox`](@ref).
 """
 struct PrecomputedNeighborhoodSearch{NDIMS, NHS, NL, PB}
     neighborhood_search :: NHS
     neighbor_lists      :: NL
     periodic_box        :: PB
 
     function PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points;
-                                                  periodic_box_min_corner = nothing,
-                                                  periodic_box_max_corner = nothing) where {
-                                                                                            NDIMS
-                                                                                            }
+                                                  periodic_box = nothing) where {NDIMS}
         nhs = GridNeighborhoodSearch{NDIMS}(search_radius, n_points,
-                                            periodic_box_min_corner = periodic_box_min_corner,
-                                            periodic_box_max_corner = periodic_box_max_corner)
+                                            periodic_box = periodic_box)
 
         neighbor_lists = Vector{Vector{Int}}()
 
         new{NDIMS, typeof(nhs),
             typeof(neighbor_lists),
-            typeof(nhs.periodic_box)}(nhs, neighbor_lists, nhs.periodic_box)
+            typeof(periodic_box)}(nhs, neighbor_lists, periodic_box)
     end
 end
 

src/nhs_trivial.jl

@@ -1,53 +1,30 @@
 @doc raw"""
-    TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint)
+    TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint, periodic_box = nothing)
 
 Trivial neighborhood search that simply loops over all points.
 
 # Arguments
 - `NDIMS`:          Number of dimensions.
 - `search_radius`:  The uniform search radius.
-- `eachpoint`:   `UnitRange` of all point indices. Usually just `1:n_points`.
+- `eachparticle`:   Iterator for all point indices. Usually just `1:n_points`.
 
 # Keywords
-- `periodic_box_min_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in negative coordinate
-                                directions.
-- `periodic_box_max_corner`:    In order to use a (rectangular) periodic domain, pass the
-                                coordinates of the domain corner in positive coordinate
-                                directions.
+- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a
+                            [`PeriodicBox`](@ref).
 """
 struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodSearch
     search_radius :: ELTYPE
     eachpoint     :: EP
     periodic_box  :: PB
 
     function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint;
-                                              periodic_box_min_corner = nothing,
-                                              periodic_box_max_corner = nothing) where {
-                                                                                        NDIMS
-                                                                                        }
-        if search_radius < eps() ||
-           (periodic_box_min_corner === nothing && periodic_box_max_corner === nothing)
-
-            # No periodicity
-            periodic_box = nothing
-        elseif periodic_box_min_corner !== nothing && periodic_box_max_corner !== nothing
-            periodic_box = PeriodicBox(periodic_box_min_corner, periodic_box_max_corner)
-        else
-            throw(ArgumentError("`periodic_box_min_corner` and `periodic_box_max_corner` " *
-                                "must either be both `nothing` or both an array or tuple"))
-        end
-
+                                              periodic_box = nothing) where {NDIMS}
         new{NDIMS, typeof(search_radius),
             typeof(eachpoint), typeof(periodic_box)}(search_radius, eachpoint, periodic_box)
     end
 end
 
-@inline function Base.ndims(neighborhood_search::TrivialNeighborhoodSearch{NDIMS}) where {
-                                                                                          NDIMS
-                                                                                          }
-    return NDIMS
-end
+@inline Base.ndims(::TrivialNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS
 
 @inline initialize!(search::TrivialNeighborhoodSearch, x, y) = search
 

test/neighborhood_search.jl

@@ -24,11 +24,11 @@
         ]
 
         periodic_boxes = [
-            ([-0.1, -0.2], [0.2, 0.4]),
+            PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [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]),
+            PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]),
+            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]),
         ]
 
         @testset verbose=true "$(names[i])" for i in eachindex(names)
@@ -40,14 +40,11 @@
 
             neighborhood_searches = [
                 TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_points,
-                                                 periodic_box_min_corner = periodic_boxes[i][1],
-                                                 periodic_box_max_corner = periodic_boxes[i][2]),
+                                                 periodic_box = periodic_boxes[i]),
                 GridNeighborhoodSearch{NDIMS}(search_radius, n_points,
-                                              periodic_box_min_corner = periodic_boxes[i][1],
-                                              periodic_box_max_corner = periodic_boxes[i][2]),
+                                              periodic_box = periodic_boxes[i]),
                 PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points,
-                                                     periodic_box_min_corner = periodic_boxes[i][1],
-                                                     periodic_box_max_corner = periodic_boxes[i][2]),
+                                                     periodic_box = periodic_boxes[i]),
             ]
             neighborhood_searches_names = [
                 "`TrivialNeighborhoodSearch`",

test/nhs_grid.jl

@@ -148,19 +148,18 @@
         ]
 
         periodic_boxes = [
-            ([-0.1, -0.2], [0.2, 0.4]),
+            PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [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]),
+            PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]),
+            PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [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])
+                                                          periodic_box = periodic_boxes[i])
 
             initialize_grid!(nhs, coords)
 
@@ -189,13 +188,18 @@
 
             # 5 x 1 cells
             nhs = GridNeighborhoodSearch{2}(1.0, size(coords, 2),
-                                            periodic_box_min_corner = [-1.5, 0.0],
-                                            periodic_box_max_corner = [2.5, 3.0])
+                                            periodic_box = PeriodicBox(min_corner = [
+                                                                           -1.5,
+                                                                           0.0,
+                                                                       ],
+                                                                       max_corner = [
+                                                                           2.5,
+                                                                           3.0,
+                                                                       ]))
 
             initialize_grid!(nhs, coords)
 
-            neighbors = [sort(unique(collect(PointNeighbors.eachneighbor(coords[:,
-                                                                                i],
+            neighbors = [sort(unique(collect(PointNeighbors.eachneighbor(coords[:, i],
                                                                          nhs))))
                          for i in 1:2]