Skip to content

Fix bug due to rounding errors in periodicity computation #146

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
efaulhaber wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Fix bug due to rounding errors in periodicity computation #146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 0 additions & 29 deletions src/cell_lists/full_grid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,35 +81,6 @@ function FullGridCellList(; min_corner, max_corner,
return FullGridCellList(cells, linear_indices, min_corner, max_corner)
end

@inline function cell_coords(coords, periodic_box::Nothing, cell_list::FullGridCellList,
cell_size)
(; min_corner) = cell_list

# Subtract `min_corner` to offset coordinates so that the min corner of the grid
# corresponds to the (1, 1, 1) cell.
return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1
end

@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list::FullGridCellList,
cell_size)
# Subtract `periodic_box.min_corner` to offset coordinates so that the min corner
# of the grid corresponds to the (0, 0, 0) cell.
offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner

# Add 2, so that the min corner will be the (2, 2, 2)-cell.
# With this, we still have one padding layer in each direction around the periodic box,
# just like without using a periodic box.
# This is not needed for finding neighbor cells, but to make the bounds check
# work the same way as without a periodic box.
return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 2
end

@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells,
cell_list::FullGridCellList)
# 2-based modulo to match the indexing of the periodic box explained above.
return mod.(cell_index .- 2, n_cells) .+ 2
end

function Base.empty!(cell_list::FullGridCellList)
(; cells) = cell_list

Expand Down
42 changes: 21 additions & 21 deletions src/nhs_grid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -579,37 +579,37 @@ end
end

@inline function periodic_cell_index(cell_index, neighborhood_search)
(; n_cells, periodic_box, cell_list) = neighborhood_search
(; n_cells, periodic_box) = neighborhood_search

periodic_cell_index(cell_index, periodic_box, n_cells, cell_list)
periodic_cell_index(cell_index, periodic_box, n_cells)
end

@inline periodic_cell_index(cell_index, ::Nothing, n_cells, cell_list) = cell_index

@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells, cell_list)
# 1-based modulo
return mod1.(cell_index, n_cells)
@inline periodic_cell_index(cell_index, ::Nothing, n_cells) = cell_index

@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells)
# With the `DictionaryCellList`, the grid is conceptually infinite,
# so we can use any offset for the modulo operation.
# With the `FullGridCellList`, we need 2-based modulo, so that the min corner
# will be the (2, 2, 2)-cell.
# With this, we still have one padding layer in each direction around the periodic box,
# just like without using a periodic box.
# This is not needed for finding neighbor cells, but to make the bounds check
# work the same way as without a periodic box (cells in bounds are 2:end-1).
return mod.(cell_index .- 2, n_cells) .+ 2
end

@inline function cell_coords(coords, neighborhood_search)
(; periodic_box, cell_list, cell_size) = neighborhood_search
# Compute cell coordinates ignoring periodicity at first
nonperiodic_cell_coords_ = nonperiodic_cell_coords(coords, neighborhood_search)

return cell_coords(coords, periodic_box, cell_list, cell_size)
# Now return (without periodicity) or convert to periodic cell index (with periodicity)
return periodic_cell_index(nonperiodic_cell_coords_, neighborhood_search)
end

@inline function cell_coords(coords, periodic_box::Nothing, cell_list, cell_size)
return Tuple(floor_to_int.(coords ./ cell_size))
end
@inline function nonperiodic_cell_coords(coords, neighborhood_search)
(; cell_size) = neighborhood_search

@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list, cell_size)
# Subtract `min_corner` to offset coordinates so that the min corner of the periodic
# box corresponds to the (0, 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

# Add one for 1-based indexing. The min corner will be the (1, 1, 1)-cell.
return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 1
return Tuple(floor_to_int.(coords ./ cell_size))
end
Comment on lines +609 to 613

function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, n_points;
Expand Down
Loading