Upgrade fields to store boundary conditions

Project.toml

@@ -1,6 +1,6 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.22.0"
+version = "0.23.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/model_setup/boundary_conditions.md

@@ -75,18 +75,19 @@ white_noise_T_bc = FluxBoundaryCondition(Q)
 You can also specify the boundary condition via a function. For z boundary conditions the function will be called with
 the signature
 ```
-f(i, j, grid, t, U, C, params)
+f(i, j, grid, time, iteration, U, C, params)
 ```
-where `i, j` is the grid index, `grid` is `model.grid`, `t` is the `model.clock.time`, `U` is the named tuple
-`model.velocities`, `C` is the named tuple `C.tracers`, and `params` is the user-defined `model.parameters`. The
-signature is similar for x and y boundary conditions expect that `i, j` is replaced with `j, k` and `i, k` respectively.
+where `i, j` is the grid index, `grid` is `model.grid`, `time` is the `model.clock.time`, `iteration` is the
+`model.clock.iteration`, `U` is the named tuple `model.velocities`, `C` is the named tuple `C.tracers`, and `params` is
+the user-defined `model.parameters`. The signature is similar for x and y boundary conditions expect that `i, j` is
+replaced with `j, k` and `i, k` respectively.
 
 We can add a fourth example now:
 4. A spatially varying and time-dependent heating representing perhaps a localized source of heating modulated by a
    diurnal cycle.
 ```@example
 using Oceananigans # hide
-@inline Q(i, j, grid, t, U, C, params) = @inbounds exp(-(grid.xC[i]^2 + grid.yC[j]^2)) * sin(2π*t)
+@inline Q(i, j, grid, t, I, U, C, params) = @inbounds exp(-(grid.xC[i]^2 + grid.yC[j]^2)) * sin(2π*t)
 localized_heating_bc = FluxBoundaryCondition(Q)
 ```
 

src/BoundaryConditions/BoundaryConditions.jl

@@ -9,24 +9,20 @@ export
     FieldBoundaryConditions, UVelocityBoundaryConditions, VVelocityBoundaryConditions,
     WVelocityBoundaryConditions, TracerBoundaryConditions, PressureBoundaryConditions,
     DiffusivityBoundaryConditions,
-    SolutionBoundaryConditions, TendenciesBoundaryConditions, PressureBoundaryConditions,
-    DiffusivityBoundaryConditions, DiffusivitiesBoundaryConditions,
-    ModelBoundaryConditions,
     BoundaryFunction,
     apply_z_bcs!, apply_y_bcs!,
     fill_halo_regions!, zero_halo_regions!
 
 using CUDAnative
 
+using Oceananigans: Cell, Face
 using Oceananigans.Architectures
 using Oceananigans.Grids
-using Oceananigans.Fields
 
 include("boundary_condition_types.jl")
 include("boundary_condition.jl")
 include("coordinate_boundary_conditions.jl")
 include("field_boundary_conditions.jl")
-include("solution_and_model_boundary_conditions.jl")
 include("boundary_function.jl")
 include("show_boundary_conditions.jl")
 

src/BoundaryConditions/apply_flux_bcs.jl

@@ -7,39 +7,34 @@ using Oceananigans.Utils: @loop_xy, @loop_xz, launch_config
 ##### Inhomogeneous Value and Gradient boundary conditions are handled by filling halos.
 #####
 
-# Avoid some computation / memory accesses for Value, Gradient, Periodic, NoPenetration,
-# and no-flux boundary conditions --- every boundary condition that does *not* prescribe
-# a non-trivial flux.
-const NotFluxBC = Union{VBC, GBC, PBC, NPBC, NFBC}
-apply_z_bcs!(Gc, arch, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing
-apply_y_bcs!(Gc, arch, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing
-
 """
-    apply_z_bcs!(Gc, arch, grid, bottom_bc, top_bc, args...)
+    apply_z_bcs!(Gc, arch, grid, args...)
 
 Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
 the source term `Gc` at the top and bottom.
 """
-function apply_z_bcs!(Gc, arch, grid, bottom_bc, top_bc, args...)
-    @launch(device(arch), config=launch_config(grid, :xy),
-            _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...))
+function apply_z_bcs!(Gc, c, arch, args...)
+    @launch(device(arch), config=launch_config(Gc.grid, :xy),
+            _apply_z_bcs!(Gc.data, Gc.grid, c.boundary_conditions.z.bottom,
+                          c.boundary_conditions.z.top, args...))
     return nothing
 end
 
 """
-    apply_y_bcs!(Gc, arch, grid, left_bc, right_bc, args...)
+    apply_y_bcs!(Gc, arch, grid, args...)
 
 Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
 the source term `Gc` at the left and right.
 """
-function apply_y_bcs!(Gc, arch, grid, left_bc, right_bc, args...)
-    @launch(device(arch), config=launch_config(grid, :xz),
-            _apply_y_bcs!(Gc, grid, left_bc, right_bc, args...))
+function apply_y_bcs!(Gc, c, arch, args...)
+    @launch(device(arch), config=launch_config(Gc.grid, :xz),
+            _apply_y_bcs!(Gc.data, Gc.grid, c.boundary_conditions.y.left,
+                          c.boundary_conditions.y.right, args...))
     return nothing
 end
 
 """
-    _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
+    _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
 
 Apply a top and/or bottom boundary condition to variable `c`.
 """
@@ -111,5 +106,5 @@ If `bottom_bc.condition` is a function, the function must have the signature
 @inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) =
     @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / ΔzF(i, j, 1, grid)
 
-@inline apply_z_left_bc!(Gc, bottom_flux::BC{<:Flux}, i, k, grid, args...) =
+@inline apply_y_left_bc!(Gc, left_flux::BC{<:Flux}, i, k, grid, args...) =
     @inbounds Gc[i, 1, k] += getbc(left_flux, i, k, grid, args...) / Δy(i, 1, k, grid)

src/BoundaryConditions/boundary_function.jl

@@ -27,12 +27,3 @@ struct BoundaryFunction{B, X1, X2, F} <: Function
         new{B, X1, X2, typeof(func)}(func)
     end
 end
-
-@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, time, args...) where {Y, Z} =
-    bc.func(ynode(Y, j, grid), znode(Z, k, grid), time)
-
-@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, time, args...) where {X, Z} =
-    bc.func(xnode(X, i, grid), znode(Z, k, grid), time)
-
-@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, time, args...) where {X, Y} =
-    bc.func(xnode(X, i, grid), ynode(Y, j, grid), time)

src/BoundaryConditions/field_boundary_conditions.jl

@@ -14,63 +14,53 @@ Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for
 """
 FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))
 
-default_tracer_bc(::Grids.Periodic) = PeriodicBoundaryCondition()
-default_tracer_bc(::Bounded)  = NoFluxBoundaryCondition()
-default_tracer_bc(::Flat)     = PeriodicBoundaryCondition()
+DefaultBoundaryCondition(::Union{Grids.Periodic, Flat}, loc) = PeriodicBoundaryCondition()
 
-# Right now it seems all defaults are the same except for no-normal
-# flow velocity BCs (which are treated in default_bc below), but
-# that might change in the future.
-const default_velocity_bc = default_tracer_bc
-const default_pressure_bc = default_tracer_bc
-const default_diffusivity_bc = default_tracer_bc
-
-"""
-    default_bc(grid, field, dim)
-
-Returns the default boundary condition for a `field` on some `grid` and along some dimension `dim`.
-"""
-function default_bc(grid, field, dim)
-    top = topology(grid)[dim]
-
-    # No normal flow velocity boundary conditions by default.
-    top isa Bounded && field == :u && dim == 1 && return NoPenetrationBoundaryCondition()
-    top isa Bounded && field == :v && dim == 2 && return NoPenetrationBoundaryCondition()
-    top isa Bounded && field == :w && dim == 3 && return NoPenetrationBoundaryCondition()
-
-    field in (:u, :v, :w) && (field = :velocity)
-    default_field_bc = Symbol(:default_, field, :_bc)
-    return @eval $default_field_bc($top)
-end
+DefaultBoundaryCondition(::Bounded, ::Cell) = NoFluxBoundaryCondition()
+DefaultBoundaryCondition(::Bounded, ::Face) = NoPenetrationBoundaryCondition()
 
 function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir)
     if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc)
-        e = "Cannot specify $left_name or $right_name boundary conditions with $topology topology in $dir-direction."
+        e = "Cannot specify $left_name or $right_name boundary conditions with " *
+            "$topology topology in $dir-direction."
         throw(ArgumentError(e))
     end
     return true
 end
 
 """
-    FieldBoundaryConditions(grid::AbstractGrid; field_type,
-                            west=default_bc(grid, field_type, 1), east=default_bc(grid, field_type, 1),
-                            south=default_bc(grid, field_type, 2), north=default_bc(grid, field_type, 2),
-                            bottom=default_bc(grid, field_type, 3), top=default_bc(grid, field_type, 3))
+    FieldBoundaryConditions(grid, loc;
+          east = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
+          west = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
+         south = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
+         north = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
+        bottom = DefaultBoundaryCondition(topology(grid)[3], loc[3]),
+           top = DefaultBoundaryCondition(topology(grid)[3], loc[3]))
+
+Construct `FieldBoundaryConditions` for a field with location `loc` (a 3-tuple of `Face` or `Cell`)
+defined on `grid` (the grid's topology is what defined the default boundary conditions that are
+imposed).
 
-Construct `FieldBoundaryConditions` for a field with type `field_type` (choices are :u, :v, :w, :tracer,
-:pressure, :diffusivity). Specific boundary conditions can be applied along the x dimension with the
-`west` and `east` kwargs, along the y-dimension with the `south` and `north` kwargs, and along the
-z-dimension with the `bottom` and `top` kwargs. Default boundary conditions are applied depending on
-the `field_type`.
+Specific boundary conditions can be applied along the x dimension with the `west` and `east` kwargs,
+along the y-dimension with the `south` and `north` kwargs, and along the z-dimension with the `bottom`
+and `top` kwargs.
 """
-function FieldBoundaryConditions(grid::AbstractGrid; field_type,
-                                 west=default_bc(grid, field_type, 1), east=default_bc(grid, field_type, 1),
-                                 south=default_bc(grid, field_type, 2), north=default_bc(grid, field_type, 2),
-                                 bottom=default_bc(grid, field_type, 3), top=default_bc(grid, field_type, 3))
+function FieldBoundaryConditions(grid, loc;
+      east = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
+      west = DefaultBoundaryCondition(topology(grid)[1], loc[1]),
+     south = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
+     north = DefaultBoundaryCondition(topology(grid)[2], loc[2]),
+    bottom = DefaultBoundaryCondition(topology(grid)[3], loc[3]),
+       top = DefaultBoundaryCondition(topology(grid)[3], loc[3]))
+
     TX, TY, TZ = topology(grid)
-    validate_bcs(TX, west,   east, default_bc(grid, field_type, 1), :west,   :east, :x)
-    validate_bcs(TY, south, north, default_bc(grid, field_type, 2), :south, :north, :y)
-    validate_bcs(TZ, bottom,  top, default_bc(grid, field_type, 3), :bottom,  :top, :z)
+    x_default_bc = DefaultBoundaryCondition(topology(grid)[1], loc[1])
+    y_default_bc = DefaultBoundaryCondition(topology(grid)[2], loc[2])
+    z_default_bc = DefaultBoundaryCondition(topology(grid)[3], loc[3])
+
+    validate_bcs(TX, west,   east, x_default_bc, :west,   :east, :x)
+    validate_bcs(TY, south, north, y_default_bc, :south, :north, :y)
+    validate_bcs(TZ, bottom,  top, z_default_bc, :bottom,  :top, :z)
 
     x = CoordinateBoundaryConditions(west, east)
     y = CoordinateBoundaryConditions(south, north)
@@ -79,9 +69,9 @@ function FieldBoundaryConditions(grid::AbstractGrid; field_type,
     return FieldBoundaryConditions(x, y, z)
 end
 
-  UVelocityBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:u; kwargs...)
-  VVelocityBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:v; kwargs...)
-  WVelocityBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:w; kwargs...)
-     TracerBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:tracer; kwargs...)
-#   PressureBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:pressure, kwargs...)
-#DiffusivityBoundaryConditions(grid; kwargs...) = FieldBoundaryConditions(grid, field_type=:diffusivity, kwargs...)
+  UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face(), Cell(), Cell()); user_defined_bcs...)
+  VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Face(), Cell()); user_defined_bcs...)
+  WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Face()); user_defined_bcs...)
+     TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)
+   PressureBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)
+DiffusivityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)

src/BoundaryConditions/fill_halo_regions.jl

@@ -5,57 +5,28 @@
 fill_halo_regions!(::Nothing, args...) = nothing
 
 """
-    fill_halo_regions!(fields, bcs, arch, grid)
+    fill_halo_regions!(fields, arch)
 
-Fill halo regions for each field in the tuple `fields` according
-to the single instance of `FieldBoundaryConditions` in `bcs`, possibly recursing into
-`fields` if it is a nested tuple-of-tuples.
+Fill halo regions for each field in the tuple `fields` according to their boundary
+conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples.
 """
-function fill_halo_regions!(fields::NamedTuple, bcs::FieldBoundaryConditions, arch, grid, args...)
+function fill_halo_regions!(fields::Union{Tuple,NamedTuple}, arch, args...)
     for field in fields
-        fill_halo_regions!(field, bcs, arch, grid, args...)
+        fill_halo_regions!(field, arch, args...)
     end
     return nothing
 end
 
+fill_halo_regions!(field, arch, args...) =
+    fill_halo_regions!(field.data, field.boundary_conditions, arch, field.grid, args...)
+
 "Fill halo regions in x, y, and z for a given field."
 function fill_halo_regions!(c::AbstractArray, fieldbcs, arch, grid, args...)
-      fill_west_halo!(c, fieldbcs.x.left,  arch, grid, args...)
-      fill_east_halo!(c, fieldbcs.x.right, arch, grid, args...)
-
-     fill_south_halo!(c, fieldbcs.y.left,  arch, grid, args...)
-     fill_north_halo!(c, fieldbcs.y.right, arch, grid, args...)
-
-     fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
-        fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
-    return nothing
-end
-
-#####
-##### Convenience methods for filling halo regions of tupled fields
-#####
-
-function tupled_fill_halo_regions!(fields, bcs, arch, grid, args...)
-    for (field, fieldbcs) in zip(fields, bcs)
-        fill_halo_regions!(field, fieldbcs, arch, grid, args...)
-    end
+      fill_west_halo!(c, fieldbcs.x.left,   arch, grid, args...)
+      fill_east_halo!(c, fieldbcs.x.right,  arch, grid, args...)
+     fill_south_halo!(c, fieldbcs.y.left,   arch, grid, args...)
+     fill_north_halo!(c, fieldbcs.y.right,  arch, grid, args...)
+    fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
+       fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
     return nothing
 end
-
-"""
-    fill_halo_regions!(fields, bcs, arch, grid)
-
-Fill halo regions for all fields in the `Tuple` `fields` according
-to the corresponding `Tuple` of `bcs`.
-"""
-fill_halo_regions!(fields::Tuple, bcs::Tuple, arch, grid, args...) =
-    tupled_fill_halo_regions!(fields, bcs, arch, grid, args...)
-
-"""
-    fill_halo_regions!(fields, bcs, arch, grid)
-
-Fill halo regions for all fields in the `NamedTuple` `fields` according
-to the corresponding `NamedTuple` of `bcs`.
-"""
-fill_halo_regions!(fields::NamedTuple{S}, bcs::NamedTuple{S}, arch, grid, args...) where S =
-    tupled_fill_halo_regions!(fields, bcs, arch, grid, args...)

src/BoundaryConditions/show_boundary_conditions.jl

@@ -1,40 +1,32 @@
+import Oceananigans.Grids: short_show
+
 #####
 ##### BoundaryCondition
 #####
 
 Base.show(io::IO, bc::BC{C, T}) where {C, T} =
-    println(io, "BoundaryCondition: type=$C, condition=$(bc.condition)")
+    print(io, "BoundaryCondition: type=$C, condition=$(bc.condition)")
 
 #####
 ##### FieldBoundaryConditions
 #####
 
+bctype_str(::FBC)  = "Flux"
+bctype_str(::PBC)  = "Periodic"
+bctype_str(::NPBC) = "NoPenetration"
+bctype_str(::VBC)  = "Value"
+bctype_str(::GBC)  = "Gradient"
+bctype_str(::NFBC) = "NoFlux"
+
+short_show(fbcs::FieldBoundaryConditions) =
+    string("x=(west=$(bctype_str(fbcs.x.left)), east=$(bctype_str(fbcs.x.right))), ",
+           "y=(south=$(bctype_str(fbcs.y.left)), north=$(bctype_str(fbcs.y.right))), ",
+           "z=(bottom=$(bctype_str(fbcs.z.left)), top=$(bctype_str(fbcs.z.right)))")
+
 show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") =
     string("Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions", '\n',
            padding, "├── x: ", typeof(bcs.x), '\n',
            padding, "├── y: ", typeof(bcs.y), '\n',
            padding, "└── z: ", typeof(bcs.z))
 
 Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs))
-
-#####
-##### ModelBoundaryConditions
-#####
-
-function show_solution_boundary_conditions(bcs, padding)
-    stringtuple = Tuple(string(
-                  padding, "├── ", field, ": ",
-                  show_field_boundary_conditions(getproperty(bcs, field), padding * "│   "), '\n')
-                  for field in propertynames(bcs)[1:end-1])
-    return string("Oceananigans.SolutionBoundaryConditions ",
-                  "(NamedTuple{(:u, :v, :w, ...)}) with boundary conditions ", '\n', stringtuple...,
-                  padding, "└── ", propertynames(bcs)[end], ": ",
-                  show_field_boundary_conditions(bcs[end], padding * "    "))
-end
-
-Base.show(io::IO, bcs::ModelBoundaryConditions) =
-    print(io,
-          "Oceananigans.ModelBoundaryConditions (NamedTuple{(:solution, :tendency, :pressure, :diffusivities)}) with ", '\n',
-          "├── solution: ", show_solution_boundary_conditions(bcs.solution, "│   "), '\n',
-          "├── tendency: ", show_solution_boundary_conditions(bcs.tendency, "│   "), '\n',
-          "└── pressure: ", show_field_boundary_conditions(bcs.pressure, "    "))