Completely overhaul grid utils + min_Δx/y/z -> minimum_spacing + move x/y/zspacing to Grids

examples/one_dimensional_diffusion.jl

@@ -88,8 +88,8 @@ current_figure() # hide
 # Next we set-up a `Simulation` that time-steps the model forward and manages output.
 
 ## Time-scale for diffusion across a grid cell
-using Oceananigans.Grids: min_Δz
-diffusion_time_scale = min_Δz(model.grid)^2 / model.closure.κ.T
+min_Δz = minimum_spacing(:z, (Face, Face, Face), model.grid)
+diffusion_time_scale = min_Δz^2 / model.closure.κ.T
 
 simulation = Simulation(model, Δt = 0.1 * diffusion_time_scale, stop_iteration = 1000)
 

examples/tilted_bottom_boundary_layer.jl

@@ -142,9 +142,8 @@ model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure,
 # conservatively, based on the smallest grid size of our domain and set-up a 
 
 using Oceananigans.Units
-using Oceananigans.Grids: min_Δz
 
-simulation = Simulation(model, Δt = 0.5 * min_Δz(grid) / V∞, stop_time = 2days)
+simulation = Simulation(model, Δt = 0.5 * minimum_spacing(:z, (Face, Face, Face), grid) / V∞, stop_time = 2days)
 
 # We use `TimeStepWizard` to adapt our time-step and print a progress message,
 

src/Advection/Advection.jl

@@ -64,5 +64,6 @@ include("topologically_conditional_interpolation.jl")
 include("momentum_advection_operators.jl")
 include("tracer_advection_operators.jl")
 include("positivity_preserving_tracer_advection_operators.jl")
+include("cell_advection_timescale.jl")
 
 end # module

src/Advection/cell_advection_timescale.jl

@@ -0,0 +1,18 @@
+using Oceananigans.Grids: topology, min_Δx, min_Δy, min_Δz
+
+function cell_advection_timescale(grid, velocities)
+    u, v, w = velocities
+    τ = KernelFunctionOperation{Center, Center, Center}(cell_advection_timescaleᶜᶜᶜ, grid, u, v, w)
+    return minimum(τ)
+end
+
+@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
+    Δx = Δxᶠᶜᶜ(i, j, k, grid)
+    Δy = Δyᶜᶠᶜ(i, j, k, grid)
+    Δz = Δzᶜᶜᶠ(i, j, k, grid)
+
+    return @inbounds min(Δx / abs(u[i, j, k]),
+                         Δy / abs(v[i, j, k]),
+                         Δz / abs(w[i, j, k])) 
+end
+

src/Diagnostics/cfl.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Utils: cell_advection_timescale
+using Oceananigans.Advection: cell_advection_timescale
 using Oceananigans.TurbulenceClosures: cell_diffusion_timescale
 
 """

src/Fields/abstract_field.jl

@@ -53,7 +53,7 @@ Returns the size of an `AbstractField{LX, LY, LZ}` located at `LX, LY, LZ`.
 This is a 3-tuple of integers corresponding to the number of interior nodes
 of `f` along `x, y, z`.
 """
-Base.size(f::AbstractField) = size(location(f), f.grid)
+Base.size(f::AbstractField) = size(f.grid, location(f))
 Base.length(f::AbstractField) = prod(size(f))
 Base.parent(f::AbstractField) = f
 
@@ -63,7 +63,7 @@ Base.parent(f::AbstractField) = f
 Returns a 3-tuple that gives the "total" size of a field including
 both interior points and halo points.
 """
-total_size(f::AbstractField) = total_size(location(f), f.grid)
+total_size(f::AbstractField) = total_size(f.grid, location(f))
 
 interior(f::AbstractField) = f
 

src/Fields/field.jl

@@ -35,7 +35,7 @@ end
 #####
 
 function validate_field_data(loc, data, grid, indices)
-    Fx, Fy, Fz = total_size(loc, grid, indices)
+    Fx, Fy, Fz = total_size(grid, loc, indices)
 
     if size(data) != (Fx, Fy, Fz)
         LX, LY, LZ = loc    
@@ -234,18 +234,20 @@ Return an `OffsetArray` of a `view` of `parent(data)` with `indices`.
 
 If `indices === (:, :, :)`, return an `OffsetArray` of `parent(data)`.
 """
-function offset_windowed_data(data, loc, grid, indices)
+function offset_windowed_data(data, Loc, grid, indices)
     halo = halo_size(grid)
-    topo = topology(grid)
+    topo = map(instantiate, topology(grid))
+    loc = map(instantiate, Loc)
 
     if indices isa typeof(default_indices(3))
         windowed_parent = parent(data)
     else
-        parent_indices = parent_index_range.(indices, loc, topo, halo)
+        parent_indices = map(parent_index_range, indices, loc, topo, halo)
         windowed_parent = view(parent(data), parent_indices...)
     end
 
     sz = size(grid)
+
     return offset_data(windowed_parent, loc, topo, sz, halo, indices)
 end
 
@@ -345,9 +347,7 @@ function boundary_conditions(f::Field)
 end
 
 immersed_boundary_condition(f::Field) = f.boundary_conditions.immersed
-
 data(field::Field) = field.data
-
 indices(obj, i=default_indices(3)) = i
 indices(f::Field, i=default_indices(3)) = f.indices
 indices(a::SubArray, i=default_indices(ndims(a))) = a.indices
@@ -357,15 +357,20 @@ indices(a::OffsetArray, i=default_indices(ndims(a))) = indices(parent(a), i)
 interior_view_indices(field_indices, interior_indices) = Colon()
 interior_view_indices(::Colon,       interior_indices) = interior_indices
 
+instantiate(T::Type) = T()
+instantiate(t) = t
+
 function interior(a::OffsetArray,
-                  loc::Tuple,
-                  topo::Tuple,
-                  size::NTuple{N, Int},
-                  halo_size::NTuple{N, Int},
+                  Loc::Tuple,
+                  Topo::Tuple,
+                  sz::NTuple{N, Int},
+                  halo_sz::NTuple{N, Int},
                   ind::Tuple=default_indices(3)) where N
 
-    i_interior = interior_parent_indices.(loc, topo, size, halo_size)
-    i_view = interior_view_indices.(ind, i_interior)
+    loc = map(instantiate, Loc)
+    topo = map(instantiate, Topo)
+    i_interior = map(interior_parent_indices, loc, topo, sz, halo_sz)
+    i_view = map(interior_view_indices, ind, i_interior)
     return view(parent(a), i_view...)
 end
 
@@ -399,11 +404,8 @@ Base.fill!(f::Field, val) = fill!(parent(f), val)
 Base.parent(f::Field) = parent(f.data)
 Adapt.adapt_structure(to, f::Field) = Adapt.adapt(to, f.data)
 
-length_indices(N, i::Colon) = N
-length_indices(N, i::UnitRange) = length(i)
-
-total_size(f::Field) = length_indices.(total_size(location(f), f.grid), f.indices)
-Base.size(f::Field)  = length_indices.(      size(location(f), f.grid), f.indices)
+total_size(f::Field) = total_size(f.grid, location(f), f.indices)
+Base.size(f::Field)  = size(f.grid, location(f), f.indices)
 
 ==(f::Field, a) = interior(f) == a
 ==(a, f::Field) = a == interior(f)
@@ -723,3 +725,4 @@ function fill_halo_regions!(field::Field, args...; kwargs...)
 
     return nothing
 end
+

src/Grids/Grids.jl

@@ -10,7 +10,7 @@ export AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid
 export LatitudeLongitudeGrid, XRegLatLonGrid, YRegLatLonGrid, ZRegLatLonGrid
 export OrthogonalSphericalShellGrid, ConformalCubedSphereGrid
 export node, xnode, ynode, znode, xnodes, ynodes, znodes, nodes
-export xspacings, yspacings, zspacings
+export xspacings, yspacings, zspacings, xspacing, yspacing, zspacing, minimum_spacing
 export offset_data, new_data
 export on_architecture
 

src/Grids/grid_generation.jl

@@ -12,20 +12,21 @@ get_coord_face(coord::Nothing, i) = 1
 get_coord_face(coord::Function, i) = coord(i)
 get_coord_face(coord::AbstractVector, i) = CUDA.@allowscalar coord[i]
 
-lower_exterior_Δcoordᶠ(topology, Fi, Hcoord) = [Fi[end - Hcoord + i] - Fi[end - Hcoord + i - 1] for i = 1:Hcoord]
-lower_exterior_Δcoordᶠ(::Type{<:BoundedTopology}, Fi, Hcoord) = [Fi[2]  - Fi[1] for i = 1:Hcoord]
+const AT = AbstractTopology
+lower_exterior_Δcoordᶠ(::AT, Fi, Hcoord) = [Fi[end - Hcoord + i] - Fi[end - Hcoord + i - 1] for i = 1:Hcoord]
+lower_exterior_Δcoordᶠ(::BoundedTopology, Fi, Hcoord) = [Fi[2]  - Fi[1] for i = 1:Hcoord]
 
-upper_exterior_Δcoordᶠ(topology, Fi, Hcoord) = [Fi[i + 1] - Fi[i] for i = 1:Hcoord]
-upper_exterior_Δcoordᶠ(::Type{<:BoundedTopology}, Fi, Hcoord) = [Fi[end]   - Fi[end - 1] for i = 1:Hcoord]
+upper_exterior_Δcoordᶠ(::AT, Fi, Hcoord) = [Fi[i + 1] - Fi[i] for i = 1:Hcoord]
+upper_exterior_Δcoordᶠ(::BoundedTopology, Fi, Hcoord) = [Fi[end]   - Fi[end - 1] for i = 1:Hcoord]
 
-upper_interior_F(topology, coord, Δ)               = coord - Δ
-upper_interior_F(::Type{<:BoundedTopology}, coord) = coord
+upper_interior_F(::AT, coord, Δ)               = coord - Δ
+upper_interior_F(::BoundedTopology, coord) = coord
 
-total_interior_length(topology, N)                  = N
-total_interior_length(::Type{<:BoundedTopology}, N) = N + 1
+total_interior_length(::AT, N)                  = N
+total_interior_length(::BoundedTopology, N) = N + 1
 
 # generate a stretched coordinate passing the explicit coord faces as vector of functionL
-function generate_coordinate(FT, topology, N, H, coord, arch)
+function generate_coordinate(FT, topo::AT, N, H, coord, arch)
 
     # Ensure correct type for F and derived quantities
     interiorF = zeros(FT, N+1)
@@ -37,8 +38,8 @@ function generate_coordinate(FT, topology, N, H, coord, arch)
     L = interiorF[N+1] - interiorF[1]
 
     # Build halo regions
-    Δᶠ₋ = lower_exterior_Δcoordᶠ(topology, interiorF, H)
-    Δᶠ₊ = reverse(upper_exterior_Δcoordᶠ(topology, interiorF, H))
+    Δᶠ₋ = lower_exterior_Δcoordᶠ(topo, interiorF, H)
+    Δᶠ₊ = reverse(upper_exterior_Δcoordᶠ(topo, interiorF, H))
 
     c¹, cᴺ⁺¹ = interiorF[1], interiorF[N+1]
 
@@ -48,12 +49,12 @@ function generate_coordinate(FT, topology, N, H, coord, arch)
     F = vcat(F₋, interiorF, F₊)
 
     # Build cell centers, cell center spacings, and cell interface spacings
-    TC = total_length(Center, topology, N, H)
+    TC = total_length(Center(), topo, N, H)
      C = [ (F[i + 1] + F[i]) / 2 for i = 1:TC ]
     Δᶠ = [  C[i] - C[i - 1]      for i = 2:TC ]
 
     # Trim face locations for periodic domains
-    TF = total_length(Face, topology, N, H)
+    TF = total_length(Face(), topo, N, H)
     F  = F[1:TF]
 
     Δᶜ = [F[i + 1] - F[i] for i = 1:TF-1]
@@ -77,7 +78,7 @@ function generate_coordinate(FT, topology, N, H, coord, arch)
 end
 
 # generate a regular coordinate passing the domain extent (2-tuple) and number of points
-function generate_coordinate(FT, topology, N, H, coord::Tuple{<:Number, <:Number}, arch)
+function generate_coordinate(FT, topo::AT, N, H, coord::Tuple{<:Number, <:Number}, arch)
 
     @assert length(coord) == 2
 
@@ -89,13 +90,13 @@ function generate_coordinate(FT, topology, N, H, coord::Tuple{<:Number, <:Number
     Δᶠ = Δᶜ = Δ = L / N
 
     F₋ = c₁ - H * Δ
-    F₊ = F₋ + total_extent(topology, H, Δ, L)
+    F₊ = F₋ + total_extent(topo, H, Δ, L)
 
     C₋ = F₋ + Δ / 2
     C₊ = C₋ + L + Δ * (2H - 1)
 
-    TF = total_length(Face,   topology, N, H)
-    TC = total_length(Center, topology, N, H)
+    TF = total_length(Face(),   topo, N, H)
+    TC = total_length(Center(), topo, N, H)
 
     F = range(FT(F₋), FT(F₊), length = TF)
     C = range(FT(C₋), FT(C₊), length = TC)
@@ -107,6 +108,6 @@ function generate_coordinate(FT, topology, N, H, coord::Tuple{<:Number, <:Number
 end
 
 # Flat domains
-function generate_coordinate(FT, ::Type{Flat}, N, H, coord::Tuple{<:Number, <:Number}, arch)
+function generate_coordinate(FT, ::Flat, N, H, coord::Tuple{<:Number, <:Number}, arch)
     return FT(1), range(1, 1, length=N), range(1, 1, length=N), FT(1), FT(1)
 end