Skip to content

Commit

Permalink
Update docs for NHS (#368)
Browse files Browse the repository at this point in the history 
* Update docs for NHS

* Rename kwargs

---------

Co-authored-by: Niklas Neher <[email protected]>
  • Loading branch information
@efaulhaber @LasNikas
efaulhaber and LasNikas authored Feb 7, 2024
1 parent 0c934a9 commit b27aa2b
Show file tree
Hide file tree
Showing 4 changed files with 119 additions and 37 deletions.
21 changes: 12 additions & 9 deletions src/general/semidiscretization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,12 @@ the keyword argument `neighborhood_search`. A value of `nothing` means no neighb
By default, the [`GridNeighborhoodSearch`](@ref) is used.
Use [`TrivialNeighborhoodSearch`](@ref) to loop over all particles
(no neighborhood search).
- `periodic_box_min_corner`: In order to use a periodic domain, pass the coordinates
of the domain corner in negative coordinate directions.
- `periodic_box_max_corner`: In order to use a periodic domain, pass the coordinates
of the domain corner in positive coordinate directions.
- `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.
# Examples
```julia
Expand Down Expand Up @@ -95,18 +97,19 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization)
end

function create_neighborhood_search(system, neighbor, ::Val{nothing},
min_corner, max_corner)
periodic_box_min_corner, periodic_box_max_corner)
radius = compact_support(system, neighbor)
TrivialNeighborhoodSearch{ndims(system)}(radius, eachparticle(neighbor),
min_corner=min_corner, max_corner=max_corner)
periodic_box_min_corner=periodic_box_min_corner,
periodic_box_max_corner=periodic_box_max_corner)
end

function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSearch},
min_corner, max_corner)
periodic_box_min_corner, periodic_box_max_corner)
radius = compact_support(system, neighbor)
search = GridNeighborhoodSearch{ndims(system)}(radius, nparticles(neighbor),
min_corner=min_corner,
max_corner=max_corner)
periodic_box_min_corner=periodic_box_min_corner,
periodic_box_max_corner=periodic_box_max_corner)

# Initialize neighborhood search
initialize!(search, initial_coordinates(neighbor))
Expand Down
68 changes: 53 additions & 15 deletions src/neighborhood_search/grid_nhs.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
@doc raw"""
GridNeighborhoodSearch{NDIMS}(search_radius, n_particles)
GridNeighborhoodSearch{NDIMS}(search_radius, n_particles;
periodic_box_min_corner=nothing,
periodic_box_max_corner=nothing)
Simple grid-based neighborhood search with uniform search radius ``d``.
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
Expand All @@ -11,15 +13,49 @@ indexed by the cell index tuple
\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.
where ``x, y, z`` are the space coordinates and ``d`` is the search radius.
To find particles within a radius around a point, only particles in the neighboring cells
are considered.
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 that makes our implementation a lot faster (although less parallelizable).
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.
!!! 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):
```julia
semi = Semidiscretization(system1, system2, system3,
neighborhood_search=GridNeighborhoodSearch)
```
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:
```julia
semi = Semidiscretization(system1, system2,
neighborhood_search=GridNeighborhoodSearch,
periodic_box_min_corner=[0.0, -0.25],
periodic_box_max_corner=[1.0, 0.75]))
```
## References:
- M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán.
Expand All @@ -42,22 +78,24 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, PB}
cell_size :: NTuple{NDIMS, ELTYPE}

function GridNeighborhoodSearch{NDIMS}(search_radius, n_particles;
min_corner=nothing,
max_corner=nothing) where {NDIMS}
periodic_box_min_corner=nothing,
periodic_box_max_corner=nothing) 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 (min_corner === nothing && max_corner === nothing) || search_radius < eps()
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 min_corner !== nothing && max_corner !== nothing
periodic_box = PeriodicBox(min_corner, max_corner)
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
Expand All @@ -73,8 +111,8 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, PB}
"Please use no NHS for very small problems."))
end
else
throw(ArgumentError("`min_corner` and `max_corner` must either be " *
"both `nothing` or both an array or tuple"))
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,
Expand Down Expand Up @@ -333,8 +371,8 @@ function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, u)
search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs))
else
search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs),
min_corner=nhs.periodic_box.min_corner,
max_corner=nhs.periodic_box.max_corner)
periodic_box_min_corner=nhs.periodic_box.min_corner,
periodic_box_max_corner=nhs.periodic_box.max_corner)
end

# Initialize neighborhood search
Expand Down
53 changes: 46 additions & 7 deletions src/neighborhood_search/trivial_nhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,62 @@
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):
```julia
semi = Semidiscretization(system1, system2, system3,
neighborhood_search=TrivialNeighborhoodSearch)
```
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:
```julia
semi = Semidiscretization(system1, system2,
neighborhood_search=TrivialNeighborhoodSearch,
periodic_box_min_corner=[0.0, -0.25],
periodic_box_max_corner=[1.0, 0.75]))
```
"""
struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB}
search_radius :: ELTYPE
eachparticle :: EP
periodic_box :: PB

function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle;
min_corner=nothing,
max_corner=nothing) where {NDIMS}
if (min_corner === nothing && max_corner === nothing) || search_radius < eps()
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 min_corner !== nothing && max_corner !== nothing
periodic_box = PeriodicBox(min_corner, max_corner)
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("`min_corner` and `max_corner` must either be " *
"both `nothing` or both an array or tuple"))
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),
Expand Down
14 changes: 8 additions & 6 deletions test/neighborhood_search/grid_nhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,8 @@

# 3 x 6 cells
nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2),
min_corner=[-0.1, -0.2], max_corner=[0.2, 0.4])
periodic_box_min_corner=[-0.1, -0.2],
periodic_box_max_corner=[0.2, 0.4])

TrixiParticles.initialize!(nhs, coords)

Expand Down Expand Up @@ -173,8 +174,8 @@

# 3 x 6 cells
nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2),
min_corner=[-0.1, -0.2],
max_corner=[0.205, 0.43])
periodic_box_min_corner=[-0.1, -0.2],
periodic_box_max_corner=[0.205, 0.43])

TrixiParticles.initialize!(nhs, coords)

Expand Down Expand Up @@ -220,7 +221,8 @@

# 5 x 1 cells
nhs = GridNeighborhoodSearch{2}(1.0, size(coords, 2),
min_corner=[-1.5, 0.0], max_corner=[2.5, 3.0])
periodic_box_min_corner=[-1.5, 0.0],
periodic_box_max_corner=[2.5, 3.0])

TrixiParticles.initialize!(nhs, coords)

Expand All @@ -239,8 +241,8 @@

# 3 x 6 x 3 cells
nhs = GridNeighborhoodSearch{3}(0.1, size(coords, 2),
min_corner=[-0.1, -0.2, 0.05],
max_corner=[0.2, 0.4, 0.35])
periodic_box_min_corner=[-0.1, -0.2, 0.05],
periodic_box_max_corner=[0.2, 0.4, 0.35])

TrixiParticles.initialize!(nhs, coords)

Expand Down

0 comments on commit b27aa2b

Please sign in to comment.