Move neighborhood search to a new standalone package

.github/workflows/Documenter.yml

@@ -8,10 +8,12 @@ on:
     paths-ignore:
       - '.github/workflows/ci.yml'
       - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
   pull_request:
     paths-ignore:
       - '.github/workflows/ci.yml'
       - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
   workflow_dispatch:
 
 concurrency:

Project.toml

@@ -13,16 +13,15 @@ FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Morton = "2a6d852e-3fac-5a38-885c-fe708af2d09e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
-ThreadingUtilities = "8290d209-cae3-49c0-8002-c8c24d57dab5"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
@@ -35,16 +34,15 @@ DiffEqCallbacks = "2.25, 3"
 FastPow = "0.1"
 ForwardDiff = "0.10"
 JSON = "0.21"
-Morton = "0.1"
 MuladdMacro = "0.2"
 Polyester = "0.7.5"
 RecipesBase = "1"
 Reexport = "1"
 SciMLBase = "1, 2"
 StaticArrays = "1"
 StrideArrays = "0.1"
-ThreadingUtilities = "0.5"
 TimerOutputs = "0.5"
 TrixiBase = "0.1"
+PointNeighbors = "0.2"
 WriteVTK = "1"
 julia = "1.9"

docs/src/general/neighborhood_search.md

@@ -1,6 +1,44 @@
 # Neighborhood Search
 
-```@autodocs
-Modules = [TrixiParticles]
-Pages = map(file -> joinpath("neighborhood_search", file), readdir(joinpath("..", "src", "neighborhood_search")))
-```
+The neighborhood search is the most essential component for performance.
+We provide several implementations in the package
+[PointNeighbors.jl](https://github.com/trixi-framework/PointNeighbors.jl).
+See the docs of this package for an overview and a comparison of different implementations.
+
+!!! note "Usage"
+    To run a simulation with a neighborhood search implementation, just pass the type
+    to the constructor of the [`Semidiscretization`](@ref):
+    ```jldoctest semi_example; output=false, setup = :(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); system1 = fluid_system; system2 = boundary_system)
+    semi = Semidiscretization(system1, system2,
+                              neighborhood_search=GridNeighborhoodSearch)
+
+    # output
+    ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+    │ Semidiscretization                                                                               │
+    │ ══════════════════                                                                               │
+    │ #spatial dimensions: ………………………… 2                                                                │
+    │ #systems: ……………………………………………………… 2                                                                │
+    │ neighborhood search: ………………………… GridNeighborhoodSearch                                           │
+    │ total #particles: ………………………………… 636                                                              │
+    └──────────────────────────────────────────────────────────────────────────────────────────────────┘
+    ```
+    The keyword arguments `periodic_box_min_corner` and `periodic_box_max_corner` mentioned
+    in the PointNeighbors.jl docs can also be passed to the
+    [`Semidiscretization`](@ref) and will internally be forwarded to the neighborhood search.
+    See the docs of [`Semidiscretization`](@ref) for more details.
+    ```jldoctest semi_example; output = false
+    semi = Semidiscretization(system1, system2,
+                              neighborhood_search=GridNeighborhoodSearch,
+                              periodic_box_min_corner=[0.0, -0.25],
+                              periodic_box_max_corner=[1.0, 0.75])
+
+    # output
+    ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+    │ Semidiscretization                                                                               │
+    │ ══════════════════                                                                               │
+    │ #spatial dimensions: ………………………… 2                                                                │
+    │ #systems: ……………………………………………………… 2                                                                │
+    │ neighborhood search: ………………………… GridNeighborhoodSearch                                           │
+    │ total #particles: ………………………………… 636                                                              │
+    └──────────────────────────────────────────────────────────────────────────────────────────────────┘
+    ```

examples/n_body/n_body_system.jl

@@ -31,9 +31,12 @@ function TrixiParticles.write_v0!(v0, system::NBodySystem)
 end
 
 # NHS update
-function TrixiParticles.nhs_coords(system::NBodySystem,
-                                   neighbor::NBodySystem, u)
-    return u
+function TrixiParticles.update_nhs!(neighborhood_search,
+                                    system::NBodySystem, neighbor::NBodySystem,
+                                    u_system, u_neighbor)
+    TrixiParticles.PointNeighbors.update!(neighborhood_search,
+                                          u_system, u_neighbor,
+                                          particles_moving=(true, true))
 end
 
 function TrixiParticles.compact_support(system::NBodySystem,

src/TrixiParticles.jl

@@ -11,7 +11,6 @@ using FastPow: @fastpow
 using ForwardDiff: ForwardDiff
 using JSON: JSON
 using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
-using Morton: cartesian2morton
 using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
@@ -21,16 +20,16 @@ using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
-using ThreadingUtilities
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
 using TrixiBase: trixi_include
+@reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch
+using PointNeighbors: PointNeighbors, for_particle_neighbor
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 
 # util needs to be first because of macro @trixi_timeit
 include("util.jl")
 include("callbacks/callbacks.jl")
 include("general/general.jl")
-include("neighborhood_search/neighborhood_search.jl")
 include("setups/setups.jl")
 include("schemes/schemes.jl")
 
@@ -59,7 +58,6 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
        PressureMirroring, PressureZeroing
 export BoundaryMovement
-export GridNeighborhoodSearch, TrivialNeighborhoodSearch
 export examples_dir, validation_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape

src/general/general.jl

@@ -26,3 +26,4 @@ include("system.jl")
 include("interpolation.jl")
 include("file_system.jl")
 include("custom_quantities.jl")
+include("neighborhood_search.jl")

src/general/initial_condition.jl

@@ -291,11 +291,10 @@ function find_too_close_particles(coords1, coords2, max_distance)
     result = Int[]
 
     nhs = GridNeighborhoodSearch{NDIMS}(max_distance, size(coords2, 2))
-    TrixiParticles.initialize!(nhs, coords2)
+    PointNeighbors.initialize!(nhs, coords1, coords2)
 
     # We are modifying the vector `result`, so this cannot be parallel
-    TrixiParticles.for_particle_neighbor(coords1, coords2, nhs,
-                                         parallel=false) do particle, _, _, _
+    for_particle_neighbor(coords1, coords2, nhs, parallel=false) do particle, _, _, _
         if !(particle in result)
             append!(result, particle)
         end

src/general/interpolation.jl

@@ -485,7 +485,7 @@ end
 function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length)
     if isapprox(smoothing_length, ref_system.smoothing_length)
         # Update existing NHS
-        update_nhs(u_ode, semi)
+        update_nhs!(semi, u_ode)
         neighborhood_searches = semi.neighborhood_searches[system_indices(ref_system, semi)]
     else
         ref_smoothing_kernel = ref_system.smoothing_kernel
@@ -494,7 +494,8 @@ function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length
             u = wrap_u(u_ode, system, semi)
             system_coords = current_coordinates(u, system)
             old_nhs = get_neighborhood_search(ref_system, system, semi)
-            nhs = copy_neighborhood_search(old_nhs, search_radius, system_coords)
+            nhs = PointNeighbors.copy_neighborhood_search(old_nhs, search_radius,
+                                                          system_coords, system_coords)
             return nhs
         end
     end
@@ -534,13 +535,15 @@ end
         system_coords = current_coordinates(u, system)
 
         # This is basically `for_particle_neighbor` unrolled
-        for particle in eachneighbor(point_coords, nhs)
+        for particle in PointNeighbors.eachneighbor(point_coords, nhs)
             coords = extract_svector(system_coords, Val(ndims(system)), particle)
 
             pos_diff = point_coords - coords
             distance2 = dot(pos_diff, pos_diff)
-            pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2,
-                                                            search_radius, periodic_box)
+            pos_diff, distance2 = PointNeighbors.compute_periodic_distance(pos_diff,
+                                                                           distance2,
+                                                                           search_radius,
+                                                                           periodic_box)
             if distance2 > search_radius^2
                 continue
             end

src/general/neighborhood_search.jl

@@ -0,0 +1,11 @@
+# Loop over all pairs of particles and neighbors within the kernel cutoff.
+# `f(particle, neighbor, pos_diff, distance)` is called for every particle-neighbor pair.
+# By default, loop over `each_moving_particle(system)`.
+function PointNeighbors.for_particle_neighbor(f, system, neighbor_system,
+                                              system_coords, neighbor_coords,
+                                              neighborhood_search;
+                                              particles=each_moving_particle(system),
+                                              parallel=true)
+    for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search,
+                          particles=particles, parallel=parallel)
+end