Bug in determination of indices of parent array by parent_index_range

[siddharthabishnu]

Consider the following MWE:

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(3, 3, 1), x=(0, 1), y=(0, 1), z=(0, 1));

julia> model = HydrostaticFreeSurfaceModel(; grid,
                                           momentum_advection = nothing,
                                           free_surface = ExplicitFreeSurface(; gravitational_acceleration = 10),
                                           closure = nothing,
                                           tracers = nothing,
                                           buoyancy = nothing)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 3×3×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×1 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── free surface: ExplicitFreeSurface with gravitational acceleration 10.0 m s⁻²
├── advection scheme:
│   └── momentum: Nothing
└── coriolis: Nothing

julia> Hc, Nc, Nz = grid.Hx, grid.Nx, grid.Nz
(2, 4, 1)

julia> model.free_surface.η[1:Hc, 1:Nc, Nz+1:Nz+1] .= 1
ERROR: BoundsError: attempt to access 9×9×1 Array{Float64, 3} at index [4:6, 4:6, 3:3]
Stacktrace:
 [1] throw_boundserror(A::Array{Float64, 3}, I::Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] view
   @ ./subarray.jl:184 [inlined]
 [4] offset_windowed_data(data::OffsetArrays.OffsetArray{…}, Loc::Tuple{…}, grid::RectilinearGrid{…}, indices::Tuple{…})
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:248
 [5] view(f::Field{…}, i::UnitRange{…}, j::UnitRange{…}, k::UnitRange{…})
   @ Oceananigans.Fields /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/main/src/Fields/field.jl:316
 [6] maybeview
   @ ./views.jl:148 [inlined]
 [7] dotview(::Field{…}, ::UnitRange{…}, ::UnitRange{…}, ::UnitRange{…})
   @ Base.Broadcast ./broadcast.jl:1244
 [8] top-level scope
   @ REPL[40]:1
Some type information was truncated. Use `show(err)` to see complete types.

The error occurs because the parent array of model.free_surface.η takes the z-halo into consideration, which does not exist for the ssh field. Generally, when dealing with 1D or 2D fields, the parent array should consider halos only in the dimensions along which the field is defined, rather than in all dimensions as currently implemented. PR #3573 aims to close this issue.

[glwagner]

Can you find another way to say this? I can't make sense of it --- how can an array "consider" something?

[navidcy]

What do you mean exactly with "the parent array should consider halos"?

[navidcy]

oh I just saw @glwagner was asking the same..

[navidcy]

Btw I have no idea what's happening in the code you posted above @siddharthabishnu. Why grid.Nx seems to be 4 and grid.Hx seems to be 2 ???

When I run this on main I get:

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(3, 3, 1), x=(0, 1), y=(0, 1), z=(0, 1));

julia> model = HydrostaticFreeSurfaceModel(; grid,
                                           momentum_advection = nothing,
                                           free_surface = ExplicitFreeSurface(; gravitational_acceleration = 10),
                                           closure = nothing,
                                           tracers = nothing,
                                           buoyancy = nothing)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 3×3×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×1 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── free surface: ExplicitFreeSurface with gravitational acceleration 10.0 m s⁻²
├── advection scheme: 
│   └── momentum: Nothing
└── coriolis: Nothing

julia> Hc, Nc, Nz = grid.Hx, grid.Nx, grid.Nz
(3, 3, 1)

So given this ambiguity, I then don't even know what model.free_surface.η[1:Hc, 1:Nc, Nz+1:Nz+1] is trying to do....

[glwagner]

I think this is a better MWE:

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(3, 3, 1), x=(0, 1), y=(0, 1), z=(0, 1));

julia> model = HydrostaticFreeSurfaceModel(; grid);

julia> view(model.free_surface.η, :, :, 2:2)
ERROR: BoundsError: attempt to access 9×9×1 Array{Float64, 3} at index [1:9, 1:9, 3:3]
Stacktrace:
 [1] throw_boundserror(A::Array{Float64, 3}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}})
   @ Base ./abstractarray.jl:734
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] view
   @ ./subarray.jl:179 [inlined]
 [4] offset_windowed_data(data::OffsetArray{Float64, 3, Array{…}}, Loc::Tuple{DataType, DataType, DataType}, grid::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetVector{…}, OffsetVector{…}, OffsetVector{…}, CPU}, indices::Tuple{Colon, Colon, UnitRange{…}})
   @ Oceananigans.Fields ~/Projects/Oceananigans.jl/src/Fields/field.jl:248
 [5] view(f::Field{Center, Center, Face, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, i::Function, j::Function, k::UnitRange{Int64})
   @ Oceananigans.Fields ~/Projects/Oceananigans.jl/src/Fields/field.jl:316
 [6] top-level scope
   @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

We should be able to call view(model.free_surface.η, :, :, 2:2) since the free surface is a windowed field with vertical indices 2:2 (the index slice corresponding to the top surface of the model...)

[glwagner]

We don't even need to make a model (that makes the MWE take longer than it should):

julia> grid = RectilinearGrid(size=(1, 2, 3), x=(0, 1), y=(0, 1), z=(0, 1))
1×2×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×2×3 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=1.0
├── Periodic y ∈ [0.0, 1.0) regularly spaced with Δy=0.5
└── Bounded  z ∈ [0.0, 1.0] regularly spaced with Δz=0.333333

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

julia> view(windowed, :, :, 2:2)
ERROR: BoundsError: attempt to access 3×6×1 Array{Float64, 3} at index [1:3, 1:6, 5:5]
Stacktrace:
 [1] throw_boundserror(A::Array{Float64, 3}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}})
   @ Base ./abstractarray.jl:734
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] view
   @ ./subarray.jl:179 [inlined]
 [4] offset_windowed_data(data::OffsetArray{Float64, 3, Array{…}}, Loc::Tuple{DataType, DataType, DataType}, grid::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, Float64, OffsetVector{…}, OffsetVector{…}, OffsetVector{…}, CPU}, indices::Tuple{Colon, Colon, UnitRange{…}})
   @ Oceananigans.Fields ~/Projects/Oceananigans.jl/src/Fields/field.jl:248
 [5] view(f::Field{Center, Center, Center, Nothing, RectilinearGrid{…}, Tuple{…}, OffsetArray{…}, Float64, FieldBoundaryConditions{…}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{…}}, i::Function, j::Function, k::UnitRange{Int64})
   @ Oceananigans.Fields ~/Projects/Oceananigans.jl/src/Fields/field.jl:316
 [6] top-level scope
   @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

[glwagner]

The bug is here:

Oceananigans.jl/src/Grids/grid_utils.jl

Line 226 in ed77780

parent_index_range(index::UnitRange, loc, topo, halo) = index .+ interior_parent_offset(loc, topo, halo)

I think it should just be

parent_index_range(index::UnitRange, loc, topo, halo) = UnitRange(1, length(index))

🤦

[siddharthabishnu]

@navidcy, the title's been fixed and the issue identified with @glwagner. It's been resolved in commit bbfc656. Regarding your earlier query, the output of Hc, Nc, Nz = grid.Hx, grid.Nx, grid.Nz should be (3, 3, 1) rather than (2, 4, 1). Initially, I used a ConformalCubedSphereGrid with panel_size = (4, 4, 1) and horizontal_direction_halo = 2, but later switched to a simpler RectilinearGrid with size = (3, 3, 1) and overlooked updating the line you referenced in the MWE.

[navidcy]

@siddharthabishnu, just copy-paste the input and output from the REPL directly. Don't edit the code snippet otherwise there is the danger of making a typo/mistake.

[navidcy]

btw, this -> #3260 seems related or duplicate of this issue?

[navidcy]

It's been resolved in commit bbfc656.

Did you test that the issue is resolved with this edit or you assume that @glwagner was right. Because many tests seems to have broken in the PR after this change.

[siddharthabishnu]

@siddharthabishnu, just copy-paste the input and output from the REPL directly. Don't edit the code snippet otherwise there is the danger of making a typo/mistake.

Agreed. I learnt that the hard way. I won't do it in the future.

[siddharthabishnu]

It's been resolved in commit bbfc656.

Did you test that the issue is resolved with this edit or you assume that @glwagner was right. Because many tests seems to have broken in the PR after this change.

@navidcy, I concur that the issue remains unresolved. Initially, I thought it was addressed during the manual bounds-checking that @glwagner and I conducted together. However, since other tests are failing, we cannot proceed with merging this PR at this time.

[navidcy]

Well sure, we'll merge when it's sorted out.