Adding reciprocal metric in the z-direction

src/Grids/grid_generation.jl

@@ -88,7 +88,9 @@ function generate_coordinate(FT, topo::AT, N, H, node_generator, coordinate_name
     C = OffsetArray(on_architecture(arch, C.parent), C.offsets...)
 
     if coordinate_name == :z
-        return L, StaticVerticalDiscretization(F, C, Δᶠ, Δᶜ)
+        Δ⁻¹ᶜ = OffsetArray(1 ./ parent(Δᶜ), -H)
+        Δ⁻¹ᶠ = OffsetArray(1 ./ parent(Δᶠ), -H - 1)
+        return L, StaticVerticalDiscretization(F, C, Δᶠ, Δᶜ, Δ⁻¹ᶠ, Δ⁻¹ᶜ)
     else
         return L, F, C, Δᶠ, Δᶜ
     end
@@ -107,7 +109,7 @@ function generate_coordinate(FT, topo::AT, N, H, node_interval::Tuple{<:Number,
     L = c₂ - c₁
 
     # Convert to get the correct type also when using single precision
-    Δᶠ = Δᶜ = Δ = L / N
+    Δᶠ = Δᶜ = Δ = FT(L / N)
 
     F₋ = c₁ - H * Δ
     F₊ = F₋ + total_extent(topo, H, Δ, L)
@@ -128,16 +130,18 @@ function generate_coordinate(FT, topo::AT, N, H, node_interval::Tuple{<:Number,
     C = OffsetArray(C, -H)
 
     if coordinate_name == :z
-        return FT(L), StaticVerticalDiscretization(F, C, FT(Δᶠ), FT(Δᶜ))
+        Δ⁻¹ᶜ = 1 / Δᶜ
+        Δ⁻¹ᶠ = 1 / Δᶠ
+        return FT(L), StaticVerticalDiscretization(F, C, Δᶠ, Δᶜ, Δ⁻¹ᶠ, Δ⁻¹ᶜ)
     else
-        return FT(L), F, C, FT(Δᶠ), FT(Δᶜ)
+        return FT(L), F, C, Δᶠ, Δᶜ
     end
 end
 
 # Flat domains
 function generate_coordinate(FT, ::Flat, N, H, c::Number, coordinate_name, arch)
     if coordinate_name == :z
-        return FT(1), StaticVerticalDiscretization(range(FT(c), FT(c), length=N), range(FT(c), FT(c), length=N), FT(1), FT(1))
+        return FT(1), StaticVerticalDiscretization(range(FT(c), FT(c), length=N), range(FT(c), FT(c), length=N), FT(1), FT(1), FT(1), FT(1))
     else
         return FT(1), range(FT(c), FT(c), length=N), range(FT(c), FT(c), length=N), FT(1), FT(1)
     end
@@ -148,8 +152,8 @@ end
 #     FT(1), c, c, FT(1), FT(1)
 function generate_coordinate(FT, ::Flat, N, H, ::Nothing, coordinate_name, arch)
     if coordinate_name == :z
-        return FT(1), StaticVerticalDiscretization(nothing, nothing, FT(1), FT(1))
+        return FT(1), StaticVerticalDiscretization(nothing, nothing, FT(1), FT(1), FT(1), FT(1))
     else
         return FT(1), nothing, nothing, FT(1), FT(1)
     end
-end
+end

src/Grids/vertical_discretization.jl

@@ -23,10 +23,12 @@ Fields
 - `Δᶠ::F`: Face-centered grid spacing.
 """
 struct StaticVerticalDiscretization{C, D, E, F} <: AbstractVerticalCoordinate
-    cᵃᵃᶠ :: C
-    cᵃᵃᶜ :: D
-    Δᵃᵃᶠ :: E
-    Δᵃᵃᶜ :: F
+      cᵃᵃᶠ :: C
+      cᵃᵃᶜ :: D
+      Δᵃᵃᶠ :: E
+      Δᵃᵃᶜ :: F
+    Δ⁻¹ᵃᵃᶠ :: E
+    Δ⁻¹ᵃᵃᶜ :: F
 end
 
 # Summaries
@@ -36,17 +38,23 @@ const AbstractStaticGrid  = AbstractUnderlyingGrid{<:Any, <:Any, <:Any, <:Any, <
 coordinate_summary(topo, z::StaticVerticalDiscretization, name) = coordinate_summary(topo, z.Δᵃᵃᶜ, name)
 
 struct MutableVerticalDiscretization{C, D, E, F, H, CC, FC, CF, FF} <: AbstractVerticalCoordinate
-    cᵃᵃᶠ :: C
-    cᵃᵃᶜ :: D
-    Δᵃᵃᶠ :: E
-    Δᵃᵃᶜ :: F
-      ηⁿ :: H
-    σᶜᶜⁿ :: CC
-    σᶠᶜⁿ :: FC
-    σᶜᶠⁿ :: CF
-    σᶠᶠⁿ :: FF
-    σᶜᶜ⁻ :: CC
-    ∂t_σ :: CC
+      cᵃᵃᶠ :: C
+      cᵃᵃᶜ :: D
+      Δᵃᵃᶠ :: E
+      Δᵃᵃᶜ :: F
+    Δ⁻¹ᵃᵃᶠ :: E
+    Δ⁻¹ᵃᵃᶜ :: F
+        ηⁿ :: H
+      σᶜᶜⁿ :: CC
+      σᶠᶜⁿ :: FC
+      σᶜᶠⁿ :: CF
+      σᶠᶠⁿ :: FF
+    σ⁻¹ᶜᶜⁿ :: CC
+    σ⁻¹ᶠᶜⁿ :: FC
+    σ⁻¹ᶜᶠⁿ :: CF
+    σ⁻¹ᶠᶠⁿ :: FF
+      σᶜᶜ⁻ :: CC
+      ∂t_σ :: CC
 end
 
 ####
@@ -68,7 +76,7 @@ or an `AbstractArray`. A `MutableVerticalDiscretization` defines a vertical coor
 following certain rules. Examples of `MutableVerticalDiscretization`s are free-surface following coordinates,
 or sigma coordinates.
 """
-MutableVerticalDiscretization(r) = MutableVerticalDiscretization(r, r, (nothing for i in 1:9)...)
+MutableVerticalDiscretization(r) = MutableVerticalDiscretization(r, r, (nothing for i in 1:15)...)
 
 coordinate_summary(::Bounded, z::RegularMutableVerticalDiscretization, name) =
     @sprintf("regularly spaced with mutable Δr=%s", prettysummary(z.Δᵃᵃᶜ))
@@ -137,41 +145,58 @@ end
 
 Adapt.adapt_structure(to, coord::StaticVerticalDiscretization) =
     StaticVerticalDiscretization(Adapt.adapt(to, coord.cᵃᵃᶠ),
-                                 Adapt.adapt(to, coord.cᵃᵃᶜ),
-                                 Adapt.adapt(to, coord.Δᵃᵃᶠ),
-                                 Adapt.adapt(to, coord.Δᵃᵃᶜ))
+                             Adapt.adapt(to, coord.cᵃᵃᶜ),
+                             Adapt.adapt(to, coord.Δᵃᵃᶠ),
+                             Adapt.adapt(to, coord.Δᵃᵃᶜ),
+                             Adapt.adapt(to, coord.Δ⁻¹ᵃᵃᶠ),
+                             Adapt.adapt(to, coord.Δ⁻¹ᵃᵃᶜ))
 
 on_architecture(arch, coord::StaticVerticalDiscretization) =
     StaticVerticalDiscretization(on_architecture(arch, coord.cᵃᵃᶠ),
-                                 on_architecture(arch, coord.cᵃᵃᶜ),
-                                 on_architecture(arch, coord.Δᵃᵃᶠ),
-                                 on_architecture(arch, coord.Δᵃᵃᶜ))
+                             on_architecture(arch, coord.cᵃᵃᶜ),
+                             on_architecture(arch, coord.Δᵃᵃᶠ),
+                             on_architecture(arch, coord.Δᵃᵃᶜ),
+                             on_architecture(arch, coord.Δ⁻¹ᵃᵃᶠ),
+                             on_architecture(arch, coord.Δ⁻¹ᵃᵃᶜ))
 
 Adapt.adapt_structure(to, coord::MutableVerticalDiscretization) =
     MutableVerticalDiscretization(Adapt.adapt(to, coord.cᵃᵃᶠ),
-                                  Adapt.adapt(to, coord.cᵃᵃᶜ),
-                                  Adapt.adapt(to, coord.Δᵃᵃᶠ),
-                                  Adapt.adapt(to, coord.Δᵃᵃᶜ),
-                                  Adapt.adapt(to, coord.ηⁿ),
-                                  Adapt.adapt(to, coord.σᶜᶜⁿ),
-                                  Adapt.adapt(to, coord.σᶠᶜⁿ),
-                                  Adapt.adapt(to, coord.σᶜᶠⁿ),
-                                  Adapt.adapt(to, coord.σᶠᶠⁿ),
-                                  Adapt.adapt(to, coord.σᶜᶜ⁻),
-                                  Adapt.adapt(to, coord.∂t_σ))
+                              Adapt.adapt(to, coord.cᵃᵃᶜ),
+                              Adapt.adapt(to, coord.Δᵃᵃᶠ),
+                              Adapt.adapt(to, coord.Δᵃᵃᶜ),
+                              Adapt.adapt(to, coord.Δ⁻¹ᵃᵃᶠ),
+                              Adapt.adapt(to, coord.Δ⁻¹ᵃᵃᶜ),
+                              Adapt.adapt(to, coord.ηⁿ),
+                              Adapt.adapt(to, coord.σᶜᶜⁿ),
+                              Adapt.adapt(to, coord.σᶠᶜⁿ),
+                              Adapt.adapt(to, coord.σᶜᶠⁿ),
+                              Adapt.adapt(to, coord.σᶠᶠⁿ),
+                              Adapt.adapt(to, coord.σ⁻¹ᶜᶜⁿ),
+                              Adapt.adapt(to, coord.σ⁻¹ᶠᶜⁿ),
+                              Adapt.adapt(to, coord.σ⁻¹ᶜᶠⁿ),
+                              Adapt.adapt(to, coord.σ⁻¹ᶠᶠⁿ),
+                              Adapt.adapt(to, coord.σᶜᶜ⁻),
+                              Adapt.adapt(to, coord.∂t_σ))
 
 on_architecture(arch, coord::MutableVerticalDiscretization) =
     MutableVerticalDiscretization(on_architecture(arch, coord.cᵃᵃᶠ),
-                                  on_architecture(arch, coord.cᵃᵃᶜ),
-                                  on_architecture(arch, coord.Δᵃᵃᶠ),
-                                  on_architecture(arch, coord.Δᵃᵃᶜ),
-                                  on_architecture(arch, coord.ηⁿ),
-                                  on_architecture(arch, coord.σᶜᶜⁿ),
-                                  on_architecture(arch, coord.σᶠᶜⁿ),
-                                  on_architecture(arch, coord.σᶜᶠⁿ),
-                                  on_architecture(arch, coord.σᶠᶠⁿ),
-                                  on_architecture(arch, coord.σᶜᶜ⁻),
-                                  on_architecture(arch, coord.∂t_σ))
+                              on_architecture(arch, coord.cᵃᵃᶜ),
+                              on_architecture(arch, coord.Δᵃᵃᶠ),
+                              on_architecture(arch, coord.Δᵃᵃᶜ),
+                              on_architecture(arch, coord.Δ⁻¹ᵃᵃᶠ),
+                              on_architecture(arch, coord.Δ⁻¹ᵃᵃᶜ),
+                              on_architecture(arch, coord.ηⁿ),
+                              on_architecture(arch, coord.σᶜᶜⁿ),
+                              on_architecture(arch, coord.σᶠᶜⁿ),
+                              on_architecture(arch, coord.σᶜᶠⁿ),
+                              on_architecture(arch, coord.σᶠᶠⁿ),
+                              on_architecture(arch, coord.σ⁻¹ᶜᶜⁿ),
+                              on_architecture(arch, coord.σ⁻¹ᶠᶜⁿ),
+                              on_architecture(arch, coord.σ⁻¹ᶜᶠⁿ),
+                              on_architecture(arch, coord.σ⁻¹ᶠᶠⁿ),
+                              on_architecture(arch, coord.σᶜᶜ⁻),
+                              on_architecture(arch, coord.∂t_σ))
+
 
 #####
 ##### Nodes and spacings (common to every grid)...

src/ImmersedBoundaries/immersed_grid_metrics.jl

@@ -25,6 +25,11 @@ import Oceananigans.Operators: intrinsic_vector, extrinsic_vector
 @inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid)
 @inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid)
 
+@inline Δr⁻¹ᵃᵃᶠ(i, j, k, ibg::IBG) = Δr⁻¹ᵃᵃᶠ(i, j, k, ibg.underlying_grid)
+@inline Δr⁻¹ᵃᵃᶜ(i, j, k, ibg::IBG) = Δr⁻¹ᵃᵃᶜ(i, j, k, ibg.underlying_grid)
+@inline Δz⁻¹ᵃᵃᶜ(i, j, k, ibg::IBG) = Δz⁻¹ᵃᵃᶜ(i, j, k, ibg.underlying_grid)
+@inline Δz⁻¹ᵃᵃᶠ(i, j, k, ibg::IBG) = Δz⁻¹ᵃᵃᶠ(i, j, k, ibg.underlying_grid)
+
 # 1D Horizontal spacings
 
 @inline Δxᶠᵃᵃ(i, j, k, ibg::RGIBG) = Δxᶠᵃᵃ(i, j, k, ibg.underlying_grid)

src/ImmersedBoundaries/mutable_immersed_grid.jl

@@ -29,8 +29,9 @@ const MutableGridOfSomeKind = Union{MutableImmersedGrid, AbstractMutableGrid}
 @inline column_depthᶠᶠᵃ(i, j, grid::MutableGridOfSomeKind) = column_depthᶠᶠᵃ(i, j, 1, grid, grid.z.ηⁿ)
 
 # Fallbacks
-@inline σⁿ(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = σⁿ(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
-@inline σ⁻(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = σ⁻(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline σⁿ(i,   j, k, ibg::IBG, ℓx, ℓy, ℓz) = σⁿ(i,   j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline σ⁻(i,   j, k, ibg::IBG, ℓx, ℓy, ℓz) = σ⁻(i,   j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
+@inline σ⁻¹ⁿ(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = σ⁻¹ⁿ(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz)
 
 @inline ∂t_σ(i, j, k, ibg::IBG) = ∂t_σ(i, j, k, ibg.underlying_grid)
 
@@ -43,6 +44,9 @@ for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ)
     zspacing = Symbol(:Δz, LX, LY, LZ)
     rspacing = Symbol(:Δr, LX, LY, LZ)
 
+    rcp_zspacing = Symbol(:Δz⁻¹, LX, LY, LZ)
+    rcp_rspacing = Symbol(:Δr⁻¹, LX, LY, LZ)
+
     ℓx = superscript_location(LX)
     ℓy = superscript_location(LY)
     ℓz = superscript_location(LZ)
@@ -54,5 +58,9 @@ for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ)
         @inline $zspacing(i, j, k, grid::IMRG)  = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
         @inline $zspacing(i, j, k, grid::IMLLG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
         @inline $zspacing(i, j, k, grid::IMOSG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+
+        @inline $rcp_zspacing(i, j, k, grid::IMRG)  = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+        @inline $rcp_zspacing(i, j, k, grid::IMLLG) = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+        @inline $rcp_zspacing(i, j, k, grid::IMOSG) = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
     end
 end

src/Models/HydrostaticFreeSurfaceModels/z_star_vertical_spacing.jl

@@ -27,7 +27,7 @@ function update_grid!(model, grid::MutableGridOfSomeKind, ::ZStar; parameters =
     η     = model.free_surface.η
 
     launch!(architecture(grid), grid, parameters, _update_grid_scaling!,
-            σᶜᶜⁿ, σᶠᶜⁿ, σᶜᶠⁿ, σᶠᶠⁿ, σᶜᶜ⁻, ηⁿ, grid, η)
+            σᶜᶜⁿ, σᶠᶜⁿ, σᶜᶠⁿ, σᶠᶠⁿ, σ⁻¹ᶜᶜⁿ, σ⁻¹ᶠᶜⁿ, σ⁻¹ᶜᶠⁿ, σ⁻¹ᶠᶠⁿ, σᶜᶜ⁻, ηⁿ, grid, η)
 
     # the barotropic velocities are retrieved from the free surface model for a
     # SplitExplicitFreeSurface and are calculated for other free surface models
@@ -41,7 +41,7 @@ function update_grid!(model, grid::MutableGridOfSomeKind, ::ZStar; parameters =
     return nothing
 end
 
-@kernel function _update_grid_scaling!(σᶜᶜⁿ, σᶠᶜⁿ, σᶜᶠⁿ, σᶠᶠⁿ, σᶜᶜ⁻, ηⁿ, grid, η)
+@kernel function _update_grid_scaling!(σᶜᶜⁿ, σᶠᶜⁿ, σᶜᶠⁿ, σᶠᶠⁿ, σ⁻¹ᶜᶜⁿ, σ⁻¹ᶠᶜⁿ, σ⁻¹ᶜᶠⁿ, σ⁻¹ᶠᶠⁿ, σᶜᶜ⁻, ηⁿ, grid, η)
     i, j = @index(Global, NTuple)
     k_top = size(grid, 3) + 1
 
@@ -70,6 +70,12 @@ end
         σᶜᶠⁿ[i, j, 1] = σᶜᶠ
         σᶠᶠⁿ[i, j, 1] = σᶠᶠ
 
+        # update current scaling reciprocal 
+        σ⁻¹ᶜᶜⁿ[i, j, 1] = 1 / σᶜᶜ
+        σ⁻¹ᶠᶜⁿ[i, j, 1] = 1 / σᶠᶜ
+        σ⁻¹ᶜᶠⁿ[i, j, 1] = 1 / σᶜᶠ
+        σ⁻¹ᶠᶠⁿ[i, j, 1] = 1 / σᶠᶠ
+
         # Update η in the grid
         ηⁿ[i, j, 1] = η[i, j, k_top]
     end
@@ -153,6 +159,6 @@ end
     i, j = @index(Global, NTuple)
 
     @unroll for k in -Hz+1:Nz+Hz
-        tracer[i, j, k] /= σⁿ(i, j, k, grid, Center(), Center(), Center())
+        tracer[i, j, k] *= σ⁻¹ⁿ(i, j, k, grid, Center(), Center(), Center())
     end
 end

src/Operators/Operators.jl

@@ -92,8 +92,8 @@ include("difference_operators.jl")
 include("interpolation_operators.jl")
 include("interpolation_utils.jl")
 
-include("reciprocal_metric_operators.jl")
 include("spacings_and_areas_and_volumes.jl")
+include("reciprocal_metric_operators.jl")
 include("products_between_fields_and_grid_metrics.jl")
 
 include("derivative_operators.jl")

src/Operators/reciprocal_metric_operators.jl

@@ -1,15 +1,46 @@
 # Reciprocal metrics operators; equal to the inverse of the metric operators
 for L1 in (:ᶜ, :ᶠ, :ᵃ), L2 in (:ᶜ, :ᶠ, :ᵃ), L3 in (:ᶜ, :ᶠ, :ᵃ)
-    for func in (:Δ, :A)
-        for dir in (:x, :y, :λ, :φ, :z, :r)
-            rcp_metric = Symbol(func, dir, :⁻¹, L1, L2, L3)
-            metric = Symbol(func, dir, L1, L2, L3)
-            @eval @inline $rcp_metric(i, j, k, grid) = 1 / $metric(i, j, k, grid)
-        end
+    for dir in (:x, :y)
+        rcp_spacing = Symbol(:Δ, dir, :⁻¹, L1, L2, L3)
+        spacing = Symbol(:Δ, dir, L1, L2, L3)
+        @eval @inline $rcp_spacing(i, j, k, grid) = 1 / $spacing(i, j, k, grid)
+    end
+
+    for dir in (:x, :y, :z)
+        rcp_area = Symbol(:A, dir, :⁻¹, L1, L2, L3)
+        area = Symbol(:A, dir, L1, L2, L3)
+        @eval @inline $rcp_area(i, j, k, grid) = 1 / $area(i, j, k, grid)
     end
 
-    rcp_volume = Symbol(:V⁻¹, L1, L2, L3)
-    volume = Symbol(:V, L1, L2, L3)
-    @eval @inline $rcp_volume(i, j, k, grid) = 1 / $volume(i, j, k, grid)
+    rcp_vol = Symbol(:V⁻¹, L1, L2, L3)
+    vol = Symbol(:V, L1, L2, L3)
+    @eval @inline $rcp_vol(i, j, k, grid) = 1 / $vol(i, j, k, grid)
 end
 
+@inline Δr⁻¹ᵃᵃᶜ(i, j, k, grid) = getspacing(k, grid.z.Δ⁻¹ᵃᵃᶜ)
+@inline Δr⁻¹ᵃᵃᶠ(i, j, k, grid) = getspacing(k, grid.z.Δ⁻¹ᵃᵃᶠ)
+
+@inline Δz⁻¹ᵃᵃᶜ(i, j, k, grid) = getspacing(k, grid.z.Δ⁻¹ᵃᵃᶜ)
+@inline Δz⁻¹ᵃᵃᶠ(i, j, k, grid) = getspacing(k, grid.z.Δ⁻¹ᵃᵃᶠ)
+
+for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ)
+    Δz⁻¹ᵃᵃˡ = Symbol(:Δz⁻¹, :ᵃ, :ᵃ, L1)
+    Δr⁻¹ᵃᵃˡ = Symbol(:Δr⁻¹, :ᵃ, :ᵃ, L1)
+    Δz⁻¹ˡᵃˡ = Symbol(:Δz⁻¹, L2, :ᵃ, L1)
+    Δr⁻¹ˡᵃˡ = Symbol(:Δr⁻¹, L2, :ᵃ, L1)
+    Δz⁻¹ᵃˡˡ = Symbol(:Δz⁻¹, :ᵃ, L2, L1)
+    Δr⁻¹ᵃˡˡ = Symbol(:Δr⁻¹, :ᵃ, L2, L1)
+
+    @eval @inline $Δz⁻¹ˡᵃˡ(i, j, k, grid) = $Δz⁻¹ᵃᵃˡ(i, j, k, grid)
+    @eval @inline $Δz⁻¹ᵃˡˡ(i, j, k, grid) = $Δz⁻¹ᵃᵃˡ(i, j, k, grid)
+    @eval @inline $Δr⁻¹ˡᵃˡ(i, j, k, grid) = $Δr⁻¹ᵃᵃˡ(i, j, k, grid)
+    @eval @inline $Δr⁻¹ᵃˡˡ(i, j, k, grid) = $Δr⁻¹ᵃᵃˡ(i, j, k, grid)
+
+    for L3 in (:ᶜ, :ᶠ)
+        Δz⁻¹ˡˡˡ = Symbol(:Δz⁻¹, L2, L3, L1)
+        Δr⁻¹ˡˡˡ = Symbol(:Δr⁻¹, L2, L3, L1)
+
+        @eval @inline $Δz⁻¹ˡˡˡ(i, j, k, grid) = $Δz⁻¹ˡᵃˡ(i, j, k, grid)
+        @eval @inline $Δr⁻¹ˡˡˡ(i, j, k, grid) = $Δr⁻¹ˡᵃˡ(i, j, k, grid)
+    end
+end

src/Operators/time_variable_grid_operators.jl

@@ -15,6 +15,11 @@ const AMG = AbstractMutableGrid
 @inline σⁿ(i, j, k, grid::AMG, ::C, ::F, ℓz) = @inbounds grid.z.σᶜᶠⁿ[i, j, 1]
 @inline σⁿ(i, j, k, grid::AMG, ::F, ::F, ℓz) = @inbounds grid.z.σᶠᶠⁿ[i, j, 1]
 
+@inline σ⁻¹ⁿ(i, j, k, grid::AMG, ::C, ::C, ℓz) = @inbounds grid.z.σ⁻¹ᶜᶜⁿ[i, j, 1]
+@inline σ⁻¹ⁿ(i, j, k, grid::AMG, ::F, ::C, ℓz) = @inbounds grid.z.σ⁻¹ᶠᶜⁿ[i, j, 1]
+@inline σ⁻¹ⁿ(i, j, k, grid::AMG, ::C, ::F, ℓz) = @inbounds grid.z.σ⁻¹ᶜᶠⁿ[i, j, 1]
+@inline σ⁻¹ⁿ(i, j, k, grid::AMG, ::F, ::F, ℓz) = @inbounds grid.z.σ⁻¹ᶠᶠⁿ[i, j, 1]
+
 # σ⁻ is needed only at centers
 @inline σ⁻(i, j, k, grid::AMG, ::C, ::C, ℓz) = @inbounds grid.z.σᶜᶜ⁻[i, j, 1]
 
@@ -35,6 +40,9 @@ for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ)
     zspacing = Symbol(:Δz, LX, LY, LZ)
     rspacing = Symbol(:Δr, LX, LY, LZ)
 
+    rcp_zspacing = Symbol(:Δz, :⁻¹, LX, LY, LZ)
+    rcp_rspacing = Symbol(:Δr, :⁻¹, LX, LY, LZ)
+
     ℓx = superscript_location(LX)
     ℓy = superscript_location(LY)
     ℓz = superscript_location(LZ)
@@ -43,6 +51,10 @@ for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ)
         @inline $zspacing(i, j, k, grid::MRG)  = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
         @inline $zspacing(i, j, k, grid::MLLG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
         @inline $zspacing(i, j, k, grid::MOSG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+
+        @inline $rcp_zspacing(i, j, k, grid::MRG)  = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+        @inline $rcp_zspacing(i, j, k, grid::MLLG) = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
+        @inline $rcp_zspacing(i, j, k, grid::MOSG) = $rcp_rspacing(i, j, k, grid) * σ⁻¹ⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz())
     end
 end