GitHub

Boundary conditions

Boundary conditions are intimately related to the grid topology, and only need to be considered in directions with Bounded topology or across immersed boundaries. In Bounded directions, tracer and momentum fluxes are conservative or "zero flux" by default. Non-default boundary conditions are therefore required to specify non-zero fluxes of tracers and momentum across Bounded directions, and across immersed boundaries when using ImmersedBoundaryGrid.

See Numerical implementation of boundary conditions for more details.

Example: no-slip conditions on every boundary

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(16, 16, 16), x=(0, 2π), y=(0, 1), z=(0, 1), topology=(Periodic, Bounded, Bounded))
16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 6.28319) regularly spaced with Δx=0.392699
├── Bounded  y ∈ [0.0, 1.0]     regularly spaced with Δy=0.0625
└── Bounded  z ∈ [0.0, 1.0]     regularly spaced with Δz=0.0625

julia> no_slip_bc = ValueBoundaryCondition(0.0)
ValueBoundaryCondition: 0.0

A "no-slip" BoundaryCondition specifies that velocity components tangential to Bounded directions decay to 0 at the boundary, leading to a viscous loss of momentum.

julia> no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc);

julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=no_slip_field_bcs, v=no_slip_field_bcs, w=no_slip_field_bcs))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: Nothing

julia> model.velocities.u.boundary_conditions
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: ValueBoundaryCondition: 0.0
├── north: ValueBoundaryCondition: 0.0
├── bottom: ValueBoundaryCondition: 0.0
├── top: ValueBoundaryCondition: 0.0
└── immersed: FluxBoundaryCondition: Nothing

Boundary conditions are passed to FieldBoundaryCondition to build boundary conditions for each field individually, and then onto the model constructor (here NonhydrotaticModel) via the keyword argument boundary_conditions. The model constructor then "interprets" the input and builds appropriate boundary conditions for the grid topology, given the user-specified no_slip default boundary condition for Bounded directions. In the above example, note that the west and east boundary conditions are PeriodicBoundaryCondition because the x-topology of the grid is Periodic.

Example: specifying boundary conditions on individual boundaries

To specify no-slip boundary conditions on every Bounded direction except the surface, we write

julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing));

julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs));

julia> model.velocities.u.boundary_conditions
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: ValueBoundaryCondition: 0.0
├── north: ValueBoundaryCondition: 0.0
├── bottom: ValueBoundaryCondition: 0.0
├── top: FluxBoundaryCondition: Nothing
└── immersed: FluxBoundaryCondition: Nothing

julia> model.velocities.v.boundary_conditions
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: OpenBoundaryCondition{Nothing}: Nothing
├── north: OpenBoundaryCondition{Nothing}: Nothing
├── bottom: ValueBoundaryCondition: 0.0
├── top: FluxBoundaryCondition: Nothing
└── immersed: FluxBoundaryCondition: Nothing

Now both u and v have FluxBoundaryCondition(nothing) at the top boundary, which is Oceananigans lingo for "no-flux boundary condition".

Boundary condition classifications

There are three primary boundary condition classifications:

  1. FluxBoundaryCondition specifies fluxes directly.

    Some applications of FluxBoundaryCondition are:

    • surface momentum fluxes due to wind, or "wind stress";
    • linear or quadratic bottom drag;
    • surface temperature fluxes due to heating or cooling;
    • surface salinity fluxes due to precipitation and evaporation;
    • relaxation boundary conditions that restores a field to some boundary distribution over a given time-scale.

  2. ValueBoundaryCondition (Dirichlet) specifies the value of a field on the given boundary, which when used in combination with a turbulence closure results in a flux across the boundary.

    Note: Do not use ValueBoundaryCondition on a wall-normal velocity component (see the note below about ImpenetrableBoundaryCondition).

    Some applications of ValueBoundaryCondition are:

    • no-slip boundary condition for wall-tangential velocity components via ValueBoundaryCondition(0);
    • surface temperature distribution, where heat fluxes in and out of the domain at a rate controlled by the near-surface temperature gradient and the temperature diffusivity;
    • constant velocity tangential to a boundary as in a driven-cavity flow (for example), where the top boundary is moving. Momentum will flux into the domain do the difference between the top boundary velocity and the interior velocity, and the prescribed viscosity.

  3. GradientBoundaryCondition (Neumann) specifies the gradient of a field on a boundary. For example, if there is a known diffusivity, we can express FluxBoundaryCondition(flux) using GradientBoundaryCondition(-flux / diffusivity) (aka "Neumann" boundary condition).

In addition to these primary boundary conditions, ImpenetrableBoundaryCondition applies to velocity components in wall-normal directions.

`ImpenetrableBoundaryCondition`

ImpenetrableBoundaryCondition is internally enforced for fields created inside the model constructor. As a result, ImpenetrableBoundaryCondition is only used for additional velocity components that are not evolved by a model, such as a velocity component used for (AdvectiveForcing)[@ref].

Finally, note that Periodic boundary conditions are internally enforced for Periodic directions, and DefaultBoundaryConditions may exist before boundary conditions are "materialized" by a model.

Default boundary conditions

The default boundary condition in Bounded directions is no-flux, or FluxBoundaryCondition(nothing). The default boundary condition can be changed by passing a positional argument to FieldBoundaryConditions, as in

julia> no_slip_bc = ValueBoundaryCondition(0.0)
ValueBoundaryCondition: 0.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)
├── top: FluxBoundaryCondition: Nothing
└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0)

Boundary condition structures

Oceananigans uses a hierarchical structure to express boundary conditions:

  1. Each boundary of each field has one BoundaryCondition
  2. Each field has seven BoundaryCondition (west, east, south, north, bottom, top and immersed)
  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 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 further illustrations of boundary condition specification.

Creating individual boundary conditions with BoundaryCondition

Boundary conditions may be specified with constants, functions, or arrays. In this section we illustrate usage of the different BoundaryCondition constructors.

1. Constant Value (Dirchlet) boundary condition

julia> constant_T_bc = ValueBoundaryCondition(20.0)
ValueBoundaryCondition: 20.0

A constant Value boundary condition can be used to specify constant tracer (such as temperature), or a constant tangential velocity component at a boundary. Note that boundary conditions on the normal velocity component must use the Open boundary condition type.

Finally, note that ValueBoundaryCondition(condition) is an alias for BoundaryCondition(Value, condition).

2. Constant Flux boundary condition

julia> ρ₀ = 1027;  # Reference density [kg/m³]

julia> τₓ = 0.08;  # Wind stress [N/m²]

julia> wind_stress_bc = FluxBoundaryCondition(-τₓ/ρ₀)
FluxBoundaryCondition: -7.78968e-5

A constant Flux boundary condition can be imposed on tracers and tangential velocity components that can be used, for example, to specify cooling, heating, evaporation, or wind stress at the ocean surface.

The flux convention in Oceananigans

Oceananigans uses the convention that positive fluxes produce transport in the positive direction (east, north, and up for $x$, $y$, $z$). This means, for example, that a negative flux of momentum or velocity at a top boundary, such as in the above example, produces currents in the positive direction, because it prescribes a downwards flux of momentum into the domain from the top. Likewise, a positive temperature flux at the top boundary causes cooling, because it transports heat upwards, out of the domain. Conversely, a positive flux at a bottom boundary acts to increase the interior values of a quantity.

3. Spatially- and temporally-varying flux

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)
FluxBoundaryCondition: ContinuousBoundaryFunction surface_flux at (Nothing, Nothing, Nothing)
Boundary condition functions

By default, a function boundary condition is called with the signature

f(ξ, η, t)

where t is time and ξ, η are spatial coordinates that vary along the boundary:

  • f(y, z, t) on x-boundaries;
  • f(x, z, t) on y-boundaries;
  • f(x, y, t) on z-boundaries.

Alternative function signatures are specified by keyword arguments to BoundaryCondition, as illustrated in subsequent examples.

4. Spatially- and temporally-varying flux with parameters

Boundary condition functions may be 'parameterized',

julia> @inline wind_stress(x, y, t, p) = - p.τ * cos(p.k * x) * cos(p.ω * t); # function with parameters

julia> top_u_bc = FluxBoundaryCondition(wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4))
FluxBoundaryCondition: ContinuousBoundaryFunction wind_stress at (Nothing, Nothing, Nothing)
Boundary condition functions with parameters

The keyword argument parameters above specifies that wind_stress is called with the signature wind_stress(x, y, t, parameters). In principle, parameters is arbitrary. However, relatively simple objects such as floating point numbers or NamedTuples must be used when running on the GPU.

5. 'Field-dependent' boundary conditions

Boundary conditions may also depend on model fields. For example, a linear drag boundary condition is implemented with

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)
FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing)

field_dependencies specifies the name of the dependent fields either with a Symbol or Tuple of Symbols.

6. 'Field-dependent' boundary conditions with parameters

When boundary conditions depends on fields and parameters, their functions take the form

julia> @inline quadratic_drag(x, y, t, u, v, drag_coeff) = - drag_coeff * u * sqrt(u^2 + v^2)
quadratic_drag (generic function with 1 method)

julia> u_bottom_bc = FluxBoundaryCondition(quadratic_drag, field_dependencies=(:u, :v), parameters=1e-3)
FluxBoundaryCondition: ContinuousBoundaryFunction quadratic_drag at (Nothing, Nothing, Nothing)

Put differently, ξ, η, t come first in the function signature, followed by field dependencies, followed by parameters is !isnothing(parameters).

7. Discrete-form boundary condition with parameters

Discrete field data may also be accessed directly from boundary condition functions using the discrete_form. For example:

@inline filtered_drag(i, j, grid, clock, model_fields) =
    @inbounds - 0.05 * (model_fields.u[i-1, j, 1] + 2 * model_fields.u[i, j, 1] + model_fields.u[i-1, j, 1])

u_bottom_bc = FluxBoundaryCondition(filtered_drag, discrete_form=true)

# output
FluxBoundaryCondition: DiscreteBoundaryFunction with filtered_drag
The 'discrete form' for boundary condition functions

The argument discrete_form=true indicates to BoundaryCondition that filtered_drag uses the 'discrete form'. Boundary condition functions that use the 'discrete form' are called with the signature

f(i, j, grid, clock, model_fields)

where i, j are grid indices that vary along the boundary, grid is model.grid, clock is the model.clock, and model_fields is a NamedTuple containing u, v, w and the fields in model.tracers. The signature is similar for $x$ and $y$ boundary conditions expect that i, j is replaced with j, k and i, k respectively.

8. Discrete-form boundary condition with parameters

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 = FluxBoundaryCondition(linear_drag, discrete_form=true, parameters=Cd)
FluxBoundaryCondition: DiscreteBoundaryFunction linear_drag with parameters 0.2
Inlining and avoiding bounds-checking in boundary condition functions

Boundary condition functions should be decorated with @inline when running on CPUs for performance reasons. On the GPU, all functions are force-inlined by default. In addition, the annotation @inbounds should be used when accessing the elements of an array in a boundary condition function (such as model_fields.u[i, j, 1] in the above example). Using @inbounds will avoid a relatively expensive check that the index i, j, 1 is 'in bounds'.

9. A random, spatially-varying, constant-in-time temperature flux specified by an array

julia> Nx = Ny = 16;  # Number of grid points.

julia> Q = randn(Nx, Ny); # temperature flux

julia> white_noise_T_bc = FluxBoundaryCondition(Q)
FluxBoundaryCondition: 16×16 Matrix{Float64}

When running on the GPU, Q must be converted to a CuArray.

Building boundary conditions on a field

To create a set of FieldBoundaryConditions for a temperature field, we write

julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0),
                                       bottom = GradientBoundaryCondition(0.01))
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: GradientBoundaryCondition: 0.01
├── top: ValueBoundaryCondition: 20.0
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

If the grid is, e.g., horizontally-periodic, then each horizontal DefaultBoundaryCondition is converted to PeriodicBoundaryCondition inside the model's constructor, before assigning the boundary conditions to temperature T.

In general, boundary condition defaults are inferred from the field location and topology(grid).

Specifying model boundary conditions

To specify non-default boundary conditions, a named tuple of FieldBoundaryConditions objects is passed to the keyword argument boundary_conditions in the NonhydrostaticModel constructor. The keys of boundary_conditions indicate the field to which the boundary condition is applied. Below, non-default boundary conditions are imposed on the $u$-velocity and tracer $c$.

julia> topology = (Periodic, Periodic, Bounded);

julia> grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology);

julia> u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(+0.1),
                                       bottom = ValueBoundaryCondition(-0.1));

julia> c_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0),
                                       bottom = GradientBoundaryCondition(0.01));

julia> model = NonhydrostaticModel(grid=grid, boundary_conditions=(u=u_bcs, c=c_bcs), tracers=:c)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: c
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: Nothing

julia> model.velocities.u
16×16×16 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Value, top: Value, immersed: ZeroFlux
└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19
    └── max=0.0, min=0.0, mean=0.0

julia> model.tracers.c
16×16×16 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Gradient, top: Value, immersed: ZeroFlux
└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19
    └── max=0.0, min=0.0, mean=0.0

Notice that the specified non-default boundary conditions have been applied at top and bottom of both model.velocities.u and model.tracers.c.

Immersed boundary conditions

Immersed boundary conditions are supported experimentally. A no-slip boundary condition is specified with

# Generate a simple ImmersedBoundaryGrid
hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2)
underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded))
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill))

# Create a no-slip boundary condition for velocity fields.
# Note that the no-slip boundary condition is _only_ applied on immersed boundaries.
velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0))
model = NonhydrostaticModel(; grid, boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs))

# Insepct the boundary condition on the vertical velocity:
model.velocities.w.boundary_conditions.immersed

# output
ImmersedBoundaryCondition:
├── west: ValueBoundaryCondition: 0.0
├── east: ValueBoundaryCondition: 0.0
├── south: ValueBoundaryCondition: 0.0
├── north: ValueBoundaryCondition: 0.0
├── bottom: Nothing
└── top: Nothing
`NonhydrostaticModel` on `ImmersedBoundaryGrid`

The pressure solver for NonhydrostaticModel is approximate, and is unable to produce a velocity field that is simultaneously divergence-free while also satisfying impenetrability on the immersed boundary. As a result, simulated dynamics with NonhydrostaticModel can exhibit egregiously unphysical errors and should be interpreted with caution.

An ImmersedBoundaryCondition encapsulates boundary conditions on each potential boundary-facet of a boundary-adjacent cell. Boundary conditions on specific faces of immersed-boundary-adjacent cells may also be specified by manually building an ImmersedBoundaryCondition:

bottom_drag_bc = ImmersedBoundaryCondition(bottom=ValueBoundaryCondition(0))

# output
ImmersedBoundaryCondition:
├── west: Nothing
├── east: Nothing
├── south: Nothing
├── north: Nothing
├── bottom: ValueBoundaryCondition: 0
└── top: Nothing

The ImmersedBoundaryCondition may then be incorporated into the boundary conditions for a Field by prescribing it to the immersed boundary label,

velocity_bcs = FieldBoundaryConditions(immersed=bottom_drag_bc)

# output
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Value, top=Nothing
`ImmersedBoundaryCondition`

ImmersedBoundaryCondition is experimental. Therefore, one should use it only when a finer level of control over the boundary conditions at the immersed boundary is required, and the user is familiar with the implementation of boundary conditions on staggered grids. For all other cases , using the immersed argument of FieldBoundaryConditions is preferred.

A boundary condition that depends on the fields may be prescribed using the immersed keyword argument in FieldBoundaryConditions. We illustrate field-dependent boundary conditions with an example that imposes linear bottom drag on u on both the bottom facets of cells adjacent to an immersed boundary, and the bottom boundary of the underlying grid.

First we create the boundary condition for the grid's bottom:

@inline linear_drag(x, y, t, u) = - 0.2 * u
drag_u = FluxBoundaryCondition(linear_drag, field_dependencies=:u)

# output
FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing)

Next, we create the immersed boundary condition by adding the argument z to linear_drag and imposing drag only on "bottom" facets of cells that neighbor immersed cells:

@inline immersed_linear_drag(x, y, z, t, u) = - 0.2 * u
immersed_drag_u = FluxBoundaryCondition(immersed_linear_drag, field_dependencies=:u)

u_immersed_bc = ImmersedBoundaryCondition(bottom = immersed_drag_u)

# output
ImmersedBoundaryCondition:
├── west: Nothing
├── east: Nothing
├── south: Nothing
├── north: Nothing
├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction immersed_linear_drag at (Nothing, Nothing, Nothing)
└── top: Nothing

Finally, we combine the two:

u_bcs = FieldBoundaryConditions(bottom = drag_u, immersed = u_immersed_bc)

# output
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Flux, top=Nothing
Positional argument requirements

Note the difference between the arguments required for the function within the bottom boundary condition versus the arguments for the function within the immersed boundary condition. E.g., x, y, t in linear_drag() versus x, y, z, t in immersed_linear_drag().