Fix spacings methods for `ImmersedBoundaryGrid`s

[tomchor]

The spacings functions aren't working for ImmersedBoundaryGrids on main. This PR fixes that.

On main:

julia> using Oceananigans

julia> grid_base = RectilinearGrid(size=(4, 5, 6), extent=(1, 1, 1));

julia> ibg = ImmersedBoundaryGrid(grid_base, GridFittedBottom((x, y) -> 1))
4×5×6 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: GridFittedBottom(min(h)=1.00e+00, max(h)=1.00e+00)
├── underlying_grid: 4×5×6 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.25
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.2
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.166667

julia> xspacings(ibg, Center())
ERROR: MethodError: no method matching xspacings(::ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded, RectilinearGrid{Float64, Periodic, Periodic, Bounded, 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}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, GridFittedBottom{OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, ::Center)
Closest candidates are:
  xspacings(::Any, ::Any, ::Any, ::Any; with_halos) at ~/repos/Oceananigans.jl/src/Grids/grid_utils.jl:363
  xspacings(::Any, ::ImmersedBoundaryGrid) at ~/repos/Oceananigans.jl/src/ImmersedBoundaries/immersed_grid_metrics.jl:37
  xspacings(::Oceananigans.Grids.XRegRectilinearGrid, ::Center; with_halos) at ~/repos/Oceananigans.jl/src/Grids/rectilinear_grid.jl:473
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[18]:1

On this branch:

julia> xspacings(ibg, Center())
0.25

julia> yspacings(ibg, Center())
0.2

julia> zspacings(ibg, Center())
0.16666666666666666

[tomchor]

Note: I couldn't find any tests for this, so it'd probably be a good idea to add some here before merging.

[navidcy]

hm... so we get this error although we added a method.

Is it because we defined, e.g., xspacing inside Oceananigans.Grids and then added a new method in Oceananigans.ImmersedBoundaries but in the unit tests we only called: using Oceananigans.Grids: xspacing?

Grid initialization: Error During Test at /var/lib/buildkite-agent/builds/tartarus-7/clima/oceananigans/test/test_grids.jl:203
Test threw exception
Expression: minimum_xspacing(grid) ≈ FT(π / 3)
MethodError: xspacing(::Int64, ::Int64, ::Int64, ::ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, ::Center, ::Center, ::Center) is ambiguous. Candidates:
xspacing(i, j, k, grid::ImmersedBoundaryGrid, args...) in Oceananigans.ImmersedBoundaries at /var/lib/buildkite-agent/builds/tartarus-7/clima/oceananigans/src/ImmersedBoundaries/immersed_grid_metrics.jl:40
xspacing(i, j, k, grid, ::Center, ::Center, ::Center) in Oceananigans.Operators at /var/lib/buildkite-agent/builds/tartarus-7/clima/oceananigans/src/Operators/spacings_and_areas_and_volumes.jl:264
Possible fix, define
xspacing(::Any, ::Any, ::Any, ::ImmersedBoundaryGrid, ::Center, ::Center, ::Center)
Stacktrace:
[1] getindex
@ ~/builds/tartarus-7/clima/oceananigans/src/AbstractOperations/kernel_function_operation.jl:73 [inlined]
[2] getindex
@ ~/builds/tartarus-7/clima/oceananigans/src/AbstractOperations/conditional_operations.jl:120 [inlined]
[3] getindex
@ ./subarray.jl:282 [inlined]
[4] _getindex
@ ./abstractarray.jl:1291 [inlined]
[5] getindex
@ ./abstractarray.jl:1241 [inlined]
[6] map!(f::typeof(identity), dest::SubArray{Float32, 3, Array{Float32, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}, false}, A::SubArray{Float32, 3, ConditionalOperation{Center, Center, Center, KernelFunctionOperation{Center, Center, Center, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Float32, typeof(xspacing), Tuple{Center, Center, Center}}, typeof(identity), ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Oceananigans.ImmersedBoundaries.NotImmersed{typeof(Oceananigans.AbstractOperations.truefunc)}, Float64, Float32}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}, false})
@ Base ./abstractarray.jl:2926
[7] mapfirst!(f::typeof(identity), R::SubArray{Float32, 3, Array{Float32, 3}, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}, false}, A::ConditionalOperation{Center, Center, Center, KernelFunctionOperation{Center, Center, Center, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Float32, typeof(xspacing), Tuple{Center, Center, Center}}, typeof(identity), ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Oceananigans.ImmersedBoundaries.NotImmersed{typeof(Oceananigans.AbstractOperations.truefunc)}, Float64, Float32})
@ Base ./reducedim.jl:274
[8] initialize_reduced_field!(#unused#::typeof(minimum!), f::Function, r::Field{Nothing, Nothing, Nothing, Nothing, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Tuple{Colon, Colon, Colon}, OffsetArray{Float32, 3, Array{Float32, 3}}, Float32, FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, c::ConditionalOperation{Center, Center, Center, KernelFunctionOperation{Center, Center, Center, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Float32, typeof(xspacing), Tuple{Center, Center, Center}}, typeof(identity), ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Oceananigans.ImmersedBoundaries.NotImmersed{typeof(Oceananigans.AbstractOperations.truefunc)}, Float64, Float32})
@ Oceananigans.Fields ~/builds/tartarus-7/clima/oceananigans/src/Fields/field.jl:556
[9] minimum(f::Function, c::KernelFunctionOperation{Center, Center, Center, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Float32, typeof(xspacing), Tuple{Center, Center, Center}}; condition::Nothing, mask::Float64, dims::Function)
@ Oceananigans.Fields ~/builds/tartarus-7/clima/oceananigans/src/Fields/field.jl:647
[10] minimum
@ ~/builds/tartarus-7/clima/oceananigans/src/Fields/field.jl:637 [inlined]
[11] #minimum#45
@ ~/builds/tartarus-7/clima/oceananigans/src/Fields/field.jl:657 [inlined]
[12] minimum(c::KernelFunctionOperation{Center, Center, Center, ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, Float32, typeof(xspacing), Tuple{Center, Center, Center}})
@ Oceananigans.Fields ~/builds/tartarus-7/clima/oceananigans/src/Fields/field.jl:657
[13] minimum_spacing(dir::Symbol, grid::ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU}, ℓx::Center, ℓy::Center, ℓz::Center)
@ Oceananigans.Grids ~/builds/tartarus-7/clima/oceananigans/src/Grids/grid_utils.jl:408
[14] minimum_xspacing(grid::ImmersedBoundaryGrid{Float32, Periodic, Periodic, Bounded, RectilinearGrid{Float32, Periodic, Periodic, Bounded, Float32, Float32, Float32, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, CPU}, GridFittedBottom{OffsetMatrix{Float32, Matrix{Float32}}, Oceananigans.ImmersedBoundaries.CenterImmersedCondition}, Nothing, CPU})
@ Oceananigans.Grids ~/builds/tartarus-7/clima/oceananigans/src/Grids/grid_utils.jl:429
[15] macro expansion
@ /storage5/buildkite-agent/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:464 [inlined]
[16] test_regular_rectilinear_xnode_ynode_znode_and_spacings(arch::CPU, FT::Type{Float32})
@ Main ~/builds/tartarus-7/clima/oceananigans/test/test_grids.jl:203


<br class="Apple-interchange-newline" style="caret-color: rgb(0, 0, 0); color: rgb(0, 0, 0); font-style: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none;">

[navidcy]

OK, I think I found the culprit.

[tomchor]

I'm going to try and put some some into this, but I'm a bit unsure of what to do. First, to summarize what the current state on main is. Taking a RectilinearGrid and zspacing at the Center as an example.

---

Current behavior for RectilinearGrid

Without immersed boundaries, this is pretty straightforward. The call to zspacings() dispatches to

https://github.com/CliMA/Oceananigans.jl/blob/56fc69bcb08623424afe25e39ae616c29208095e/src/Grids/rectilinear_grid.jl#L477-L478

which directly accesses the grid's Δzᵃᵃᶜ property.

The call to zspacing() goes through some metaprogramming here but ultimately also directly accesses the grid's Δzᵃᵃᶜ property here:

https://github.com/CliMA/Oceananigans.jl/blob/56fc69bcb08623424afe25e39ae616c29208095e/src/Operators/spacings_and_areas_and_volumes.jl#L104-L112

---

Current behavior for ImmersedBoundaryGrid

For the same grid wrapped around an ImmersedBoundaryGrid on main, when calling zspacings() we get the error at the top comment (MethodError: no method matching; which this PR aims to fix). The call to zspacing() works the same way as with RectilinearGrid, with the exception of an extra dispatch here.

---

Summary of the PR as of this post

This PR is originally trying to fix the issue with zspacings(::IBG) simply by adding

zspacings(grid::IBG, args...; kwargs...) = zspacings(grid.underlying_grid, args...; kwargs...)

If I understand correctly, @glwagner's point is that this is wrong in cases where cells have a fraction of "wet" volume and a fraction of "immersed solid" volume. (For now the two such cases in sight are PartialCellBottom and CutCellBottom (for which there's open PR https://github.com/CliMA/Oceananigans.jl/pull/3146).)

Again, if I understand correctly, the simple solution would be to redirect zspacings() calls to

1. Create a KernelFunctionOperation using zspacing(), which we know has the correct behavior for all cases
2. Calculate and collect the values in a Array
3. Return to user

The main challenge here (imo) is that

• if we wanna keep the user-interface simple (e.g. return a float or a 1D array from calls to zspacings(::RectilinearGrid) instead of a Field or a 3D array), then we lose consistency (because in the general case a call to zspacings(::ImmersedBoundaryGrid) must return a 3D array (or Field, or AbstractOp...)).
• if we wanna keep the user-interface consistent (i.e. always return the same type of object) then we lose on simplicity (e.g. a call to zspacings(::RectilinearGrid) would always return something 3D, even if the grid is regular).

@navidcy @glwagner is this a fair assessment of the situation? Feel free to edit the text above if not.

[tomchor]

My take is that I don't think there's a way around the choices I just mentioned in the previous post.

To make it ever harder, based on @glwagner's comment I think we also need to discuss what exactly (regardless of format) zspacings should return to the user. Taking PartialGridBottom as an example, if we're to mimic the internals of the code, then it zspacings(ibm_grid::IBG) should return

1. zspacings(ibm_grid.underlying_grid) if z is above the immersed boundary
2. zspacing(ibm_grid) if the grid is at the immersed boundary (i.e. if the cell is partially wet, partially immersed)
3. zspacings(ibm_grid.underlying_grid) if z is fully below/inside the immersed boundary/solid

afaik this is exactly how the code internals work. However, since zspacings (and I guess also zspacing?) is a user-facing function, it's not clear if this is the best approach. For example, as a user, I wonder why we don't include the space occupied by the solid object when the cell is in the immersed boundary, but we do include it fully when the cell is completely inside/below the immersed boundary? In other words, why don't we do:

1. zspacings(ibm_grid.underlying_grid) if z is above the immersed boundary
2. zspacing(ibm_grid) if the grid is at the immersed boundary (i.e. if the cell is partially wet, partially immersed)
3. zero

?

I think this latter option is more consistent (i.e. a solid never counts as "wet space") and it makes stuff like offline integrals easier (since sum(quantity * Δz) / sum(Δz) will return the correct integral (ignoring x and y)). This is also more in line with the xgcm (also SGRID I think) conventions that I'm trying to implement in https://github.com/CliMA/Oceananigans.jl/pull/2652. I'm proposing the name active_zspacing() there:

https://github.com/CliMA/Oceananigans.jl/blob/b8b80598fb7250495f0d1925bfb406e15b3305b2/src/Grids/grid_utils.jl#L482-L496

(The main difference there is that the space outside the domain (i.e. halos) also doesn't count as wet volume)

[tomchor]

I'm tempted to change this PR to fix the issue at hand only for GridFittedBoundarys (since ultimately this part is an easy/fast fix), merge this PR to fix main asap, and then create an issue discussing how to move forward with the other issues.

@glwagner @navidcy what do you guys think?

[navidcy]

Why don't we deal with the full issue in this PR rather than merging and opening a new one. There is no rush, right?

[glwagner]

@tomchor FYI if you paste permalinks to code directly in the text, eg:

https://github.com/CliMA/Oceananigans.jl/blob/56fc69bcb08623424afe25e39ae616c29208095e/src/Operators/spacings_and_areas_and_volumes.jl#L247-L269

then github will generate a view into the code. This will save some time creating screenshots, etc. You have to paste the whole link for this to work, rather than appending a link onto some text (eg linking to the word "here")

[glwagner]

I'm tempted to change this PR to fix the issue at hand only for GridFittedBoundarys (since ultimately this part is an easy/fast fix), merge this PR to fix main asap, and then create an issue discussing how to move forward with the other issues.

@glwagner @navidcy what do you guys think?

I'm hesitant because I think we create more work for ourselves when we write code that may have to be revisited in the future. Better to do it right the first time I think...

[glwagner]

I think the right way forward is to implement something that works for all grids. Then we can implement the grid-specific versions --- which should be viewed as conveniences or optimizations rather than necessities --- as time allows. I think this is a better and more efficient approach then implementing convenience versions first, and figuring out the general version later. You might find that the convenience versions aren't really necessary because things are convenient enough...

[tomchor]

My suggestion of opening a separate issue comes from the fact that I don't have much time to work on this PR (as you guys could probably see by the month-delay in answering the last comments). And in the meantime there's a bug in the code that could be fixed rn with one extra commit. But fair enough. Let's solve this here.

I think the right way forward is to implement something that works for all grids. Then we can implement the grid-specific versions --- which should be viewed as conveniences or optimizations rather than necessities --- as time allows. I think this is a better and more efficient approach then implementing convenience versions first, and figuring out the general version later. You might find that the convenience versions aren't really necessary because things are convenient enough...

If I understand what you're proposing, the way forward would be to implement something that always returns the 3D field/array, since that's the only way that works for all grids no matter what. Correct?

[tomchor]

@tomchor FYI if you paste permalinks to code directly in the text, eg: then github will generate a view into the code. This will save some time creating screenshots, etc. You have to paste the whole link for this to work, rather than appending a link onto some text (eg linking to the word "here")

Thanks! I know about that feature actually. The pieces of code I pasted as hyperlink were ones I considered to be less important for the argument and didn't want it taking up space in the comment.

[navidcy]

But it's possible to write code that works for all grids by referencing the kernel functions directly (rather than trying to extract the spacings from the grid properties). That's how we get spacings inside the code, so if we use the same method in the user interface then all is good.

@glwagner, what do you mean here by "kernel functions"?

[glwagner]

But it's possible to write code that works for all grids by referencing the kernel functions directly (rather than trying to extract the spacings from the grid properties). That's how we get spacings inside the code, so if we use the same method in the user interface then all is good.

@glwagner, what do you mean here by "kernel functions"?

By "kernel function" I mean "a function meant to be used inside a kernel". In our convention these functions have the first four arguments i, j, k, grid. An example is Δxᶠᶜᶜ(i, j, k, grid), which gets the x-spacing at face, center, center.

Annoyingly they are all defined via metaprogramming, eg:

https://github.com/CliMA/Oceananigans.jl/blob/f4d7a94faeaa00cc48d0f92a26027a90d917ec62/src/Operators/spacings_and_areas_and_volumes.jl#L93

the 2D spacings though are defined

https://github.com/CliMA/Oceananigans.jl/blob/f4d7a94faeaa00cc48d0f92a26027a90d917ec62/src/Operators/spacings_and_areas_and_volumes.jl#L138