Numerical implementation of boundary conditions

We adopt a mixed approach for implementing boundary conditions that uses both halo regions and "direct" imposition of boundary conditions, depending on the condition prescribed.

We illustrate how boundary conditions are implemented by considering the tracer equation

where is the diffusive flux of and is an arbitrary source term.

See Model setup: boundary conditions for how to create and use these boundary conditions in Oceananigans.

Gradient boundary conditions

Users impose gradient boundary conditions by prescribing the gradient of a field across an external boundary . The prescribed gradient may be a constant, discrete array of values, or an arbitrary function. The gradient boundary condition is enforced by setting the value of halo points located outside the domain interior such that

where is the vector normal to .

Across the bottom boundary in , for example, this requires that

where   are the heights of the finite volume at and   and  . This prescription implies that the -derivative of across the boundary at   is

as prescribed by the user.

Gradient boundary conditions are represented by the Gradient type.

Value boundary conditions

Users impose value boundary conditions by prescribing , the value of on the external boundary . The value may be a constant, array of discrete values, or an arbitrary function. To enforce a value boundary condition, the gradient associated with the difference between and at boundary-adjacent nodes is diagnosed and used to set the value of the halo point located outside the boundary.

At the bottom boundary in , for example, this means that the gradient of is determined by

which is then used to set the halo point via linear extrapolation.

Value boundary conditions are represented by the Value type.

Flux boundary conditions

Users impose flux boundary conditions by prescribing the flux   of across the external boundary . The flux   may be a constant, array of discrete values, or an arbitrary function. To explain how flux boundary conditions are imposed in Oceananigans.jl, we note that the average of the tracer conservation equation over a finite volume yields

where the surface integral over averages the flux of across the six faces of the finite volume. The right-hand-side of above is denoted as .

An external boundary of a finite volume is associated with a no-penetration condition such that    , where is the vector normal to . Furthermore, the closures currently available in Oceananigans.jl have the property that  . Thus setting     on the external boundary implies that the total flux of across the external boundary is

Oceananigans.jl exploits this fact to define an algorithm that prescribes fluxes across external boundaries :

1. Impose a constant gradient     across external boundaries by using halo points (similar to ), which ensures that the evaluation of in boundary-adjacent cells does not include fluxes across the external boundary, and;

2. Add the prescribed flux to the boundary-adjacent volumes prior to calculating :       , where   denotes values of in boundary-adjacent volumes,   is the flux prescribed along the boundary, is the volume of the boundary-adjacent cell, and is the area of the external boundary of the boundary-adjacent cell. The factor is $-$1 and $+$1 on "left" and "right" boundaries, and accounts for the fact that a positive flux on a left boundary where    implies an "inward" flux of that increases interior values of , whereas a positive flux on a right boundary where   implies an "outward" flux that decreases interior values of .

Flux boundary conditions are represented by the Flux type.

Open boundary conditions

Open boundary conditions directly specify the value of the halo points. Typically this is used to impose no penetration boundary conditions, i.e. setting wall normal velocity components to zero on the boundary.

The nuance here is that open boundaries behave differently for fields on face points in the boundary direction due to the staggered grid. For example, the u-component of velocity lies on (Face, Center, Center) points so for open west or east boundaries the point specified by the boundary condition is the point lying on the boundary, whereas for a tracer on (Center, Center, Center) points the open boundary condition specifies a point outside of the domain (hence the difference with Value boundary conditions).

The other important detail is that open (including no-penetration) boundary conditions are the only conditions used on wall normal velocities when the domain is not periodic. This means that their value affects the pressure calculation for nonhydrostatic models as it is involved in calculating the divergence in the boundary adjacent center point (as described in the fractional step method documentation). Usually boundary points are filled for the predictor velocity (i.e. before the pressure is calculated), and on the corrected field (i.e. after the pressure correction is applied), but for open boundaries this would result in the boundary adjacent center point becoming divergent so open boundaries are only filled for the predictor velocity and stay the same after the pressure correction (so the boundary point is filled with the final corrected velocity at the predictor step).

The restriction arises as the boundary condition is specifying the wall normal velocity,  , which leads to the pressure boundary condition

implying that there is a pressure gradient across the boundary. Since we solve the pressure poisson equation ( ) using the method described by Schumann and Sweet (1988) we have to move inhomogeneous boundary conditions on the pressure to the right hand side. In order to do this we define a new field where

This moves the boundary condition to the right hand side as becomes

Given the boundary condition on pressure given above, we can define a new modified predictor velocity which is equal to the predictor velocity within the domain but shares boundary conditions with the corrected field,

The modified pressure Poisson equation becomes   which can easily be solved.

Perhaps a more intuitive way to consider this is to recall that the corrector step projects to the space of divergence free velocity by applying

but we have changed to and to so for   the modified predictor velocity must equal the corrected velocity on the boundary.

For simple open boundary conditions such as no penetration or a straight forward prescription of a known velocity at this is simple to implement as we just set the boundary condition on the predictor velocity and don't change it after the correction. But some open boundary methods calculate the boundary value based on the interior solution. As a simple example, if we wanted to set the wall normal velocity gradient to zero at the west boundary then we would set the boundary point to

but we then pressure correct the interior so a new error is introduced as

This is preferred to a divergent interior solution as open boundary conditions (except no penetration) are typically already unphysical and only used in an attempt to allow information to enter or exit the domain.

Open boundary conditions are represented by the Open type.

Open boundary condition "schemes"

Except for trivial cases (i.e. no-penetration) the velocity on the boundary point has to be approximated as it is outside the computed domain. There is insufficient information to step the full equation of motion as gradients across the boundary cannot be computed and simply prescribing a boundary normal velocity is unphysical and reflects energy leaving the domain (Orlanski, 1976).

Perturbation advection

A common method for allowing information to exit is to perform a one-sided advection operation. For example, in the Orlanski (1976) boundary condition fields are advected out of the domain by a locally determined phase speed. We can show that this is the first-order approximation of the full equations of motion in the predictor velocity step. Consider a right boundary normal to the u velocity component (the east boundary):

let    with   where is an externally determined "background" wall-normal flow in the proximity of the boundary, and assume that the stress tensor gradient is small,

then, taking only first-order terms:

Now consider the dominant forcing to be relaxation to the background state (we will explain why later) so that   , and recall that when we compute the boundary normal point the interior domain is already stepped, we can take an upwind backward Euler step:

where is the boundary point and and are interpolated to the boundary point. Of course if   (i.e. it is directed into the domain), then we need which is outside of the domain, so this scheme cannot be used. We therefore disregard this term, taking   and  , and rearranging we recover:

For numerical stability, we also need to apply a CFL-like constraint and set  .

Previously we considered the dominant forcing to be relaxation to so that we can take the backward Euler step. And because it is useful to prescribe relaxation on the boundary both to damp numerical oscillations, and to prevent inconsistency when is directed into the domain.