Should vertical spacings `Δzᵃᵃᶜ` and `Δzᵃᵃᶠ` depend on horizontal location?

[glwagner]

They don't right now, even though right now the lateral areas Ax and Ay do depend on horizontal location:

https://github.com/CliMA/Oceananigans.jl/blob/24e766481cebbc8f61099b386623d175218acedb/src/Operators/spacings_and_areas_and_volumes.jl#L151

Though our underlying grids are "extruded" in the vertical (put another way, we make the thin shell approximation for spherical shell grids), vertical spacings may in principle need to depend on horizontal location to accurately represent bathymetry.

Another way to see this is that the vertical spacings for GridFittedBottom (implemented in #2023) are

https://github.com/CliMA/Oceananigans.jl/blob/8eadc493e04002066448323573edcfde046b9c30/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl#L56-L62

thus conditioning on the function is_immersed. But is_immersed is evaluated at (Center, Center, Center). But more general operators such as Δzᶠᶜᶜ may be needed for correct bathymetry representation.

Note that we might need this note only to compute depths correctly, but also for correct diagnostics.

cc @sandreza @jm-c

[sandreza]

For the use case that I was thinking we would need to distinguish, in particular, Δzᶠᶜᶜ from Δzᶜᶠᶜ

[francispoulin]

In dealing with a stretched grid, don't we need these to depenend on the horizontal grid?

[christophernhill]

@glwagner this

https://github.com/CliMA/Oceananigans.jl/blob/24e766481cebbc8f61099b386623d175218acedb/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_lateral_areas.jl#L9

looks like it needs some Δzᶠᶜᶜ etc... to me? Something like

ifelse( ( is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary) | is_immersed(i-1, j, k, ibg.grid, ibg.immersed_boundary) ), 
                                              zero(eltype(ibg.grid)), 
                                              Δzᵃᵃᶠ(i, j, k, ibg.grid)) 

should work while things are step, I think?

Lopping and/or sigma would be more fancy and might be better precomputed?

[glwagner]

I think you're right @christophernhill, we probably need this for correct implicit free surface too...

[glwagner]

In dealing with a stretched grid, don't we need these to depenend on the horizontal grid?

Not in general; without bathymetry we would only need this for three-dimensional curvilinearity. If we have only horizontal curvilinearity (as arises in a thin approximation to the spherical shell), then vertical spacing are independent of horizontal location.

But as noted in this issue, vertical spacing do depend on horizontal location with an immersed boundary.

[christophernhill]

I think you're right @christophernhill, we probably need this for correct implicit free surface too...

@glwagner (plus @sandreza and @jm-c ) I have put some first steps here

https://github.com/CliMA/Oceananigans.jl/tree/cnh/stagger-immeresed-bc-operators

its based on the global-lat-lon branch so it should have latest immersed_boundaries e.g.

$ git diff global-lat-lon..

diff --git a/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl b/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl
index 7b955541..e9e16770 100644
--- a/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl
+++ b/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl
@@ -57,6 +57,18 @@ const GFBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:GridFit
                                              zero(eltype(ibg.grid)),
                                              Δzᵃᵃᶜ(i, j, k, ibg.grid))

+@inline Δzᶠᶜᶜ(i, j, k, ibg::GFBIBG) = ifelse(  is_immersed(i-1, j  , k, ibg.grid, ibg.immersed_boundary)
+                                             | is_immersed(i  , j  , k, ibg.grid, ibg.immersed_boundary)
+                                             ,
+                                             zero(eltype(ibg.grid)),
+                                             Δzᵃᵃᶜ(i, j, k, ibg.grid))
+
+@inline Δzᶜᶠᶜ(i, j, k, ibg::GFBIBG) = ifelse(  is_immersed(i  , j-1, k, ibg.grid, ibg.immersed_boundary)
+                                             | is_immersed(i  , j  , k, ibg.grid, ibg.immersed_boundary)
+                                             ,
+                                             zero(eltype(ibg.grid)),
+                                             Δzᵃᵃᶜ(i, j, k, ibg.grid))
+
 @inline Δzᵃᵃᶠ(i, j, k, ibg::GFBIBG) = ifelse(is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary),
                                              zero(eltype(ibg.grid)),
                                              Δzᵃᵃᶠ(i, j, k, ibg.grid))
diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl
index ecc1e434..677818c0 100644
--- a/src/ImmersedBoundaries/immersed_grid_metrics.jl
+++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl
@@ -11,6 +11,8 @@ for metric in (

                :Δzᵃᵃᶠ,
                :Δzᵃᵃᶜ,
+               :Δzᶠᶜᶜ,
+               :Δzᶜᶠᶜ,

                :Azᶜᶜᵃ,
                :Azᶜᶠᵃ,
diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_lateral_areas.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_lateral_areas.jl
index d2eb51c6..85be87b5 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_lateral_areas.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_lateral_areas.jl
@@ -6,8 +6,8 @@
         ∫ᶻ_A.yᶜᶠᶜ[i, j, 1] = 0

         @unroll for k in 1:grid.Nz
-            ∫ᶻ_A.xᶠᶜᶜ[i, j, 1] += Δyᶠᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid)
-            ∫ᶻ_A.yᶜᶠᶜ[i, j, 1] += Δxᶜᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid)
+            ∫ᶻ_A.xᶠᶜᶜ[i, j, 1] += Δyᶠᶜᵃ(i, j, k, grid) * Δzᶠᶜᶜ(i, j, k, grid)
+            ∫ᶻ_A.yᶜᶠᶜ[i, j, 1] += Δxᶜᶠᵃ(i, j, k, grid) * Δzᶜᶠᶜ(i, j, k, grid)
         end
     end
 end
diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl
index 8bd00859..00a445eb 100644
--- a/src/Operators/Operators.jl
+++ b/src/Operators/Operators.jl
@@ -1,7 +1,7 @@
 module Operators

 export Δx, Δy, Δz, Ax, Ay, Az, volume
-export ΔzF, ΔzC, Δzᵃᵃᶜ, Δzᵃᵃᶠ
+export ΔzF, ΔzC, Δzᵃᵃᶜ, Δzᵃᵃᶠ, Δzᶠᶜᶜ, Δzᶜᶠᶜ
 export Δxᶜᶜᵃ, Δxᶠᶜᵃ, Δxᶠᶠᵃ, Δxᶜᶠᵃ
 export Δyᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶠᶠᵃ, Δyᶜᶠᵃ
 export Axᵃᵃᶜ, Axᵃᵃᶠ, Axᶜᶜᶜ, Axᶠᶜᶜ, Axᶠᶠᶜ, Axᶠᶜᶠ, Axᶜᶠᶜ, Axᶜᶜᶠ
diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl
index 8b7aaeb7..9988b320 100644
--- a/src/Operators/spacings_and_areas_and_volumes.jl
+++ b/src/Operators/spacings_and_areas_and_volumes.jl
@@ -46,6 +46,9 @@ The operators in this file fall into three categories:
 @inline Δzᵃᵃᶜ(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
 @inline Δzᵃᵃᶜ(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶜ[k]

+@inline Δzᶠᶜᶜ(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
+@inline Δzᶜᶠᶜ(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶜ[k]
+
 #####
 ##### "Spacings" in Flat directions for rectilinear grids.
 ##### Here we dispatch all spacings to `one`. This abuse of notation
diff --git a/test/runtests.jl b/test/runtests.jl
index 246192f6..d9d520f0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -131,6 +131,7 @@ group = get(ENV, "TEST_GROUP", :all) |> Symbol
             include("test_vertical_vorticity_field.jl")
             include("test_implicit_free_surface_solver.jl")
             include("test_hydrostatic_free_surface_immersed_boundaries_apply_surf_bc.jl")
+            include("test_hydrostatic_free_surface_immersed_boundaries_vertical_integrals.jl")
         end
     end

diff --git a/test/test_hydrostatic_free_surface_immersed_boundaries_vertical_integrals.jl b/test/test_hydrostatic_free_surface_immersed_boundaries_vertical_integrals.jl
new file mode 100644
index 00000000..5241929b
--- /dev/null
+++ b/test/test_hydrostatic_free_surface_immersed_boundaries_vertical_integrals.jl
@@ -0,0 +1,56 @@
+using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
+using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization
+
+@testset "Immersed boundaries with hydrostatic free surface models" begin
+    @info "Testing immersed boundaries vertical integrals"
+
+    for arch in archs
+        Nx = 5
+        Ny = 5
+
+        # A spherical domain
+        underlying_grid =
+        RegularRectilinearGrid(size=(Nx, Ny, 3), extent=(Nx, Ny, 3), topology=(Periodic,Periodic,Bounded))
+
+        B = [-3. for i=1:Nx, j=1:Ny ]
+        B[2:Nx-1,2:Ny-1] .= [-2. for i=2:Nx-1, j=2:Ny-1 ]
+        B[3:Nx-2,3:Ny-2] .= [-1. for i=3:Nx-2, j=3:Ny-2 ]
+        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(B))
+
+        free_surface = ImplicitFreeSurface(gravitational_acceleration=0.1)
+
+        model = HydrostaticFreeSurfaceModel(grid = grid,
+                                           architecture = arch,
+                                           #free_surface = ExplicitFreeSurface(),
+                                           #free_surface = ImplicitFreeSurface(maximum_iterations=10),
+                                           free_surface = ImplicitFreeSurface(),
+                                           momentum_advection = nothing,
+                                           tracer_advection = WENO5(),
+                                           coriolis = nothing,
+                                           buoyancy = nothing,
+                                           tracers = nothing,
+                                           closure = nothing)
+
+        x_ref = [0.0  0.0  0.0  0.0  0.0  0.0  0.0
+                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
+                 0.0  3.0  2.0  2.0  2.0  2.0  0.0
+                 0.0  3.0  2.0  1.0  1.0  2.0  0.0
+                 0.0  3.0  2.0  2.0  2.0  2.0  0.0
+                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
+                 0.0  0.0  0.0  0.0  0.0  0.0  0.0]'
+
+        y_ref = [0.0  0.0  0.0  0.0  0.0  0.0  0.0
+                 0.0  3.0  3.0  3.0  3.0  3.0  0.0
+                 0.0  3.0  2.0  2.0  2.0  3.0  0.0
+                 0.0  3.0  2.0  1.0  2.0  3.0  0.0
+                 0.0  3.0  2.0  1.0  2.0  3.0  0.0
+                 0.0  3.0  2.0  2.0  2.0  3.0  0.0
+                 0.0  0.0  0.0  0.0  0.0  0.0  0.0]'
+
+        fs=model.free_surface
+        xok=parent(fs.implicit_step_solver.vertically_integrated_lateral_areas.xᶠᶜᶜ.data[:,:,1])-x_ref == zeros(7,7)
+        yok=parent(fs.implicit_step_solver.vertically_integrated_lateral_areas.yᶜᶠᶜ.data[:,:,1])-y_ref == zeros(7,7)
+        @test (xok & yok)
+    end
+end
+

[glwagner]

Do we want to define these metrics next to the other ones in the Operators module? We should use the existing function solid_node I think.

Edited to say that I think this is a good start and we can update all the operators in a future PR. It's probably best done after #2050.