Interpolation over periodic dimension returns strange result?

[navidcy]

In an attempt to debug #3100 I ran onto this

using Oceananigans

grid = RectilinearGrid(CPU(), size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))

buoyancy = SeawaterBuoyancy(gravitational_acceleration=1, equation_of_state = LinearEquationOfState(thermal_expansion=1, haline_contraction=1))

model = NonhydrostaticModel(; grid, buoyancy, tracers = (:T, :S))

set!(model.velocities.u, 1)

model.velocities.u

gives

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

But then, if I define a new u interpolated on cell centers:

new_u = @at (Center, Center, Center) u
u_new = Field(new_u)
compute!(u_new)
interior(u_new)

I get

4×4×4 view(::Array{Float64, 3}, 4:7, 4:7, 4:7) with eltype Float64:
[:, :, 1] =
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 0.5  0.5  0.5  0.5

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

[:, :, 3] =
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 0.5  0.5  0.5  0.5

[:, :, 4] =
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 0.5  0.5  0.5  0.5

I would expect that the interpolated u_new still would have values 1 everywhere. No?

If I change the topology to Bounded then the interpolated u_new is all 1s.

(The above behaviour occurs with either Julia v1.8 or v1.9.)

[simone-silvestri]

compute! sets the halo region of the final computed field after computation (u_new in this case), not of the operands (u in this case). You need to call fill_halo_regions!(u) before interpolating.

I wonder if we should fill also the operand's halos automatically before computing or we should leave it to the user...

[navidcy]

Thanks! I often forget this!

but then, how do you explain that this test passes without gilling the halos??

https://github.com/CliMA/Oceananigans.jl/blob/4551a78b1f3fe4bb3b238676c128dc751be9b934/test/test_computed_field.jl#L372

[simone-silvestri]

that test is explicitly not testing the points near the boundary

https://github.com/CliMA/Oceananigans.jl/blob/4551a78b1f3fe4bb3b238676c128dc751be9b934/test/test_computed_field.jl#L67-L68

[navidcy]

🙈

[glwagner]

Huh. Seems like we do need to fill the halo regions of the operand for the result of compute! to be correct. Should we re-open this issue?

[navidcy]

Hm... Perhaps we should! Feel free to reopen and change the title to better represent the issue?

[glwagner]

You're getting a wrong result here right? I'm not sure I'm interpreting this correctly.

[glwagner]

I would expect that the interpolated u_new still would have values 1 everywhere. No?

Yes, u_new should have value 1 everywhere here. This is a bug.

I think @simone-silvestri accurately explains why we get this result, but just because we understand something doesn't mean it is correct?

[glwagner]

The issue is a little tricky. Typically we expect abstract operations to be computed during time-stepping. In that case, the halos should be correctly filled.

However, @navidcy expects that abstract operations should be correct at any time and does not expect to have to call fill halo regions. Thus for compute! to be more generally useful to users I think we do want this behavior.

The problem is that fill halo regions can be expensive eg for distributed models. Therefore to both serve expected user behavior and provide a performant interface we perhaps have to add a flag to compute! like fill_halo_regions=false so that computation for output does not trigger extra calls to fill halo regions.

Note @navidcy you can also use the simpler and more transparent

parent(model.velocities.u) .= 1

or just fill!(model.velocities.u, 1).

I think your result would be correct then. But still if we are setting to functions then we need set!.

[navidcy]

Sure! I was just doing what the test does to reproduce why CI fails... :) Let's simplify the test also!

[glwagner]

maybe the kwarg is fill_operand_halo_regions=true. But this requires some extensive changes to compute! since we have to propagate this kwarg throughout the expression tree.

[jagoosw]

I've just come across an issue that I think is probably the same as this:

using Oceananigans

using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Fields: interpolate!

grid = RectilinearGrid(topology = (Periodic, Flat, Flat), size = (8,), extent = (1,))

c = CenterField(grid)
u = Field{Face, Nothing, Nothing}(grid)

c[1:8, 1, 1] = [1:8;]

fill_halo_regions!(c)

interpolate!(u, c)

u[:, 1, 1]

returns

 5.5
 6.5
 7.5
 1.5
 1.5
 2.5
 3.5
 4.5
 5.5
 6.5
 7.5
 1.5
 1.5
 2.5

so I think the i=1 entry is incorrect and should be 4.5. I think the issue is coming from this: https://github.com/CliMA/Oceananigans.jl/blob/00f028bb37f13692e24921588aeb8a9150f6dd55/src/Fields/interpolate.jl#L253-L256 as both -0.5 and 0.5 truncate to 0 and then the assumption is that $i^+$ will be to the right which isn't true in this edge case.

This is fixed if we change to:

i⁺ = i⁻ + 1 * Int(sign(fractional_idx))

and now u is:

 5.5
 6.5
 7.5
 4.5
 1.5
 2.5
 3.5
 4.5
 5.5
 6.5
 7.5
 4.5
 1.5
 2.5