interpolate doesn't work with Nothing locations

[glwagner]

The function interpolate is currently mostly used for Lagrangian particle advection:

Oceananigans.jl/src/Fields/interpolate.jl

Lines 186 to 198 in d0b7ec8

	@inline function interpolate(field, LX, LY, LZ, grid, x, y, z)
	i, j, k = fractional_indices(x, y, z, (LX, LY, LZ), grid)
	# We use mod and trunc as CUDA.modf is not defined.
	# For why we use Base.unsafe_trunc instead of trunc see:
	# https://github.com/CliMA/Oceananigans.jl/issues/828
	# https://github.com/CliMA/Oceananigans.jl/pull/997
	ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i)
	η, j = mod(j, 1), Base.unsafe_trunc(Int, j)
	ζ, k = mod(k, 1), Base.unsafe_trunc(Int, k)
	return _interpolate(field, ξ, η, ζ, i+1, j+1, k+1)
	end

However, it's also potentially useful for many other things. Right now though it doesn't work (or make sense) for fields with Nothing locations:

source_grid = RectilinearGrid(size=(2, 2), x=(0, 1), y=(0, 1), topology=(Periodic, Periodic, Flat))
s = Field{Center, Center, Nothing}(source_grid)
set!(s, (x, y) -> x + y)
loc = Tuple(L() for L in location(s))
x = y = 0.67
z = 0
interpolate(s, loc..., source_grid, x, y, z)

gives

julia> interpolate(s, loc..., source_grid, x, y, z)
ERROR: BoundsError: attempt to access Tuple{Float64, Float64} at index [3]
Stacktrace:
 [1] getindex(t::Tuple, i::Int64)
   @ Base ./tuple.jl:29
 [2] fractional_z_index
   @ ~/.julia/packages/Oceananigans/ususO/src/Fields/interpolate.jl:110 [inlined]
 [3] fractional_indices
   @ ~/.julia/packages/Oceananigans/ususO/src/Fields/interpolate.jl:133 [inlined]
 [4] interpolate(field::Field{Center, Center, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Periodic, Flat, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}}, LX::Center, LY::Center, LZ::Nothing, grid::RectilinearGrid{Float64, Periodic, Periodic, Flat, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, CPU}, x::Float64, y::Float64, z::Int64)
   @ Oceananigans.Fields ~/.julia/packages/Oceananigans/ususO/src/Fields/interpolate.jl:187
 [5] top-level scope
   @ REPL[19]:1

It doesn't really make sense to have to pass a z location into interpolate. I think we should refactor this function to read something like

interpolate(to_node, from_field, from_location, from_grid)

Then to_node can be a tuple of appropriate length for the from_location and from_grid topology.

[glwagner]

Functions like fractional_z_index assume 3D location / node:

Oceananigans.jl/src/Fields/interpolate.jl

Lines 109 to 113 in d0b7ec8

	@inline function fractional_z_index(z::FT, locs, grid::ZRegGrid) where FT
	z₀ = @inbounds node(1, 1, 1, grid, locs...)[3]
	Δz = @inbounds zspacings(grid, locs...)
	return convert(FT, (z - z₀) / Δz)
	end

[glwagner]

Oceananigans.jl/src/Fields/interpolate.jl

Line 111 in d0b7ec8

Δz = @inbounds zspacings(grid, locs...)

Shouldn't this be zspacing? How do we know that this will always be a scalar appropriate for using in division? I don't it's documented or guaranteed anywhere that "ZReg" also implies that zspacings returns a scalar. It's confusing to begin with that this pluralized function returns a single value...

[glwagner]

@simone-silvestri @jagoosw

[glwagner]

Is there some kind of performance benefit of

Oceananigans.jl/src/Fields/interpolate.jl

Lines 63 to 67 in d0b7ec8

	@inline function fractional_x_index(x::FT, locs, grid::XRegRectilinearGrid) where FT
	x₀ = @inbounds node(1, 1, 1, grid, locs...)[1]
	Δx = @inbounds xspacings(grid, locs...)
	return convert(FT, (x - x₀) / Δx)
	end

compared to

Oceananigans.jl/src/Fields/interpolate.jl

Lines 75 to 81 in d0b7ec8

	@inline function fractional_x_index(x, locs, grid)
	loc = @inbounds locs[1]
	Tx = topology(grid, 1)()
	L = length(loc, Tx, grid.Nx)
	xn = @inbounds nodes(grid, locs)[1]
	return fractional_index(L, x, xn) - 1
	end

?

Why do we have both?

[jagoosw]

I'm not sure if we've tested but I've assumed there is a performance benefit to the simpler version for regularly spaced grids rather than using the binary search.

It would probably be sensible to change the differentiation between the methods to just fractional_index though.

[glwagner]

I'm not sure if we've tested but I've assumed there is a performance benefit to the simpler version for regularly spaced grids rather than using the binary search.

@simone-silvestri have you ever benchmarked this?

[simone-silvestri]

I have never tried benchmarking this. maybe the gain in performance is negligible. I guess it will depend on the number of particles and the size of the grid.