Merge pull request #1641 from CliMA/ali/output-readers

docs/src/model_setup/boundary_conditions.md

@@ -283,13 +283,13 @@ IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0)
 
 julia> model.velocities.u
 Field located at (Face, Center, Center)
-├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (18, 18, 18)
+├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (16, 16, 16)
 ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=16, Ny=16, Nz=16)
 └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=Value, top=Value)
 
 julia> model.tracers.T
 Field located at (Center, Center, Center)
-├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (18, 18, 18)
+├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (16, 16, 16)
 ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=16, Ny=16, Nz=16)
 └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=Gradient, top=Value)
 ```

docs/src/model_setup/tracers.md

@@ -26,13 +26,13 @@ whose fields can be accessed via `model.tracers.T` and `model.tracers.S`.
 ```jldoctest tracers
 julia> model.tracers.T
 Field located at (Center, Center, Center)
-├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (66, 66, 66)
+├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (64, 64, 64)
 ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=64, Ny=64, Nz=64)
 └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux)
 
 julia> model.tracers.S
 Field located at (Center, Center, Center)
-├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (66, 66, 66)
+├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (64, 64, 64)
 ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=64, Ny=64, Nz=64)
 └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux)
 

src/AbstractOperations/AbstractOperations.jl

@@ -26,7 +26,7 @@ import Oceananigans.Fields: data, compute_at!
 ##### Basic functionality
 #####
 
-abstract type AbstractOperation{X, Y, Z, A, G, T} <: AbstractField{X, Y, Z, A, G, T} end
+abstract type AbstractOperation{X, Y, Z, A, G, T} <: AbstractField{X, Y, Z, A, G, T, 3} end
 
 const AF = AbstractField
 

src/BuoyancyModels/buoyancy_field.jl

@@ -11,14 +11,14 @@ import Oceananigans.Fields: compute!, compute_at!
 
 import Oceananigans: short_show
 
-struct BuoyancyField{B, S, A, D, G, T, C} <: AbstractDataField{Center, Center, Center, A, G, T}
+struct BuoyancyField{B, S, A, D, G, T, C} <: AbstractDataField{Center, Center, Center, A, G, T, 3}
             data :: D
     architecture :: A
             grid :: G
         buoyancy :: B
          tracers :: C
           status :: S
-    
+
     """
         BuoyancyField(data, grid, buoyancy, tracers)
 

src/Fields/Fields.jl

@@ -1,15 +1,14 @@
 module Fields
 
-export
-    Face, Center,
-    AbstractField, Field,
-    CenterField, XFaceField, YFaceField, ZFaceField,
-    ReducedField, AveragedField, ComputedField, KernelComputedField, BackgroundField,
-    interior, data,
-    xnode, ynode, znode, location,
-    set!, compute!, @compute,
-    VelocityFields, TracerFields, tracernames, PressureFields, TendencyFields,
-    interpolate, FieldSlicer
+export Face, Center
+export AbstractField, AbstractDataField, Field
+export CenterField, XFaceField, YFaceField, ZFaceField
+export ReducedField, AveragedField, ComputedField, KernelComputedField, BackgroundField
+export interior, data
+export xnode, ynode, znode, location
+export set!, compute!, @compute
+export VelocityFields, TracerFields, tracernames, PressureFields, TendencyFields
+export interpolate, FieldSlicer
 
 using Oceananigans.Architectures
 using Oceananigans.Grids
@@ -36,4 +35,8 @@ include("show_fields.jl")
 include("broadcasting_abstract_fields.jl")
 include("mapreduce_abstract_fields.jl")
 
+# Fallback: cannot infer boundary conditions.
+boundary_conditions(field) = nothing
+boundary_conditions(f::Union{Field, ReducedField, ComputedField, KernelComputedField}) = f.boundary_conditions
+
 end

src/Fields/abstract_field.jl

@@ -21,20 +21,21 @@ const ArchOrNothing = Union{AbstractArchitecture, Nothing}
 const GridOrNothing = Union{AbstractGrid, Nothing}
 
 """
-    AbstractField{X, Y, Z, A, G, T}
+    AbstractField{X, Y, Z, A, G, T, N}
 
 Abstract supertype for fields located at `(X, Y, Z)` on architecture `A`
-and defined on a grid `G` with eltype `T`.
+and defined on a grid `G` with eltype `T` and `N` dimensions.
 """
-abstract type AbstractField{X, Y, Z, A <: ArchOrNothing, G <: GridOrNothing, T} <: AbstractArray{T, 3} end
+abstract type AbstractField{X, Y, Z, A <: ArchOrNothing, G <: GridOrNothing, T, N} <: AbstractArray{T, N} end
 
 """
-    AbstractDataField{X, Y, Z, A, G}
+    AbstractDataField{X, Y, Z, A, G, T, N}
 
 Abstract supertype for fields with concrete data in settable underlying arrays,
-located at `(X, Y, Z)` on architecture `A` and defined on a grid `G` with eltype `T`.
+located at `(X, Y, Z)` on architecture `A` and defined on a grid `G` with eltype `T`
+and `N` dimensions.
 """
-abstract type AbstractDataField{X, Y, Z, A, G, T} <: AbstractField{X, Y, Z, A, G, T} end
+abstract type AbstractDataField{X, Y, Z, A, G, T, N} <: AbstractField{X, Y, Z, A, G, T, N} end
 
 Base.IndexStyle(::AbstractField) = IndexCartesian()
 
@@ -186,7 +187,7 @@ Base.fill!(f::AbstractDataField, val) = fill!(parent(f), val)
 # Don't use axes(f) to checkbounds; use axes(f.data)
 Base.checkbounds(f::AbstractField, I...) = Base.checkbounds(f.data, I...)
 
-@propagate_inbounds Base.getindex(f::AbstractDataField, i, j, k) = f.data[i, j, k]
+@propagate_inbounds Base.getindex(f::AbstractDataField, inds...) = getindex(f.data, inds...)
 
 # Linear indexing
 @propagate_inbounds Base.getindex(f::AbstractDataField, i::Int)  = parent(f)[i]
@@ -223,8 +224,6 @@ znodes(ψ::AbstractField) = znodes(location(ψ, 3), ψ.grid)
 
 nodes(ψ::AbstractField; kwargs...) = nodes(location(ψ), ψ.grid; kwargs...)
 
-Base.iterate(f::AbstractDataField, state=1) = iterate(f.data, state)
-
 #####
 ##### fill_halo_regions!
 #####
@@ -235,53 +234,49 @@ fill_halo_regions!(field::AbstractField, arch, args...) = fill_halo_regions!(fie
 ##### Field reductions
 #####
 
+const AbstractGPUDataField = AbstractDataField{X, Y, Z, GPU} where {X, Y, Z}
+
 """
     minimum(field::AbstractDataField; dims=:)
-
 Compute the minimum value of an Oceananigans `field` over the given dimensions (not including halo points).
 By default all dimensions are included.
 """
-minimum(field::AbstractDataField; dims=:) = minimum(interior_copy(field); dims=dims)
+minimum(field::AbstractGPUDataField; dims=:) = minimum(interior_copy(field); dims=dims)
 
 """
     minimum(f, field::AbstractDataField; dims=:)
-
 Returns the smallest result of calling the function `f` on each element of an Oceananigans `field`
 (not including halo points) over the given dimensions. By default all dimensions are included.
 """
-minimum(f, field::AbstractDataField; dims=:) = minimum(f, interior_copy(field); dims=dims)
+minimum(f, field::AbstractGPUDataField; dims=:) = minimum(f, interior_copy(field); dims=dims)
 
 """
     maximum(field::AbstractDataField; dims=:)
-
 Compute the maximum value of an Oceananigans `field` over the given dimensions (not including halo points).
 By default all dimensions are included.
 """
-maximum(field::AbstractDataField; dims=:) = maximum(interior_copy(field); dims=dims)
+maximum(field::AbstractGPUDataField; dims=:) = maximum(interior_copy(field); dims=dims)
 
 """
     maximum(f, field::AbstractDataField; dims=:)
-
 Returns the largest result of calling the function `f` on each element of an Oceananigans `field`
 (not including halo points) over the given dimensions. By default all dimensions are included.
 """
-maximum(f, field::AbstractDataField; dims=:) = maximum(f, interior_copy(field); dims=dims)
+maximum(f, field::AbstractGPUDataField; dims=:) = maximum(f, interior_copy(field); dims=dims)
 
 """
     mean(field::AbstractDataField; dims=:)
-
 Compute the mean of an Oceananigans `field` over the given dimensions (not including halo points).
 By default all dimensions are included.
 """
-mean(field::AbstractDataField; dims=:) = mean(interior_copy(field); dims=dims)
+mean(field::AbstractGPUDataField; dims=:) = mean(interior_copy(field); dims=dims)
 
 """
     mean(f::Function, field::AbstractDataField; dims=:)
-
 Apply the function `f` to each element of an Oceananigans `field` and take the mean over dimensions `dims`
 (not including halo points). By default all dimensions are included.
 """
-mean(f::Function, field::AbstractDataField; dims=:) = mean(f, interior_copy(field); dims=dims)
+mean(f::Function, field::AbstractGPUDataField; dims=:) = mean(f, interior_copy(field); dims=dims)
 
 # Risky to use these without tests. Docs would also be nice.
 Statistics.norm(a::AbstractField) = sqrt(mapreduce(x -> x * x, +, interior(a)))

src/Fields/computed_field.jl

@@ -4,7 +4,7 @@ using KernelAbstractions: @kernel, @index, Event
 using Oceananigans.Grids
 using Oceananigans.BoundaryConditions: fill_halo_regions!
 
-struct ComputedField{X, Y, Z, S, O, A, D, G, T, C} <: AbstractDataField{X, Y, Z, A, G, T}
+struct ComputedField{X, Y, Z, S, O, A, D, G, T, C} <: AbstractDataField{X, Y, Z, A, G, T, 3}
                    data :: D
            architecture :: A
                    grid :: G
@@ -76,7 +76,7 @@ function ComputedField(LX, LY, LZ, operand, arch, grid;
                        data = nothing,
                        recompute_safely = true,
                        boundary_conditions = ComputedFieldBoundaryConditions(grid, (LX, LY, LZ)))
-    
+
     # Architecturanigans
     operand_arch = architecture(operand)
     arch = isnothing(operand_arch) ? arch : operand_arch
@@ -105,7 +105,7 @@ function compute!(comp::ComputedField{LX, LY, LZ}, time=nothing) where {LX, LY,
                                       include_right_boundaries=true,
                                       location=(LX, LY, LZ))
 
-    compute_kernel! = _compute!(device(arch), workgroup, worksize) 
+    compute_kernel! = _compute!(device(arch), workgroup, worksize)
 
     event = compute_kernel!(comp.data, comp.operand; dependencies=Event(device(arch)))
 

src/Fields/field.jl

@@ -1,6 +1,6 @@
 using Adapt
 
-struct Field{X, Y, Z, A, D, G, T, B} <: AbstractDataField{X, Y, Z, A, G, T}
+struct Field{X, Y, Z, A, D, G, T, B} <: AbstractDataField{X, Y, Z, A, G, T, 3}
                    data :: D
            architecture :: A
                    grid :: G
@@ -29,7 +29,7 @@ julia> using Oceananigans, Oceananigans.Grids
 
 julia> ω = Field(Face, Face, Center, CPU(), RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
 Field located at (Face, Face, Center)
-├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (3, 3, 3)
+├── data: OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}}, size: (1, 1, 1)
 ├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1)
 └── boundary conditions: x=(west=Periodic, east=Periodic), y=(south=Periodic, north=Periodic), z=(bottom=ZeroFlux, top=ZeroFlux)
 ```

src/Fields/function_field.jl

@@ -1,5 +1,5 @@
 
-struct FunctionField{X, Y, Z, C, P, F, G, T} <: AbstractField{X, Y, Z, Nothing, G, T}
+struct FunctionField{X, Y, Z, C, P, F, G, T} <: AbstractField{X, Y, Z, Nothing, G, T, 3}
           func :: F
           grid :: G
          clock :: C

src/Fields/kernel_computed_field.jl

@@ -2,7 +2,7 @@ using Oceananigans: AbstractModel
 using Oceananigans.Grids
 using Oceananigans.Utils: tupleit
 
-struct KernelComputedField{X, Y, Z, A, S, D, G, T, K, B, F, P} <: AbstractDataField{X, Y, Z, A, G, T}
+struct KernelComputedField{X, Y, Z, A, S, D, G, T, K, B, F, P} <: AbstractDataField{X, Y, Z, A, G, T, 3}
                      data :: D
              architecture :: A
                      grid :: G
@@ -34,18 +34,18 @@ struct KernelComputedField{X, Y, Z, A, S, D, G, T, K, B, F, P} <: AbstractDataFi
 end
 
 """
-    KernelComputedField(X, Y, Z, kernel, model; 
-                        boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)), 
-                        computed_dependencies = (), 
-                        parameters = nothing, 
+    KernelComputedField(X, Y, Z, kernel, model;
+                        boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)),
+                        computed_dependencies = (),
+                        parameters = nothing,
                         data = nothing,
                         recompute_safely = true)
 
 Builds a `KernelComputedField` at `X, Y, Z` computed with `kernel` and `model.architecture` and `model.grid`, with `boundary_conditions`.
 
 `computed_dependencies` are an iterable of `AbstractField`s or other objects on which `compute!` is called prior to launching `kernel`.
 
-`data` is a three-dimensional `OffsetArray` of scratch space where the kernel computation is stored. 
+`data` is a three-dimensional `OffsetArray` of scratch space where the kernel computation is stored.
 
 If `data=nothing` (the default) then additional memory will be allocated to store the `data` of `KernelComputedField`.
 
@@ -57,10 +57,10 @@ Otherwise, `kernel` is launched with the function signature
 
 `kernel(data, grid, computed_dependencies..., parameters)`
 
-`recompute_safely` (default: `true`) determines whether the `KernelComputedField` is "recomputed" if embedded in the expression 
-tree of another operation. 
-    - If `recompute_safely=true`, the `KernelComputedField` is always recomputed. 
-    - If `recompute_safely=false`, the `KernelComputedField` will not be recomputed if its status is up-to-date. 
+`recompute_safely` (default: `true`) determines whether the `KernelComputedField` is "recomputed" if embedded in the expression
+tree of another operation.
+    - If `recompute_safely=true`, the `KernelComputedField` is always recomputed.
+    - If `recompute_safely=false`, the `KernelComputedField` will not be recomputed if its status is up-to-date.
 
 Example
 =======
@@ -111,7 +111,7 @@ function compute!(kcf::KernelComputedField{X, Y, Z}) where {X, Y, Z}
         tuple(kcf.data, kcf.grid, kcf.computed_dependencies...) :
         tuple(kcf.data, kcf.grid, kcf.computed_dependencies..., kcf.parameters)
 
-    event = launch!(arch, kcf.grid, :xyz, kcf.kernel, args...;   
+    event = launch!(arch, kcf.grid, :xyz, kcf.kernel, args...;
                     location=(X, Y, Z), include_right_boundaries=true)
 
     wait(device(arch), event)

src/Fields/reduced_field.jl

@@ -2,7 +2,11 @@ using Adapt
 
 import Oceananigans.BoundaryConditions: fill_halo_regions!
 
-abstract type AbstractReducedField{X, Y, Z, A, G, T, N} <: AbstractDataField{X, Y, Z, A, G, T} end
+#####
+##### AbstractReducedField stuff
+#####
+
+abstract type AbstractReducedField{X, Y, Z, A, G, T, N} <: AbstractDataField{X, Y, Z, A, G, T, 3} end
 
 const ARF = AbstractReducedField
 
@@ -16,7 +20,7 @@ function validate_reduced_dims(dims)
 
     # Check type
     dims isa DimsType || error("Reduced dims must be an integer or tuple of integers.")
-    
+
     # Check length
     length(dims) > 3  && error("Models are 3-dimensional. Cannot reduce over 4+ dimensions.")
 

src/Fields/show_fields.jl

@@ -14,14 +14,14 @@ short_show(field::AbstractField) = string(typeof(field).name.wrapper, " located
 short_show(field::AveragedField) = string("AveragedField over dims=$(field.dims) located at ", show_location(field), " of ", short_show(field.operand))
 short_show(field::ComputedField) = string("ComputedField located at ", show_location(field), " of ", short_show(field.operand))
 
-Base.show(io::IO, field::AbstractField{X, Y, Z, A, G}) where {X, Y, Z, A, G} =
+Base.show(io::IO, field::AbstractField{X, Y, Z, A}) where {X, Y, Z, A} =
     print(io, "$(short_show(field))\n",
-          "├── architecture: $A",
-          "└── grid: $G")
+          "├── architecture: $A\n",
+          "└── grid: $(short_show(field.grid))")
 
 Base.show(io::IO, field::Field) =
     print(io, "$(short_show(field))\n",
-          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
+          "├── data: $(typeof(field.data)), size: $(size(field))\n",
           "├── grid: $(short_show(field.grid))\n",
           "└── boundary conditions: $(short_show(field.boundary_conditions))")
 
@@ -30,26 +30,26 @@ show_status(status) = "time=$(status.time)"
 
 Base.show(io::IO, field::AveragedField) =
     print(io, "$(short_show(field))\n",
-          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
-          "├── grid: $(short_show(field.grid))", '\n',
-          "├── dims: $(field.dims)", '\n',
-          "├── operand: $(short_show(field.operand))", '\n',
-          "└── status: ", show_status(field.status), '\n')
+          "├── data: $(typeof(field.data)), size: $(size(field))\n",
+          "├── grid: $(short_show(field.grid))\n",
+          "├── dims: $(field.dims)\n",
+          "├── operand: $(short_show(field.operand))\n",
+          "└── status: ", show_status(field.status))
 
 Base.show(io::IO, field::ComputedField) =
     print(io, "$(short_show(field))\n",
-          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
-          "├── grid: $(short_show(field.grid))", '\n',
-          "├── operand: $(short_show(field.operand))", '\n',
-          "└── status: ", show_status(field.status), '\n')
+          "├── data: $(typeof(field.data)), size: $(size(field))\n",
+          "├── grid: $(short_show(field.grid))\n",
+          "├── operand: $(short_show(field.operand))\n",
+          "└── status: $(show_status(field.status))")
 
 Base.show(io::IO, field::KernelComputedField) =
     print(io, "$(short_show(field))\n",
-          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
-          "├── grid: $(short_show(field.grid))", '\n',
-          "├── computed_dependencies: $(Tuple(short_show(d) for d in field.computed_dependencies))", '\n',
-          "├── kernel: $(short_show(field.kernel))", '\n',
-          "└── status: ", show_status(field.status), '\n')
+          "├── data: $(typeof(field.data)), size: $(size(field))\n",
+          "├── grid: $(short_show(field.grid))\n",
+          "├── computed_dependencies: $(Tuple(short_show(d) for d in field.computed_dependencies))"\n,
+          "├── kernel: $(short_show(field.kernel))\n",
+          "└── status: $(show_status(field.status))")
 
 short_show(array::OffsetArray{T, D, A}) where {T, D, A} = string("OffsetArray{$T, $D, $A}")
 

src/Fields/zero_field.jl

@@ -1,3 +1,3 @@
-struct ZeroField <: AbstractField{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing} end
+struct ZeroField <: AbstractField{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 3} end
 
 @inline Base.getindex(::ZeroField, i, j, k) = 0