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

Manifest.toml

@@ -318,9 +318,9 @@ version = "4.15.0"
 
 [[deps.LLVMExtra_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg", "TOML"]
-git-tree-sha1 = "771bfe376249626d3ca12bcd58ba243d3f961576"
+git-tree-sha1 = "7718cf44439c676bc0ec66a87099f41015a522d6"
 uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
-version = "0.0.16+0"
+version = "0.0.16+2"
 
 [[deps.LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]

docs/src/appendix/library.md

@@ -37,7 +37,7 @@ Modules = [Oceananigans.BoundaryConditions]
 Private = false
 ```
 
-## BuoyancyModels
+## Buoyancy models
 
 ```@autodocs
 Modules = [Oceananigans.BuoyancyModels]
@@ -127,7 +127,6 @@ Private = false
 ```@autodocs
 Modules = [Oceananigans.Models.HydrostaticFreeSurfaceModels]
 Private = false
-]
 ```
 
 ### Shallow-water models
@@ -149,23 +148,20 @@ Private = false
 ```@autodocs
 Modules = [Oceananigans.Operators]
 Private = false
-]
 ```
 
 ## Output readers
 
 ```@autodocs
 Modules = [Oceananigans.OutputReaders]
 Private = false
-]
 ```
 
 ## Output writers
 
 ```@autodocs
 Modules = [Oceananigans.OutputWriters]
 Private = false
-]
 ```
 
 ## Simulations

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_zspacing(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_zspacing(grid) / V∞, stop_time = 2days)
 
 # We use `TimeStepWizard` to adapt our time-step and print a progress message,
 

src/AbstractOperations/at.jl

@@ -77,16 +77,19 @@ end
 interpolate_index(::Colon, ::Colon, args...)       = Colon()
 interpolate_index(::Colon, b::UnitRange, args...)  = b
 
+instantiate(T::Type) = T()
+instantiate(t) = t
+
 function interpolate_index(a::UnitRange, ::Colon, loc, new_loc)  
-    a = corrected_index(a, loc, new_loc)
+    a = corrected_index(a, instantiate(loc), instantiate(new_loc))
 
     # Abstract operations that require an interpolation of a sliced fields are not supported!
     first(a) > last(a) && throw(ArgumentError("Cannot interpolate a sliced field from $loc to $(new_loc)!"))
     return a
 end
 
 function interpolate_index(a::UnitRange, b::UnitRange, loc, new_loc)   
-    a = corrected_index(a, loc, new_loc)
+    a = corrected_index(a, instantiate(loc), instantiate(new_loc))
 
     # Abstract operations that require an interpolation of a sliced fields are not supported!
     first(a) > last(a) && throw(ArgumentError("Cannot interpolate a sliced field from $loc to $(new_loc)!"))
@@ -100,7 +103,7 @@ end
 
 # Windowed fields interpolated from `Center`s to `Face`s lose the first index.
 # Viceverse, windowed fields interpolated from `Face`s to `Center`s lose the last index
-corrected_index(a, ::Type{Face},   ::Type{Face})   = UnitRange(first(a),   last(a))
-corrected_index(a, ::Type{Center}, ::Type{Center}) = UnitRange(first(a),   last(a))
-corrected_index(a, ::Type{Face},   ::Type{Center}) = UnitRange(first(a),   last(a)-1)
-corrected_index(a, ::Type{Center}, ::Type{Face})   = UnitRange(first(a)+1, last(a))
+corrected_index(a, ::Face,   ::Face)   = UnitRange(first(a),   last(a)  )
+corrected_index(a, ::Center, ::Center) = UnitRange(first(a),   last(a)  )
+corrected_index(a, ::Face,   ::Center) = UnitRange(first(a),   last(a)-1)
+corrected_index(a, ::Center, ::Face)   = UnitRange(first(a)+1, last(a)  )

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,19 @@
+using Oceananigans.AbstractOperations: KernelFunctionOperation
+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/CubedSpheres/conformal_cubed_sphere_grid.jl

@@ -229,10 +229,11 @@ const OSSG = OrthogonalSphericalShellGrid
 ##### Grid utils
 #####
 
-Base.size(grid::ConformalCubedSphereGrid)      = (size(grid.faces[1])..., length(grid.faces))
-Base.size(loc, grid::ConformalCubedSphereGrid) = size(loc, grid.faces[1])
-Base.size(grid::ConformalCubedSphereGrid, i)   = size(grid)[i]
-halo_size(ccsg::ConformalCubedSphereGrid)      = halo_size(first(ccsg.faces)) # hack
+Base.size(grid::ConformalCubedSphereGrid)             = (size(grid.faces[1])..., length(grid.faces))
+Base.size(grid::ConformalCubedSphereGrid, loc::Tuple) = size(grid.faces[1], loc)
+Base.size(grid::ConformalCubedSphereGrid, i::Int)     = size(grid)[i]
+
+halo_size(ccsg::ConformalCubedSphereGrid) = halo_size(first(ccsg.faces)) # hack
 
 Base.eltype(grid::ConformalCubedSphereGrid{FT}) where FT = FT
 

src/CubedSpheres/cubed_sphere_utils.jl

@@ -11,95 +11,98 @@ using Oceananigans.Grids:
 ##### Viewing halos
 #####
 
+instantiate(T::Type) = T()
+instantiate(t) = t
+
 west_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
-                      view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
-                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
-                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))
+    include_corners ? view(f.data, left_halo_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx, f.grid.Hx), :, :) :
+                      view(f.data, left_halo_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx, f.grid.Hx),
+                                   interior_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny),
+                                   interior_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz))
 
 east_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
-                      view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
-                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
-                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))
+    include_corners ? view(f.data, right_halo_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx, f.grid.Hx), :, :) :
+                      view(f.data, right_halo_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx, f.grid.Hx),
+                                   interior_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny),
+                                   interior_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz))
 
 south_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, :, left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
-                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
-                                   left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
-                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))
+    include_corners ? view(f.data, :, left_halo_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny, f.grid.Hy), :) :
+                      view(f.data, interior_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx),
+                                   left_halo_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny, f.grid.Hy),
+                                   interior_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz))
 
 north_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, :, right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
-                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
-                                   right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
-                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))
+    include_corners ? view(f.data, :, right_halo_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny, f.grid.Hy), :) :
+                      view(f.data, interior_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx),
+                                   right_halo_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny, f.grid.Hy),
+                                   interior_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz))
 
 bottom_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, :, :, left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
-                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
-                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
-                                   left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))
+    include_corners ? view(f.data, :, :, left_halo_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz, f.grid.Hz)) :
+                      view(f.data, interior_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx),
+                                   interior_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny),
+                                   left_halo_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz, f.grid.Hz))
 
 top_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
-    include_corners ? view(f.data, :, :, right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
-                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
-                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
-                                   right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))
+    include_corners ? view(f.data, :, :, right_halo_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz, f.grid.Hz)) :
+                      view(f.data, interior_indices(instantiate(LX), instantiate(topology(f, 1)), f.grid.Nx),
+                                   interior_indices(instantiate(LY), instantiate(topology(f, 2)), f.grid.Ny),
+                                   right_halo_indices(instantiate(LZ), instantiate(topology(f, 3)), f.grid.Nz, f.grid.Hz))
 
 underlying_west_halo(f, grid, location, topo=topology(grid, 1)) =
-    view(f.parent, underlying_left_halo_indices(location, topo, grid.Nx, grid.Hx), :, :)
+    view(f.parent, underlying_left_halo_indices(instantiate(location), instantiate(topo), grid.Nx, grid.Hx), :, :)
 
 underlying_east_halo(f, grid, location, topo=topology(grid, 1)) =
-    view(f.parent, underlying_right_halo_indices(location, topo, grid.Nx, grid.Hx), :, :)
+    view(f.parent, underlying_right_halo_indices(instantiate(location), instantiate(topo), grid.Nx, grid.Hx), :, :)
 
 underlying_south_halo(f, grid, location, topo=topology(grid, 2)) =
-    view(f.parent, :, underlying_left_halo_indices(location, topo, grid.Ny, grid.Hy), :)
+    view(f.parent, :, underlying_left_halo_indices(instantiate(location), instantiate(topo), grid.Ny, grid.Hy), :)
 
 underlying_north_halo(f, grid, location, topo=topology(grid, 2)) =
-    view(f.parent, :, underlying_right_halo_indices(location, topo, grid.Ny, grid.Hy), :)
+    view(f.parent, :, underlying_right_halo_indices(instantiate(location), instantiate(topo), grid.Ny, grid.Hy), :)
 
 underlying_bottom_halo(f, grid, location, topo=topology(grid, 3)) =
-    view(f.parent, :, :, underlying_left_halo_indices(location, topo, grid.Nz, grid.Hz))
+    view(f.parent, :, :, underlying_left_halo_indices(instantiate(location), instantiate(topo), grid.Nz, grid.Hz))
 
 underlying_top_halo(f, grid, location, topo=topology(grid, 3)) =
-    view(f.parent, :, :, underlying_right_halo_indices(location, topo, grid.Nz, grid.Hz))
+    view(f.parent, :, :, underlying_right_halo_indices(instantiate(location), instantiate(topo), grid.Nz, grid.Hz))
 
 #####
 ##### Viewing boundary grid points (used to fill other halos)
 #####
 
 left_boundary_indices(loc, topo, N, H) = 1:H
-left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty
+left_boundary_indices(::Nothing, topo, N, H) = 1:0 # empty
 
 right_boundary_indices(loc, topo, N, H) = N-H+1:N
-right_boundary_indices(::Type{Face}, ::Type{Bounded}, N, H) = N-H:N+1
-right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty
+right_boundary_indices(::Face, ::Bounded, N, H) = N-H:N+1
+right_boundary_indices(::Nothing, topo, N, H) = 1:0 # empty
 
 underlying_left_boundary_indices(loc, topo, N, H) = 1+H:2H
-underlying_left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty
+underlying_left_boundary_indices(::Nothing, topo, N, H) = 1:0 # empty
 
 underlying_right_boundary_indices(loc, topo, N, H) = N+1:N+H
-underlying_right_boundary_indices(::Type{Face}, ::Type{Bounded}, N, H) = N+2:N+H+1
-underlying_right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty
+underlying_right_boundary_indices(::Face, ::Bounded, N, H) = N+2:N+H+1
+underlying_right_boundary_indices(::Nothing, topo, N, H) = 1:0 # empty
 
 underlying_west_boundary(f, grid, location, topo=topology(grid, 1)) =
-    view(f.parent, underlying_left_boundary_indices(location, topo, grid.Nx, grid.Hx), :, :)
+    view(f.parent, underlying_left_boundary_indices(instantiate(location), instantiate(topo), grid.Nx, grid.Hx), :, :)
 
 underlying_east_boundary(f, grid, location, topo=topo = topology(grid, 1)) =
-    view(f.parent, underlying_right_boundary_indices(location, topo, grid.Nx, grid.Hx), :, :)
+    view(f.parent, underlying_right_boundary_indices(instantiate(location), instantiate(topo), grid.Nx, grid.Hx), :, :)
 
 underlying_south_boundary(f, grid, location, topo=topology(grid, 2)) =
-    view(f.parent, :, underlying_left_boundary_indices(location, topo, grid.Ny, grid.Hy), :)
+    view(f.parent, :, underlying_left_boundary_indices(instantiate(location), instantiate(topo), grid.Ny, grid.Hy), :)
 
 underlying_north_boundary(f, grid, location, topo=topology(grid, 2)) =
-    view(f.parent, :, underlying_right_boundary_indices(location, topo, grid.Ny, grid.Hy), :)
+    view(f.parent, :, underlying_right_boundary_indices(instantiate(location), instantiate(topo), grid.Ny, grid.Hy), :)
 
 underlying_bottom_boundary(f, grid, location, topo=topology(grid, 3)) =
-    view(f.parent, :, :, underlying_left_boundary_indices(location, topo, grid.Nz, grid.Hz))
+    view(f.parent, :, :, underlying_left_boundary_indices(instantiate(location), instantiate(topo), grid.Nz, grid.Hz))
 
 underlying_top_boundary(f, grid, location, topo=topology(grid, 3)) =
-    view(f.parent, :, :, underlying_right_boundary_indices(location, topo, grid.Nz, grid.Hz))
+    view(f.parent, :, :, underlying_right_boundary_indices(instantiate(location), instantiate(topo), grid.Nz, grid.Hz))
 
 #####
 ##### Convenience functions

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/Distributed/distributed_fields.jl

@@ -12,11 +12,12 @@ function Field((LX, LY, LZ)::Tuple, grid::DistributedGrid, data, old_bcs, indice
     indices = validate_indices(indices, (LX, LY, LZ), grid)
     new_bcs = inject_halo_communication_boundary_conditions(old_bcs, arch.local_rank, arch.connectivity)
     buffers = FieldBoundaryBuffers(grid, data, new_bcs)
+
     return Field{LX, LY, LZ}(grid, data, new_bcs, indices, op, status, buffers)
 end
 
 const DistributedField      = Field{<:Any, <:Any, <:Any, <:Any, <:DistributedGrid}
 const DistributedFieldTuple = NamedTuple{S, <:NTuple{N, DistributedField}} where {S, N}
 
 # TODO: make sure the definition of architecture is consistent
-architecture(f::DistributedField) = child_architecture(architecture(f.grid))
+architecture(f::DistributedField) = child_architecture(architecture(f.grid))