Merge pull request #33 from trixi-framework/dev

Merge breaking changes from `dev`

Project.toml

@@ -1,7 +1,7 @@
 name = "PointNeighbors"
 uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 authors = ["Erik Faulhaber <erik.faulhaber@uni-koeln.de>"]
-version = "0.2.4-dev"
+version = "0.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/PointNeighbors.jl

@@ -11,10 +11,11 @@ include("neighborhood_search.jl")
 include("nhs_trivial.jl")
 include("cell_lists/cell_lists.jl")
 include("nhs_grid.jl")
-include("nhs_neighbor_lists.jl")
+include("nhs_precomputed.jl")
 
-export for_particle_neighbor, foreach_neighbor
+export foreach_point_neighbor, foreach_neighbor
 export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch
 export initialize!, update!, initialize_grid!, update_grid!
+export PeriodicBox, copy_neighborhood_search
 
 end # module PointNeighbors

src/cell_lists/dictionary.jl

@@ -16,13 +16,13 @@ function Base.empty!(cell_list::DictionaryCellList)
     return cell_list
 end
 
-function push_cell!(cell_list::DictionaryCellList, cell, particle)
+function push_cell!(cell_list::DictionaryCellList, cell, point)
     (; hashtable) = cell_list
 
     if haskey(hashtable, cell)
-        append!(hashtable[cell], particle)
+        append!(hashtable[cell], point)
     else
-        hashtable[cell] = [particle]
+        hashtable[cell] = [point]
     end
 
     return cell_list

src/neighborhood_search.jl

@@ -15,7 +15,7 @@ See also [`update!`](@ref).
 @inline initialize!(search::AbstractNeighborhoodSearch, x, y) = search
 
 """
-    update!(search::AbstractNeighborhoodSearch, x, y; particles_moving = (true, true))
+    update!(search::AbstractNeighborhoodSearch, x, y; points_moving = (true, true))
 
 Update an already initialized neighborhood search with the two coordinate arrays `x` and `y`.
 
@@ -29,70 +29,112 @@ If incremental updates are not possible for an implementation, `update!` will fa
 to a regular `initialize!`.
 
 Some neighborhood searches might not need to update when only `x` changed since the last
-`update!` or `initialize!` and `y` did not change. Pass `particles_moving = (true, false)`
+`update!` or `initialize!` and `y` did not change. Pass `points_moving = (true, false)`
 in this case to avoid unnecessary updates.
-The first flag in `particles_moving` indicates if points in `x` are moving.
+The first flag in `points_moving` indicates if points in `x` are moving.
 The second flag indicates if points in `y` are moving.
 
 See also [`initialize!`](@ref).
 """
 @inline function update!(search::AbstractNeighborhoodSearch, x, y;
-                         particles_moving = (true, true))
+                         points_moving = (true, true))
     return search
 end
 
+"""
+    copy_neighborhood_search(search::AbstractNeighborhoodSearch, search_radius, n_points;
+                             eachpoint = 1:n_points)
+
+Create a new **uninitialized** neighborhood search of the same type and with the same
+configuration options as `search`, but with a different search radius and number of points.
+
+The [`TrivialNeighborhoodSearch`](@ref) also requires an iterator `eachpoint`, which most
+of the time will be `1:n_points`. If the `TrivialNeighborhoodSearch` is never going to be
+used, the keyword argument `eachpoint` can be ignored.
+
+This is useful when a simulation code requires multiple neighborhood searches of the same
+kind. One can then just pass an empty neighborhood search as a template and use
+this function inside the simulation code to generate similar neighborhood searches with
+different search radii and different numbers of points.
+```jldoctest; filter = r"GridNeighborhoodSearch{2,.*"
+# Template
+nhs = GridNeighborhoodSearch{2}()
+
+# Inside the simulation code, generate similar neighborhood searches
+nhs1 = copy_neighborhood_search(nhs, 1.0, 100)
+
+# output
+GridNeighborhoodSearch{2, Float64, ...}(...)
+```
+"""
+@inline function copy_neighborhood_search(search::AbstractNeighborhoodSearch,
+                                          search_radius, n_points; eachpoint = 1:n_points)
+    return nothing
+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
 
 # The type annotation is to make Julia specialize on the type of the function.
 # Otherwise, unspecialized code will cause a lot of allocations
 # and heavily impact performance.
 # See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
-function for_particle_neighbor(f::T, system_coords, neighbor_coords, neighborhood_search;
-                               particles = axes(system_coords, 2),
-                               parallel = true) where {T}
-    for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particles,
-                          Val(parallel))
+function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborhood_search;
+                                points = axes(system_coords, 2),
+                                parallel = true) where {T}
+    foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search, points,
+                           Val(parallel))
 end
 
-@inline function for_particle_neighbor(f, system_coords, neighbor_coords,
-                                       neighborhood_search, particles, parallel::Val{true})
-    @threaded for particle in particles
-        foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particle)
+@inline function foreach_point_neighbor(f, system_coords, neighbor_coords,
+                                        neighborhood_search, points, parallel::Val{true})
+    @threaded for point in points
+        foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, point)
     end
 
     return nothing
 end
 
-@inline function for_particle_neighbor(f, system_coords, neighbor_coords,
-                                       neighborhood_search, particles, parallel::Val{false})
-    for particle in particles
-        foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particle)
+@inline function foreach_point_neighbor(f, system_coords, neighbor_coords,
+                                        neighborhood_search, points, parallel::Val{false})
+    for point in points
+        foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, point)
     end
 
     return nothing
 end
 
 @inline function foreach_neighbor(f, system_coords, neighbor_system_coords,
-                                  neighborhood_search, particle;
+                                  neighborhood_search, point;
                                   search_radius = neighborhood_search.search_radius)
     (; periodic_box) = neighborhood_search
 
-    particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)),
-                                      particle)
-    for neighbor in eachneighbor(particle_coords, neighborhood_search)
+    point_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), point)
+    for neighbor in eachneighbor(point_coords, neighborhood_search)
         neighbor_coords = extract_svector(neighbor_system_coords,
                                           Val(ndims(neighborhood_search)), neighbor)
 
-        pos_diff = particle_coords - neighbor_coords
+        pos_diff = point_coords - neighbor_coords
         distance2 = dot(pos_diff, pos_diff)
 
         pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius,
@@ -102,8 +144,8 @@ end
             distance = sqrt(distance2)
 
             # Inline to avoid loss of performance
-            # compared to not using `for_particle_neighbor`.
-            @inline f(particle, neighbor, pos_diff, distance)
+            # compared to not using `foreach_point_neighbor`.
+            @inline f(point, neighbor, pos_diff, distance)
         end
     end
 end
@@ -113,8 +155,7 @@ end
     return pos_diff, distance2
 end
 
-@inline function compute_periodic_distance(pos_diff, distance2, search_radius,
-                                           periodic_box)
+@inline function compute_periodic_distance(pos_diff, distance2, search_radius, periodic_box)
     if distance2 > search_radius^2
         # Use periodic `pos_diff`
         pos_diff -= periodic_box.size .* round.(pos_diff ./ periodic_box.size)