(0.91.7) Open boundary conditions for NonhydrostaticModel

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.91.6"
+version = "0.91.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/model_setup/boundary_conditions.md

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

docs/src/numerical_implementation/boundary_conditions.md

@@ -123,3 +123,104 @@ boundaries ``\partial \Omega_b``:
    values of ``c``.
 
 Flux boundary conditions are represented by the [`Flux`](@ref) 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 on 
+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](@ref finite_volume). 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, where as 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](@ref time_stepping) 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 arrises as the boundary condition is specifying the wall normal velocity, 
+``\hat{\boldsymbol{n}}\cdot\boldsymbol{u}``, which leads to the pressure boundary condition
+```math
+    \begin{equation}
+    \label{eq:pressure_boundary_condition}
+    \Delta t \, \hat{\boldsymbol{n}}\cdot\boldsymbol{\nabla}p^{n+1}\big |_{\partial\Omega} = \left[\Delta t \, \hat{\boldsymbol{n}}\cdot\boldsymbol{u}^\star - \hat{\boldsymbol{n}}\cdot\boldsymbol{u}^{n+1}\right],
+    \end{equation}
+```
+implying that there is a pressure gradient across the boundary. Since we solve the pressure poisson 
+equation (``\nabla^2p^{n+1}=\frac{\boldsymbol{\nabla}\cdot\boldsymbol{u}^\star}{\Delta t}``)
+using the method described by [Schumann88](@citet) we have to move inhomogeneus boundary conditions
+on the pressure to the right hand side. In order to do this we define a new field ``\phi`` where
+```math
+    \begin{equation}
+    \label{eq:modified_pressure_field}
+    \phi = p^{n+1} \quad \text{inside} \quad \Omega \quad \text{but} \quad \boldsymbol{\nabla} \cdot \boldsymbol{\nabla} \phi \, \big |_{\partial\Omega} = 0.
+    \end{equation}
+```
+This moves the boundary condition to the right hand side as ``\phi`` becomes
+```math
+    \begin{equation}
+    \label{eq:modified_pressure_poisson}
+    \boldsymbol{\nabla}^2\phi^{n+1} = \boldsymbol{\nabla}\cdot\left[\frac{\boldsymbol{u}^\star}{\Delta t} - \delta\left(\boldsymbol{x} - \boldsymbol{x}_\Omega\right)\boldsymbol{\nabla}p\right].
+    \end{equation}
+```
+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,
+```math
+    \begin{equation}
+    \label{eq:quasi_predictor_velocity}
+    \tilde{\boldsymbol{u}}^\star:=\boldsymbol{u}^\star + \delta\left(\boldsymbol{x} - \boldsymbol{x}_\Omega\right)(\boldsymbol{u}^{n+1} - \boldsymbol{u}^\star).
+    \end{equation}
+```
+The modified pressure poisson equation becomes ``\nabla^2p^{n+1}=\frac{\boldsymbol{\nabla}\cdot\tilde{\boldsymbol{u}}^\star}{\Delta t}``
+which can easily be solved. 
+
+Perhaps a more intuitive way to consider this is to recall that the corrector step projects ``\boldsymbol{u}^\star``
+to the space of divergenece free velocity by applying
+```math
+    \begin{equation}
+    \label{eq:pressure_correction_step}
+    \boldsymbol{u}^{n+1} = \boldsymbol{u}^\star - \Delta t\boldsymbol{\nabla}p^{n+1},
+    \end{equation}
+```
+but we have changed ``p^{n+1}`` to ``\phi`` and ``\boldsymbol{u}^\star`` to ``\tilde{\boldsymbol{u}}^\star``
+so for ``\boldsymbol{\nabla}\phi \big |_{\partial\Omega} = 0`` 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 ``t^{n+1}`` 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 veloicty gradient to zero at the west boundary then we would set the boundary 
+point to
+```math
+    \begin{equation}
+    \label{eq:zero_wall_normal_velocity_gradient}
+    u^\star_{1jk} \approx u^\star_{3jk} + (u^\star_{2jk} - u^\star_{jk4}) / 2 + \mathcal{O}(\Delta x^2),
+    \end{equation}
+```
+but we then pressure correct the interior so a new ``\mathcal{O}(\Delta t)`` error is introduced as
+```math
+    \begin{equation}
+    \begin{align}
+    u^{n+1}_{1jk} &\approx u^{n+1}_{3jk} + (u^{n+1}_{2jk} - u^{n+1}_{jk4}) / 2 + \mathcal{O}(\Delta x^2),\\
+    &= u^\star_{1jk} - \Delta t \left(\boldsymbol{\nabla}p^{n+1}_{3jk} + (\boldsymbol{\nabla}p^{n+1}_{2jk} - \boldsymbol{\nabla}p^{n+1}_{4jk}) / 2\right) + \mathcal{O}(\Delta x^2),\\
+    &\approx u^\star_{1jk} + \mathcal{O}(\Delta x^2) + \mathcal{O}(\Delta t).
+    \end{align}
+    \end{equation}
+```
+This is prefered to a divergent interior solution as open boundary conditions (except no penetration)
+are typlically 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`](@ref) type.

docs/src/numerical_implementation/finite_volume.md

@@ -1,4 +1,4 @@
-# Finite volume method on a staggered grid
+# [Finite volume method on a staggered grid](@id finite_volume)
 
 The `Oceananigans.jl` staggered grid is defined by a rectilinear array of cuboids of horizontal dimensions 
 ``\Delta x_{i, j, k}, \Delta y_{i, j, k}`` and vertical dimension 

src/BoundaryConditions/BoundaryConditions.jl

@@ -10,14 +10,16 @@ export
     apply_x_bcs!, apply_y_bcs!, apply_z_bcs!,
     fill_halo_regions!
 
-using CUDA
+using CUDA, Adapt
 using KernelAbstractions: @index, @kernel
 
 using Oceananigans.Architectures: CPU, GPU, device
 using Oceananigans.Utils: work_layout, launch!
 using Oceananigans.Operators: Ax, Ay, Az, volume
 using Oceananigans.Grids
 
+import Adapt: adapt_structure
+
 include("boundary_condition_classifications.jl")
 include("boundary_condition.jl")
 include("discrete_boundary_function.jl")
@@ -34,4 +36,7 @@ include("fill_halo_regions_nothing.jl")
 
 include("apply_flux_bcs.jl")
 
+include("update_boundary_conditions.jl")
+
+include("flat_extrapolation_open_boundary_matching_scheme.jl")
 end # module

src/BoundaryConditions/boundary_condition.jl

@@ -1,4 +1,3 @@
-import Adapt
 import Oceananigans.Architectures: on_architecture
 
 """
@@ -12,21 +11,12 @@ struct BoundaryCondition{C<:AbstractBoundaryConditionClassification, T}
 end
 
 """
-    BoundaryCondition(Classification::DataType, condition)
-
-Construct a boundary condition of type `BC` with a number or array as a `condition`.
-
-Boundary condition types include `Periodic`, `Flux`, `Value`, `Gradient`, and `Open`.
-"""
-BoundaryCondition(Classification::DataType, condition) = BoundaryCondition(Classification(), condition)
-
-"""
-    BoundaryCondition(Classification::DataType, condition::Function;
+    BoundaryCondition(classification::AbstractBoundaryConditionClassification, condition::Function;
                       parameters = nothing,
                       discrete_form = false,
                       field_dependencies=())
 
-Construct a boundary condition of type `Classification` with a function boundary `condition`.
+Construct a boundary condition of type `classification` with a function boundary `condition`.
 
 By default, the function boudnary `condition` is assumed to have the 'continuous form'
 `condition(ξ, η, t)`, where `t` is time and `ξ` and `η` vary along the boundary.
@@ -50,7 +40,7 @@ where `i`, and `j` are indices that vary along the boundary. If `discrete_form =
 condition(i, j, grid, clock, model_fields, parameters)
 ```
 """
-function BoundaryCondition(Classification::DataType, condition::Function;
+function BoundaryCondition(classification::AbstractBoundaryConditionClassification, condition::Function;
                            parameters = nothing,
                            discrete_form = false,
                            field_dependencies=())
@@ -63,17 +53,16 @@ function BoundaryCondition(Classification::DataType, condition::Function;
         condition = ContinuousBoundaryFunction(condition, parameters, field_dependencies)
     end
 
-    return BoundaryCondition(Classification(), condition)
+    return BoundaryCondition(classification, condition)
 end
 
 # Adapt boundary condition struct to be GPU friendly and passable to GPU kernels.
-Adapt.adapt_structure(to, b::BoundaryCondition{Classification}) where Classification =
-    BoundaryCondition(Classification(), Adapt.adapt(to, b.condition))
-
+Adapt.adapt_structure(to, b::BoundaryCondition) =
+    BoundaryCondition(Adapt.adapt(to, b.classification), Adapt.adapt(to, b.condition))
 
 # Adapt boundary condition struct to be GPU friendly and passable to GPU kernels.
-on_architecture(to, b::BoundaryCondition{Classification}) where Classification =
-    BoundaryCondition(Classification(), on_architecture(to, b.condition))
+on_architecture(to, b::BoundaryCondition) =
+    BoundaryCondition(on_architecture(to, b.classification), on_architecture(to, b.condition))
 
 #####
 ##### Some abbreviations to make life easier.
@@ -91,18 +80,18 @@ const MCBC = BoundaryCondition{<:MultiRegionCommunication}
 const DCBC = BoundaryCondition{<:DistributedCommunication}
 
 # More readable BC constructors for the public API.
-                PeriodicBoundaryCondition() = BoundaryCondition(Periodic,                 nothing)
-                  NoFluxBoundaryCondition() = BoundaryCondition(Flux,                     nothing)
-            ImpenetrableBoundaryCondition() = BoundaryCondition(Open,                     nothing)
-MultiRegionCommunicationBoundaryCondition() = BoundaryCondition(MultiRegionCommunication, nothing)
-DistributedCommunicationBoundaryCondition() = BoundaryCondition(DistributedCommunication, nothing)
-
-                    FluxBoundaryCondition(val; kwargs...) = BoundaryCondition(Flux, val; kwargs...)
-                   ValueBoundaryCondition(val; kwargs...) = BoundaryCondition(Value, val; kwargs...)
-                GradientBoundaryCondition(val; kwargs...) = BoundaryCondition(Gradient, val; kwargs...)
-                    OpenBoundaryCondition(val; kwargs...) = BoundaryCondition(Open, val; kwargs...)
-MultiRegionCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(MultiRegionCommunication, val; kwargs...)
-DistributedCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(DistributedCommunication, val; kwargs...)
+                PeriodicBoundaryCondition() = BoundaryCondition(Periodic(),                 nothing)
+                  NoFluxBoundaryCondition() = BoundaryCondition(Flux(),                     nothing)
+            ImpenetrableBoundaryCondition() = BoundaryCondition(Open(), nothing)
+MultiRegionCommunicationBoundaryCondition() = BoundaryCondition(MultiRegionCommunication(), nothing)
+DistributedCommunicationBoundaryCondition() = BoundaryCondition(DistributedCommunication(), nothing)
+
+                    FluxBoundaryCondition(val; kwargs...) = BoundaryCondition(Flux(), val; kwargs...)
+                   ValueBoundaryCondition(val; kwargs...) = BoundaryCondition(Value(), val; kwargs...)
+                GradientBoundaryCondition(val; kwargs...) = BoundaryCondition(Gradient(), val; kwargs...)
+                    OpenBoundaryCondition(val; kwargs...) = BoundaryCondition(Open(nothing), val; kwargs...)
+MultiRegionCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(MultiRegionCommunication(), val; kwargs...)
+DistributedCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(DistributedCommunication(), val; kwargs...)
 
 # Support for various types of boundary conditions.
 #
@@ -125,9 +114,6 @@ DistributedCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(Di
 const NumberRef = Base.RefValue{<:Number}
 @inline getbc(bc::BC{<:Any, <:NumberRef}, args...) = bc.condition[]
 
-Adapt.adapt_structure(to, bc::BoundaryCondition) = BoundaryCondition(Adapt.adapt(to, bc.classification),
-                                                                     Adapt.adapt(to, bc.condition))
-
 #####
 ##### Validation with topology
 #####

src/BoundaryConditions/boundary_condition_classifications.jl

@@ -59,7 +59,15 @@ For fields located at `Faces`, `Open` also specifies field value _on_ the bounda
 Open boundary conditions are used to specify the component of a velocity field normal to a boundary
 and can also be used to describe nested or linked simulation domains.
 """
-struct Open <: AbstractBoundaryConditionClassification end
+struct Open{MS} <: AbstractBoundaryConditionClassification
+    matching_scheme::MS
+end
+
+Open() = Open(nothing)
+
+(open::Open)() = open
+
+Adapt.adapt_structure(to, open::Open) = Open(adapt(to, open.matching_scheme))
 
 """
     struct MultiRegionCommunication <: AbstractBoundaryConditionClassification

src/BoundaryConditions/continuous_boundary_function.jl

@@ -88,7 +88,7 @@ function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoun
                                                                              boundary_func.field_dependencies,
                                                                              indices, interps)
 
-    return BoundaryCondition(C, regularized_boundary_func)
+    return BoundaryCondition(bc.classification, regularized_boundary_func)
 end
 
 @inline domain_boundary_indices(::LeftBoundary, N) = 1, 1

src/BoundaryConditions/field_boundary_conditions.jl

@@ -57,7 +57,8 @@ function FieldBoundaryConditions(indices::Tuple, west, east, south, north, botto
 end
 
 FieldBoundaryConditions(indices::Tuple, bcs::FieldBoundaryConditions) =
-    FieldBoundaryConditions(indices, (getproperty(bcs, side) for side in fieldnames(FieldBoundaryConditions))...)
+    FieldBoundaryConditions(indices, (getproperty(bcs, side) for side in propertynames(bcs))...)
+
 
 window_boundary_conditions(::Colon,     left, right) = left, right
 window_boundary_conditions(::UnitRange, left, right) = nothing, nothing

src/BoundaryConditions/fill_halo_regions.jl

@@ -45,10 +45,14 @@ end
 const MaybeTupledData = Union{OffsetArray, NTuple{<:Any, OffsetArray}}
 
 "Fill halo regions in ``x``, ``y``, and ``z`` for a given field's data."
-function fill_halo_regions!(c::MaybeTupledData, boundary_conditions, indices, loc, grid, args...; kwargs...)
-
+function fill_halo_regions!(c::MaybeTupledData, boundary_conditions, indices, loc, grid, args...; 
+                            fill_boundary_normal_velocities = true, kwargs...)
     arch = architecture(grid)
 
+    if fill_boundary_normal_velocities
+        fill_open_boundary_regions!(c, boundary_conditions, indices, loc, grid, args...; kwargs...)
+    end
+
     fill_halos!, bcs = permute_boundary_conditions(boundary_conditions)
     number_of_tasks  = length(fill_halos!)
 
@@ -208,7 +212,6 @@ end
     _fill_bottom_halo!(i, j, grid, c, bottom_bc, loc, args...)
        _fill_top_halo!(i, j, grid, c, top_bc,    loc, args...)
 end
-
 #####
 ##### Single-sided fill_halo! kernels
 #####