(0.94.0) Overhaul spacings functions to return KernelFunctionOperations

[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]

tests are a good idea!

[navidcy]

I pushed some tests. Given that they pass then this lgtm! :)

[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.

[glwagner]

This is wrong in general, are we sure we want this?

Except for the "GridFitted" immersed boundary methods, immersed boundaries in general may change the metrics of the underlying grid. I think we want this user-facing method to return the correct spacings, ie when using PartialCellBottom, and the in-development CutCellBottom.

[glwagner]

We can proceed by dispatching this method only when valid, ie for GridFittedBoundary and GridFittedBottom. We'll need a different implementation for PartialCellBottom.

I think the right approach is to use KernelFunctionOperation for xspacing. Then every spacings method can simply collect the output of xspacing. This is generically valid and would presumably reduce the amount of code we have to write by a bit.

[navidcy]

oh good point!

so, for PartialCellBottom only the zspacings is not correct, right? the horizontal spacings don't change.

[glwagner]

Correct

[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.

[tomchor]

Agree that KernelFunctionOperations would reduce the amount of code. It would also return an AbstractOperation instead of the current options (which are either a float or an array), which I think has extra computational costs, no? Also I think we'd lose a bit of the simplicity on the user side.

For example the following line of code (taken from this doc example would become way less readable since zspacings() would now return an AbstractOperation, which would need to be first converted into a 3D Field, computed, and transformed into a 1D Array before plotting:

lines!(ax, zspacings(grid, Center()), znodes(grid, Center()))

Finally, following the same logic, would xnodes etc. also become KernelFunctionOperations?

[tomchor]

Just to clarify, I'm happy to make the changes once we agree on the way forward. I just wanna make sure we're all on the same page and also get an idea of the size of the PR.

[glwagner]

Hmm right I was a little confused. zspacing is actually the kernel function itself (not an operation or array).

For zspacings we do either have to (1) represent it with an operation or (2) compute the metrics. There's no way around this in general: with partial cells, we compute the metrics on the fly (this computation seems to have no cost on the GPU). For cut cells, computing the metrics on the fly is required in fact, because precomputing all the metrics (3D arrays) would greatly increase our memory footprint, and also probably incur a big cost since we are memory-bound.

We could implement optimized zspacings for non-immersed grids or GridFittedBottom. But is there a point? I don't think the costs would be noticeable compared to, for example, the computational cost of displaying a plot.

It also seems simpler to compute zspacings from zspacing, since more generic code might be implemented in that case. Note that for partial cells though, zspacings has to return a 3D array. For many other grids its just a 1D array.

I think it's a good idea to think carefully about the implementation before doing anything hasty.

[glwagner]

which would need to be first converted into a 3D Field, computed, and transformed into a 1D Array before plotting:

We have 1D fields. Also operations like "tranforming a field into an array" are not expensive. For example, we might do something like this O(100) times in a single time step.

I suggest we put performance to the side right now and think first about a clean, simple implementation. We can optimize performance later if needed; these functions neede primarily for plotting aren't performance critical anyways

[glwagner]

For example the following line of code (taken from this doc example would become way less readable

I'm not suggesting to change the user API. I'm suggesting to change the source code.

[navidcy]

@glwagner raised a valid point! I'm retracting my approval until we fixed this.

[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

Oceananigans.jl/src/Grids/rectilinear_grid.jl

Lines 477 to 478 in 56fc69b

	@inline xspacings(grid::RectilinearGrid, ℓx::Center; with_halos=false) = with_halos ? grid.Δxᶜᵃᵃ : view(grid.Δxᶜᵃᵃ, interior_indices(ℓx, topology(grid, 1)(), size(grid, 1)))
	@inline xspacings(grid::XRegRectilinearGrid, ℓx::Center; with_halos=false) = grid.Δxᶜᵃᵃ

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:

Oceananigans.jl/src/Operators/spacings_and_areas_and_volumes.jl

Lines 104 to 112 in 56fc69b

	@inline Δxᶠᵃᵃ(i, j, k, grid::RG) = @inbounds grid.Δxᶠᵃᵃ[i]
	@inline Δxᶜᵃᵃ(i, j, k, grid::RG) = @inbounds grid.Δxᶜᵃᵃ[i]
	@inline Δyᵃᶠᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶠᵃ[j]
	@inline Δyᵃᶜᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j]
	@inline Δxᶠᵃᵃ(i, j, k, grid::RGX) = @inbounds grid.Δxᶠᵃᵃ
	@inline Δxᶜᵃᵃ(i, j, k, grid::RGX) = @inbounds grid.Δxᶜᵃᵃ
	@inline Δyᵃᶠᵃ(i, j, k, grid::RGY) = @inbounds grid.Δyᵃᶠᵃ
	@inline Δyᵃᶜᵃ(i, j, k, grid::RGY) = @inbounds grid.Δyᵃᶜᵃ

---

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 #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 #2652. I'm proposing the name active_zspacing() there:

Oceananigans.jl/src/Grids/grid_utils.jl

Lines 482 to 496 in b8b8059

	for dir in (:x, :y, :z)
	active_spacing = Symbol(:active_, dir, :spacing)
	spacing = Symbol(dir, :spacing)
	boundary_spacing = Symbol(active_spacing, :_at_boundary)
	@eval begin
	function $active_spacing(i, j, k, grid, ℓx, ℓy, ℓz)
	Δξ = ifelse(inactive_node(i, j, k, grid, ℓx, ℓy, ℓz),
	0,
	ifelse(boundary_node(i, j, k, grid, ℓx, ℓy, ℓz),
	$boundary_spacing(i, j, k, grid, ℓx, ℓy, ℓz),
	$spacing(i, j, k, grid, ℓx, ℓy, ℓz)))
	return Δξ
	end
	end
	end

(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:

Oceananigans.jl/src/Operators/spacings_and_areas_and_volumes.jl

Lines 247 to 269 in 56fc69b

	for LX in (:Center, :Face)
	for LY in (:Center, :Face)
	for LZ in (:Center, :Face)
	LXe = @eval $LX
	LYe = @eval $LY
	LZe = @eval $LZ
	volume_function = Symbol(:V, location_code(LXe, LYe, LZe))
	@eval begin
	@inline volume(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $volume_function(i, j, k, grid)
	end
	for op in (:Δ, :A), dir in (:x, :y, :z)
	func = Symbol(op, dir)
	metric = Symbol(op, dir, location_code(LXe, LYe, LZe))
	@eval begin
	@inline $func(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $metric(i, j, k, grid)
	end
	end
	end
	end
	end

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:

Oceananigans.jl/src/Operators/spacings_and_areas_and_volumes.jl

Line 93 in f4d7a94

@inline $x_spacing_3D(i, j, k, grid) = $x_spacing_2D(i, j, k, grid)

the 2D spacings though are defined

Oceananigans.jl/src/Operators/spacings_and_areas_and_volumes.jl

Line 138 in f4d7a94

@inline Δxᶠᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j])

[ali-ramadhan]

I fixed the tests, doctests, and plotting examples in the docs! I ended up defining the spacing functions needed and used by other modules, and it should be easy to add more in the future as needed.

All pipelines should pass 🤞 And this PR should be ready for review again.

[tomchor]

I can't approve my own PR so I'm giving verbal approval here. Thanks for reviving this PR @ali-ramadhan! Also extra kudos for removing so many trailing white spaces :)

[ali-ramadhan]

@tomchor I'll approve this PR in your place haha

@glwagner Just let me know if you want to have a final look at this PR or if I should merge.

[glwagner]

where did these go?

[glwagner]

Great work @ali-ramadhan. I'm a little confused where some spacing functions went. But if tests pass I'm assuming everything is in order. It would only be nice to know so I know where to look when I need to find them.

[glwagner]

As for white space, while perhaps it is laudable to remove them it does make reviewing the PR a lot more difficult. I suggest putting in a separate PR to remove them throughout the code base so that we don't have to remove them PR-by-PR in the future (at least, it will reduce the need for that).

[navidcy]

Remember to register v0.94.0 when this is merged ;)

[glwagner]

I found the spacings definitions! They were just moved around in the file.

[glwagner]

@tomchor I'll leave it to you to delete the branch if you like

[ali-ramadhan]

Thanks for merging! And yeah the definitions should have just moved around.

[ali-ramadhan]

As for white space, while perhaps it is laudable to remove them it does make reviewing the PR a lot more difficult.

Yeah sorry about that. I have set my editor to trim whitespace. It can be disabled when looking at the changes

but dedicating a PR to it is a good idea.