Issue with out-of-bounds access and windowed field indexing

[siddharthabishnu]

Consider the following MWE:

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1), halo = (2,2,2))

julia> my_field = CenterField(grid)

julia> set!(my_field, 1)
2×3×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×2×2 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 6×7×8 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, -1:6) with eltype Float64 with indices -1:4×-1:5×-1:6
    └── max=1.0, min=1.0, mean=1.0

julia> my_field[-2, :, :]
7×8 OffsetArray(::Matrix{Float64}, -1:5, -1:6) with eltype Float64 with indices -1:5×-1:6:
 2.09002e-95  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> my_field[:, -2, :]
6×8 OffsetArray(::Matrix{Float64}, -1:4, -1:6) with eltype Float64 with indices -1:4×-1:6:
 4.06893e233  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 7.49511e247  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 2.09002e-95  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 3.83945e151  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 5.33788e223  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 2.09002e-95  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> my_field[:, :, -2]
6×7 OffsetArray(::Matrix{Float64}, -1:4, -1:5) with eltype Float64 with indices -1:4×-1:5:
   0.0  NaN      0.0  2.42941e-314  5.13924e151   3.5e-323     4.06893e233
   0.0    0.0    0.0  2.66157e-314  5.6593e-314   4.0e-323     7.49511e247
 NaN      0.0  NaN    1.4545e-320   2.98408e-314  6.09933e6    2.09002e-95
   0.0    0.0    0.0  8.46763e165   1.66e-321     4.06893e233  3.83945e151
   0.0  NaN      0.0  6.09933e6     2.75228e-318  7.49089e247  5.33788e223
   0.0    0.0    0.0  9.85509e165   3.0e-323      6.09933e6    2.09002e-95

Even by accessing array elements indices outside the interior and halo ranges, we don’t get out-of-bounds errors. Instead, we obtain junk values for these indices. This behavior applies to windowed fields as well. But why is that the case? What is the reasoning behind this design choice, since I am assuming it’s on purpose?

Secondly, when accessing the interior elements of a windowed field, the indices in the windowed dimension must start from 1 instead of their actual values, which also results in junk values for these elements. Is this also intentional, or should I create a pull request to address it?

julia> my_windowed_field = CenterField(grid, indices=(:, :, 5:6))
2×3×2 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×2×2 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 5:6)
└── data: 6×7×2 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, 5:6) with eltype Float64 with indices -1:4×-1:5×5:6
    └── max=0.0, min=0.0, mean=0.0

julia> set!(my_windowed_field, 1)
2×3×2 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×2×2 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
├── indices: (:, :, 5:6)
└── data: 6×7×2 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, 5:6) with eltype Float64 with indices -1:4×-1:5×5:6
    └── max=1.0, min=1.0, mean=1.0

julia> interior(my_windowed_field, :, :, 5:6)
2×3×2 view(::Array{Float64, 3}, 3:4, 3:5, 5:6) with eltype Float64:
[:, :, 1] =
 1.52533e-52  2.3146e-152  4.56636e257
 2.74012e-57  1.53521e-51  1.77459e248

[:, :, 2] =
 4.71658e-90  2.12781e160  2.90979e-14
 6.58984e-66  4.31647e-61  2.70148e-56

julia> interior(my_windowed_field, :, :, 1:2)
2×3×2 view(::Array{Float64, 3}, 3:4, 3:5, 1:2) with eltype Float64:
[:, :, 1] =
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 2] =
 1.0  1.0  1.0
 1.0  1.0  1.0

cc @glwagner @simone-silvestri @navidcy

[glwagner]

Will replacing an @inbounds with @propagate_inbounds fix the problem? Does this happen with --check-bounds=yes?

[simone-silvestri]

Indeed, I think the fact you do not have an out-of-bounds access error is connected to the @inbounds flag. With regards to the second problem, I think it might be more consistent to have the correct "displaced" indices in interior refer to the valid data chunk instead of the index starting from 1. We need a PR to address this.

[glwagner]

The user should always be working with "grid" indices. So yeah that should be fixed. I'm surprised that view doesn't give us something with grid indices though already using OffsetArray.

[glwagner]

Also just to note, this is our design, so its not really about being "consistent" or not. We just have to implement the design we choose, for the moment we choose that the user works with grid indices.

[siddharthabishnu]

@glwagner and @simone-silvestri, thank you for the clarification. Indeed, this does not happen with --check-bounds=yes, as shown below:

julia> my_field[-2, :, :]
ERROR: BoundsError: attempt to access 6×7×8 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, -1:6) with eltype Float64 with indices -1:4×-1:5×-1:6 at index [-2, OffsetArrays.IdOffsetRange(values=-1:5, indices=-1:5), OffsetArrays.IdOffsetRange(values=-1:6, indices=-1:6)]
Stacktrace:
 [1] throw_boundserror(A::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, I::Tuple{Int64, Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:888 [inlined]
 [4] getindex
   @ ./abstractarray.jl:1291 [inlined]
 [5] getindex(::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, ::Int64, ::Function, ::Function)
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:408
 [6] top-level scope
   @ REPL[23]:1
Some type information was truncated. Use `show(err)` to see complete types.
julia> my_field[:, -2, :]
ERROR: BoundsError: attempt to access 6×7×8 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, -1:6) with eltype Float64 with indices -1:4×-1:5×-1:6 at index [OffsetArrays.IdOffsetRange(values=-1:4, indices=-1:4), -2, OffsetArrays.IdOffsetRange(values=-1:6, indices=-1:6)]
Stacktrace:
 [1] throw_boundserror(A::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, I::Tuple{Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Int64, Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:888 [inlined]
 [4] getindex
   @ ./abstractarray.jl:1291 [inlined]
 [5] getindex(::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, ::Function, ::Int64, ::Function)
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:408
 [6] top-level scope
   @ REPL[24]:1
Some type information was truncated. Use `show(err)` to see complete types.
julia> my_field[:, :, -2]
ERROR: BoundsError: attempt to access 6×7×8 OffsetArray(::Array{Float64, 3}, -1:4, -1:5, -1:6) with eltype Float64 with indices -1:4×-1:5×-1:6 at index [OffsetArrays.IdOffsetRange(values=-1:4, indices=-1:4), OffsetArrays.IdOffsetRange(values=-1:5, indices=-1:5), -2]
Stacktrace:
 [1] throw_boundserror(A::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, I::Tuple{Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Int64})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:888 [inlined]
 [4] getindex
   @ ./abstractarray.jl:1291 [inlined]
 [5] getindex(::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArrays.OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, ::Function, ::Function, ::Int64)
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:408
 [6] top-level scope
   @ REPL[25]:1
Some type information was truncated. Use `show(err)` to see complete types.
julia> interior(my_windowed_field, :, :, 5:6)
ERROR: BoundsError: attempt to access 2×3×2 view(::Array{Float64, 3}, 3:4, 3:5, :) with eltype Float64 at index [1:2, 1:3, 5:6]
Stacktrace:
 [1] throw_boundserror(A::SubArray{Float64, 3, Array{Float64, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] view(::SubArray{Float64, 3, Array{Float64, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, ::Function, ::Function, ::UnitRange{Int64})
   @ Base ./subarray.jl:184
 [4] interior(::Field{…}, ::Function, ::Function, ::UnitRange{…})
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:395
 [5] top-level scope
   @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.

I will create a PR to address the second issue, as suggested by @simone-silvestri.