Major refactor of FieldBoundaryConditions (#1843)

* Reworking FieldBoundaryConditions

* Refactors FieldBoundaryConditions

* Fixes many tests

* More test fixes

* Update docs/src/model_setup/boundary_conditions.md

Co-authored-by: Navid C. Constantinou <[email protected]>

* Nuke AuxiliaryFieldBoundaryConditions; just use special FieldBoundaryConditions

* No more x, y, z, right, left...

* Banish DiffusivityBoundaryConditions

* Banish CoordinateBoundaryConditions

* Fixes boundary condition test with ContinuousBoundaryFunction

* Adds better default bcs for VelocityFields and TracerFields

* New syntax!

* Iterate better over boundary conditions

* Fixes bug in field_boundary_conditions

* Fix TracerFields implementation

* Were not ready for zero flux boundary conditions...

* Fixes bottom-top mix up

* Fixes broken ref to Nx in Poisson solver test

* Default for ZFaceField was changed so now broadcasting tests need to be changed too

* Fix regularize_boundary_condition for CubedSpheres

* West vs east

* Fix issues with TKEBasedVerticalDiffusivity

* Fiddling with AMD boundary conditions

* Adds regularize_boundary_conditions for regular bcs on ConformalCubedSphereGrid

* Fixes bug in determining tracer diffusivity bcs with AMD

* Adds location to FieldBoundaryConditions for broadcasting test

* Properly inject CubedSphere exchange bcs during regularization

* Adds a method for nested NamedTuple boundary conditions

* Fixes VelocityFields constructor to use valid regularized boundary conditions

* Updates to new syntax for stratified couette flow

* Elide tendency fields when velocities are prescribed

* Import regularize_boundary_conditions into field_tuples

* Try 2 with poisson solver test refactor

* Simplifies default calculation

* Fixes ommission of tracer tendencies when using prescribed velocities

* Import regularize_field_boundary_conditions into poisson solver tests

* Eta is prognostic, not w

* Adds missing comma in show methods

* Fix doctests

* Fix doctests

* Fix doctests

* Fix docstring in Field constructor

* Use FluxBoundaryCondition

* Where bc?

* Ancient issues with library.md

* Comments out library to see what happens

* Deletes misleading and impertinent comments

* Bump to 0.59.0

* Comments out appendix in docs

* Cleans up BoundaryConditions code + docstrings

* Docs despration

* update BoundaryConditions pages

* can't document call syntax; see JuliaDocs/Documenter.jl#228

* back to including everything

* improves some docstrings

* Updates validation experiments

* Fixes typos in TKEBasedVerticalMixing validation

Co-authored-by: Navid C. Constantinou <[email protected]>

Project.toml

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

docs/make.jl

@@ -116,7 +116,6 @@ numerical_pages = [
     "Large eddy simulation" => "numerical_implementation/large_eddy_simulation.md"
 ]
 
-
 appendix_pages = [
     "Staggered grid" => "appendix/staggered_grid.md",
     "Fractional step method" => "appendix/fractional_step.md",
@@ -161,7 +160,7 @@ makedocs(bib,
    doctest = true,
     strict = true,
      clean = true,
- checkdocs = :none  # Should fix our docstring so we can use checkdocs=:exports with strict=true.
+ checkdocs = :none # Should fix our docstring so we can use checkdocs=:exports with strict=true.
 )
 
 deploydocs(
@@ -173,4 +172,4 @@ deploydocs(
 @info "Cleaning up temporary .jld2 and .nc files created by doctests..."
 for file in vcat(glob("docs/*.jld2"), glob("docs/*.nc"))
     rm(file)
-end
+end

docs/make_example.jl

@@ -53,9 +53,9 @@ end
 #####
 
 example_pages = [
-                 #"Internal wave"                    => "generated/internal_wave.md",
-                 #"Geostrophic adjustment"            => "generated/geostrophic_adjustment.md",
-                 "Shallow water Bickley jet"    => "generated/shallow_water_Bickley_jet.md"
+                 #"Internal wave"            => "generated/internal_wave.md",
+                 #"Geostrophic adjustment"   => "generated/geostrophic_adjustment.md",
+                 "Shallow water Bickley jet" => "generated/shallow_water_Bickley_jet.md"
                  ]
 
 pages = [

docs/src/appendix/benchmarks.md

@@ -10,7 +10,7 @@ format the benchmark results.
 
 This is a benchmark of a simple "static ocean" configuration. The time stepping and Poisson
 solver still takes the same amount of time whether the ocean is static or active, so it is
-quite indicative of actual performance. It tests the performance of a bare-bones
+indicative of actual performance. It tests the performance of a bare-bones
 horizontally-periodic model with `topology = (Periodic, Periodic, Bounded)`.
 
 ```

docs/src/appendix/fractional_step.md

@@ -35,12 +35,3 @@ methods, which is that "*while the velocity can be reliably computed to second-o
 in time and space, the pressure is typically only first-order accurate in the ``L_\infty``-norm.*" 
 The numerical boundary conditions must be carefully accounted for to ensure the second-order 
 accuracy promised by the fractional step methods.
-
-We are currently investigating whether our projection method is indeed second-order accurate 
-in both velocity and pressure. However, it may not matter too much for simulating high Reynolds 
-number geophysical fluids as [Brown01](@cite) conclude that "*Quite often, semi-implicit projection 
-methods are applied to problems in which the viscosity is small. Since the predicted first-order 
-errors in the pressure are scaled by ``\nu``, it is not clear whether the improved pressure-update 
-formula is beneficial in such situations. [...] Finally, in some applications of projection methods, 
-second-order accuracy in the pressure may not be relevant or in some cases even possible due 
-to the treatment of other terms in the equations.*"

docs/src/appendix/library.md

@@ -36,15 +36,17 @@ Pages   = [
     "BoundaryConditions/BoundaryConditions.jl",
     "BoundaryConditions/boundary_condition_classifications.jl",
     "BoundaryConditions/boundary_condition.jl",
-    "BoundaryConditions/coordinate_boundary_conditions.jl",
+    "BoundaryConditions/discrete_boundary_function.jl",
+    "BoundaryConditions/continuous_boundary_function.jl",
     "BoundaryConditions/field_boundary_conditions.jl",
-    "BoundaryConditions/boundary_function.jl",
-    "BoundaryConditions/parameterized_boundary_condition.jl",
     "BoundaryConditions/show_boundary_conditions.jl",
     "BoundaryConditions/fill_halo_regions.jl",
+    "BoundaryConditions/fill_halo_regions_value_gradient.jl",
+    "BoundaryConditions/fill_halo_regions_open.jl",
+    "BoundaryConditions/fill_halo_regions_periodic.jl",
+    "BoundaryConditions/fill_halo_regions_flux.jl",
+    "BoundaryConditions/fill_halo_regions_nothing.jl",
     "BoundaryConditions/apply_flux_bcs.jl",
-    "BoundaryConditions/apply_value_gradient_bcs.jl",
-    "BoundaryConditions/apply_normal_flow_bcs.jl"
 ]
 ```
 
@@ -107,7 +109,6 @@ Pages   = [
     "Fields/field.jl",
     "Fields/set!.jl",
     "Fields/show_fields.jl",
-    "Fields/field_utils.jl"
 ]
 ```
 

docs/src/model_setup/boundary_conditions.md

@@ -45,10 +45,26 @@ Here's a short list of useful tips for defining and prescribing boundary conditi
    If there is a known `diffusivity`, you can express `FluxBoundaryCondition(flux)` using `GradientBoundaryCondition(-flux / diffusivity)` (aka "Neumann" boundary condition).
    But when `diffusivity` is not known or is variable (as for large eddy simulation, for example), it's convenient and more straightforward to apply `FluxBoundaryCondition`.
 
+## Default boundary conditions
+
+By default, periodic boundary conditions are applied on all fields along periodic dimensions. Otherwise tracers
+get no-flux boundary conditions and velocities get free-slip and no normal flow boundary conditions.
+
+## Boundary condition structures
+
+Oceananigans uses a hierarchical structure to express boundary conditions:
+
+1. Each boundary has one [`BoundaryCondition`](@ref)
+2. Each field has seven [`BoundaryCondition`](@ref) (`west`, `east`, `south`, `north`, `bottom`, `top` and
+   and an additional experimental condition for `immersed` boundaries)
+3. A set of `FieldBoundaryConditions`, up to one for each field, are grouped into a `NamedTuple` and passed
+   to the model constructor.
+
 ## Specifying boundary conditions for a model
 
 Boundary conditions are defined at model construction time by passing a `NamedTuple` of `FieldBoundaryConditions`
-specifying non-default boundary conditions for velocities (``u``, ``v``, ``w``) and tracers.
+specifying non-default boundary conditions for fields such as velocities and tracers.
+
 Fields for which boundary conditions are not specified are assigned a default boundary conditions.
 
 A few illustrations are provided below. See the examples for
@@ -59,7 +75,6 @@ further illustrations of boundary condition specification.
 ```@meta
 DocTestSetup = quote
    using Oceananigans
-   using Oceananigans.BoundaryConditions
 
    using Random
    Random.seed!(1234)
@@ -73,7 +88,7 @@ In this section we illustrate usage of the different [`BoundaryCondition`](@ref)
 
 ```jldoctest
 julia> constant_T_bc = ValueBoundaryCondition(20.0)
-BoundaryCondition: type=Value, condition=20.0
+BoundaryCondition: classification=Value, condition=20.0
 ```
 
 A constant [`Value`](@ref) boundary condition can be used to specify constant tracer (such as temperature),
@@ -90,7 +105,7 @@ julia> ρ₀ = 1027;  # Reference density [kg/m³]
 julia> τₓ = 0.08;  # Wind stress [N/m²]
 
 julia> wind_stress_bc = FluxBoundaryCondition(-τₓ/ρ₀)
-BoundaryCondition: type=Flux, condition=-7.789678675754625e-5
+BoundaryCondition: classification=Flux, condition=-7.789678675754625e-5
 ```
 
 A constant [`Flux`](@ref) boundary condition can be imposed on tracers and tangential velocity components
@@ -115,7 +130,7 @@ Boundary conditions may be specified by functions,
 julia> @inline surface_flux(x, y, t) = cos(2π * x) * cos(t);
 
 julia> top_tracer_bc = FluxBoundaryCondition(surface_flux)
-BoundaryCondition: type=Flux, condition=surface_flux(x, y, t) in Main at none:1
+BoundaryCondition: classification=Flux, condition=surface_flux(x, y, t) in Main at none:1
 ```
 
 !!! info "Boundary condition functions"
@@ -137,8 +152,8 @@ Boundary condition functions may be 'parameterized',
 ```jldoctest
 julia> @inline wind_stress(x, y, t, p) = - p.τ * cos(p.k * x) * cos(p.ω * t); # function with parameters
 
-julia> top_u_bc = BoundaryCondition(Flux, wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4))
-BoundaryCondition: type=Flux, condition=wind_stress(x, y, t, p) in Main at none:1
+julia> top_u_bc = FluxBoundaryCondition(wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4))
+BoundaryCondition: classification=Flux, condition=wind_stress(x, y, t, p) in Main at none:1
 ```
 
 !!! info "Boundary condition functions with parameters"
@@ -157,7 +172,7 @@ julia> @inline linear_drag(x, y, t, u) = - 0.2 * u
 linear_drag (generic function with 1 method)
 
 julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, field_dependencies=:u)
-BoundaryCondition: type=Flux, condition=linear_drag(x, y, t, u) in Main at none:1
+BoundaryCondition: classification=Flux, condition=linear_drag(x, y, t, u) in Main at none:1
 ```
 
 `field_dependencies` specifies the name of the dependent fields either with a `Symbol` or `Tuple` of `Symbol`s.
@@ -171,7 +186,7 @@ julia> @inline quadratic_drag(x, y, t, u, v, drag_coeff) = - drag_coeff * u * sq
 quadratic_drag (generic function with 1 method)
 
 julia> u_bottom_bc = FluxBoundaryCondition(quadratic_drag, field_dependencies=(:u, :v), parameters=1e-3)
-BoundaryCondition: type=Flux, condition=quadratic_drag(x, y, t, u, v, drag_coeff) in Main at none:1
+BoundaryCondition: classification=Flux, condition=quadratic_drag(x, y, t, u, v, drag_coeff) in Main at none:1
 ```
 
 Put differently, `ξ, η, t` come first in the function signature, followed by field dependencies,
@@ -189,7 +204,7 @@ using the `discrete_form`. For example:
 u_bottom_bc = FluxBoundaryCondition(filtered_drag, discrete_form=true)
 
 # output
-BoundaryCondition: type=Flux, condition=filtered_drag(i, j, grid, clock, model_fields) in Main at none:1
+BoundaryCondition: classification=Flux, condition=filtered_drag(i, j, grid, clock, model_fields) in Main at none:1
 ```
 
 !!! info "The 'discrete form' for boundary condition functions"
@@ -212,8 +227,8 @@ julia> Cd = 0.2;  # drag coefficient
 
 julia> @inline linear_drag(i, j, grid, clock, model_fields, Cd) = @inbounds - Cd * model_fields.u[i, j, 1];
 
-julia> u_bottom_bc = BoundaryCondition(Flux, linear_drag, discrete_form=true, parameters=Cd)
-BoundaryCondition: type=Flux, condition=linear_drag(i, j, grid, clock, model_fields, Cd) in Main at none:1
+julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, discrete_form=true, parameters=Cd)
+BoundaryCondition: classification=Flux, condition=linear_drag(i, j, grid, clock, model_fields, Cd) in Main at none:1
 ```
 
 !!! info "Inlining and avoiding bounds-checking in boundary condition functions"
@@ -231,36 +246,34 @@ julia> Nx = Ny = 16;  # Number of grid points.
 julia> Q = randn(Nx, Ny); # temperature flux
 
 julia> white_noise_T_bc = FluxBoundaryCondition(Q)
-BoundaryCondition: type=Flux, condition=16×16 Matrix{Float64}
+BoundaryCondition: classification=Flux, condition=16×16 Matrix{Float64}
 ```
 
 When running on the GPU, `Q` must be converted to a `CuArray`.
 
 ## Building boundary conditions on a field
 
-To create, for example, a set of horizontally-periodic field boundary conditions, write
+To create a set of [`FieldBoundaryConditions`](@ref) for a temperature field,
+we write
 
 ```jldoctest
-julia> topology = (Periodic, Periodic, Bounded);
-
-julia> grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology);
-
-julia> T_bcs = TracerBoundaryConditions(grid,    top = ValueBoundaryCondition(20),
-                                              bottom = GradientBoundaryCondition(0.01))
-Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
-├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}
-├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}
-└── z: CoordinateBoundaryConditions{BoundaryCondition{Gradient, Float64}, BoundaryCondition{Value, Int64}}
+julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20),
+                                       bottom = GradientBoundaryCondition(0.01))
+Oceananigans.FieldBoundaryConditions, with boundary conditions
+├── west: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
+├── east: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
+├── south: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
+├── north: Oceananigans.BoundaryConditions.DefaultPrognosticFieldBoundaryCondition
+├── bottom: BoundaryCondition{Gradient, Float64}
+├── top: BoundaryCondition{Value, Int64}
+└── immersed: BoundaryCondition{Flux, Nothing}
 ```
 
-`T_bcs` is a [`FieldBoundaryConditions`](@ref) object for temperature `T` appropriate
-for horizontally periodic grid topologies.
-The default `Periodic` boundary conditions in ``x`` and ``y`` are inferred from the `topology` of `grid`.
+If the grid is, e.g., horizontally-periodic, then each horizontal `DefaultPrognosticFieldBoundaryCondition`
+is converted to `PeriodicBoundaryCondition` inside the model's constructor, before assigning the
+boundary conditions to temperature `T`.
 
-For ``u``, ``v``, and ``w``, use the
-`UVelocityBoundaryConditions`
-`VVelocityBoundaryConditions`, and
-`WVelocityBoundaryConditions` constructors, respectively.
+In general, boundary condition defaults are inferred from the field location and `topology(grid)`.
 
 ## Specifying model boundary conditions
 
@@ -274,11 +287,11 @@ julia> topology = (Periodic, Periodic, Bounded);
 
 julia> grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology);
 
-julia> u_bcs = UVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(+0.1),
-                                              bottom = ValueBoundaryCondition(-0.1));
+julia> u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(+0.1),
+                                       bottom = ValueBoundaryCondition(-0.1));
 
-julia> T_bcs = TracerBoundaryConditions(grid, top = ValueBoundaryCondition(20),
-                                           bottom = GradientBoundaryCondition(0.01));
+julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20),
+                                       bottom = GradientBoundaryCondition(0.01));
 
 julia> model = IncompressibleModel(grid=grid, boundary_conditions=(u=u_bcs, T=T_bcs))
 IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0)
@@ -292,13 +305,13 @@ julia> model.velocities.u
 Field located at (Face, Center, Center)
 ├── 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)
+└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Value, top=Value, immersed=ZeroFlux
 
 julia> model.tracers.T
 Field located at (Center, Center, Center)
 ├── 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)
+└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=Gradient, top=Value, immersed=ZeroFlux
 ```
 
 Notice that the specified non-default boundary conditions have been applied at

docs/src/model_setup/tracers.md

@@ -28,14 +28,13 @@ julia> model.tracers.T
 Field located at (Center, Center, Center)
 ├── 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)
+└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux
 
 julia> model.tracers.S
 Field located at (Center, Center, Center)
 ├── 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)
-
+└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux
 ```
 
 Any number of arbitrary tracers can be appended to this list and passed to a model constructor. For example, to evolve

docs/src/simulation_tips.md

@@ -73,19 +73,26 @@ in GPU simulations by following a few simple principles.
 ### Variables that need to be used in GPU computations need to be defined as constants
 
 Any global variable that needs to be accessed by the GPU needs to be a constant or the simulation
-will crash. This includes any variables used in forcing functions and boundary conditions. For
-example, if you define a boundary condition like the example below and run your simulation on a GPU
-you'll get an error.
+will crash. This includes any variables that are referenced as global variables in functions
+used for forcing of boundary conditions. For example,
 
 ```julia
-dTdz = 0.01 # K m⁻¹
-T_bcs = TracerBoundaryConditions(grid,
-                                 bottom = GradientBoundaryCondition(dTdz))
+T₀ = 20 # ᵒC
+surface_temperature(x, y, t) = T₀ * sin(2π / 86400 * t)
+T_bcs = FieldBoundaryConditions(bottom = GradientBoundaryCondition(surface_temperature))
 ```
 
-However, if you define `dTdz` as a constant by replacing the first line with `const dTdz = 0.01`,
-then (provided everything else is done properly) your run will be successful.
+will throw an error if run on the GPU (and will run more slowly than it should on the CPU).
+Replacing the first line above with
 
+```julia
+const T₀ = 20 # ᵒC
+```
+
+fixes the issue by indicating to the compiler that `T₀` will not change.
+
+Note that the _literal_ `2π / 86400` is not an issue -- it's only the
+_variable_ `T₀` that must be declared `const`.
 
 ### Complex diagnostics using `ComputedField`s may not work on GPUs
 

examples/convecting_plankton.jl

@@ -12,7 +12,7 @@
 # The phytoplankton in our model are advected, diffuse, grow, and die according to
 #
 # ```math
-# ∂_t P + \boldsymbol{u ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, ,
+# ∂_t P + \boldsymbol{v ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, ,
 # ```
 #
 # where ``\boldsymbol{v}`` is the turbulent velocity field, ``κ`` is an isotropic diffusivity,
@@ -92,7 +92,7 @@ buoyancy_gradient_bc = GradientBoundaryCondition(N²)
 # In summary, the buoyancy boundary conditions impose a destabilizing flux
 # at the top and a stable buoyancy gradient at the bottom:
 
-buoyancy_bcs = TracerBoundaryConditions(grid, top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc)
+buoyancy_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc)
 
 # ## Phytoplankton dynamics: light-dependent growth and uniform mortality
 #

examples/eady_turbulence.jl

@@ -167,8 +167,8 @@ cᴰ = 1e-4 # quadratic drag coefficient
 drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
 drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=cᴰ)
 
-u_bcs = UVelocityBoundaryConditions(grid, bottom = drag_bc_u)
-v_bcs = VVelocityBoundaryConditions(grid, bottom = drag_bc_v)
+u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
+v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)
 
 # ## Turbulence closures
 #