static_column_depth interface

ext/OceananigansMakieExt.jl

@@ -38,7 +38,7 @@ end
 convert_arguments(pl::Type{<:AbstractPlot}, f::Field) =
     convert_arguments(pl, convert_field_argument(f)...)
 
-function convert_arguments(pl::Type{<:AbstractPlot}, fop::AbstractOperation)
+function convert_arguments(pl::Type{<:AbstractPlot}, op::AbstractOperation)
     f = Field(op)
     compute!(f)
     return convert_arguments(pl, f)

src/Grids/Grids.jl

@@ -18,6 +18,7 @@ export xnodes, ynodes, znodes, λnodes, φnodes
 export spacings
 export xspacings, yspacings, zspacings, xspacing, yspacing, zspacing
 export minimum_xspacing, minimum_yspacing, minimum_zspacing
+export static_column_depthᶜᶜᵃ, static_column_depthᶠᶜᵃ, static_column_depthᶜᶠᵃ, static_column_depthᶠᶠᵃ
 export offset_data, new_data
 export on_architecture
 

src/Grids/grid_utils.jl

@@ -314,6 +314,15 @@ coordinate_summary(topo, Δ::Union{AbstractVector, AbstractMatrix}, name) =
              name, prettysummary(minimum(parent(Δ))),
              name, prettysummary(maximum(parent(Δ))))
 
+#####
+##### Static column depth
+#####
+
+@inline static_column_depthᶜᶜᵃ(i, j, grid) = grid.Lz
+@inline static_column_depthᶜᶠᵃ(i, j, grid) = grid.Lz
+@inline static_column_depthᶠᶜᵃ(i, j, grid) = grid.Lz
+@inline static_column_depthᶠᶠᵃ(i, j, grid) = grid.Lz
+
 #####
 ##### Spherical geometry
 #####

src/ImmersedBoundaries/ImmersedBoundaries.jl

@@ -17,15 +17,18 @@ import Base: show, summary
 import Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z,
                            x_domain, y_domain, z_domain
 
-import Oceananigans.Grids: architecture, on_architecture, with_halo, inflate_halo_size_one_dimension,
+import Oceananigans.Grids: architecture, with_halo, inflate_halo_size_one_dimension,
                            xnode, ynode, znode, λnode, φnode, node,
                            ξnode, ηnode, rnode,
                            ξname, ηname, rname, node_names,
                            xnodes, ynodes, znodes, λnodes, φnodes, nodes,
                            ξnodes, ηnodes, rnodes,
+                           static_column_depthᶜᶜᵃ, static_column_depthᶠᶜᵃ, static_column_depthᶜᶠᵃ, static_column_depthᶠᶠᵃ,
                            inactive_cell
 
 
+import Oceananigans.Architectures: on_architecture
+
 import Oceananigans.Fields: fractional_x_index, fractional_y_index, fractional_z_index
 
 """
@@ -92,10 +95,6 @@ with_halo(halo, ibg::ImmersedBoundaryGrid) =
 inflate_halo_size_one_dimension(req_H, old_H, _, ::IBG)            = max(req_H + 1, old_H)
 inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, ::IBG) = 0
 
-# Defining the bottom
-@inline z_bottom(i, j, grid) = znode(i, j, 1, grid, c, c, f)
-@inline z_bottom(i, j, ibg::IBG) = error("The function `bottom` has not been defined for $(summary(ibg))!")
-
 function Base.summary(grid::ImmersedBoundaryGrid)
     FT = eltype(grid)
     TX, TY, TZ = topology(grid)

src/ImmersedBoundaries/abstract_grid_fitted_boundary.jl

@@ -18,21 +18,3 @@ const AGFB = AbstractGridFittedBoundary
     @inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, <:Any, Flat, Flat},  ib::AGFB) = _immersed_cell(i, 1, 1, grid, ib)
     @inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, Flat, Flat},   ib::AGFB) = _immersed_cell(1, 1, 1, grid, ib)
 end
-
-function clamp_bottom_height!(bottom_field, grid)
-    launch!(architecture(grid), grid, :xy, _clamp_bottom_height!, bottom_field, grid)
-    return nothing
-end
-
-const c = Center()
-const f = Face()
-
-@kernel function _clamp_bottom_height!(z, grid)
-    i, j = @index(Global, NTuple)
-    Nz = size(grid, 3)
-    zmin = znode(i, j, 1,    grid, c, c, f)
-    zmax = znode(i, j, Nz+1, grid, c, c, f)
-    @inbounds z[i, j, 1] = clamp(z[i, j, 1], zmin, zmax)
-end
-
-

src/ImmersedBoundaries/grid_fitted_bottom.jl

@@ -15,17 +15,17 @@ abstract type AbstractGridFittedBottom{H} <: AbstractGridFittedBoundary end
 
 # To enable comparison with PartialCellBottom in the limiting case that
 # fractional cell height is 1.0.
-struct CenterImmersedCondition end
 struct InterfaceImmersedCondition end
-
-Base.summary(::CenterImmersedCondition) = "CenterImmersedCondition"
-Base.summary(::InterfaceImmersedCondition) = "InterfaceImmersedCondition"
+struct CenterImmersedCondition end
 
 struct GridFittedBottom{H, I} <: AbstractGridFittedBottom{H}
     bottom_height :: H
     immersed_condition :: I
 end
 
+Base.summary(::CenterImmersedCondition) = "CenterImmersedCondition"
+Base.summary(::InterfaceImmersedCondition) = "InterfaceImmersedCondition"
+
 const GFBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:GridFittedBottom}
 
 """
@@ -36,6 +36,7 @@ Return a bottom immersed boundary.
 Keyword Arguments
 =================
 
+
 * `bottom_height`: an array or function that gives the height of the
                    bottom in absolute ``z`` coordinates.
 
@@ -72,50 +73,77 @@ Base.summary(ib::GridFittedBottom{<:Function}) = @sprintf("GridFittedBottom(%s)"
 function Base.show(io::IO, ib::GridFittedBottom)
     print(io, summary(ib), '\n')
     print(io, "├── bottom_height: ", prettysummary(ib.bottom_height), '\n')
-    print(io, "└── immersed_condition: ", summary(ib.immersed_condition))
 end
 
-@inline z_bottom(i, j, ibg::GFBIBG) = @inbounds ibg.immersed_boundary.bottom_height[i, j, 1]
+on_architecture(arch, ib::GridFittedBottom) = GridFittedBottom(on_architecture(arch, ib.bottom_height), ib.immersed_condition)
+
+function on_architecture(arch, ib::GridFittedBottom{<:Field})
+    architecture(ib.bottom_height) == arch && return ib
+    arch_grid = on_architecture(arch, ib.bottom_height.grid)
+    new_bottom_height = Field{Center, Center, Nothing}(arch_grid)
+    set!(new_bottom_height, ib.bottom_height)
+    fill_halo_regions!(new_bottom_height)
+    return GridFittedBottom(new_bottom_height, ib.immersed_condition)
+end
+
+Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom_height), adapt(to, ib.immersed_condition))
 
 """
     ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
 
 Return a grid with `GridFittedBottom` immersed boundary (`ib`).
 
-Computes `ib.bottom_height` and wraps it in a Field.
+Computes `ib.bottom_height` and wraps it in a Field. `ib.bottom_height` is the z-coordinate of top-most interface
+of the last ``immersed`` cell in the column.
 """
 function ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
     bottom_field = Field{Center, Center, Nothing}(grid)
     set!(bottom_field, ib.bottom_height)
-    @apply_regionally clamp_bottom_height!(bottom_field, grid)
+    @apply_regionally compute_numerical_bottom_height!(bottom_field, grid, ib)
     fill_halo_regions!(bottom_field)
-    new_ib = GridFittedBottom(bottom_field, ib.immersed_condition)
+    new_ib = GridFittedBottom(bottom_field)
     TX, TY, TZ = topology(grid)
     return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib)
 end
 
-@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition})
-    z = znode(i, j, k+1, underlying_grid, c, c, f)
-    h = @inbounds ib.bottom_height[i, j, 1]
-    return z ≤ h
+compute_numerical_bottom_height!(bottom_field, grid, ib) = 
+    launch!(architecture(grid), grid, :xy, _compute_numerical_bottom_height!, bottom_field, grid, ib)
+
+@kernel function _compute_numerical_bottom_height!(bottom_field, grid, ib::GridFittedBottom)
+    i, j = @index(Global, NTuple)
+    zb = @inbounds bottom_field[i, j, 1]
+    @inbounds bottom_field[i, j, 1] = znode(i, j, 1, grid, c, c, f)
+    condition = ib.immersed_condition
+    for k in 1:grid.Nz
+        z⁺ = znode(i, j, k+1, grid, c, c, f)
+        z  = znode(i, j, k,   grid, c, c, c)
+        bottom_cell = ifelse(condition isa CenterImmersedCondition, z ≤ zb, z⁺ ≤ zb)
+        @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, z⁺, zb)
+    end
 end
 
-@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:CenterImmersedCondition})
-    z = znode(i, j, k, underlying_grid, c, c, c)
-    h = @inbounds ib.bottom_height[i, j, 1]
-    return z ≤ h
+@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom)
+    z  = znode(i, j, k, underlying_grid, c, c, c)
+    zb = @inbounds ib.bottom_height[i, j, 1]
+    return z ≤ zb
 end
 
-on_architecture(arch, ib::GridFittedBottom) = GridFittedBottom(ib.bottom_height, ib.immersed_condition)
+#####
+##### Static column depth
+#####
 
-function on_architecture(arch, ib::GridFittedBottom{<:Field})
-    architecture(ib.bottom_height) == arch && return ib
-    arch_grid = on_architecture(arch, ib.bottom_height.grid)
-    new_bottom_height = Field{Center, Center, Nothing}(arch_grid)
-    set!(new_bottom_height, ib.bottom_height)
-    fill_halo_regions!(new_bottom_height)
-    return GridFittedBottom(new_bottom_height, ib.immersed_condition)
-end
+const AGFBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:AbstractGridFittedBottom}
+
+@inline static_column_depthᶜᶜᵃ(i, j, ibg::AGFBIBG) = @inbounds znode(i, j, ibg.Nz+1, ibg, c, c, f) - ibg.immersed_boundary.bottom_height[i, j, 1] 
+@inline static_column_depthᶜᶠᵃ(i, j, ibg::AGFBIBG) = min(static_column_depthᶜᶜᵃ(i, j-1, ibg), static_column_depthᶜᶜᵃ(i, j, ibg))
+@inline static_column_depthᶠᶜᵃ(i, j, ibg::AGFBIBG) = min(static_column_depthᶜᶜᵃ(i-1, j, ibg), static_column_depthᶜᶜᵃ(i, j, ibg))
+@inline static_column_depthᶠᶠᵃ(i, j, ibg::AGFBIBG) = min(static_column_depthᶠᶜᵃ(i, j-1, ibg), static_column_depthᶠᶜᵃ(i, j, ibg))
+
+# Make sure column_height works for horizontally-Flat topologies.
+XFlatAGFIBG = ImmersedBoundaryGrid{<:Any, <:Flat, <:Any, <:Any, <:Any, <:AbstractGridFittedBottom}
+YFlatAGFIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:AbstractGridFittedBottom}
 
-Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom_height),
-                                                                             ib.immersed_condition)
+@inline static_column_depthᶠᶜᵃ(i, j, ibg::XFlatAGFIBG) = static_column_depthᶜᶜᵃ(i, j, ibg)
+@inline static_column_depthᶜᶠᵃ(i, j, ibg::YFlatAGFIBG) = static_column_depthᶜᶜᵃ(i, j, ibg)
+@inline static_column_depthᶠᶠᵃ(i, j, ibg::XFlatAGFIBG) = static_column_depthᶜᶠᵃ(i, j, ibg)
+@inline static_column_depthᶠᶠᵃ(i, j, ibg::YFlatAGFIBG) = static_column_depthᶠᶜᵃ(i, j, ibg)

src/ImmersedBoundaries/partial_cell_bottom.jl

@@ -63,13 +63,40 @@ end
 function ImmersedBoundaryGrid(grid, ib::PartialCellBottom)
     bottom_field = Field{Center, Center, Nothing}(grid)
     set!(bottom_field, ib.bottom_height)
-    @apply_regionally clamp_bottom_height!(bottom_field, grid)
+    @apply_regionally compute_numerical_bottom_height!(bottom_field, grid, ib)
     fill_halo_regions!(bottom_field)
     new_ib = PartialCellBottom(bottom_field, ib.minimum_fractional_cell_height)
     TX, TY, TZ = topology(grid)
     return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib)
 end
 
+@kernel function _compute_numerical_bottom_height!(bottom_field, grid, ib::PartialCellBottom)
+    i, j = @index(Global, NTuple)
+
+    # Save analytical bottom height
+    zb = @inbounds bottom_field[i, j, 1]
+
+    # Cap bottom height at Lz and at rnode(i, j, grid.Nz+1, grid, c, c, f)
+    domain_bottom = znode(i, j, 1, grid, c, c, f)
+    domain_top    = znode(i, j, grid.Nz+1, grid, c, c, f)
+    @inbounds bottom_field[i, j, 1] = clamp(zb, domain_bottom, domain_top)
+
+    ϵ  = ib.minimum_fractional_cell_height
+    for k in 1:grid.Nz
+        z⁻ = znode(i, j, k,   grid, c, c, f)
+        z⁺ = znode(i, j, k+1, grid, c, c, f)
+        Δz = Δzᶜᶜᶜ(i, j, k, grid)
+        bottom_cell = (z⁻ ≤ zb) & (z⁺ ≥ zb)
+
+        # If the size of the bottom cell is less than ϵ Δz, use 
+        # a simple `GridFittedBottom` approach where the bottom 
+        # height is the top interface of cell k.
+        capped_zb   = ifelse(zb < z⁻ + Δz * (1 - ϵ), zb, z⁺)
+
+        @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, capped_zb, bottom_field[i, j, 1])
+    end
+end
+
 function on_architecture(arch, ib::PartialCellBottom{<:Field})
     architecture(ib.bottom_height) == arch && return ib
     arch_grid = on_architecture(arch, ib.bottom_height.grid)
@@ -106,13 +133,12 @@ Criterion is zb ≥ z - ϵ Δz
 
 """
 @inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom)
-    # Face node below current cell
-    z  = znode(i, j, k, underlying_grid, c, c, f)
-    zb = @inbounds ib.bottom_height[i, j, 1]
+    z⁻ = znode(i, j, k, underlying_grid, c, c, f)
     ϵ  = ib.minimum_fractional_cell_height
-    # z + Δz is equal to the face above the current cell
     Δz = Δzᶜᶜᶜ(i, j, k, underlying_grid)
-    return (z + Δz * (1 - ϵ)) ≤ zb
+    z★ = z⁻ + Δz * (1 - ϵ)
+    zb = @inbounds ib.bottom_height[i, j, 1]
+    return z★ < zb
 end
 
 @inline function bottom_cell(i, j, k, ibg::PCBIBG)
@@ -129,15 +155,14 @@ end
     # Get node at face above and defining nodes on c,c,f
     z = znode(i, j, k+1, underlying_grid, c, c, f)
 
-    # Get bottom height and fractional Δz parameter
-    h = @inbounds ib.bottom_height[i, j, 1]
-    ϵ = ibg.immersed_boundary.minimum_fractional_cell_height
+    # 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 = Δzᶜᶜᶜ(i, j, k, ibg.underlying_grid)
-    partial_Δz = max(ϵ * full_Δz, z - h)
+    full_Δz    = Δzᶜᶜᶜ(i, j, k, ibg.underlying_grid)
+    partial_Δz = z - zb
 
     return ifelse(at_the_bottom, partial_Δz, full_Δz)
 end
@@ -175,6 +200,3 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial
 @inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg)
 @inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg)
 @inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg)
-
-@inline z_bottom(i, j, ibg::PCBIBG) = @inbounds ibg.immersed_boundary.bottom_height[i, j, 1]
-

src/Models/HydrostaticFreeSurfaceModels/distributed_split_explicit_free_surface.jl

@@ -10,27 +10,13 @@ function SplitExplicitAuxiliaryFields(grid::DistributedGrid)
     Gᵁ = Field((Face,   Center, Nothing), grid)
     Gⱽ = Field((Center, Face,   Nothing), grid)
 
-    Hᶠᶜ = Field((Face,   Center, Nothing), grid)
-    Hᶜᶠ = Field((Center, Face,   Nothing), grid)
-
-    calculate_column_height!(Hᶠᶜ, (Face, Center, Center))
-    calculate_column_height!(Hᶜᶠ, (Center, Face, Center))
-
-    fill_halo_regions!((Hᶠᶜ, Hᶜᶠ))
-
     # In a non-parallel grid we calculate only the interior
     kernel_size    = augmented_kernel_size(grid)
     kernel_offsets = augmented_kernel_offsets(grid)
 
     kernel_parameters = KernelParameters(kernel_size, kernel_offsets)
 
-    return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ, kernel_parameters)
-end
-
-"""Integrate z at locations `location` and set! `height`` with the result"""
-@inline function calculate_column_height!(height, location)
-    dz = GridMetricOperation(location, Δz, height.grid)
-    return sum!(height, dz)
+    return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, kernel_parameters)
 end
 
 @inline function augmented_kernel_size(grid::DistributedGrid)

src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl

@@ -209,10 +209,6 @@ Base.@kwdef struct SplitExplicitAuxiliaryFields{𝒞ℱ, ℱ𝒞, 𝒦}
     Gᵁ :: ℱ𝒞
     "Vertically-integrated slow barotropic forcing function for `V` (`ReducedField` over ``z``)"
     Gⱽ :: 𝒞ℱ
-    "Depth at `(Face, Center)` (`ReducedField` over ``z``)"
-    Hᶠᶜ :: ℱ𝒞
-    "Depth at `(Center, Face)` (`ReducedField` over ``z``)"
-    Hᶜᶠ :: 𝒞ℱ
     "kernel size for barotropic time stepping"
     kernel_parameters :: 𝒦
 end
@@ -227,20 +223,9 @@ function SplitExplicitAuxiliaryFields(grid::AbstractGrid)
     Gᵁ = Field((Face,   Center, Nothing), grid)
     Gⱽ = Field((Center, Face,   Nothing), grid)
 
-    Hᶠᶜ = Field((Face,   Center, Nothing), grid)
-    Hᶜᶠ = Field((Center, Face,   Nothing), grid)
-
-    dz = GridMetricOperation((Face, Center, Center), Δz, grid)
-    sum!(Hᶠᶜ, dz)
-
-    dz = GridMetricOperation((Center, Face, Center), Δz, grid)
-    sum!(Hᶜᶠ, dz)
-
-    fill_halo_regions!((Hᶠᶜ, Hᶜᶠ))
-
     kernel_parameters = :xy
 
-    return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ, kernel_parameters)
+    return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, kernel_parameters)
 end
 
 """