zspacing does not adjust for immersed boundaries

[tomchor]

I noticed that zspacing is blind to whether or not the grid has an immersed boundary. For example:

julia> grid = RectilinearGrid(size=(3, 3, 3), extent=(1, 1, 1));

julia> bottom(x, y) = - 1/2
bottom (generic function with 1 method)

julia> grid = ImmersedBoundaryGrid(grid, PartialCellBottom(bottom))
3×3×3 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── immersed_boundary: PartialCellBottom(mean(zb)=-0.5, min(zb)=-0.5, max(zb)=-0.5, ϵ=0.2)
├── underlying_grid: 3×3×3 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.333333
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.333333
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=0.333333

julia> @compute dz = Field(zspacings(grid, Center(), Center(), Center()));

julia> interior(dz)
3×3×3 view(::Array{Float64, 3}, 4:6, 4:6, 4:6) with eltype Float64:
[:, :, 1] =
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333

[:, :, 2] =
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333

[:, :, 3] =
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333

I would expect dz to be zero in at k=1 (inside the immersed boundary), 1/6 at k=2 and 1/3 at k=3. @glwagner suggested that this is a bug.

[glwagner]

which underlying function is being called for zspacings?

[tomchor]

Seems to be Δz:

Oceananigans.jl/src/AbstractOperations/grid_metrics.jl

Lines 156 to 160 in 4e38f12

	function zspacings(grid, ℓx, ℓy, ℓz)
	LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
	Δz_op = KernelFunctionOperation{LX, LY, LZ}(Δz, grid, ℓx, ℓy, ℓz)
	return Δz_op
	end

[simone-silvestri]

Hmmm, I guess it's a bug introduced with zstar, what if you try with

function rspacings(grid, ℓx, ℓy, ℓz)
    LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
    Δx_op = KernelFunctionOperation{LX, LY, LZ}(Δr, grid, ℓx, ℓy, ℓz)
    return Δr_op
end

[glwagner]

I think this bug was introduced by the transition to z*, since

Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl

Lines 154 to 205 in 4e38f12

	@inline function Δrᶜᶜᶜ(i, j, k, ibg::PCBIBG)
	underlying_grid = ibg.underlying_grid
	ib = ibg.immersed_boundary
	# Get node at face above and defining nodes on c,c,f
	z = rnode(i, j, k+1, underlying_grid, c, c, f)
	# Get bottom z-coordinate and fractional Δz parameter
	zb = @inbounds ib.bottom_height[i, j, 1]
	# Are we in a bottom cell?
	at_the_bottom = bottom_cell(i, j, k, ibg)
	full_Δz = Δrᶜᶜᶜ(i, j, k, ibg.underlying_grid)
	partial_Δz = z - zb
	return ifelse(at_the_bottom, partial_Δz, full_Δz)
	end
	@inline function Δrᶜᶜᶠ(i, j, k, ibg::PCBIBG)
	just_above_bottom = bottom_cell(i, j, k-1, ibg)
	zc = rnode(i, j, k, ibg.underlying_grid, c, c, c)
	zf = rnode(i, j, k, ibg.underlying_grid, c, c, f)
	full_Δz = Δrᶜᶜᶠ(i, j, k, ibg.underlying_grid)
	partial_Δz = zc - zf + Δrᶜᶜᶜ(i, j, k-1, ibg) / 2
	Δz = ifelse(just_above_bottom, partial_Δz, full_Δz)
	return Δz
	end
	@inline Δrᶠᶜᶜ(i, j, k, ibg::PCBIBG) = min(Δrᶜᶜᶜ(i-1, j, k, ibg), Δrᶜᶜᶜ(i, j, k, ibg))
	@inline Δrᶜᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δrᶜᶜᶜ(i, j-1, k, ibg), Δrᶜᶜᶜ(i, j, k, ibg))
	@inline Δrᶠᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δrᶠᶜᶜ(i, j-1, k, ibg), Δrᶠᶜᶜ(i, j, k, ibg))
	@inline Δrᶠᶜᶠ(i, j, k, ibg::PCBIBG) = min(Δrᶜᶜᶠ(i-1, j, k, ibg), Δrᶜᶜᶠ(i, j, k, ibg))
	@inline Δrᶜᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δrᶜᶜᶠ(i, j-1, k, ibg), Δrᶜᶜᶠ(i, j, k, ibg))
	@inline Δrᶠᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δrᶠᶜᶠ(i, j-1, k, ibg), Δrᶠᶜᶠ(i, j, k, ibg))
	# Make sure Δz works for horizontally-Flat topologies.
	# (There's no point in using z-Flat with PartialCellBottom).
	XFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Flat, <:Any, <:Any, <:Any, <:PartialCellBottom}
	YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:PartialCellBottom}
	@inline Δrᶠᶜᶜ(i, j, k, ibg::XFlatPCBIBG) = Δrᶜᶜᶜ(i, j, k, ibg)
	@inline Δrᶠᶜᶠ(i, j, k, ibg::XFlatPCBIBG) = Δrᶜᶜᶠ(i, j, k, ibg)
	@inline Δrᶜᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δrᶜᶜᶜ(i, j, k, ibg)
	@inline Δrᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δrᶜᶜᶠ(i, j, k, ibg)
	@inline Δrᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δrᶜᶠᶜ(i, j, k, ibg)
	@inline Δrᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δrᶠᶜᶜ(i, j, k, ibg)

[glwagner]

I also think that it's warranted to clean up the names in partial_cell_bottom which mix the concepts of dr and dz and are a bit confusing.

We need to define dz correctly right? Including the time-fluctuating metrics?

[tomchor]

Here's a doubt that I have. In the first example, would we also expect a different result for zspacing() if we had used a GridFittedBottom? i.e., would we expect dz to be zero in regions inside the immersed boundary?

[simone-silvestri]

No, the dz is not affected by a grid fitted bottom (because the bottom is fitted), so the dz will not be zero inside the immersed boundary

[glwagner]

the abstraction for the currently implemented immersed boundaries has two components:

• definition of immersed_cell
• possibly, changes to grid metrics, right now just for PartialCellBottom

if we do have cut cells in the future, we'll need some additional things I think.

[tomchor]

Hmmm, I guess it's a bug introduced with zstar, what if you try with

function rspacings(grid, ℓx, ℓy, ℓz)
LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
Δx_op = KernelFunctionOperation{LX, LY, LZ}(Δr, grid, ℓx, ℓy, ℓz)
return Δr_op
end

This seems to work. Although your version has a small typo so I'm pasting a working version here for reference:

import Oceananigans: rspacings
using Oceananigans.Operators: Δr
function rspacings(grid, ℓx, ℓy, ℓz)
    LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
    Δr_op = KernelFunctionOperation{LX, LY, LZ}(Δr, grid, ℓx, ℓy, ℓz)
    return Δr_op
end

Output:

3×3×3 view(::Array{Float64, 3}, 4:6, 4:6, 4:6) with eltype Float64:
[:, :, 1] =
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333

[:, :, 2] =
 0.166667  0.166667  0.166667
 0.166667  0.166667  0.166667
 0.166667  0.166667  0.166667

[:, :, 3] =
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333
 0.333333  0.333333  0.333333

[glwagner]

I don't think we want rspacings though. That gives the reference spacings but not the actual spacings

[tomchor]

I'm happy to create a PR that fixes this and adds tests, but can I get some guidance on the best to do so? I'm a bit out of the loop on the current system of defining vertical coordinates and their spacings, so I'm unsure on what the best approach is here.

[glwagner]

Ok, here is a summary of how things work with a time-dependent grid:

• zspacings returns the internode distances on the current grid
• rspacings returns the internode distances on a "resting" or "reference" grid

For a static grid, zspacings and rspacings are identical.