`BoundsError` using WENO and stretched grids

[tomchor]

I get a BoundsError when running the following MWE using the latest version of Oceananigans:

using Oceananigans

Nx=Ny=Nz=10

z_faces(k) = k/Nz
grid = RectilinearGrid(topology=(Bounded, Bounded, Bounded),
                       size=(Nx, Ny, Nz),
                       x=(0,1), y=(0,1), 
                       z=z_faces,
                       halo=(3,3,3),
                       )

advection = WENO(grid=grid, order=7)

The error I get is:

ERROR: LoadError: BoundsError: attempt to access 19-element OffsetArray(::Vector{Float64}, -3:15) with eltype Float64 with indices -3:15 at index [-4]
Stacktrace:
  [1] throw_boundserror(A::OffsetArrays.OffsetVector{Float64, Vector{Float64}}, I::Tuple{Int64})
    @ Base ./abstractarray.jl:651
  [2] checkbounds
    @ ./abstractarray.jl:616 [inlined]
  [3] getindex
    @ ~/.julia/packages/OffsetArrays/80Lkc/src/OffsetArrays.jl:435 [inlined]
  [4] #1
    @ ./none:0 [inlined]
  [5] MappingRF
    @ ./reduce.jl:93 [inlined]
  [6] FilteringRF
    @ ./reduce.jl:105 [inlined]
  [7] _foldl_impl(op::Base.FilteringRF{Oceananigans.Advection.var"#2#4"{Int64, Int64}, Base.MappingRF{Oceananigans.Advection.var"#1#3"{Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64, typeof(-)}, Base.BottomRF{typeof(Base.mul_prod)}}}, init::Base._InitialValue, itr::UnitRange{Int64})
    @ Base ./reduce.jl:58
  [8] foldl_impl
    @ ./reduce.jl:48 [inlined]
  [9] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
 [10] #mapfoldl#214
    @ ./reduce.jl:160 [inlined]
 [11] mapfoldl
    @ ./reduce.jl:160 [inlined]
 [12] #mapreduce#218
    @ ./reduce.jl:287 [inlined]
 [13] mapreduce
    @ ./reduce.jl:287 [inlined]
 [14] #prod#225
    @ ./reduce.jl:582 [inlined]
 [15] prod
    @ ./reduce.jl:582 [inlined]
 [16] num_prod
    @ ~/.julia/packages/Oceananigans/W63bs/src/Advection/reconstruction_coefficients.jl:31 [inlined]
 [17] #18
    @ ./none:0 [inlined]
 [18] MappingRF
    @ ./reduce.jl:93 [inlined]
 [19] (::Base.FilteringRF{Oceananigans.Advection.var"#19#23"{Int64}, Base.MappingRF{Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}, Base.BottomRF{typeof(Base.add_sum)}}})(acc::Float64, x::Int64)
    @ Base ./reduce.jl:105
 [20] _foldl_impl(op::Base.FilteringRF{Oceananigans.Advection.var"#19#23"{Int64}, Base.MappingRF{Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}, Base.BottomRF{typeof(Base.add_sum)}}}, init::Base._InitialValue, itr::UnitRange{Int64})
    @ Base ./reduce.jl:62
 [21] foldl_impl(op::Base.FilteringRF{Oceananigans.Advection.var"#19#23"{Int64}, Base.MappingRF{Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}, Base.BottomRF{typeof(Base.add_sum)}}}, nt::Base._InitialValue, itr::UnitRange{Int64})
    @ Base ./reduce.jl:48
 [22] mapfoldl_impl(f::typeof(identity), op::typeof(Base.add_sum), nt::Base._InitialValue, itr::Base.Generator{Base.Iterators.Filter{Oceananigans.Advection.var"#19#23"{Int64}, UnitRange{Int64}}, Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}})
    @ Base ./reduce.jl:44
 [23] mapfoldl(f::Function, op::Function, itr::Base.Generator{Base.Iterators.Filter{Oceananigans.Advection.var"#19#23"{Int64}, UnitRange{Int64}}, Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}}; init::Base._InitialValue)
    @ Base ./reduce.jl:160
 [24] mapfoldl
    @ ./reduce.jl:160 [inlined]
 [25] mapreduce(f::Function, op::Function, itr::Base.Generator{Base.Iterators.Filter{Oceananigans.Advection.var"#19#23"{Int64}, UnitRange{Int64}}, Oceananigans.Advection.var"#18#22"{Int64, typeof(-), Int64, Nothing, Int64, Int64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Int64}}; kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./reduce.jl:287
 [26] mapreduce
    @ ./reduce.jl:287 [inlined]
 [27] #sum#221
    @ ./reduce.jl:501 [inlined]
 [28] sum
    @ ./reduce.jl:501 [inlined]
 [29] #sum#222
    @ ./reduce.jl:528 [inlined]
 [30] sum
    @ ./reduce.jl:528 [inlined]
 [31] #stencil_coefficients#17
    @ ~/.julia/packages/Oceananigans/W63bs/src/Advection/reconstruction_coefficients.jl:61 [inlined]
 [32] create_reconstruction_coefficients(FT::Type, r::Int64, cpu_coord::OffsetArrays.OffsetVector{Float64, Vector{Float64}}, arch::CPU, N::Int64; order::Int64)
    @ Oceananigans.Advection ~/.julia/packages/Oceananigans/W63bs/src/Advection/reconstruction_coefficients.jl:278
 [33] #calc_reconstruction_coefficients#35
    @ ~/.julia/packages/Oceananigans/W63bs/src/Advection/reconstruction_coefficients.jl:268 [inlined]
 [34] top-level scope
    @ none:1
 [35] eval
    @ ./boot.jl:360 [inlined]
 [36] #compute_reconstruction_coefficients#26
    @ ~/.julia/packages/Oceananigans/W63bs/src/Advection/reconstruction_coefficients.jl:227 [inlined]
 [37] WENO(FT::DataType; order::Int64, grid::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, CPU}, zweno::Bool, vector_invariant::Nothing, bounds::Nothing)
    @ Oceananigans.Advection ~/.julia/packages/Oceananigans/W63bs/src/Advection/weno_reconstruction.jl:133
 [38] top-level scope
    @ ~/repos/convddvs/simulations/mwe.jl:12
 [39] include(fname::String)
    @ Base.MainInclude ./client.jl:444
 [40] top-level scope
    @ REPL[2]:1
 [41] top-level scope
    @ ~/.julia/packages/CUDA/DfvRa/src/initialization.jl:52

Some useful notes:

• When I increase the halo to 4 on the z direction the error disappears. Interestingly this code runs even with halo=(1,1,4). (Is the WENO order decreasing as you approach the boundary? Should we throw a warning about this?)
• This behavior is the same if decrease the WENO order and the halo size in the z direction, as expected.
• When I run the same code with the stretched grid (e.g. with z=(0,1)) the error disappears, which seems to point to a bug.

If this isn't a bug and is instead expected behavior, maybe a more useful error would be helpful here?

[tomchor]

CC @simone-silvestri

[simone-silvestri]

The out-of-bounds error comes from the precomputation of stretched coefficients:

The order is not preserved in Bounded directions, so halos larger than 1 are not necessary when time-stepping. This is why also with (1, 1, 1) would be fine if directions are not stretched.

If the direction is stretched, the stretched coefficients are precomputed on every grid point. The error you see there comes from the coefficients near the boundary (which, indeed are not used) that require halo cells to be computed. Note that a 7th order WENO would require 4 halos in Periodic directions

We could remove the pre-computation of boundary coefficients in case of Bounded directions...

[tomchor]

We could remove the pre-computation of boundary coefficients in case of Bounded directions...

As long as it doesn't make the code slower, I'm okay with that. I don't really know enough to have an informed opinion about how to fix it and will defer to whatever you guys think is best :)

[glwagner]

We need halos larger than 1 in Bounded directions, because we should not use short-circuiting logic in GPU kernels.

[navidcy]

Can't we use something like inflate_grid_halo_size inside the WENO constructors?

[glwagner]

still an issue?

[tomchor]

still an issue?

Yes. I just tested this on main

[simone-silvestri]

I do not really like the idea to inflate the grid inside the advection scheme though, if you want you can issue an error message

[glwagner]

Agree --- the issue here is just that we would need halo (4, 4, 4)?

[simone-silvestri]

yep

[simone-silvestri]

The out-of-bounds error comes from the precomputation of stretched coefficients:

The order is not preserved in Bounded directions, so halos larger than 1 are not necessary when time-stepping. This is why also with (1, 1, 1) would be fine if directions are not stretched.

If the direction is stretched, the stretched coefficients are precomputed on every grid point. The error you see there comes from the coefficients near the boundary (which, indeed are not used) that require halo cells to be computed. Note that a 7th order WENO would require 4 halos in Periodic directions

We could remove the pre-computation of boundary coefficients in case of Bounded directions...

Edit: this is wrong, high order halos are always required (also in Bounded directions) because we use ifelse statement to discriminate between low and high order. ifelse executes both branches so we need the halos for the branch that we eventually discard

[glwagner]

Right so just to clarify we should throw an error in the WENO constructor when the halo isn't big enough. This is fine:

using Oceananigans

H = 4

grid = RectilinearGrid(topology=(Bounded, Bounded, Bounded),
                       size = (10, 10, 10),
                       x = (0, 1),
                       y =(0, 1),
                       z = k -> k,
                       halo = (H, H, H))

advection = WENO(grid=grid, order=7)