Validate location input to FieldBoundaryConditions

docs/src/model_setup/boundary_conditions.md

@@ -153,18 +153,18 @@ The default boundary condition can be changed by passing a positional argument t
 as in
 
 ```jldoctest
-julia> no_slip_bc = ValueBoundaryCondition(0.0)
-ValueBoundaryCondition: 0.0
+julia> no_slip_bc = ValueBoundaryCondition(0)
+ValueBoundaryCondition: 0
 
 julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing))
 Oceananigans.FieldBoundaryConditions, with boundary conditions
-├── west: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
-├── east: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
-├── south: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
-├── north: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
-├── bottom: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
+├── west: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
+├── east: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
+├── south: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
+├── north: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
+├── bottom: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
 ├── top: FluxBoundaryCondition: Nothing
-└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)
+└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0)
 ```
 
 ## Boundary condition structures

docs/src/model_setup/forcing_functions.md

@@ -21,7 +21,7 @@ By default, momentum and tracer forcing functions are assumed to be functions of
 u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t)
 
 grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
-model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,))
+model = NonhydrostaticModel(; grid, forcing=(; u=u_forcing))
 
 model.forcing.u
 
@@ -59,16 +59,15 @@ _(i)_ a single scalar parameter `s`, and _(ii)_ a `NamedTuple` of parameters, `p
 ```jldoctest parameterized_forcing
 # Forcing that depends on a scalar parameter `s`
 u_forcing_func(x, y, z, t, s) = s * z
-
 u_forcing = Forcing(u_forcing_func, parameters=0.1)
 
 # Forcing that depends on a `NamedTuple` of parameters `p`
 T_forcing_func(x, y, z, t, p) = - p.μ * exp(z / p.λ) * cos(p.k * x) * sin(p.ω * t)
-
 T_forcing = Forcing(T_forcing_func, parameters=(μ=1, λ=0.5, k=2π, ω=4π))
 
 grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
-model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing, T=T_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
+forcing = (u=u_forcing, T=T_forcing)
+model = NonhydrostaticModel(; grid, forcing, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
 
 model.forcing.T
 
@@ -113,7 +112,8 @@ S_forcing_func(x, y, z, t, S, μ) = - μ * S
 S_forcing = Forcing(S_forcing_func, parameters=0.01, field_dependencies=:S)
 
 grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
-model = NonhydrostaticModel(grid=grid, forcing=(w=w_forcing, S=S_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
+forcing = (w=w_forcing, S=S_forcing)
+model = NonhydrostaticModel(; grid, forcing, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
 
 model.forcing.w
 
@@ -182,17 +182,16 @@ using Oceananigans.Operators: ∂zᶠᶜᶠ, ℑxzᶠᵃᶜ
 function u_forcing_func(i, j, k, grid, clock, model_fields, ε)
     # The vertical derivative of buoyancy, interpolated to the u-velocity location:
     N² = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶠᶜᶠ, model_fields.b)
-
-    # Set to zero in unstable stratification where N² < 0:
-    N² = max(N², zero(typeof(N²)))
-
-    return @inbounds - ε * sqrt(N²) * model_fields.u[i, j, k]
+    N²⁺ = max(0, N²) # limit buoyancy frequency to strictly positive values
+    u_ijk = @inbounds model_fields.u[i, j, k]
+    return - ε * sqrt(N²⁺) * u_ijk
 end
 
 u_forcing = Forcing(u_forcing_func, discrete_form=true, parameters=1e-3)
 
 grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
-model = NonhydrostaticModel(grid=grid, tracers=:b, buoyancy=BuoyancyTracer(), forcing=(u=u_forcing, b=b_forcing))
+forcing = (u=u_forcing, b=b_forcing)
+model = NonhydrostaticModel(; grid, forcing, tracers=:b, buoyancy=BuoyancyTracer())
 
 model.forcing.b
 
@@ -227,7 +226,7 @@ of the velocity field are damped to zero everywhere on a time-scale of 1000 seco
 damping = Relaxation(rate = 1/1000)
 
 grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
-model = NonhydrostaticModel(grid=grid, forcing=(u=damping, v=damping, w=damping))
+model = NonhydrostaticModel(; grid, forcing=(u=damping, v=damping, w=damping))
 
 model.forcing.w
 
@@ -250,17 +249,18 @@ velocity fields to zero and restores temperature to a linear gradient in the bot
 ```jldoctest sponge_layer
 grid = RectilinearGrid(size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(-1, 0))
 
-        damping_rate = 1/100 # relax fields on a 100 second time-scale
+damping_rate = 1/100 # relax fields on a 100 second time-scale
 temperature_gradient = 0.001 # ⁰C m⁻¹
- surface_temperature = 20    # ⁰C (at z=0)
+surface_temperature = 20    # ⁰C (at z=0)
 
 target_temperature = LinearTarget{:z}(intercept=surface_temperature, gradient=temperature_gradient)
-       bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)
+bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)
 
 uvw_sponge = Relaxation(rate=damping_rate, mask=bottom_mask)
-  T_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=target_temperature)
+T_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=target_temperature)
 
-model = NonhydrostaticModel(grid=grid, forcing=(u=uvw_sponge, v=uvw_sponge, w=uvw_sponge, T=T_sponge), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
+forcing=(u=uvw_sponge, v=uvw_sponge, w=uvw_sponge, T=T_sponge)
+model = NonhydrostaticModel(; grid, forcing, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
 
 model.forcing.u
 
@@ -324,7 +324,7 @@ grid = RectilinearGrid(size=(32, 32, 32), x=(-10, 10), y=(-10, 10), z=(-4, 4),
                        topology=(Periodic, Periodic, Bounded))
 
 no_penetration = ImpenetrableBoundaryCondition()
-slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face),
+slip_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Face()),
                                    top=no_penetration, bottom=no_penetration)
 
 w_slip = ZFaceField(grid, boundary_conditions=slip_bcs)

src/BoundaryConditions/field_boundary_conditions.jl

@@ -138,7 +138,7 @@ FieldBoundaryConditions(default_bounded_bc::BoundaryCondition = NoFluxBoundaryCo
     FieldBoundaryConditions(west, east, south, north, bottom, top, immersed)
 
 """
-    FieldBoundaryConditions(grid, location, indices=(:, :, :);
+    FieldBoundaryConditions(grid, location::Tuple, indices=(:, :, :);
                             west     = default_auxiliary_bc(grid, boundary, loc),
                             east     = default_auxiliary_bc(grid, boundary, loc),
                             south    = default_auxiliary_bc(grid, boundary, loc),
@@ -148,7 +148,8 @@ FieldBoundaryConditions(default_bounded_bc::BoundaryCondition = NoFluxBoundaryCo
                             immersed = NoFluxBoundaryCondition())
 
 Return boundary conditions for auxiliary fields (fields whose values are
-derived from a model's prognostic fields) on `grid` and at `location`.
+derived from a model's prognostic fields) on `grid` and at `location`,
+which is a 3-tuple of either `Center()`, `Face()`, or `nothing`.
 
 Keyword arguments
 =================
@@ -171,14 +172,24 @@ and the topology in the boundary-normal direction is used:
 - `nothing` for `Bounded` directions and `Face`-located fields
 - `nothing` for `Flat` directions and/or `Nothing`-located fields
 """
-function FieldBoundaryConditions(grid::AbstractGrid, loc, indices=(:, :, :);
-                                 west     = default_auxiliary_bc(grid, Val(:west),   loc),
-                                 east     = default_auxiliary_bc(grid, Val(:east),   loc),
-                                 south    = default_auxiliary_bc(grid, Val(:south),  loc),
-                                 north    = default_auxiliary_bc(grid, Val(:north),  loc),
-                                 bottom   = default_auxiliary_bc(grid, Val(:bottom), loc),
-                                 top      = default_auxiliary_bc(grid, Val(:top),    loc),
-                                 immersed = DefaultBoundaryCondition())
+function FieldBoundaryConditions(grid::AbstractGrid, loc::Tuple, indices=(:, :, :); kwargs...)
+
+    for ℓ in loc
+        if !(ℓ isa Union{Nothing, Face, Center})
+            msg = string("Location $ℓ in $loc is not a valid location!", '\n',
+                         "Locations must be Center(), Face(), or nothing.")
+            throw(ArgumentError(msg))
+        end
+    end
+
+    # Build defaults _after_ validating the location
+    west     = get(kwargs, :west,     default_auxiliary_bc(grid, Val(:west),   loc))
+    east     = get(kwargs, :east,     default_auxiliary_bc(grid, Val(:east),   loc))
+    south    = get(kwargs, :south,    default_auxiliary_bc(grid, Val(:south),  loc))
+    north    = get(kwargs, :north,    default_auxiliary_bc(grid, Val(:north),  loc))
+    bottom   = get(kwargs, :bottom,   default_auxiliary_bc(grid, Val(:bottom), loc))
+    top      = get(kwargs, :top,      default_auxiliary_bc(grid, Val(:top),    loc))
+    immersed = get(kwargs, :immersed, DefaultBoundaryCondition())
 
     bcs = FieldBoundaryConditions(indices, west, east, south, north, bottom, top, immersed)
     return regularize_field_boundary_conditions(bcs, grid, loc)

src/MultiRegion/multi_region_boundary_conditions.jl

@@ -100,8 +100,11 @@ function (::MultiRegionFillHalo{<:WestAndEast})(c, westbc, eastbc, loc, grid, bu
     device_copy_to!(westdst, westsrc)
     device_copy_to!(eastdst, eastsrc)
 
-    view(parent(c), 1:H, :, :) .= westdst
-    view(parent(c), N+H+1:N+2H, :, :) .= eastdst
+    west_halo = view(parent(c), 1:H, :, :)
+    west_halo .= westdst
+
+    east_halo = view(parent(c), N+H+1:N+2H, :, :)
+    east_halo .= eastdst
 
     return nothing
 end
@@ -119,8 +122,11 @@ function (::MultiRegionFillHalo{<:SouthAndNorth})(c, southbc, northbc, loc, grid
     device_copy_to!(southdst, southsrc)
     device_copy_to!(northdst, northsrc)
 
-    view(parent(c), :, 1:H, :) .= southdst
-    view(parent(c), :, N+H+1:N+2H, :) .= northdst
+    south_halo = view(parent(c), :, 1:H, :)
+    south_halo .= southdst
+
+    north_halo = view(parent(c), :, N+H+1:N+2H, :)
+    north_halo .= northdst
 
     return nothing
 end
@@ -138,8 +144,8 @@ function (::MultiRegionFillHalo{<:West})(c, bc, loc, grid, buffers)
 
     device_copy_to!(dst, src)
 
-    p  = view(parent(c), 1:H, :, :)
-    p .= dst
+    west_halo = view(parent(c), 1:H, :, :)
+    west_halo .= dst
 
     return nothing
 end
@@ -153,8 +159,8 @@ function (::MultiRegionFillHalo{<:East})(c, bc, loc, grid, buffers)
 
     device_copy_to!(dst, src)
 
-    p  = view(parent(c), N+H+1:N+2H, :, :)
-    p .= dst
+    east_halo = view(parent(c), N+H+1:N+2H, :, :)
+    east_halo .= dst
 
     return nothing
 end
@@ -168,8 +174,8 @@ function (::MultiRegionFillHalo{<:South})(c, bc, loc, grid, buffers)
 
     device_copy_to!(dst, src)
 
-    p  = view(parent(c), :, 1:H, :)
-    p .= dst
+    south_halo = view(parent(c), :, 1:H, :)
+    south_halo .= dst
 
     return nothing
 end
@@ -183,8 +189,8 @@ function (::MultiRegionFillHalo{<:North})(c, bc, loc, grid, buffers)
 
     device_copy_to!(dst, src)
 
-    p = view(parent(c), :, N+H+1:N+2H, :)
-    p .= dst
+    north_halo = view(parent(c), :, N+H+1:N+2H, :)
+    north_halo .= dst
 
     return nothing
 end
@@ -194,15 +200,15 @@ end
 #####
 
 @inline getregion(fc::FieldBoundaryConditions, i) =
-            FieldBoundaryConditions(_getregion(fc.west, i),
-                                    _getregion(fc.east, i),
-                                    _getregion(fc.south, i),
-                                    _getregion(fc.north, i),
-                                    _getregion(fc.bottom, i),
-                                    _getregion(fc.top, i),
-                                    fc.immersed,
-                                    fc.kernels,
-                                    _getregion(fc.ordered_bcs, i))
+    FieldBoundaryConditions(_getregion(fc.west, i),
+                            _getregion(fc.east, i),
+                            _getregion(fc.south, i),
+                            _getregion(fc.north, i),
+                            _getregion(fc.bottom, i),
+                            _getregion(fc.top, i),
+                            fc.immersed,
+                            fc.kernels,
+                            _getregion(fc.ordered_bcs, i))
 
 @inline getregion(bc::BoundaryCondition, i) = BoundaryCondition(bc.classification, _getregion(bc.condition, i))
 
@@ -214,27 +220,27 @@ end
                                                 cf.field_dependencies_interp)
 
 @inline getregion(df::DiscreteBoundaryFunction, i) =
-            DiscreteBoundaryFunction(_getregion(df.func, i), _getregion(df.parameters, i))
+    DiscreteBoundaryFunction(_getregion(df.func, i), _getregion(df.parameters, i))
 
 @inline _getregion(fc::FieldBoundaryConditions, i) =
-            FieldBoundaryConditions(getregion(fc.west, i),
-                                    getregion(fc.east, i),
-                                    getregion(fc.south, i),
-                                    getregion(fc.north, i),
-                                    getregion(fc.bottom, i),
-                                    getregion(fc.top, i),
-                                    fc.immersed,
-                                    fc.kernels,
-                                    getregion(fc.ordered_bcs, i))
+    FieldBoundaryConditions(getregion(fc.west, i),
+                            getregion(fc.east, i),
+                            getregion(fc.south, i),
+                            getregion(fc.north, i),
+                            getregion(fc.bottom, i),
+                            getregion(fc.top, i),
+                            fc.immersed,
+                            fc.kernels,
+                            getregion(fc.ordered_bcs, i))
 
 @inline _getregion(bc::BoundaryCondition, i) = BoundaryCondition(bc.classification, getregion(bc.condition, i))
 
 @inline _getregion(cf::ContinuousBoundaryFunction{X, Y, Z, I}, i) where {X, Y, Z, I} =
-            ContinuousBoundaryFunction{X, Y, Z, I}(cf.func,
-                                                getregion(cf.parameters, i),
-                                                cf.field_dependencies,
-                                                cf.field_dependencies_indices,
-                                                cf.field_dependencies_interp)
+    ContinuousBoundaryFunction{X, Y, Z, I}(cf.func,
+                                           getregion(cf.parameters, i),
+                                           cf.field_dependencies,
+                                           cf.field_dependencies_indices,
+                                           cf.field_dependencies_interp)
 
 @inline _getregion(df::DiscreteBoundaryFunction, i) = DiscreteBoundaryFunction(getregion(df.func, i), getregion(df.parameters, i))
 

src/MultiRegion/multi_region_field.jl

@@ -193,7 +193,7 @@ function inject_regional_bcs(grid, connectivity, loc, indices;
     return FieldBoundaryConditions(indices, west, east, south, north, bottom, top, immersed)
 end
 
-FieldBoundaryConditions(mrg::MultiRegionGrids, loc, indices; kwargs...) =
+FieldBoundaryConditions(mrg::MultiRegionGrids, loc::Tuple, indices::MultiRegionObject; kwargs...) =
     construct_regionally(inject_regional_bcs, mrg, mrg.connectivity, Reference(loc), indices; kwargs...)
 
 function Base.show(io::IO, field::MultiRegionField)

test/test_boundary_conditions.jl

@@ -258,6 +258,10 @@ end
 
         grid = bbb_grid
 
+        @test_throws ArgumentError FieldBoundaryConditions(grid, (Center, Center(), Center()))
+        @test_throws ArgumentError FieldBoundaryConditions(grid, (Center(), Face, Center()))
+        @test_throws ArgumentError FieldBoundaryConditions(grid, (Center(), Center(), Nothing))
+
         T_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()),
                                         east = ValueBoundaryCondition(simple_bc),
                                         west = ValueBoundaryCondition(simple_bc),

test/test_broadcasting.jl

@@ -41,7 +41,9 @@ include("dependencies_for_runtests.jl")
 
         a2 = CenterField(three_point_grid)
 
-        b2_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), top=OpenBoundaryCondition(0), bottom=OpenBoundaryCondition(0))
+
+        loc = (Center(), Center(), Face())
+        b2_bcs = FieldBoundaryConditions(grid, loc, top=OpenBoundaryCondition(0), bottom=OpenBoundaryCondition(0))
         b2 = ZFaceField(three_point_grid, boundary_conditions=b2_bcs)
 
         b2 .= 1

test/test_forcings.jl

@@ -163,7 +163,8 @@ function advective_and_multiple_forcing(arch)
     constant_slip = AdvectiveForcing(w=1)
     zero_slip = AdvectiveForcing(w=0)
     no_penetration = ImpenetrableBoundaryCondition()
-    slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), top=no_penetration, bottom=no_penetration)
+    loc = (Center(), Center(), Face())
+    slip_bcs = FieldBoundaryConditions(grid, loc, top=no_penetration, bottom=no_penetration)
     slip_velocity = ZFaceField(grid, boundary_conditions=slip_bcs)
     set!(slip_velocity, 1)
     velocity_field_slip = AdvectiveForcing(w=slip_velocity)