diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 941fef735..4cc31a611 100644 --- a/.github/workflows/Documenter.yml +++ b/.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: diff --git a/Project.toml b/Project.toml index f47a5b2c6..20b40306c 100644 --- a/Project.toml +++ b/Project.toml @@ -13,8 +13,8 @@ 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" @@ -22,7 +22,6 @@ 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,7 +34,6 @@ 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" @@ -43,8 +41,8 @@ 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" diff --git a/docs/src/general/neighborhood_search.md b/docs/src/general/neighborhood_search.md index 5a1512b6d..e611f5350 100644 --- a/docs/src/general/neighborhood_search.md +++ b/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 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘ + ``` diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 068eae953..96b6cdebf 100644 --- a/examples/n_body/n_body_system.jl +++ b/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, diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index a7ca8e274..c95451a25 100644 --- a/src/TrixiParticles.jl +++ b/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 diff --git a/src/general/general.jl b/src/general/general.jl index 6651c03c5..4c9531683 100644 --- a/src/general/general.jl +++ b/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") diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 9ede75648..8e94b1615 100644 --- a/src/general/initial_condition.jl +++ b/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 diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index fb917f91a..80116a605 100644 --- a/src/general/interpolation.jl +++ b/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 diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl new file mode 100644 index 000000000..0713f153d --- /dev/null +++ b/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 diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 36eb54766..e699f4bf0 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -13,8 +13,8 @@ the keyword argument `neighborhood_search`. A value of `nothing` means no neighb # Keywords - `neighborhood_search`: The type of neighborhood search to be used in the simulation. - By default, the [`GridNeighborhoodSearch`](@ref) is used. - Use [`TrivialNeighborhoodSearch`](@ref) or `nothing` to loop + By default, the `GridNeighborhoodSearch` is used. + Use `TrivialNeighborhoodSearch` or `nothing` to loop over all particles (no neighborhood search). - `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the coordinates of the domain corner in negative coordinate @@ -140,7 +140,8 @@ function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSear periodic_box_max_corner=periodic_box_max_corner, threaded_nhs_update=threaded_nhs_update) # Initialize neighborhood search - initialize!(search, initial_coordinates(neighbor)) + PointNeighbors.initialize!(search, initial_coordinates(system), + initial_coordinates(neighbor)) return search end @@ -437,7 +438,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) end # Update NHS - @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi) + @trixi_timeit timer() "update nhs" update_nhs!(semi, u_ode) # Second update step. # This is used to calculate density and pressure of the fluid systems @@ -467,14 +468,16 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) end end -function update_nhs(u_ode, semi) +function update_nhs!(semi, u_ode) # Update NHS for each pair of systems foreach_system(semi) do system + u_system = wrap_u(u_ode, system, semi) + foreach_system(semi) do neighbor u_neighbor = wrap_u(u_ode, neighbor, semi) neighborhood_search = get_neighborhood_search(system, neighbor, semi) - update!(neighborhood_search, nhs_coords(system, neighbor, u_neighbor)) + update_nhs!(neighborhood_search, system, neighbor, u_system, u_neighbor) end end end @@ -604,82 +607,120 @@ end end # NHS updates -# To prevent hard to spot errors there is not default version - -function nhs_coords(system::DEMSystem, neighbor::DEMSystem, u) - return current_coordinates(u, neighbor) -end - -function nhs_coords(system::BoundaryDEMSystem, - neighbor::Union{BoundaryDEMSystem, DEMSystem}, u) - return nothing -end -function nhs_coords(system::DEMSystem, neighbor::BoundaryDEMSystem, u) - return nothing -end - -function nhs_coords(system::FluidSystem, - neighbor::FluidSystem, u) - return current_coordinates(u, neighbor) -end - -function nhs_coords(system::FluidSystem, - neighbor::TotalLagrangianSPHSystem, u) - return current_coordinates(u, neighbor) -end - -function nhs_coords(system::FluidSystem, - neighbor::BoundarySPHSystem, u) - if neighbor.ismoving[] - return current_coordinates(u, neighbor) - end - - # Don't update - return nothing -end - -function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::FluidSystem, u) - return current_coordinates(u, neighbor) -end - -function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::TotalLagrangianSPHSystem, u) - # Don't update - return nothing -end - -function nhs_coords(system::TotalLagrangianSPHSystem, - neighbor::BoundarySPHSystem, u) - if neighbor.ismoving[] - return current_coordinates(u, neighbor) - end - - # Don't update - return nothing -end - -function nhs_coords(system::BoundarySPHSystem, - neighbor::FluidSystem, u) - # Don't update - return nothing -end - -function nhs_coords(system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, - neighbor::FluidSystem, u) - return current_coordinates(u, neighbor) -end - -function nhs_coords(system::BoundarySPHSystem, - neighbor::TotalLagrangianSPHSystem, u) - # Don't update - return nothing -end - -function nhs_coords(system::BoundarySPHSystem, - neighbor::BoundarySPHSystem, u) - # Don't update - return nothing +# To prevent hard-to-find bugs, there is not default version +function update_nhs!(neighborhood_search, + system::FluidSystem, + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, + u_system, u_neighbor) + # The current coordinates of fluids and solids change over time + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + +function update_nhs!(neighborhood_search, + system::FluidSystem, neighbor::BoundarySPHSystem, + u_system, u_neighbor) + # Boundary coordinates only change over time when `neighbor.ismoving[]` + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, neighbor.ismoving[])) +end + +function update_nhs!(neighborhood_search, + system::TotalLagrangianSPHSystem, neighbor::FluidSystem, + u_system, u_neighbor) + # The current coordinates of fluids and solids change over time + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + +function update_nhs!(neighborhood_search, + system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem, + u_system, u_neighbor) + # Don't update. Neighborhood search works on the initial coordinates, which don't change. + return neighborhood_search +end + +function update_nhs!(neighborhood_search, + system::TotalLagrangianSPHSystem, neighbor::BoundarySPHSystem, + u_system, u_neighbor) + # The current coordinates of solids change over time. + # Boundary coordinates only change over time when `neighbor.ismoving[]`. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, neighbor.ismoving[])) +end + +function update_nhs!(neighborhood_search, + system::BoundarySPHSystem, + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem, + BoundarySPHSystem}, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + +function update_nhs!(neighborhood_search, + system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, + neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, + u_system, u_neighbor) + # Depending on the density calculator of the boundary model, this NHS is used for + # - kernel summation (`SummationDensity`) + # - continuity equation (`ContinuityDensity`) + # - pressure extrapolation (`AdamiPressureExtrapolation`) + # + # Boundary coordinates only change over time when `neighbor.ismoving[]`. + # The current coordinates of fluids and solids change over time. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(system.ismoving[], true)) +end + +function update_nhs!(neighborhood_search, + system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, + neighbor::BoundarySPHSystem, + u_system, u_neighbor) + # `system` coordinates only change over time when `system.ismoving[]`. + # `neighbor` coordinates only change over time when `neighbor.ismoving[]`. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(system.ismoving[], neighbor.ismoving[])) +end + +function update_nhs!(neighborhood_search, + system::DEMSystem, neighbor::DEMSystem, + u_system, u_neighbor) + # Both coordinates change over time + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + +function update_nhs!(neighborhood_search, + system::DEMSystem, neighbor::BoundaryDEMSystem, + u_system, u_neighbor) + # DEM coordinates change over time, the boundary coordinates don't + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, false)) +end + +function update_nhs!(neighborhood_search, + system::BoundaryDEMSystem, + neighbor::Union{DEMSystem, BoundaryDEMSystem}, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search end function check_configuration(systems) diff --git a/src/neighborhood_search/grid_nhs.jl b/src/neighborhood_search/grid_nhs.jl deleted file mode 100644 index 82526a6b4..000000000 --- a/src/neighborhood_search/grid_nhs.jl +++ /dev/null @@ -1,423 +0,0 @@ -@doc raw""" - GridNeighborhoodSearch{NDIMS}(search_radius, n_particles; periodic_box_min_corner=nothing, - periodic_box_max_corner=nothing, threaded_nhs_update=true) - -Simple grid-based neighborhood search with uniform search radius. -The domain is divided into a regular grid. -For each (non-empty) grid cell, a list of particles in this cell is stored. -Instead of representing a finite domain by an array of cells, a potentially infinite domain -is represented by storing cell lists in a hash table (using Julia's `Dict` data structure), -indexed by the cell index tuple -```math -\left( \left\lfloor \frac{x}{d} \right\rfloor, \left\lfloor \frac{y}{d} \right\rfloor \right) \quad \text{or} \quad -\left( \left\lfloor \frac{x}{d} \right\rfloor, \left\lfloor \frac{y}{d} \right\rfloor, \left\lfloor \frac{z}{d} \right\rfloor \right), -``` -where ``x, y, z`` are the space coordinates and ``d`` is the search radius. - -To find particles within the search radius around a point, only particles in the neighboring -cells are considered. - -See also (Chalela et al., 2021), (Ihmsen et al. 2011, Section 4.4). - -As opposed to (Ihmsen et al. 2011), we do not sort the particles in any way, -since not sorting makes our implementation a lot faster (although less parallelizable). - -# Arguments -- `NDIMS`: Number of dimensions. -- `search_radius`: The uniform search radius. -- `n_particles`: Total number of particles. - -# 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 particle ordering changes. - -!!! warning "Internal use only" - Please note that this constructor is intended for internal use only. It is *not* part of - the public API of TrixiParticles.jl, and it thus can altered (or be removed) at any time - without it being considered a breaking change. - - To run a simulation with this neighborhood search, just pass the type to the constructor - of [`Semidiscretization`](@ref): - ```jldoctest semi_example; output=false, setup = :(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` explained - above can also be passed to the [`Semidiscretization`](@ref) and will internally be - forwarded to the neighborhood search: - ```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 │ - └──────────────────────────────────────────────────────────────────────────────────────────────────┘ - ``` - -## References -- M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán. - "GriSPy: A Python package for fixed-radius nearest neighbors search". - In: Astronomy and Computing 34 (2021). - [doi: 10.1016/j.ascom.2020.100443](https://doi.org/10.1016/j.ascom.2020.100443) -- Markus Ihmsen, Nadir Akinci, Markus Becker, Matthias Teschner. - "A Parallel SPH Implementation on Multi-Core CPUs". - In: Computer Graphics Forum 30.1 (2011), pages 99–112. - [doi: 10.1111/J.1467-8659.2010.01832.X](https://doi.org/10.1111/J.1467-8659.2010.01832.X) -""" -struct GridNeighborhoodSearch{NDIMS, ELTYPE, PB} - hashtable :: Dict{NTuple{NDIMS, Int}, Vector{Int}} - search_radius :: ELTYPE - empty_vector :: Vector{Int} # Just an empty vector (used in `eachneighbor`) - cell_buffer :: Array{NTuple{NDIMS, Int}, 2} # Multithreaded buffer for `update!` - cell_buffer_indices :: Vector{Int} # Store which entries of `cell_buffer` are initialized - periodic_box :: PB - n_cells :: NTuple{NDIMS, Int} - cell_size :: NTuple{NDIMS, ELTYPE} - threaded_nhs_update :: Bool - - function GridNeighborhoodSearch{NDIMS}(search_radius, n_particles; - periodic_box_min_corner=nothing, - periodic_box_max_corner=nothing, - threaded_nhs_update=true) where {NDIMS} - ELTYPE = typeof(search_radius) - - hashtable = Dict{NTuple{NDIMS, Int}, Vector{Int}}() - empty_vector = Int[] - cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_particles, Threads.nthreads()) - cell_buffer_indices = zeros(Int, Threads.nthreads()) - - if search_radius < eps() || - (periodic_box_min_corner === nothing && periodic_box_max_corner === nothing) - - # 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) - - # 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. - # Allow small tolerance to avoid inefficient larger cells due to machine - # rounding errors. - n_cells = Tuple(floor.(Int, (periodic_box.size .+ 10eps()) / search_radius)) - cell_size = Tuple(periodic_box.size ./ n_cells) - - if any(i -> i < 3, n_cells) - throw(ArgumentError("the `GridNeighborhoodSearch` needs at least 3 cells " * - "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(periodic_box)}(hashtable, search_radius, empty_vector, - cell_buffer, cell_buffer_indices, - periodic_box, n_cells, cell_size, threaded_nhs_update) - end -end - -@inline Base.ndims(neighborhood_search::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS - -@inline function nparticles(neighborhood_search::GridNeighborhoodSearch) - return size(neighborhood_search.cell_buffer, 1) -end - -function initialize!(neighborhood_search::GridNeighborhoodSearch, ::Nothing) - # No particle coordinates function -> don't initialize. - return neighborhood_search -end - -function initialize!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, - x::AbstractArray) where {NDIMS} - initialize!(neighborhood_search, i -> extract_svector(x, Val(NDIMS), i)) -end - -function initialize!(neighborhood_search::GridNeighborhoodSearch, coords_fun) - (; hashtable) = neighborhood_search - - empty!(hashtable) - - # This is needed to prevent lagging on macOS ARM. - # See https://github.com/JuliaSIMD/Polyester.jl/issues/89 - ThreadingUtilities.sleep_all_tasks() - - for particle in 1:nparticles(neighborhood_search) - # Get cell index of the particle's cell - cell = cell_coords(coords_fun(particle), neighborhood_search) - - # Add particle to corresponding cell or create cell if it does not exist - if haskey(hashtable, cell) - append!(hashtable[cell], particle) - else - hashtable[cell] = [particle] - end - end - - return neighborhood_search -end - -function update!(neighborhood_search::GridNeighborhoodSearch, ::Nothing) - # No particle coordinates function -> don't update. - return neighborhood_search -end - -function update!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, - x::AbstractArray) where {NDIMS} - update!(neighborhood_search, i -> extract_svector(x, Val(NDIMS), i)) -end - -# Modify the existing hash table by moving particles into their new cells -function update!(neighborhood_search::GridNeighborhoodSearch, coords_fun) - (; hashtable, cell_buffer, cell_buffer_indices, threaded_nhs_update) = neighborhood_search - - # Reset `cell_buffer` by moving all pointers to the beginning. - cell_buffer_indices .= 0 - - # Find all cells containing particles that now belong to another cell. - # `collect` the keyset to be able to loop over it with `@threaded`. - mark_changed_cell!(neighborhood_search, hashtable, coords_fun, Val(threaded_nhs_update)) - - # This is needed to prevent lagging on macOS ARM. - # See https://github.com/JuliaSIMD/Polyester.jl/issues/89 - ThreadingUtilities.sleep_all_tasks() - - # Iterate over all marked cells and move the particles into their new cells. - for thread in 1:Threads.nthreads() - # Only the entries `1:cell_buffer_indices[thread]` are initialized for `thread`. - for i in 1:cell_buffer_indices[thread] - cell = cell_buffer[i, thread] - particles = hashtable[cell] - - # Find all particles whose coordinates do not match this cell - moved_particle_indices = (i for i in eachindex(particles) - if cell_coords(coords_fun(particles[i]), - neighborhood_search) != cell) - - # Add moved particles to new cell - for i in moved_particle_indices - particle = particles[i] - new_cell_coords = cell_coords(coords_fun(particle), neighborhood_search) - - # Add particle to corresponding cell or create cell if it does not exist - if haskey(hashtable, new_cell_coords) - append!(hashtable[new_cell_coords], particle) - else - hashtable[new_cell_coords] = [particle] - end - end - - # Remove moved particles from this cell or delete the cell if it is now empty - if count(_ -> true, moved_particle_indices) == length(particles) - delete!(hashtable, cell) - else - deleteat!(particles, moved_particle_indices) - end - end - end - - return neighborhood_search -end - -@inline function mark_changed_cell!(neighborhood_search, hashtable, coords_fun, - threaded_nhs_update::Val{true}) - @threaded for cell in collect(keys(hashtable)) - mark_changed_cell!(neighborhood_search, cell, coords_fun) - end -end - -@inline function mark_changed_cell!(neighborhood_search, hashtable, coords_fun, - threaded_nhs_update::Val{false}) - for cell in collect(keys(hashtable)) - mark_changed_cell!(neighborhood_search, cell, coords_fun) - end -end - -# Use this function barrier and unpack inside to avoid passing closures to Polyester.jl -# with `@batch` (`@threaded`). -# Otherwise, `@threaded` does not work here with Julia ARM on macOS. -# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. -@inline function mark_changed_cell!(neighborhood_search, cell, coords_fun) - (; hashtable, cell_buffer, cell_buffer_indices) = neighborhood_search - - for particle in hashtable[cell] - if cell_coords(coords_fun(particle), neighborhood_search) != cell - # Mark this cell and continue with the next one. - # - # `cell_buffer` is preallocated, - # but only the entries 1:i are used for this thread. - i = cell_buffer_indices[Threads.threadid()] += 1 - cell_buffer[i, Threads.threadid()] = cell - break - end - end -end - -# 1D -@inline function eachneighbor(coords, neighborhood_search::GridNeighborhoodSearch{1}) - cell = cell_coords(coords, neighborhood_search) - x = cell[1] - # Generator of all neighboring cells to consider - neighboring_cells = ((x + i) for i in -1:1) - - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) - for cell in neighboring_cells) -end - -# 2D -@inline function eachneighbor(coords, neighborhood_search::GridNeighborhoodSearch{2}) - cell = cell_coords(coords, neighborhood_search) - x, y = cell - # Generator of all neighboring cells to consider - neighboring_cells = ((x + i, y + j) for i in -1:1, j in -1:1) - - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) - for cell in neighboring_cells) -end - -# 3D -@inline function eachneighbor(coords, neighborhood_search::GridNeighborhoodSearch{3}) - cell = cell_coords(coords, neighborhood_search) - x, y, z = cell - # Generator of all neighboring cells to consider - neighboring_cells = ((x + i, y + j, z + k) for i in -1:1, j in -1:1, k in -1:1) - - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) - for cell in neighboring_cells) -end - -@inline function particles_in_cell(cell_index, neighborhood_search) - (; hashtable, empty_vector) = neighborhood_search - - # Return an empty vector when `cell_index` is not a key of `hashtable` and - # reuse the empty vector to avoid allocations - return get(hashtable, periodic_cell_index(cell_index, neighborhood_search), - empty_vector) -end - -@inline function periodic_cell_index(cell_index, neighborhood_search) - (; n_cells, periodic_box) = neighborhood_search - - periodic_cell_index(cell_index, periodic_box, n_cells) -end - -@inline periodic_cell_index(cell_index, ::Nothing, n_cells) = cell_index - -@inline function periodic_cell_index(cell_index, periodic_box, n_cells) - return rem.(cell_index, n_cells, RoundDown) -end - -@inline function cell_coords(coords, neighborhood_search) - (; periodic_box, cell_size) = neighborhood_search - - return cell_coords(coords, periodic_box, cell_size) -end - -@inline function cell_coords(coords, periodic_box::Nothing, cell_size) - return Tuple(floor_to_int.(coords ./ cell_size)) -end - -@inline function cell_coords(coords, periodic_box, cell_size) - # Subtract `min_corner` to offset coordinates so that the min corner of the periodic - # box corresponds to the (0, 0) cell of the NHS. - # This way, there are no partial cells in the domain if the domain size is an integer - # multiple of the cell size (which is required, see the constructor). - offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner - - return Tuple(floor_to_int.(offset_coords ./ cell_size)) -end - -# When particles end up with coordinates so big that the cell coordinates -# exceed the range of Int, then `floor(Int, i)` will fail with an InexactError. -# In this case, we can just use typemax(Int), since we can assume that particles -# that far away will not interact with anything, anyway. -# This usually indicates an instability, but we don't want the simulation to crash, -# since adaptive time integration methods may detect the instability and reject the -# time step. -# If we threw an error here, we would prevent the time integration method from -# retrying with a smaller time step, and we would thus crash perfectly fine simulations. -@inline function floor_to_int(i) - if isnan(i) || i > typemax(Int) - return typemax(Int) - elseif i < typemin(Int) - return typemin(Int) - end - - return floor(Int, i) -end - -# Sorting only really makes sense in longer simulations where particles -# end up very unordered. -# WARNING: This is currently unmaintained. -function z_index_sort!(coordinates, system) - (; mass, pressure, neighborhood_search) = system - - perm = sortperm(eachparticle(system), - by=(i -> cell_z_index(extract_svector(coordinates, system, i), - neighborhood_search))) - - permute!(mass, perm) - permute!(pressure, perm) - Base.permutecols!!(u, perm) - - return nothing -end - -@inline function cell_z_index(coords, neighborhood_search) - cell = cell_coords(coords, neighborhood_search.search_radius) .+ 1 - - return cartesian2morton(SVector(cell)) -end - -# Create a copy of a neighborhood search but with a different search radius -function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, u) - if nhs.periodic_box === nothing - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs)) - else - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs), - periodic_box_min_corner=nhs.periodic_box.min_corner, - periodic_box_max_corner=nhs.periodic_box.max_corner) - end - - # Initialize neighborhood search - initialize!(search, u) - - return search -end - -# Create a copy of a neighborhood search but with a different search radius -function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, u) - return nhs -end diff --git a/src/neighborhood_search/neighborhood_search.jl b/src/neighborhood_search/neighborhood_search.jl deleted file mode 100644 index e9062c5db..000000000 --- a/src/neighborhood_search/neighborhood_search.jl +++ /dev/null @@ -1,110 +0,0 @@ -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) - end -end - -# 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)`. -@inline function 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 - -@inline function for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search; - particles=axes(system_coords, 2), parallel=true) - for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particles, - Val(parallel)) -end - -@inline function for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, particles, parallel::Val{true}) - @threaded for particle in particles - for_particle_neighbor_inner(f, system_coords, neighbor_coords, neighborhood_search, - particle) - end - - return nothing -end - -@inline function for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, particles, parallel::Val{false}) - for particle in particles - for_particle_neighbor_inner(f, system_coords, neighbor_coords, neighborhood_search, - particle) - end - - return nothing -end - -# Use this function barrier and unpack inside to avoid passing closures to Polyester.jl -# with `@batch` (`@threaded`). -# Otherwise, `@threaded` does not work here with Julia ARM on macOS. -# See https://github.com/JuliaSIMD/Polyester.jl/issues/88. -@inline function for_particle_neighbor_inner(f, system_coords, neighbor_system_coords, - neighborhood_search, particle) - (; 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) - neighbor_coords = extract_svector(neighbor_system_coords, - Val(ndims(neighborhood_search)), neighbor) - - pos_diff = particle_coords - neighbor_coords - distance2 = dot(pos_diff, pos_diff) - - pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, - periodic_box) - - if distance2 <= search_radius^2 - distance = sqrt(distance2) - - # Inline to avoid loss of performance - # compared to not using `for_particle_neighbor`. - @inline f(particle, neighbor, pos_diff, distance) - end - end -end - -@inline function compute_periodic_distance(pos_diff, distance2, search_radius, - periodic_box::Nothing) - return pos_diff, distance2 -end - -@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) - distance2 = dot(pos_diff, pos_diff) - end - - return pos_diff, distance2 -end - -@inline function periodic_coords(coords, periodic_box) - (; min_corner, size) = periodic_box - - # Move coordinates into the periodic box - box_offset = floor.((coords .- min_corner) ./ size) - - return coords - box_offset .* size -end - -@inline function periodic_coords(coords, periodic_box::Nothing) - return coords -end - -include("trivial_nhs.jl") -include("grid_nhs.jl") diff --git a/src/neighborhood_search/trivial_nhs.jl b/src/neighborhood_search/trivial_nhs.jl deleted file mode 100644 index 3902f1d3f..000000000 --- a/src/neighborhood_search/trivial_nhs.jl +++ /dev/null @@ -1,97 +0,0 @@ -@doc raw""" - TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle) - -Trivial neighborhood search that simply loops over all particles. -The search radius still needs to be passed in order to sort out particles outside the -search radius in the internal function `for_particle_neighbor`, but it's not used in the -internal function `eachneighbor`. - -# Arguments -- `NDIMS`: Number of dimensions. -- `search_radius`: The uniform search radius. -- `eachparticle`: `UnitRange` of all particle indices. Usually just `1:n_particles`. - -# 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. - -!!! warning "Internal use only" - Please note that this constructor is intended for internal use only. It is *not* part of - the public API of TrixiParticles.jl, and it thus can altered (or be removed) at any time - without it being considered a breaking change. - - To run a simulation with this neighborhood search, just pass the type to the constructor - of [`Semidiscretization`](@ref): - ```jldoctest semi_example; output=false, setup = :(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=TrivialNeighborhoodSearch) - - # output - ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ - │ Semidiscretization │ - │ ══════════════════ │ - │ #spatial dimensions: ………………………… 2 │ - │ #systems: ……………………………………………………… 2 │ - │ neighborhood search: ………………………… TrivialNeighborhoodSearch │ - │ total #particles: ………………………………… 636 │ - └──────────────────────────────────────────────────────────────────────────────────────────────────┘ - ``` - The keyword arguments `periodic_box_min_corner` and `periodic_box_max_corner` explained - above can also be passed to the [`Semidiscretization`](@ref) and will internally be - forwarded to the neighborhood search: - ```jldoctest semi_example; output = false - semi = Semidiscretization(system1, system2, - neighborhood_search=TrivialNeighborhoodSearch, - 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: ………………………… TrivialNeighborhoodSearch │ - │ total #particles: ………………………………… 636 │ - └──────────────────────────────────────────────────────────────────────────────────────────────────┘ - ``` -""" -struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} - search_radius :: ELTYPE - eachparticle :: EP - periodic_box :: PB - - function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle; - 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 - - new{NDIMS, typeof(search_radius), - typeof(eachparticle), typeof(periodic_box)}(search_radius, eachparticle, - periodic_box) - end -end - -@inline function Base.ndims(neighborhood_search::TrivialNeighborhoodSearch{NDIMS}) where { - NDIMS - } - return NDIMS -end - -@inline initialize!(search::TrivialNeighborhoodSearch, coords_fun) = search -@inline update!(search::TrivialNeighborhoodSearch, coords_fun) = search -@inline eachneighbor(coords, search::TrivialNeighborhoodSearch) = search.eachparticle diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 81581494e..bd1658e32 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -85,7 +85,8 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix # Reset the collection when the iteration is 0 pvd = paraview_collection(collection_file; append=iter > 0) - points = periodic_coords(current_coordinates(u, system), periodic_box) + points = PointNeighbors.periodic_coords(current_coordinates(u, system), + periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] if abs(maximum(points)) > max_coordinates || abs(minimum(points)) > max_coordinates diff --git a/test/count_allocations.jl b/test/count_allocations.jl index 6a87e9868..3b5110ba1 100644 --- a/test/count_allocations.jl +++ b/test/count_allocations.jl @@ -27,7 +27,11 @@ end end # No update -@inline TrixiParticles.update!(search::NoUpdateNeighborhoodSearch, coords_fun) = search +@inline function TrixiParticles.PointNeighbors.update!(search::NoUpdateNeighborhoodSearch, + x, y; + particles_moving=(true, true)) + return search +end # Count allocations of one call to the right-hand side (`kick!` + `drift!`) function count_rhs_allocations(sol, semi) diff --git a/test/neighborhood_search/grid_nhs.jl b/test/neighborhood_search/grid_nhs.jl deleted file mode 100644 index d686efc5b..000000000 --- a/test/neighborhood_search/grid_nhs.jl +++ /dev/null @@ -1,272 +0,0 @@ -@testset verbose=true "GridNeighborhoodSearch" begin - @testset "Coordinate Limits" begin - # Test the threshold for very large and very small coordinates. - coords1 = [Inf, -Inf] - coords2 = [NaN, 0] - coords3 = [typemax(Int) + 1.0, -typemax(Int) - 1.0] - - @test TrixiParticles.cell_coords(coords1, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) - @test TrixiParticles.cell_coords(coords2, nothing, (1.0, 1.0)) == (typemax(Int), 0) - @test TrixiParticles.cell_coords(coords3, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) - end - - @testset "Rectangular Point Cloud 2D" begin - #### Setup - # Rectangular filled with equidistant spaced particles - # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35) - range = -0.25:0.1:0.35 - coordinates1 = hcat(collect.(Iterators.product(range, range))...) - nparticles = size(coordinates1, 2) - - particle_position1 = [0.05, 0.05] - particle_spacing = 0.1 - radius = particle_spacing - - # Create neighborhood search - nhs1 = GridNeighborhoodSearch{2}(radius, nparticles) - - coords_fun(i) = coordinates1[:, i] - TrixiParticles.initialize!(nhs1, coords_fun) - - # Get each neighbor for `particle_position1` - neighbors1 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs1))) - - # Move particles - coordinates2 = coordinates1 .+ [1.4, -3.5] - - # Update neighborhood search - coords_fun2(i) = coordinates2[:, i] - TrixiParticles.update!(nhs1, coords_fun2) - - # Get each neighbor for updated NHS - neighbors2 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs1))) - - # Change position - particle_position2 = particle_position1 .+ [1.4, -3.5] - - # Get each neighbor for `particle_position2` - neighbors3 = sort(collect(TrixiParticles.eachneighbor(particle_position2, nhs1))) - - # Double search radius - nhs2 = GridNeighborhoodSearch{2}(2 * radius, size(coordinates1, 2)) - TrixiParticles.initialize!(nhs2, coords_fun) - - # Get each neighbor in double search radius - neighbors4 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs2))) - - # Move particles - coordinates2 = coordinates1 .+ [0.4, -0.4] - - # Update neighborhood search - TrixiParticles.update!(nhs2, coords_fun2) - - # Get each neighbor in double search radius - neighbors5 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs2))) - - #### Verification - @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33] - - @test neighbors2 == Int[] - - @test neighbors3 == [17, 18, 19, 24, 25, 26, 31, 32, 33] - - @test neighbors4 == [9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, - 26, 27, 28, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, - 48, 49] - - @test neighbors5 == [36, 37, 38, 43, 44, 45] - end - - @testset "Rectangular Point Cloud 3D" begin - #### Setup - # Rectangular filled with equidistant spaced particles - # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y) = (0.35, 0.35, 0.35) - range = -0.25:0.1:0.35 - coordinates1 = hcat(collect.(Iterators.product(range, range, range))...) - nparticles = size(coordinates1, 2) - - particle_position1 = [0.05, 0.05, 0.05] - particle_spacing = 0.1 - radius = particle_spacing - - # Create neighborhood search - nhs1 = GridNeighborhoodSearch{3}(radius, nparticles) - - coords_fun(i) = coordinates1[:, i] - TrixiParticles.initialize!(nhs1, coords_fun) - - # Get each neighbor for `particle_position1` - neighbors1 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs1))) - - # Move particles - coordinates2 = coordinates1 .+ [1.4, -3.5, 0.8] - - # Update neighborhood search - coords_fun2(i) = coordinates2[:, i] - TrixiParticles.update!(nhs1, coords_fun2) - - # Get each neighbor for updated NHS - neighbors2 = sort(collect(TrixiParticles.eachneighbor(particle_position1, nhs1))) - - # Change position - particle_position2 = particle_position1 .+ [1.4, -3.5, 0.8] - - # Get each neighbor for `particle_position2` - neighbors3 = sort(collect(TrixiParticles.eachneighbor(particle_position2, nhs1))) - - #### Verification - @test neighbors1 == - [115, 116, 117, 122, 123, 124, 129, 130, 131, 164, 165, 166, 171, 172, - 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] - - @test neighbors2 == Int[] - - @test neighbors3 == - [115, 116, 117, 122, 123, 124, 129, 130, 131, 164, 165, 166, 171, 172, - 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] - end - - @testset verbose=true "Periodicity 2D" begin - @testset "Simple Example" begin - coords = [-0.08 0.0 0.18 0.1 -0.08 - -0.12 -0.05 -0.09 0.15 0.39] - - # 3 x 6 cells - nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2), - periodic_box_min_corner=[-0.1, -0.2], - periodic_box_max_corner=[0.2, 0.4]) - - TrixiParticles.initialize!(nhs, coords) - - neighbors = [sort(collect(TrixiParticles.eachneighbor(coords[:, i], nhs))) - for i in 1:5] - - # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells - @test neighbors[1] == [1, 2, 3, 5] - @test neighbors[2] == [1, 2, 3] - @test neighbors[3] == [1, 2, 3] - @test neighbors[4] == [4] - @test neighbors[5] == [1, 5] - - neighbors_loop = [Int[] for _ in axes(coords, 2)] - - TrixiParticles.for_particle_neighbor(nothing, nothing, - coords, coords, nhs, - particles=axes(coords, 2)) do particle, - neighbor, - pos_diff, - distance - append!(neighbors_loop[particle], neighbor) - end - - @test sort(neighbors_loop[1]) == [1, 3, 5] - @test sort(neighbors_loop[2]) == [2] - @test sort(neighbors_loop[3]) == [1, 3] - @test sort(neighbors_loop[4]) == [4] - @test sort(neighbors_loop[5]) == [1, 5] - end - - @testset "Rounding Up Cell Sizes" begin - coords = [-0.08 0.0 0.18 0.1 -0.08 - -0.12 -0.05 -0.09 0.15 0.42] - - # 3 x 6 cells - nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2), - periodic_box_min_corner=[-0.1, -0.2], - periodic_box_max_corner=[0.205, 0.43]) - - TrixiParticles.initialize!(nhs, coords) - - neighbors = [sort(collect(TrixiParticles.eachneighbor(coords[:, i], nhs))) - for i in 1:5] - - # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells - @test neighbors[1] == [1, 2, 3, 5] - @test neighbors[2] == [1, 2, 3] - @test neighbors[3] == [1, 2, 3] - @test neighbors[4] == [4] - @test neighbors[5] == [1, 5] - - neighbors_loop = [Int[] for _ in axes(coords, 2)] - - TrixiParticles.for_particle_neighbor(nothing, nothing, - coords, coords, nhs, - particles=axes(coords, 2)) do particle, - neighbor, - pos_diff, - distance - append!(neighbors_loop[particle], neighbor) - end - - @test sort(neighbors_loop[1]) == [1, 3, 5] - @test sort(neighbors_loop[2]) == [2] - @test sort(neighbors_loop[3]) == [1, 3] - @test sort(neighbors_loop[4]) == [4] - @test sort(neighbors_loop[5]) == [1, 5] - end - - @testset "Offset Domain Triggering Split Cells" begin - # This used to trigger a "split cell bug", where the left and right boundary - # cells were only partially contained in the domain. - # The left particle was placed inside a ghost cells, which causes it to not - # see the right particle, even though it is within the search distance. - # The domain size is an integer multiple of the cell size, but the NHS did not - # offset the grid based on the domain position. - # See https://github.com/trixi-framework/TrixiParticles.jl/pull/211 for a more - # detailed explanation. - coords = [-1.4 1.9 - 0.0 0.0] - - # 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]) - - TrixiParticles.initialize!(nhs, coords) - - neighbors = [sort(unique(collect(TrixiParticles.eachneighbor(coords[:, i], nhs)))) - for i in 1:2] - - @test neighbors[1] == [1, 2] - @test neighbors[2] == [1, 2] - end - end - - @testset verbose=true "Periodicity 3D" begin - coords = [-0.08 0.0 0.18 0.1 -0.08 - -0.12 -0.05 -0.09 0.15 0.39 - 0.14 0.34 0.12 0.06 0.13] - - # 3 x 6 x 3 cells - nhs = GridNeighborhoodSearch{3}(0.1, size(coords, 2), - periodic_box_min_corner=[-0.1, -0.2, 0.05], - periodic_box_max_corner=[0.2, 0.4, 0.35]) - - TrixiParticles.initialize!(nhs, coords) - - neighbors = [sort(collect(TrixiParticles.eachneighbor(coords[:, i], nhs))) - for i in 1:5] - - # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells - @test neighbors[1] == [1, 2, 3, 5] - @test neighbors[2] == [1, 2, 3] - @test neighbors[3] == [1, 2, 3] - @test neighbors[4] == [4] - @test neighbors[5] == [1, 5] - - neighbors_loop = [Int[] for _ in axes(coords, 2)] - - TrixiParticles.for_particle_neighbor(coords, coords, - nhs) do particle, neighbor, pos_diff, distance - append!(neighbors_loop[particle], neighbor) - end - - @test sort(neighbors_loop[1]) == [1, 3, 5] - @test sort(neighbors_loop[2]) == [2] - @test sort(neighbors_loop[3]) == [1, 3] - @test sort(neighbors_loop[4]) == [4] - @test sort(neighbors_loop[5]) == [1, 5] - end -end diff --git a/test/neighborhood_search/neighborhood_search.jl b/test/neighborhood_search/neighborhood_search.jl deleted file mode 100644 index 0c7270a7e..000000000 --- a/test/neighborhood_search/neighborhood_search.jl +++ /dev/null @@ -1,2 +0,0 @@ -include("trivial_nhs.jl") -include("grid_nhs.jl") diff --git a/test/neighborhood_search/trivial_nhs.jl b/test/neighborhood_search/trivial_nhs.jl deleted file mode 100644 index aa0b21846..000000000 --- a/test/neighborhood_search/trivial_nhs.jl +++ /dev/null @@ -1,10 +0,0 @@ -@testset verbose=true "TrivialNeighborhoodSearch" begin - # Setup with 5 particles - nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(1.0, Base.OneTo(5)) - - # Get each neighbor for arbitrary coordinates - neighbors = collect(TrixiParticles.eachneighbor([1.0, 2.0], nhs)) - - #### Verification - @test neighbors == [1, 2, 3, 4, 5] -end diff --git a/test/schemes/solid/total_lagrangian_sph/rhs.jl b/test/schemes/solid/total_lagrangian_sph/rhs.jl index 6a4f6bce3..dfc9fc234 100644 --- a/test/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/test/schemes/solid/total_lagrangian_sph/rhs.jl @@ -85,7 +85,7 @@ each_moving_particle end TrixiParticles.eachparticle(::Val{:mock_system_interact}) = eachparticle - TrixiParticles.eachneighbor(_, ::Val{:nhs}) = eachneighbor + TrixiParticles.PointNeighbors.eachneighbor(_, ::Val{:nhs}) = eachneighbor function Base.getproperty(::Val{:nhs}, f::Symbol) if f === :search_radius diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index 812493fa8..369022273 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -153,7 +153,7 @@ return Val(Symbol("mock_" * string(f))) end - TrixiParticles.eachneighbor(_, ::Val{:nhs}) = neighbors + TrixiParticles.PointNeighbors.eachneighbor(_, ::Val{:nhs}) = neighbors function Base.getproperty(::Val{:nhs}, f::Symbol) if f === :search_radius diff --git a/test/unittest.jl b/test/unittest.jl index 01308bd44..5c05924e1 100644 --- a/test/unittest.jl +++ b/test/unittest.jl @@ -3,7 +3,6 @@ @testset verbose=true "Unit Tests" begin include("callbacks/callbacks.jl") include("general/general.jl") - include("neighborhood_search/neighborhood_search.jl") include("setups/setups.jl") include("systems/systems.jl") include("schemes/schemes.jl")