From 625463f2c384be9e1a98503a6d2d1ca519042f90 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 24 Feb 2026 16:22:04 +0100 Subject: [PATCH 1/6] Implement NHS handler --- src/general/neighborhood_search.jl | 60 +++++++++++++++++++++++++++--- src/general/semidiscretization.jl | 13 +++---- 2 files changed, 59 insertions(+), 14 deletions(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 962567a87c..2526d5a94f 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -104,12 +104,54 @@ function create_neighborhood_search(neighborhood_search, system::TotalLagrangian end # === Neighborhood search lookup === -@inline function get_neighborhood_search(system, semi) - (; neighborhood_searches) = semi +abstract type AbstractNHSHandler end - system_index = system_indices(system, semi) +struct PairsNHSHandler{NHS} <: AbstractNHSHandler + neighborhood_searches::NHS +end + +function create_neighborhood_search_handler(::Type{<:PairsNHSHandler}, + neighborhood_search, systems) + # Create a tuple of n neighborhood searches for each of the n systems. + # We will need one neighborhood search for each pair of systems. + return Tuple(Tuple(create_neighborhood_search(neighborhood_search, + system, neighbor) + for neighbor in systems) + for system in systems) +end + +function get_neighborhood_search(handler::PairsNHSHandler, system_index, neighbor_index, _) + return handler.neighborhood_searches[system_index][neighbor_index] +end - return neighborhood_searches[system_index][system_index] +struct GridNHSHandler{SI, NHS} <: AbstractNHSHandler + search_radii :: SI + neighborhood_searches :: NHS +end + +function create_neighborhood_search_handler(::Type{<:GridNHSHandler}, nhs, systems) + # Find a list of search radii that will be requested for each system + search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] for system in systems) + search_radii = Tuple.(unique.(search_radii)) + + # For each system, create a neighborhood search for each unique search radius + return Tuple(Tuple(copy_neighborhood_search(nhs, search_radius, nparticles(system)) + for search_radius in search_radii[system_indices(system, systems)]) + for system in systems) +end + +function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor_index, + search_radius) + # Find the index of the search radius in the list of search radii for this system + search_radii = handler.search_radii[neighbor_index] + idx = searchsortedfirst(search_radii, search_radius) + + # Return the neighborhood search at that index + return handler.neighborhood_searches[neighbor_index][idx] +end + +@inline function get_neighborhood_search(system, semi) + return get_neighborhood_search(system, system, semi) end @inline function get_neighborhood_search(system::TotalLagrangianSPHSystem, semi) @@ -124,7 +166,10 @@ end system_index = system_indices(system, semi) neighbor_index = system_indices(neighbor_system, semi) - return neighborhood_searches[system_index][neighbor_index] + search_radius = compact_support(system, neighbor_system) + + return get_neighborhood_search(neighborhood_searches, system_index, neighbor_index, + search_radius) end @inline function get_neighborhood_search(system::TotalLagrangianSPHSystem, @@ -140,7 +185,10 @@ end return system.self_interaction_nhs end - return neighborhood_searches[system_index][neighbor_index] + search_radius = compact_support(system, neighbor_system) + + return get_neighborhood_search(neighborhood_searches, system_index, neighbor_index, + search_radius) end # === Initialization === diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 5def744bbb..0c23cb96e4 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -98,12 +98,7 @@ function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) - # Create a tuple of n neighborhood searches for each of the n systems. - # We will need one neighborhood search for each pair of systems. - searches = Tuple(Tuple(create_neighborhood_search(neighborhood_search, - system, neighbor) - for neighbor in systems) - for system in systems) + searches = create_neighborhood_search_handler(PairsNHSHandler, neighborhood_search, systems) # These will be set to true inside the `UpdateCallback`. # Some techniques require the use of this callback, and this flag can be used @@ -152,10 +147,12 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) end end -@inline function system_indices(system, semi) +@inline system_indices(system, semi::Semidiscretization) = system_indices(system, semi.systems) + +@inline function system_indices(system, systems) # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict` # is ~30x slower because `hash(::System)` is very slow. - index = findfirst(==(system), semi.systems) + index = findfirst(==(system), systems) if isnothing(index) throw(ArgumentError("system is not in the semidiscretization")) From 671b64ab40813cf49e9b9441c2a24fce22aa4516 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 24 Feb 2026 16:27:33 +0100 Subject: [PATCH 2/6] Add search radius tolerance --- src/general/neighborhood_search.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 2526d5a94f..904ba26698 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -144,7 +144,7 @@ function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor search_radius) # Find the index of the search radius in the list of search radii for this system search_radii = handler.search_radii[neighbor_index] - idx = searchsortedfirst(search_radii, search_radius) + idx = searchsortedfirst(search_radii, search_radius - eps(search_radius)) # Return the neighborhood search at that index return handler.neighborhood_searches[neighbor_index][idx] From a4dcd480d6350866d757cfea24e82e38320476d8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 24 Feb 2026 16:40:42 +0100 Subject: [PATCH 3/6] Revise everything --- src/general/neighborhood_search.jl | 68 ++++++++++++++++++++++++------ 1 file changed, 56 insertions(+), 12 deletions(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 904ba26698..599e44ab11 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -124,30 +124,74 @@ function get_neighborhood_search(handler::PairsNHSHandler, system_index, neighbo return handler.neighborhood_searches[system_index][neighbor_index] end -struct GridNHSHandler{SI, NHS} <: AbstractNHSHandler - search_radii :: SI +struct VariableRadiusNHS{SR, NHS} <: PointNeighbors.AbstractNeighborhoodSearch + search_radii :: SR neighborhood_searches :: NHS end -function create_neighborhood_search_handler(::Type{<:GridNHSHandler}, nhs, systems) +function VariableSearchRadiusNHS(nhs_implementation, search_radii, n_particles) + searches = Tuple(copy_neighborhood_search(nhs_implementation, search_radius, n_particles) + for search_radius in search_radii) + return VariableRadiusNHS(search_radii, searches) +end + +@inline Base.ndims(::TrivialNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS + +@inline requires_update(::TrivialNeighborhoodSearch) = (false, true) + +@inline function initialize!(search::TrivialNeighborhoodSearch, x, y; + parallelization_backend = default_backend(x), + eachindex_y = axes(y, 2)) + return search +end + +@inline function update!(search::TrivialNeighborhoodSearch, x, y; + points_moving = (true, true), + parallelization_backend = default_backend(x), + eachindex_y = axes(y, 2)) + return search +end + +# Create a copy of a neighborhood search but with a different search radius +function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, x, y) + return nhs +end + +function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, + search_radius, n_points; eachpoint = 1:n_points) + return TrivialNeighborhoodSearch{ndims(nhs)}(; search_radius, eachpoint, + periodic_box = nhs.periodic_box) +end + +@inline function foreach_neighbor(f, neighbor_system_coords, + neighborhood_search::VariableRadiusNHS, + point, point_coords, search_radius) + # Find the index of the search radius in the list of search radii for this system + search_radii = neighborhood_search.search_radii + idx = searchsortedfirst(search_radii, search_radius - eps(search_radius)) + + nhs = neighborhood_search.neighborhood_searches[idx] + + PointNeighbors.foreach_neighbor(f, neighbor_system_coords, nhs, point, point_coords, search_radius) +end + +struct VariableNHSHandler{NHS} <: AbstractNHSHandler + neighborhood_searches::NHS +end + +function create_neighborhood_search_handler(::Type{<:VariableNHSHandler}, nhs, systems) # Find a list of search radii that will be requested for each system search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] for system in systems) search_radii = Tuple.(unique.(search_radii)) # For each system, create a neighborhood search for each unique search radius - return Tuple(Tuple(copy_neighborhood_search(nhs, search_radius, nparticles(system)) - for search_radius in search_radii[system_indices(system, systems)]) + return Tuple(VariableSearchRadiusNHS(nhs, search_radii[system_indices(system, systems)], nparticles(system)) for system in systems) end -function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor_index, +function get_neighborhood_search(handler::VariableNHSHandler, system_index, neighbor_index, search_radius) - # Find the index of the search radius in the list of search radii for this system - search_radii = handler.search_radii[neighbor_index] - idx = searchsortedfirst(search_radii, search_radius - eps(search_radius)) - - # Return the neighborhood search at that index - return handler.neighborhood_searches[neighbor_index][idx] + handler.neighborhood_searches[neighbor_index] end @inline function get_neighborhood_search(system, semi) From 314f8b6d28e1268d350e49db1ba13024d2e9a4ee Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Mon, 9 Mar 2026 18:44:14 +0100 Subject: [PATCH 4/6] VariableNHSHandler and GridNHSHandler are running correctly with the 2d dam break example. --- src/general/neighborhood_search.jl | 118 +++++++++++++++++++++-------- src/general/semidiscretization.jl | 36 +++++---- 2 files changed, 106 insertions(+), 48 deletions(-) diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 599e44ab11..e33c44f66c 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -114,85 +114,135 @@ function create_neighborhood_search_handler(::Type{<:PairsNHSHandler}, neighborhood_search, systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. - return Tuple(Tuple(create_neighborhood_search(neighborhood_search, - system, neighbor) + searches = Tuple(Tuple(create_neighborhood_search(neighborhood_search, + system, neighbor) for neighbor in systems) for system in systems) + + return PairsNHSHandler(searches) end function get_neighborhood_search(handler::PairsNHSHandler, system_index, neighbor_index, _) return handler.neighborhood_searches[system_index][neighbor_index] end -struct VariableRadiusNHS{SR, NHS} <: PointNeighbors.AbstractNeighborhoodSearch +struct VariableNHSHandler{NHS} <: AbstractNHSHandler + neighborhood_searches::NHS +end + +function create_neighborhood_search_handler(::Type{<:VariableNHSHandler}, nhs, systems) + # Find a list of search radii that will be requested for each system + search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] + for system in systems) + search_radii = Tuple.(unique.(search_radii)) + + # For each system, create a neighborhood search for each unique search radius + searches = Tuple(VariableSearchRadiusNHS(nhs, + search_radii[system_indices(system, + systems)], + nparticles(system)) + for system in systems) + + # For each system, create a neighborhood search for each unique search radius + return VariableNHSHandler(searches) +end + +function get_neighborhood_search(handler::VariableNHSHandler, system_index, neighbor_index, + search_radius) + handler.neighborhood_searches[neighbor_index] +end + +struct VariableSearchRadiusNHS{SR, NHS} <: PointNeighbors.AbstractNeighborhoodSearch search_radii :: SR neighborhood_searches :: NHS end function VariableSearchRadiusNHS(nhs_implementation, search_radii, n_particles) - searches = Tuple(copy_neighborhood_search(nhs_implementation, search_radius, n_particles) + searches = Tuple(PointNeighbors.copy_neighborhood_search(nhs_implementation, + search_radius, n_particles) for search_radius in search_radii) - return VariableRadiusNHS(search_radii, searches) + return VariableSearchRadiusNHS(search_radii, searches) end -@inline Base.ndims(::TrivialNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS +@inline Base.ndims(search::VariableSearchRadiusNHS) = ndims(first(search.neighborhood_searches)) -@inline requires_update(::TrivialNeighborhoodSearch) = (false, true) +@inline requires_update(::VariableSearchRadiusNHS) = (false, true) -@inline function initialize!(search::TrivialNeighborhoodSearch, x, y; - parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) +# TODO +@inline function initialize!(search::VariableSearchRadiusNHS, x, y; + parallelization_backend=default_backend(x), + eachindex_y=axes(y, 2)) return search end -@inline function update!(search::TrivialNeighborhoodSearch, x, y; - points_moving = (true, true), - parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) +# TODO +@inline function update!(search::VariableSearchRadiusNHS, x, y; + points_moving=(true, true), + parallelization_backend=default_backend(x), + eachindex_y=axes(y, 2)) return search end # Create a copy of a neighborhood search but with a different search radius -function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, x, y) - return nhs +function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, + search_radius, x, y) + search_radii = Tuple(PointNeighbors.search_radius(search) + for search in nhs.neighborhood_searches) + searches = Tuple(PointNeighbors.copy_neighborhood_search(search) + for search in nhs.neighborhood_searches) + + return VariableSearchRadiusNHS(search_radii, searches) end -function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, - search_radius, n_points; eachpoint = 1:n_points) +function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, + search_radius, n_points; + eachpoint=1:n_points) return TrivialNeighborhoodSearch{ndims(nhs)}(; search_radius, eachpoint, - periodic_box = nhs.periodic_box) + periodic_box=nhs.periodic_box) end @inline function foreach_neighbor(f, neighbor_system_coords, - neighborhood_search::VariableRadiusNHS, + neighborhood_search::VariableSearchRadiusNHS, point, point_coords, search_radius) # Find the index of the search radius in the list of search radii for this system search_radii = neighborhood_search.search_radii - idx = searchsortedfirst(search_radii, search_radius - eps(search_radius)) + idx = searchsortedfirst(SVector(search_radii), search_radius - eps(search_radius)) nhs = neighborhood_search.neighborhood_searches[idx] - PointNeighbors.foreach_neighbor(f, neighbor_system_coords, nhs, point, point_coords, search_radius) + PointNeighbors.foreach_neighbor(f, neighbor_system_coords, nhs, point, point_coords, + search_radius) end -struct VariableNHSHandler{NHS} <: AbstractNHSHandler +struct GridNHSHandler{SI, NHS} <: AbstractNHSHandler + search_radii::SI neighborhood_searches::NHS end -function create_neighborhood_search_handler(::Type{<:VariableNHSHandler}, nhs, systems) +function create_neighborhood_search_handler(::Type{<:GridNHSHandler}, nhs, systems) # Find a list of search radii that will be requested for each system - search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] for system in systems) + search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] + for system in systems) search_radii = Tuple.(unique.(search_radii)) # For each system, create a neighborhood search for each unique search radius - return Tuple(VariableSearchRadiusNHS(nhs, search_radii[system_indices(system, systems)], nparticles(system)) - for system in systems) + searches = Tuple(Tuple(copy_neighborhood_search(nhs, search_radius, nparticles(system)) + for search_radius in + search_radii[system_indices(system, systems)]) + for system in systems) + + return GridNHSHandler(search_radii, searches) end -function get_neighborhood_search(handler::VariableNHSHandler, system_index, neighbor_index, +function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor_index, search_radius) - handler.neighborhood_searches[neighbor_index] + # Find the index of the search radius in the list of search radii for this system + search_radii = handler.search_radii[neighbor_index] + idx = searchsortedfirst(SVector(search_radii), search_radius - eps(search_radius)) + # Return the neighborhood search at that index + return handler.neighborhood_searches[neighbor_index][idx] end +# ================================== @inline function get_neighborhood_search(system, semi) return get_neighborhood_search(system, system, semi) @@ -205,20 +255,21 @@ end end @inline function get_neighborhood_search(system, neighbor_system, semi) - (; neighborhood_searches) = semi + (; neighborhood_search_handler) = semi system_index = system_indices(system, semi) neighbor_index = system_indices(neighbor_system, semi) search_radius = compact_support(system, neighbor_system) - return get_neighborhood_search(neighborhood_searches, system_index, neighbor_index, + return get_neighborhood_search(neighborhood_search_handler, system_index, + neighbor_index, search_radius) end @inline function get_neighborhood_search(system::TotalLagrangianSPHSystem, neighbor_system::TotalLagrangianSPHSystem, semi) - (; neighborhood_searches) = semi + (; neighborhood_search_handler) = semi system_index = system_indices(system, semi) neighbor_index = system_indices(neighbor_system, semi) @@ -231,7 +282,8 @@ end search_radius = compact_support(system, neighbor_system) - return get_neighborhood_search(neighborhood_searches, system_index, neighbor_index, + return get_neighborhood_search(neighborhood_search_handler, system_index, + neighbor_index, search_radius) end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 0c23cb96e4..94ddf42af8 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -49,32 +49,34 @@ semi = Semidiscretization(fluid_system, boundary_system, └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -struct Semidiscretization{BACKEND, S, RU, RV, NS, UCU, IT} - systems :: S - ranges_u :: RU - ranges_v :: RV - neighborhood_searches :: NS - parallelization_backend :: BACKEND - update_callback_used :: UCU - integrate_tlsph :: IT # `false` if TLSPH integration is decoupled +struct Semidiscretization{BACKEND, S, RU, RV, NSH, UCU, IT} + systems :: S + ranges_u :: RU + ranges_v :: RV + neighborhood_search_handler :: NSH + parallelization_backend :: BACKEND + update_callback_used :: UCU + integrate_tlsph :: IT # `false` if TLSPH integration is decoupled # Dispatch at `systems` to distinguish this constructor from the one below when # 4 systems are passed. # This is an internal constructor only used in `test/count_allocations.jl` # and by Adapt.jl. - function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches, + function Semidiscretization(systems::Tuple, ranges_u, ranges_v, + neighborhood_search_handler, parallelization_backend::PointNeighbors.ParallelizationBackend, update_callback_used, integrate_tlsph) new{typeof(parallelization_backend), typeof(systems), typeof(ranges_u), - typeof(ranges_v), typeof(neighborhood_searches), + typeof(ranges_v), typeof(neighborhood_search_handler), typeof(update_callback_used), typeof(integrate_tlsph)}(systems, ranges_u, ranges_v, - neighborhood_searches, parallelization_backend, + neighborhood_search_handler, parallelization_backend, update_callback_used, integrate_tlsph) end end function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; + nhs_handler=GridNHSHandler, neighborhood_search=GridNeighborhoodSearch{ndims(first(systems))}(), parallelization_backend=PolyesterBackend()) systems = filter(system -> !isnothing(system), systems) @@ -98,7 +100,9 @@ function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) - searches = create_neighborhood_search_handler(PairsNHSHandler, neighborhood_search, systems) + neighborhood_search_handler = create_neighborhood_search_handler(nhs_handler, + neighborhood_search, + systems) # These will be set to true inside the `UpdateCallback`. # Some techniques require the use of this callback, and this flag can be used @@ -110,7 +114,8 @@ function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; # with this set to false. integrate_tlsph = Ref(true) - return Semidiscretization(systems, ranges_u, ranges_v, searches, + return Semidiscretization(systems, ranges_u, ranges_v, + neighborhood_search_handler, parallelization_backend, update_callback_used, integrate_tlsph) end @@ -139,7 +144,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) summary_line(io, "#spatial dimensions", ndims(semi.systems[1])) summary_line(io, "#systems", length(semi.systems)) summary_line(io, "neighborhood search", - semi.neighborhood_searches |> eltype |> eltype |> nameof) + semi.neighborhood_search_handler |> eltype |> eltype |> nameof) summary_line(io, "total #particles", sum(nparticles.(semi.systems))) summary_line(io, "eltype", eltype(semi.systems[1])) summary_line(io, "coordinates eltype", coordinates_eltype(semi.systems[1])) @@ -147,7 +152,8 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) end end -@inline system_indices(system, semi::Semidiscretization) = system_indices(system, semi.systems) +@inline system_indices(system, + semi::Semidiscretization) = system_indices(system, semi.systems) @inline function system_indices(system, systems) # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict` From 8511ea16a4b6c9a2776f3116868a37129fab7fda Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Wed, 25 Mar 2026 20:15:22 +0100 Subject: [PATCH 5/6] Implement `initialize!` and `update!` for VariabelSearchRadiusNHS. Add `adapt_structure` for VariableSearchRadiusNHS, PairsNHSHandler, VariableNHSHandler. Add documentation for nhs handlers and variable search radius nhs. Export VariableSearchRadiusNHS, PairsNHSHandler, VariableNHSHandler. Replace `create_neighborhood_search_handler()` with outer constructor. Sort the search_radii before calling `searchsortedfirst`. Add basic test for the different nhs handlers. --- src/TrixiParticles.jl | 2 + src/general/gpu.jl | 28 ++++++++ src/general/neighborhood_search.jl | 101 +++++++++++++++-------------- src/general/semidiscretization.jl | 6 +- test/general/nhs_handler.jl | 37 +++++++++++ 5 files changed, 121 insertions(+), 53 deletions(-) create mode 100644 test/general/nhs_handler.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 553a778aa9..c1351e989f 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -112,5 +112,7 @@ export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, SurfaceTensionMomentumMorris export ColorfieldSurfaceNormal export SymplecticPositionVerlet +export GridNHSHandler, PairsNHSHandler, VariableNHSHandler +export VariableSearchRadiusNHS end # module diff --git a/src/general/gpu.jl b/src/general/gpu.jl index 6b07775b4f..825f73e3dd 100644 --- a/src/general/gpu.jl +++ b/src/general/gpu.jl @@ -45,3 +45,31 @@ end function allocate(backend, ELTYPE, size) return Array{ELTYPE, length(size)}(undef, size) end + +function Adapt.adapt_structure(to, nhs::VariableSearchRadiusNHS) + (; search_radii, neighborhood_searches) = nhs + radii = Adapt.adapt_structure(to, search_radii) + searches = Adapt.adapt_structure(to, neighborhood_searches) + + return VariableSearchRadiusNHS(radii, searches) +end + +function Adapt.adapt_structure(to, nhs_handler::GridNHSHandler) + (; search_radii, neighborhood_searches) = nhs_handler + radii = Adapt.adapt_structure(to, search_radii) + searches = Adapt.adapt_structure(to, neighborhood_searches) + + return GridNHSHandler(radii, searches) +end + +function Adapt.adapt_structure(to, nhs_handler::PairsNHSHandler) + searches = Adapt.adapt_structure(to, nhs_handler.neighborhood_searches) + + return PairsNHSHandler(searches) +end + +function Adapt.adapt_structure(to, nhs_handler::VariableNHSHandler) + searches = Adapt.adapt_structure(to, nhs_handler.neighborhood_searches) + + return VariableNHSHandler(searches) +end diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index e33c44f66c..36efaa05f2 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -103,19 +103,21 @@ function create_neighborhood_search(neighborhood_search, system::TotalLagrangian nparticles(neighbor)) end -# === Neighborhood search lookup === +# === Neighborhood search handlers === +# Base type for all neighborhood search handlers. Handlers manage how neighborhood +# searches are stored and retrieved for different interacting systems. abstract type AbstractNHSHandler end +# Handler that stores a distinct neighborhood search for every possible pair of systems. +# Best suited for standard neighborhood searches where each pair has a fixed search radius. struct PairsNHSHandler{NHS} <: AbstractNHSHandler neighborhood_searches::NHS end -function create_neighborhood_search_handler(::Type{<:PairsNHSHandler}, - neighborhood_search, systems) +function PairsNHSHandler(nhs, systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. - searches = Tuple(Tuple(create_neighborhood_search(neighborhood_search, - system, neighbor) + searches = Tuple(Tuple(create_neighborhood_search(nhs, system, neighbor) for neighbor in systems) for system in systems) @@ -126,15 +128,48 @@ function get_neighborhood_search(handler::PairsNHSHandler, system_index, neighbo return handler.neighborhood_searches[system_index][neighbor_index] end +# Handler specifically designed for grid-based neighborhood searches that need to operate +# with multiple, pre-defined unique search radii per system. +struct GridNHSHandler{SI, NHS} <: AbstractNHSHandler + search_radii::SI + neighborhood_searches::NHS +end + +function GridNHSHandler(nhs::PointNeighbors.AbstractNeighborhoodSearch, systems) + # Find a list of search radii that will be requested for each system + search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] + for system in systems) + search_radii = Tuple.(sort.(unique.(search_radii))) + + # For each system, create a neighborhood search for each unique search radius + searches = Tuple(Tuple(copy_neighborhood_search(nhs, search_radius, nparticles(system)) + for search_radius in + search_radii[system_indices(system, systems)]) + for system in systems) + + return GridNHSHandler{typeof(search_radii), typeof(searches)}(search_radii, searches) +end + +function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor_index, + search_radius) + # Find the index of the search radius in the list of search radii for this system + search_radii = handler.search_radii[neighbor_index] + idx = searchsortedfirst(SVector(search_radii), search_radius - eps(search_radius)) + # Return the neighborhood search at that index + return handler.neighborhood_searches[neighbor_index][idx] +end + +# Handler that uses `VariableSearchRadiusNHS` to support different +# compact supports for the same system. struct VariableNHSHandler{NHS} <: AbstractNHSHandler neighborhood_searches::NHS end -function create_neighborhood_search_handler(::Type{<:VariableNHSHandler}, nhs, systems) +function VariableNHSHandler(nhs, systems) # Find a list of search radii that will be requested for each system search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] for system in systems) - search_radii = Tuple.(unique.(search_radii)) + search_radii = Tuple.(sort.(unique.(search_radii))) # For each system, create a neighborhood search for each unique search radius searches = Tuple(VariableSearchRadiusNHS(nhs, @@ -143,7 +178,6 @@ function create_neighborhood_search_handler(::Type{<:VariableNHSHandler}, nhs, s nparticles(system)) for system in systems) - # For each system, create a neighborhood search for each unique search radius return VariableNHSHandler(searches) end @@ -152,6 +186,9 @@ function get_neighborhood_search(handler::VariableNHSHandler, system_index, neig handler.neighborhood_searches[neighbor_index] end +# === Variable search radius neighborhood search === +# A wrapper around multiple neighborhood searches that share the same points +# but operate using different search radii. struct VariableSearchRadiusNHS{SR, NHS} <: PointNeighbors.AbstractNeighborhoodSearch search_radii :: SR neighborhood_searches :: NHS @@ -168,22 +205,25 @@ end @inline requires_update(::VariableSearchRadiusNHS) = (false, true) -# TODO @inline function initialize!(search::VariableSearchRadiusNHS, x, y; parallelization_backend=default_backend(x), eachindex_y=axes(y, 2)) + (; neighborhood_searches) = search + initialize!.(neighborhood_searches) + return search end -# TODO @inline function update!(search::VariableSearchRadiusNHS, x, y; points_moving=(true, true), parallelization_backend=default_backend(x), eachindex_y=axes(y, 2)) + (; neighborhood_searches) = search + update!.(neighborhood_searches) + return search end -# Create a copy of a neighborhood search but with a different search radius function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, search_radius, x, y) search_radii = Tuple(PointNeighbors.search_radius(search) @@ -194,13 +234,8 @@ function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, return VariableSearchRadiusNHS(search_radii, searches) end -function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, - search_radius, n_points; - eachpoint=1:n_points) - return TrivialNeighborhoodSearch{ndims(nhs)}(; search_radius, eachpoint, - periodic_box=nhs.periodic_box) -end - +# Iterate over neighbors by dynamically selecting the correct underlying neighborhood +# search based on the requested search radius. @inline function foreach_neighbor(f, neighbor_system_coords, neighborhood_search::VariableSearchRadiusNHS, point, point_coords, search_radius) @@ -214,36 +249,6 @@ end search_radius) end -struct GridNHSHandler{SI, NHS} <: AbstractNHSHandler - search_radii::SI - neighborhood_searches::NHS -end - -function create_neighborhood_search_handler(::Type{<:GridNHSHandler}, nhs, systems) - # Find a list of search radii that will be requested for each system - search_radii = Tuple([compact_support(system, neighbor) for neighbor in systems] - for system in systems) - search_radii = Tuple.(unique.(search_radii)) - - # For each system, create a neighborhood search for each unique search radius - searches = Tuple(Tuple(copy_neighborhood_search(nhs, search_radius, nparticles(system)) - for search_radius in - search_radii[system_indices(system, systems)]) - for system in systems) - - return GridNHSHandler(search_radii, searches) -end - -function get_neighborhood_search(handler::GridNHSHandler, system_index, neighbor_index, - search_radius) - # Find the index of the search radius in the list of search radii for this system - search_radii = handler.search_radii[neighbor_index] - idx = searchsortedfirst(SVector(search_radii), search_radius - eps(search_radius)) - # Return the neighborhood search at that index - return handler.neighborhood_searches[neighbor_index][idx] -end -# ================================== - @inline function get_neighborhood_search(system, semi) return get_neighborhood_search(system, system, semi) end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 94ddf42af8..20d48de608 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -76,8 +76,8 @@ struct Semidiscretization{BACKEND, S, RU, RV, NSH, UCU, IT} end function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; - nhs_handler=GridNHSHandler, neighborhood_search=GridNeighborhoodSearch{ndims(first(systems))}(), + neighborhood_search_handler=PairsNHSHandler(neighborhood_search, systems), parallelization_backend=PolyesterBackend()) systems = filter(system -> !isnothing(system), systems) @@ -100,10 +100,6 @@ function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)) - neighborhood_search_handler = create_neighborhood_search_handler(nhs_handler, - neighborhood_search, - systems) - # These will be set to true inside the `UpdateCallback`. # Some techniques require the use of this callback, and this flag can be used # to determine if the callback is used in a simulation. diff --git a/test/general/nhs_handler.jl b/test/general/nhs_handler.jl new file mode 100644 index 0000000000..e457f6726e --- /dev/null +++ b/test/general/nhs_handler.jl @@ -0,0 +1,37 @@ +@testset verbose=true "Test NHS Handlers" begin + nhs_handlers = (# PairsNHSHandler, + # GridNHSHandler, + VariableNHSHandler, + ) + + @testset verbose=true "$nhs_handler" for nhs_handler in nhs_handlers + + # Run the `dam_break_gate_2d.jl` file for 0 timesteps to instantiate all systems + sol = trixi_include(@__MODULE__, + joinpath(examples_dir(), "fsi", "dam_break_gate_2d.jl"), + tspan=(0.0, 0.0), + callbacks=OrdinaryDiffEq.CallbackSet() + ) + + # sol.prob.p is the original semidiscretization + semi = sol.prob.p + systems = semi.systems + nhs = TrixiParticles.get_neighborhood_search(first(systems), semi) + neighborhood_search_handler = nhs_handler(nhs, systems) + + new_semi = Semidiscretization(systems..., + neighborhood_search_handler=neighborhood_search_handler, + parallelization_backend=PolyesterBackend()) + + tspan_test = (0.0, 0.1) + ode = semidiscretize(new_semi, tspan_test) + + new_sol = solve(ode, RDPK3SpFSAL49(), + abstol=1e-6, + reltol=1e-4, + dtmax=1e-3, + save_everystep=false) + + @test new_sol.retcode == ReturnCode.Success + end +end \ No newline at end of file From d9224ca96f8bb1aa5f42bc4cf0622df19a80ff9c Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Sat, 28 Mar 2026 15:57:22 +0100 Subject: [PATCH 6/6] Move VariableSearchRadiusNHS to PointNeighbors.jl. Add search_radius arg to `PointNeighbors.foreach_point_neighbor`. Fix the NHS Handler unit test. --- src/TrixiParticles.jl | 3 +- src/general/neighborhood_search.jl | 68 ++---------------------------- src/general/semidiscretization.jl | 3 +- test/general/nhs_handler.jl | 35 +++++++-------- 4 files changed, 25 insertions(+), 84 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c1351e989f..1adb61a9b0 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -40,7 +40,8 @@ using TrixiBase: @trixi_timeit, timer, timeit_debug_enabled, SerialIncrementalUpdate, FullGridCellList, DictionaryCellList, SerialBackend, PolyesterBackend, ThreadsStaticBackend, - ThreadsDynamicBackend, default_backend + ThreadsDynamicBackend, default_backend, + VariableSearchRadiusNHS using PointNeighbors: PointNeighbors, foreach_point_neighbor, copy_neighborhood_search, @threaded using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 36efaa05f2..e9ae55385f 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -5,10 +5,11 @@ function PointNeighbors.foreach_point_neighbor(f, system, neighbor_system, system_coords, neighbor_coords, semi; points=eachparticle(system), - parallelization_backend=semi.parallelization_backend) + parallelization_backend=semi.parallelization_backend, + search_radius=nothing) neighborhood_search = get_neighborhood_search(system, neighbor_system, semi) foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search; - points, parallelization_backend) + points, parallelization_backend, search_radius) end # === Compact support selection === @@ -186,69 +187,6 @@ function get_neighborhood_search(handler::VariableNHSHandler, system_index, neig handler.neighborhood_searches[neighbor_index] end -# === Variable search radius neighborhood search === -# A wrapper around multiple neighborhood searches that share the same points -# but operate using different search radii. -struct VariableSearchRadiusNHS{SR, NHS} <: PointNeighbors.AbstractNeighborhoodSearch - search_radii :: SR - neighborhood_searches :: NHS -end - -function VariableSearchRadiusNHS(nhs_implementation, search_radii, n_particles) - searches = Tuple(PointNeighbors.copy_neighborhood_search(nhs_implementation, - search_radius, n_particles) - for search_radius in search_radii) - return VariableSearchRadiusNHS(search_radii, searches) -end - -@inline Base.ndims(search::VariableSearchRadiusNHS) = ndims(first(search.neighborhood_searches)) - -@inline requires_update(::VariableSearchRadiusNHS) = (false, true) - -@inline function initialize!(search::VariableSearchRadiusNHS, x, y; - parallelization_backend=default_backend(x), - eachindex_y=axes(y, 2)) - (; neighborhood_searches) = search - initialize!.(neighborhood_searches) - - return search -end - -@inline function update!(search::VariableSearchRadiusNHS, x, y; - points_moving=(true, true), - parallelization_backend=default_backend(x), - eachindex_y=axes(y, 2)) - (; neighborhood_searches) = search - update!.(neighborhood_searches) - - return search -end - -function PointNeighbors.copy_neighborhood_search(nhs::VariableSearchRadiusNHS, - search_radius, x, y) - search_radii = Tuple(PointNeighbors.search_radius(search) - for search in nhs.neighborhood_searches) - searches = Tuple(PointNeighbors.copy_neighborhood_search(search) - for search in nhs.neighborhood_searches) - - return VariableSearchRadiusNHS(search_radii, searches) -end - -# Iterate over neighbors by dynamically selecting the correct underlying neighborhood -# search based on the requested search radius. -@inline function foreach_neighbor(f, neighbor_system_coords, - neighborhood_search::VariableSearchRadiusNHS, - point, point_coords, search_radius) - # Find the index of the search radius in the list of search radii for this system - search_radii = neighborhood_search.search_radii - idx = searchsortedfirst(SVector(search_radii), search_radius - eps(search_radius)) - - nhs = neighborhood_search.neighborhood_searches[idx] - - PointNeighbors.foreach_neighbor(f, neighbor_system_coords, nhs, point, point_coords, - search_radius) -end - @inline function get_neighborhood_search(system, semi) return get_neighborhood_search(system, system, semi) end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 20d48de608..22c35940bc 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -77,7 +77,8 @@ end function Semidiscretization(systems::Union{AbstractSystem, Nothing}...; neighborhood_search=GridNeighborhoodSearch{ndims(first(systems))}(), - neighborhood_search_handler=PairsNHSHandler(neighborhood_search, systems), + neighborhood_search_handler=PairsNHSHandler(neighborhood_search, + systems), parallelization_backend=PolyesterBackend()) systems = filter(system -> !isnothing(system), systems) diff --git a/test/general/nhs_handler.jl b/test/general/nhs_handler.jl index e457f6726e..fc94c52c17 100644 --- a/test/general/nhs_handler.jl +++ b/test/general/nhs_handler.jl @@ -1,37 +1,38 @@ +using OrdinaryDiffEq + @testset verbose=true "Test NHS Handlers" begin - nhs_handlers = (# PairsNHSHandler, + nhs_handlers = ( + # PairsNHSHandler, # GridNHSHandler, - VariableNHSHandler, - ) + VariableNHSHandler,) @testset verbose=true "$nhs_handler" for nhs_handler in nhs_handlers - + # Run the `dam_break_gate_2d.jl` file for 0 timesteps to instantiate all systems sol = trixi_include(@__MODULE__, - joinpath(examples_dir(), "fsi", "dam_break_gate_2d.jl"), + joinpath(examples_dir(), "fsi", "dam_break_gate_2d.jl"), tspan=(0.0, 0.0), - callbacks=OrdinaryDiffEq.CallbackSet() - ) + callbacks=OrdinaryDiffEq.CallbackSet()) # sol.prob.p is the original semidiscretization semi = sol.prob.p systems = semi.systems - nhs = TrixiParticles.get_neighborhood_search(first(systems), semi) - neighborhood_search_handler = nhs_handler(nhs, systems) + nhs = TrixiParticles.get_neighborhood_search(first(systems), semi) + neighborhood_search_handler = nhs_handler(nhs, systems) new_semi = Semidiscretization(systems..., neighborhood_search_handler=neighborhood_search_handler, parallelization_backend=PolyesterBackend()) - tspan_test = (0.0, 0.1) + tspan_test = (0.0, 0.1) ode = semidiscretize(new_semi, tspan_test) - - new_sol = solve(ode, RDPK3SpFSAL49(), - abstol=1e-6, - reltol=1e-4, - dtmax=1e-3, - save_everystep=false) + + new_sol = Base.invokelatest(solve, ode, RDPK3SpFSAL49(); + abstol=1e-6, + reltol=1e-4, + dtmax=1e-3, + save_everystep=false) @test new_sol.retcode == ReturnCode.Success end -end \ No newline at end of file +end