trixi-framework · efaulhaber · Jul 9, 2024 · Jun 27, 2024 · Jun 27, 2024 · Jul 1, 2024
diff --git a/NEWS.md b/NEWS.md
@@ -4,7 +4,7 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](h
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions.
 
-## Version 0.2.x
+## Version 0.3.x
 
 ### Highlights
 
@@ -14,32 +14,34 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between
 
 ### Deprecated
 
-## Version 0.1.3
+## Version 0.2.0
 
-### Added
-Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
-
-## Version 0.1.2
+### Removed
+- Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
-### Added
-A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
 
-## Version 0.1.1
+## Development Cycle 0.1
 
 ### Highlights
 
 #### Discrete Element Method
 A basic implementation of the discrete element method was added.
 
-# Pre Initial Release (v0.1.0)
+#### Surface Tension and Adhesion Model
+A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
+
+#### Support for Open Boundaries
+Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
+
+## Pre Initial Release (v0.1.0)
 This section summarizes the initial features that TrixiParticles.jl was released with.
 
-## Highlights
-### EDAC
+### Highlights
+#### EDAC
 An implementation of EDAC (Entropically Damped Artificial Compressibility) was added,
 which allows for more stable simulations compared to basic WCSPH and reduces spurious pressure oscillations.
 
-### WCSPH
+#### WCSPH
 An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach.
 
 Features:
@@ -51,5 +53,5 @@ Features:
 - Density diffusion based on the models by Molteni & Colagrossi (2009), Ferrari et al. (2009) and Antuono et al. (2010).
 
 
-### TLSPH
+#### TLSPH
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.1.4-dev"
+version = "0.2.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -11,7 +11,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -35,11 +35,11 @@ DataFrames = "1.6"
 DiffEqCallbacks = "2.25, 3"
 FastPow = "0.1"
 ForwardDiff = "0.10"
-GPUArrays = "9, 10"
+GPUArraysCore = "0.1"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
-PointNeighbors = "0.2.3"
+PointNeighbors = "0.4.2"
 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.4"
 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
@@ -76,9 +76,10 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coef
 
 # ==========================================================================================
 # ==== Simulation
+# `nothing` will automatically choose the best update strategy. This is only to be able
+# to change this with `trixi_include`.
 semi = Semidiscretization(fluid_system, boundary_system,
-                          neighborhood_search=GridNeighborhoodSearch,
-                          threaded_nhs_update=true)
+                          neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing))
 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
@@ -9,7 +9,7 @@ using DataFrames: DataFrame
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect, PresetTimeCallback
 using FastPow: @fastpow
 using ForwardDiff: ForwardDiff
-using GPUArrays: AbstractGPUArray
+using GPUArraysCore: AbstractGPUArray
 using JSON: JSON
 using KernelAbstractions: KernelAbstractions, @kernel, @index
 using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
@@ -25,8 +25,11 @@ 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,
+                                ParallelUpdate, SemiParallelUpdate, SerialUpdate
+using PointNeighbors: PointNeighbors, foreach_point_neighbor, copy_neighborhood_search,
+                      @threaded
 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/gpu.jl b/src/general/gpu.jl
@@ -22,21 +22,14 @@ function Adapt.adapt_structure(to, ic::InitialCondition)
     return nothing
 end
 
-# `adapt(CuArray, ::SVector)::SVector`, but `adapt(Array, ::SVector)::Vector`.
-# We don't want to change the type of the `SVector` here.
-function Adapt.adapt_structure(to::typeof(Array), svector::SVector)
-    return svector
-end
-
-# `adapt(CuArray, ::UnitRange)::UnitRange`, but `adapt(Array, ::UnitRange)::Vector`.
-# We don't want to change the type of the `UnitRange` here.
-function Adapt.adapt_structure(to::typeof(Array), range::UnitRange)
-    return range
-end
-
 KernelAbstractions.get_backend(::PtrArray) = KernelAbstractions.CPU()
 KernelAbstractions.get_backend(system::System) = KernelAbstractions.get_backend(system.mass)
 
 function KernelAbstractions.get_backend(system::BoundarySPHSystem)
     KernelAbstractions.get_backend(system.coordinates)
 end
+
+# On GPUs, execute `f` inside a GPU kernel with KernelAbstractions.jl
+@inline function PointNeighbors.parallel_foreach(f, iterator, system::GPUSystem)
+    PointNeighbors.parallel_foreach(f, iterator, KernelAbstractions.get_backend(system))
+end