Merge branch 'dev' into release-0.2

trixi-framework · Jun 27, 2024 · 0382fc8 · 0382fc8
2 parents cb84f30 + 2326d03
commit 0382fc8
Show file tree

Hide file tree

Showing 36 changed files with 213 additions and 209 deletions.
diff --git a/Project.toml b/Project.toml
@@ -39,7 +39,7 @@ GPUArrays = "9, 10"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
-PointNeighbors = "0.2.3"
+PointNeighbors = "0.3"
 Polyester = "0.7.5"
 RecipesBase = "1"
 Reexport = "1"

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,9 +1,11 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 Documenter = "1"
 OrdinaryDiffEq = "6"
+PointNeighbors = "0.3"
 TrixiBase = "0.1"
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,6 +1,7 @@
 using Documenter
 using TrixiParticles
 using TrixiBase
+using PointNeighbors
 
 # Get TrixiParticles.jl root directory
 trixiparticles_root_dir = dirname(@__DIR__)
@@ -129,6 +130,7 @@ makedocs(sitename="TrixiParticles.jl",
                  "Time Integration" => "time_integration.md",
                  "Callbacks" => "callbacks.md",
                  "TrixiBase.jl API Reference" => "reference-trixibase.md",
+                 "PointNeighbors.jl API Reference" => "reference-pointneighbors.md",
              ],
              "Authors" => "authors.md",
              "Contributing" => "contributing.md",

diff --git a/docs/src/general/neighborhood_search.md b/docs/src/general/neighborhood_search.md
@@ -6,31 +6,30 @@ We provide several implementations in the package
 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):
+    To run a simulation with a neighborhood search implementation, pass a template of the
+    neighborhood search to the constructor of the [`Semidiscretization`](@ref).
+    A template is just an empty neighborhood search with search radius `0.0`.
+    See [`copy_neighborhood_search`](@ref) and the examples below for more details.
     ```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)
+                              neighborhood_search=PrecomputedNeighborhoodSearch{2}())
 
     # output
     ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
     │ Semidiscretization                                                                               │
     │ ══════════════════                                                                               │
     │ #spatial dimensions: ………………………… 2                                                                │
     │ #systems: ……………………………………………………… 2                                                                │
-    │ neighborhood search: ………………………… GridNeighborhoodSearch                                           │
+    │ neighborhood search: ………………………… PrecomputedNeighborhoodSearch                                    │
     │ 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.
+    The keyword argument `periodic_box` in the neighborhood search constructors can be used
+    to define a periodic domain. See the PointNeighbors.jl docs for more details.
     ```jldoctest semi_example; output = false
+    periodic_box = PeriodicBox(min_corner=[0.0, -0.25], max_corner=[1.0, 0.75])
     semi = Semidiscretization(system1, system2,
-                              neighborhood_search=GridNeighborhoodSearch,
-                              periodic_box_min_corner=[0.0, -0.25],
-                              periodic_box_max_corner=[1.0, 0.75])
+                              neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box))
 
     # output
     ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐

diff --git a/docs/src/reference-pointneighbors.md b/docs/src/reference-pointneighbors.md
@@ -0,0 +1,9 @@
+# PointNeighbors.jl API
+
+```@meta
+CurrentModule = PointNeighbors
+```
+
+```@autodocs
+Modules = [PointNeighbors]
+```
diff --git a/examples/dem/rectangular_tank_2d.jl b/examples/dem/rectangular_tank_2d.jl
@@ -30,8 +30,7 @@ boundary_system = BoundaryDEMSystem(tank.boundary, 10e7)
 # ==========================================================================================
 # ==== Simulation
 
-semi = Semidiscretization(rock_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch)
+semi = Semidiscretization(rock_system, boundary_system)
 
 tspan = (0.0, 5.0)
 ode = semidiscretize(semi, tspan)

diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -74,8 +74,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coef
 # ==========================================================================================
 # ==== Simulation
 semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch,
-                          threaded_nhs_update=true)
+                          neighborhood_search=GridNeighborhoodSearch{2}(threaded_update=true))
 ode = semidiscretize(semi, tspan, data_type=nothing)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl
@@ -56,9 +56,10 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system,
-                          periodic_box_min_corner=[0.0, -0.25],
-                          periodic_box_max_corner=[1.0, 0.75])
+periodic_box = PeriodicBox(min_corner=[0.0, -0.25], max_corner=[1.0, 0.75])
+neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box)
+
+semi = Semidiscretization(fluid_system, boundary_system; neighborhood_search)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)

diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl
@@ -2,7 +2,7 @@ using TrixiParticles
 using LinearAlgebra
 
 # The second type parameter of `System` can't be `Nothing`, or TrixiParticles will launch
-# GPU kernel for `for_particle_neighbor` loops.
+# GPU kernel for `foreach_point_neighbor` loops.
 struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, 0}
     initial_condition :: InitialCondition{ELTYPE}
     mass              :: Array{ELTYPE, 1} # [particle]
@@ -39,7 +39,7 @@ function TrixiParticles.update_nhs!(neighborhood_search,
                                     u_system, u_neighbor)
     TrixiParticles.PointNeighbors.update!(neighborhood_search,
                                           u_system, u_neighbor,
-                                          particles_moving=(true, true))
+                                          points_moving=(true, true))
 end
 
 function TrixiParticles.compact_support(system::NBodySystem,
@@ -59,10 +59,10 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system,
     neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, neighbor_system)
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff.
-    TrixiParticles.for_particle_neighbor(particle_system, neighbor_system,
-                                         system_coords, neighbor_coords,
-                                         neighborhood_search) do particle, neighbor,
-                                                                 pos_diff, distance
+    TrixiParticles.foreach_point_neighbor(particle_system, neighbor_system,
+                                          system_coords, neighbor_coords,
+                                          neighborhood_search) do particle, neighbor,
+                                                                  pos_diff, distance
         # Only consider particles with a distance > 0.
         distance < sqrt(eps()) && return
 

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -25,8 +25,9 @@ using StrideArrays: PtrArray, StaticInt
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
 using TrixiBase: trixi_include, @trixi_timeit, timer, timeit_debug_enabled,
                  disable_debug_timings, enable_debug_timings
-@reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch
-using PointNeighbors: PointNeighbors, for_particle_neighbor
+@reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch,
+                                PrecomputedNeighborhoodSearch, PeriodicBox
+using PointNeighbors: PointNeighbors, foreach_point_neighbor, copy_neighborhood_search
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 
 # `util.jl` depends on the `GPUSystem` type defined in `system.jl`

diff --git a/src/general/corrections.jl b/src/general/corrections.jl
@@ -193,10 +193,10 @@ function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi,
         neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff
-        for_particle_neighbor(system, neighbor_system, system_coords,
-                              neighbor_coords, neighborhood_search,
-                              particles=eachparticle(system)) do particle, neighbor,
-                                                                 pos_diff, distance
+        foreach_point_neighbor(system, neighbor_system, system_coords,
+                               neighbor_coords, neighborhood_search,
+                               points=eachparticle(system)) do particle, neighbor,
+                                                               pos_diff, distance
             rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
             m_b = hydrodynamic_mass(neighbor_system, neighbor)
             volume = m_b / rho_b
@@ -255,11 +255,10 @@ function compute_correction_values!(system,
         neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff
-        for_particle_neighbor(system, neighbor_system, system_coords,
-                              neighbor_coords,
-                              neighborhood_search,
-                              particles=eachparticle(system)) do particle, neighbor,
-                                                                 pos_diff, distance
+        foreach_point_neighbor(system, neighbor_system, system_coords,
+                               neighbor_coords, neighborhood_search,
+                               points=eachparticle(system)) do particle, neighbor,
+                                                               pos_diff, distance
             rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
             m_b = hydrodynamic_mass(neighbor_system, neighbor)
             volume = m_b / rho_b
@@ -361,11 +360,9 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search,
     set_zero!(corr_matrix)
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff.
-    for_particle_neighbor(system, system,
-                          coordinates, coordinates,
-                          neighborhood_search;
-                          particles=eachparticle(system)) do particle, neighbor,
-                                                             pos_diff, distance
+    foreach_point_neighbor(system, system, coordinates, coordinates, neighborhood_search;
+                           points=eachparticle(system)) do particle, neighbor,
+                                                           pos_diff, distance
         volume = mass[neighbor] / density_fun(neighbor)
 
         grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
@@ -397,12 +394,10 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
         neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
-        for_particle_neighbor(system, neighbor_system, coordinates, neighbor_coords,
-                              neighborhood_search;
-                              particles=eachparticle(system)) do particle,
-                                                                 neighbor,
-                                                                 pos_diff,
-                                                                 distance
+        foreach_point_neighbor(system, neighbor_system, coordinates, neighbor_coords,
+                               neighborhood_search;
+                               points=eachparticle(system)) do particle, neighbor,
+                                                               pos_diff, distance
             volume = hydrodynamic_mass(neighbor_system, neighbor) /
                      particle_density(v_neighbor_system, neighbor_system, neighbor)
 

diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl
@@ -62,8 +62,9 @@ function summation_density!(system, semi, u, u_ode, density;
         nhs = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff.
-        for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, nhs,
-                              particles=particles) do particle, neighbor, pos_diff, distance
+        foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, nhs,
+                               points=particles) do particle, neighbor,
+                                                    pos_diff, distance
             mass = hydrodynamic_mass(neighbor_system, neighbor)
             density[particle] += mass * smoothing_kernel(system, distance)
         end

diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
@@ -290,11 +290,12 @@ function find_too_close_particles(coords1, coords2, max_distance)
     NDIMS = size(coords1, 1)
     result = Int[]
 
-    nhs = GridNeighborhoodSearch{NDIMS}(max_distance, size(coords2, 2))
+    nhs = GridNeighborhoodSearch{NDIMS}(search_radius=max_distance,
+                                        n_points=size(coords2, 2))
     PointNeighbors.initialize!(nhs, coords1, coords2)
 
     # We are modifying the vector `result`, so this cannot be parallel
-    for_particle_neighbor(coords1, coords2, nhs, parallel=false) do particle, _, _, _
+    foreach_point_neighbor(coords1, coords2, nhs, parallel=false) do particle, _, _, _
         if !(particle in result)
             append!(result, particle)
         end

diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl
@@ -495,7 +495,9 @@ function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length
             system_coords = current_coordinates(u, system)
             old_nhs = get_neighborhood_search(ref_system, system, semi)
             nhs = PointNeighbors.copy_neighborhood_search(old_nhs, search_radius,
-                                                          system_coords, system_coords)
+                                                          nparticles(system))
+            PointNeighbors.initialize!(nhs, system_coords, system_coords)
+
             return nhs
         end
     end
@@ -534,7 +536,7 @@ end
 
         system_coords = current_coordinates(u, system)
 
-        # This is basically `for_particle_neighbor` unrolled
+        # This is basically `foreach_point_neighbor` unrolled
         for particle in PointNeighbors.eachneighbor(point_coords, nhs)
             coords = extract_svector(system_coords, Val(ndims(system)), particle)
 

diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl
@@ -1,22 +1,22 @@
 # 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)
+function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system,
+                                               system_coords, neighbor_coords,
+                                               neighborhood_search;
+                                               points=each_moving_particle(system),
+                                               parallel=true)
+    foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search;
+                           points, parallel)
 end
 
-function PointNeighbors.for_particle_neighbor(f, system::GPUSystem, neighbor_system,
-                                              system_coords, neighbor_coords,
-                                              neighborhood_search;
-                                              particles=each_moving_particle(system),
-                                              parallel=true)
-    @threaded system for particle in particles
+function PointNeighbors.foreach_point_neighbor(f, system::GPUSystem, neighbor_system,
+                                               system_coords, neighbor_coords,
+                                               neighborhood_search;
+                                               points=each_moving_particle(system),
+                                               parallel=true)
+    @threaded system for point in points
         PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords,
-                                        neighborhood_search, particle)
+                                        neighborhood_search, point)
     end
 end