Fix interpolation on a stretched grid with a singleton dimension that is not Flat

[simone-silvestri]

In main we cannot do this because of not loading the full coordinate array, connected to the current change in the interpolation for LatitudeLongitudeGrids. This is a very peculiar configuration that does not happen usually, but we hit it in ClimaOcean where a prescribed atmosphere has stretched coordinates that might have a singleton dimension and expose this issue.

MWE (NOTE: this will fail only with --check-bounds=yes!):

# Interpolation from 1-cell not-flat non-uniform grid
from_grid = LatitudeLongitudeGrid(size = (1, 1, 1), latitude = [-1, 1], longitude = [-1, 1], z = (0, 1))
to_grid   = LatitudeLongitudeGrid(size = 1, latitude = 0, longitude = 0, z = (0, 1), topology = (Flat, Flat, Bounded))

ff = CenterField(from_grid)
ft = CenterField(to_grid)

fill!(ff, 1)
interpolate!(ft, ff)

fails with:

julia> interpolate!(ft, ff)
ERROR: BoundsError: attempt to access 1-element view(OffsetArray(::Vector{Float64}, 0:2), 1:1) with eltype Float64 at index [2]
Stacktrace:
  [1] throw_boundserror(A::SubArray{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Tuple{UnitRange{Int64}}, true}, I::Tuple{Int64})
    @ Base ./abstractarray.jl:737
  [2] checkbounds
    @ ./abstractarray.jl:702 [inlined]
  [3] getindex
    @ ./subarray.jl:322 [inlined]
  [4] fractional_x_index
    @ ~/development/Oceananigans.jl/src/Fields/interpolate.jl:129 [inlined]
  [5] _fractional_indices
    @ ~/development/Oceananigans.jl/src/Fields/interpolate.jl:203 [inlined]
  [6] fractional_indices
    @ ~/development/Oceananigans.jl/src/Fields/interpolate.jl:192 [inlined]
  [7] interpolate
    @ ~/development/Oceananigans.jl/src/Fields/interpolate.jl:262 [inlined]
  [8] macro expansion
    @ ~/development/Oceananigans.jl/src/Fields/interpolate.jl:362 [inlined]
  [9] cpu__interpolate!
    @ ~/.julia/packages/KernelAbstractions/mD0Rj/src/macros.jl:295 [inlined]
 [10] cpu__interpolate!(__ctx__::KernelAbstractions.CompilerMetadata{…}, to_field::Field{…}, to_grid::LatitudeLongitudeGrid{…}, to_location::Tuple{…}, from_field::Field{…}, from_grid::Oceananigans.Grids.ZRegularLLG{…}, from_location::Tuple{…})
    @ Oceananigans.Fields ./none:0
 [11] __thread_run(tid::Int64, len::Int64, rem::Int64, obj::KernelAbstractions.Kernel{…}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{…}, args::Tuple{…}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck)
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/mD0Rj/src/cpu.jl:144
 [12] __run(obj::KernelAbstractions.Kernel{…}, ndrange::Nothing, iterspace::KernelAbstractions.NDIteration.NDRange{…}, args::Tuple{…}, dynamic::KernelAbstractions.NDIteration.NoDynamicCheck, static_threads::Bool)
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/mD0Rj/src/cpu.jl:111
 [13] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/mD0Rj/src/cpu.jl:45
 [14] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…})
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/mD0Rj/src/cpu.jl:38
 [15] _launch!(::CPU, ::LatitudeLongitudeGrid{…}, ::Oceananigans.Utils.KernelParameters{…}, ::Function, ::Field{…}, ::LatitudeLongitudeGrid{…}, ::Vararg{…}; exclude_periphery::Bool, reduced_dimensions::Tuple{}, active_cells_map::Nothing)
    @ Oceananigans.Utils ~/development/Oceananigans.jl/src/Utils/kernel_launching.jl:298
 [16] _launch!
    @ ~/development/Oceananigans.jl/src/Utils/kernel_launching.jl:275 [inlined]
 [17] launch!
    @ ~/development/Oceananigans.jl/src/Utils/kernel_launching.jl:258 [inlined]
 [18] interpolate!(to_field::Field{…}, from_field::Field{…})
    @ Oceananigans.Fields ~/development/Oceananigans.jl/src/Fields/interpolate.jl:393
 [19] top-level scope
    @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

This PR fixes the issue (and introduces a very simple test based on the above MWE)

[glwagner]

not sure but you may want to make with_halos a positional argument for this to be efficient in a kernel

something like

xnodes(grid, lx, ly, lz, with_halos)

plus

xnodes(grid, lx, ly, lz; with_halos=false) = xnodes(grid, lx, ly, lz, with_halos)

to preserve default behavior

[simone-silvestri]

sounds good, no keyword arguments in kernels is definitely preferred