Skip to content
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.

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

Kernel size calculation for sliced field halo filling seems to be wrong #3473

Open
glwagner opened this issue Feb 13, 2024 · 1 comment
Open
Labels
bug 🐞 Even a perfect program still has bugs

Comments

@glwagner
Copy link
Member

To fill the halos of sliced fields we call fill_halo_size:

@inline fill_halo_size(::OffsetArray, ::WEB, ::Tuple{<:Any, <:Colon, <:Colon}, args...) = :yz
@inline fill_halo_size(::OffsetArray, ::SNB, ::Tuple{<:Colon, <:Any, <:Colon}, args...) = :xz
@inline fill_halo_size(::OffsetArray, ::TBB, ::Tuple{<:Colon, <:Colon, <:Any}, args...) = :xy

However, that code incorrectly assumes that an indices of : means that the whole dimension should be filled. This is wrong for fields with Nothing locations, which have only 1 point in the Nothing direction, but can still have indices :.

There is also logic for Nothing locations:

# If the index is a Colon and the location is _NOT_ a `Nothing` (i.e. not a `ReducedField`),
# then fill the whole boundary, otherwise fill the size of the corresponding array
@inline whole_halo(idx, loc) = false
@inline whole_halo(idx, ::Nothing) = false
@inline whole_halo(::Colon, ::Nothing) = false
@inline whole_halo(::Colon, loc) = true
# Calculate kernel size
@inline fill_halo_size(c::OffsetArray, ::WEB, idx, bc, loc, grid) =
@inbounds (ifelse(whole_halo(idx[2], loc[2]), size(grid, 2), size(c, 2)), ifelse(whole_halo(idx[3], loc[3]), size(grid, 3), size(c, 3)))
@inline fill_halo_size(c::OffsetArray, ::SNB, idx, bc, loc, grid) =
@inbounds (ifelse(whole_halo(idx[1], loc[1]), size(grid, 1), size(c, 1)), ifelse(whole_halo(idx[3], loc[3]), size(grid, 3), size(c, 3)))
@inline fill_halo_size(c::OffsetArray, ::TBB, idx, bc, loc, grid) =
@inbounds (ifelse(whole_halo(idx[1], loc[1]), size(grid, 1), size(c, 1)), ifelse(whole_halo(idx[2], loc[2]), size(grid, 2), size(c, 2)))

But this is never used because of the incorrect assumption about the meaning of indices.

One possibility is to convert indices in Nothing directions to nothing. That would permit dispatching on the type of the indices like we're trying to do. Another possibility is to fix the logic here. Dispatch is annoying because loc is a tuple, but it can be done. A third possibility is just to avoid the symbol sizes :yz, xz and :xy. That last seems maybe like the best choice.

This bug would prevent anyone from using ImmersedBoundaryGrid. So I'm confused.

@glwagner glwagner added the bug 🐞 Even a perfect program still has bugs label Feb 13, 2024
@glwagner
Copy link
Member Author

@simone-silvestri

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐞 Even a perfect program still has bugs
Projects
None yet
Development

No branches or pull requests

1 participant