(0.98.0) Modify interface for OpenBoundaryConditions with schemes and use Tuples in boundary_mass_fluxes

src/BoundaryConditions/BoundaryConditions.jl

@@ -5,6 +5,7 @@ export
     BoundaryCondition, getbc, setbc!,
     PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition, MultiRegionCommunicationBoundaryCondition,
     FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, DistributedCommunicationBoundaryCondition,
+    FlatExtrapolation, PerturbationAdvection,
     validate_boundary_condition_topology, validate_boundary_condition_architecture,
     FieldBoundaryConditions,
     compute_x_bcs!, compute_y_bcs!, compute_z_bcs!,
@@ -40,6 +41,6 @@ include("compute_flux_bcs.jl")
 include("update_boundary_conditions.jl")
 include("polar_boundary_condition.jl")
 
-include("flat_extrapolation_open_boundary_matching_scheme.jl")
-include("perturbation_advection_open_boundary_matching_scheme.jl")
+include("flat_extrapolation_open_boundary_scheme.jl")
+include("perturbation_advection_open_boundary_scheme.jl")
 end # module

src/BoundaryConditions/boundary_condition.jl

@@ -96,7 +96,7 @@ MultiRegionCommunicationBoundaryCondition() = BoundaryCondition(MultiRegionCommu
                     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...)
+  OpenBoundaryCondition(val; scheme = nothing, kwargs...) = BoundaryCondition(Open(scheme), val; kwargs...)
 MultiRegionCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(MultiRegionCommunication(), val; kwargs...)
                   ZipperBoundaryCondition(val; kwargs...) = BoundaryCondition(Zipper(), val; kwargs...)
 DistributedCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(DistributedCommunication(), val; kwargs...)

src/BoundaryConditions/boundary_condition_classifications.jl

@@ -60,14 +60,14 @@ Open boundary conditions are used to specify the component of a velocity field n
 and can also be used to describe nested or linked simulation domains.
 """
 struct Open{MS} <: AbstractBoundaryConditionClassification
-    matching_scheme::MS
+    scheme::MS
 end
 
 Open() = Open(nothing)
 
 (open::Open)() = open
 
-Adapt.adapt_structure(to, open::Open) = Open(adapt(to, open.matching_scheme))
+Adapt.adapt_structure(to, open::Open) = Open(adapt(to, open.scheme))
 
 """
     struct MultiRegionCommunication <: AbstractBoundaryConditionClassification

src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl → src/BoundaryConditions/flat_extrapolation_open_boundary_scheme.jl

@@ -33,21 +33,17 @@ of 1/2 changes to ``Δx₋₁/(Δx₋₂ + Δx₋₃)`` instead, i.e.:
 f(xᵢ) ≈ f(xᵢ₋₂) - (f(xᵢ₋₁) - f(xᵢ₋₃))Δxᵢ₋₁/(Δxᵢ₋₂ + Δxᵢ₋₃) + O(Δx²)
 ```.
 """
-struct FlatExtrapolation{FT}
-    relaxation_timescale :: FT
+@kwdef struct FlatExtrapolation{FT}
+    relaxation_timescale :: FT = Inf
 end
 
-const FEOBC = BoundaryCondition{<:Open{<:FlatExtrapolation}}
-
-function FlatExtrapolationOpenBoundaryCondition(val = nothing; relaxation_timescale = Inf, kwargs...)
-    classification = Open(FlatExtrapolation(relaxation_timescale))
+Adapt.adapt_structure(to, fe::FlatExtrapolation) = FlatExtrapolation(adapt(to, fe.relaxation_timescale))
 
-    return BoundaryCondition(classification, val; kwargs...)
-end
+const FEOBC = BoundaryCondition{<:Open{<:FlatExtrapolation}}
 
 @inline function relax(l, m, grid, ϕ, bc, clock, model_fields)
     Δt = clock.last_stage_Δt
-    τ = bc.classification.matching_scheme.relaxation_timescale
+    τ = bc.classification.scheme.relaxation_timescale
 
     Δt̄ = min(1, Δt / τ)
     ϕₑₓₜ = getbc(bc, l, m, grid, clock, model_fields)
@@ -147,4 +143,4 @@ end
     @inbounds ϕ[i, j, k] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields)
 
     return nothing
-end
+end

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl → src/BoundaryConditions/perturbation_advection_open_boundary_scheme.jl

@@ -4,6 +4,10 @@ using Oceananigans: defaults
 """
     PerturbationAdvection
 
+Create a `PerturbationAdvection` scheme to be used with an `OpenBoundaryCondition`.
+This scheme will nudge the boundary velocity to the OpenBoundaryCondition's exterior value `val`,
+using a time-scale `inflow_timescale` for inflow and `outflow_timescale` for outflow.
+
 For cases where we assume that the internal flow is a small perturbation from
 an external prescribed or coarser flow, we can split the velocity into background
 and perturbation components.
@@ -48,31 +52,18 @@ struct PerturbationAdvection{FT}
    outflow_timescale :: FT
 end
 
-Adapt.adapt_structure(to, pe::PerturbationAdvection) =
-    PerturbationAdvection(adapt(to, pe.inflow_timescale),
-                          adapt(to, pe.outflow_timescale))
-
-"""
-    PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType;
-                                               outflow_timescale = Inf,
-                                               inflow_timescale = 0, kwargs...)
-
-Creates a `PerturbationAdvectionOpenBoundaryCondition` with a given exterior value `val`, to which
-the flow is forced with an `outflow_timescale` for outflow and `inflow_timescale` for inflow. For
-details about this method, refer to the docstring for `PerturbationAdvection`.
-"""
-function PerturbationAdvectionOpenBoundaryCondition(val, FT = defaults.FloatType;
-                                                    outflow_timescale = Inf,
-                                                    inflow_timescale = 0, kwargs...)
+function PerturbationAdvection(FT = defaults.FloatType;
+                               outflow_timescale = Inf,
+                               inflow_timescale = 0)
     inflow_timescale = convert(FT, inflow_timescale)
     outflow_timescale = convert(FT, outflow_timescale)
-    classification = Open(PerturbationAdvection(inflow_timescale, outflow_timescale))
-
-    @warn "`PerturbationAdvection` open boundaries matching scheme is experimental and un-tested/validated"
-
-    return BoundaryCondition(classification, val; kwargs...)
+    return PerturbationAdvection(inflow_timescale, outflow_timescale)
 end
 
+Adapt.adapt_structure(to, pe::PerturbationAdvection) =
+    PerturbationAdvection(adapt(to, pe.inflow_timescale),
+                          adapt(to, pe.outflow_timescale))
+
 const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}}
 
 @inline function step_right_boundary!(bc::PAOBC, l, m, boundary_indices, boundary_adjacent_indices,
@@ -87,7 +78,7 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}}
     uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ)
     U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹))
 
-    pa = bc.classification.matching_scheme
+    pa = bc.classification.scheme
     τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale)
     τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale
 
@@ -111,7 +102,7 @@ end
     uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ)
     U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹))
 
-    pa = bc.classification.matching_scheme
+    pa = bc.classification.scheme
     τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale)
     τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale
 

src/Models/NonhydrostaticModels/boundary_mass_fluxes.jl

@@ -4,7 +4,7 @@ using Oceananigans.Fields: Field, interior, XFaceField, YFaceField, ZFaceField
 using GPUArraysCore: @allowscalar
 
 const OBC  = BoundaryCondition{<:Open} # OpenBoundaryCondition
-const IOBC = BoundaryCondition{<:Open{<:Nothing}} # "Imposed-velocity" OpenBoundaryCondition (with no matching scheme)
+const IOBC = BoundaryCondition{<:Open{<:Nothing}} # "Imposed-velocity" OpenBoundaryCondition (with no scheme)
 const FIOBC = BoundaryCondition{<:Open{<:Nothing}, <:Number} # "Fixed-imposed-velocity" OpenBoundaryCondition
 const ZIOBC = BoundaryCondition{<:Open{<:Nothing}, <:Nothing} # "Zero-imposed-velocity" OpenBoundaryCondition (no-inflow)
 
@@ -79,7 +79,7 @@ initialize_boundary_mass_flux(velocity, ::Nothing, side) = NamedTuple()
 initialize_boundary_mass_flux(velocity, bc, side) = NamedTuple()
 
 needs_mass_flux_correction(::IOBC) = false
-needs_mass_flux_correction(bc::OBC) = bc.classification.matching_scheme !== nothing
+needs_mass_flux_correction(bc::OBC) = bc.classification.scheme !== nothing
 needs_mass_flux_correction(::Nothing) = false
 needs_mass_flux_correction(bc) = false
 
@@ -97,63 +97,63 @@ function initialize_boundary_mass_fluxes(velocities::NamedTuple)
     w_bcs = w.boundary_conditions
 
     boundary_fluxes = NamedTuple()
-    right_matching_scheme_boundaries = Symbol[]
-    left_matching_scheme_boundaries = Symbol[]
-    total_area_matching_scheme_boundaries = zero(eltype(u))
+    right_scheme_boundaries = Symbol[]
+    left_scheme_boundaries = Symbol[]
+    total_area_scheme_boundaries = zero(eltype(u))
 
     # Check west boundary (u velocity)
     west_flux_and_area = initialize_boundary_mass_flux(u, u_bcs.west, Val(:west))
     boundary_fluxes = merge(boundary_fluxes, west_flux_and_area)
     if needs_mass_flux_correction(u_bcs.west)
-        push!(left_matching_scheme_boundaries, :west)
-        total_area_matching_scheme_boundaries += boundary_fluxes.west_area
+        push!(left_scheme_boundaries, :west)
+        total_area_scheme_boundaries += boundary_fluxes.west_area
     end
 
     # Check east boundary (u velocity)
     east_flux_and_area = initialize_boundary_mass_flux(u, u_bcs.east, Val(:east))
     boundary_fluxes = merge(boundary_fluxes, east_flux_and_area)
     if needs_mass_flux_correction(u_bcs.east)
-        push!(right_matching_scheme_boundaries, :east)
-        total_area_matching_scheme_boundaries += boundary_fluxes.east_area
+        push!(right_scheme_boundaries, :east)
+        total_area_scheme_boundaries += boundary_fluxes.east_area
     end
 
     # Check south boundary (v velocity)
     south_flux_and_area = initialize_boundary_mass_flux(v, v_bcs.south, Val(:south))
     boundary_fluxes = merge(boundary_fluxes, south_flux_and_area)
     if needs_mass_flux_correction(v_bcs.south)
-        push!(left_matching_scheme_boundaries, :south)
-        total_area_matching_scheme_boundaries += boundary_fluxes.south_area
+        push!(left_scheme_boundaries, :south)
+        total_area_scheme_boundaries += boundary_fluxes.south_area
     end
 
     # Check north boundary (v velocity)
     north_flux_and_area = initialize_boundary_mass_flux(v, v_bcs.north, Val(:north))
     boundary_fluxes = merge(boundary_fluxes, north_flux_and_area)
     if needs_mass_flux_correction(v_bcs.north)
-        push!(right_matching_scheme_boundaries, :north)
-        total_area_matching_scheme_boundaries += boundary_fluxes.north_area
+        push!(right_scheme_boundaries, :north)
+        total_area_scheme_boundaries += boundary_fluxes.north_area
     end
 
     # Check bottom boundary (w velocity)
     bottom_flux_and_area = initialize_boundary_mass_flux(w, w_bcs.bottom, Val(:bottom))
     boundary_fluxes = merge(boundary_fluxes, bottom_flux_and_area)
     if needs_mass_flux_correction(w_bcs.bottom)
-        push!(left_matching_scheme_boundaries, :bottom)
-        total_area_matching_scheme_boundaries += boundary_fluxes.bottom_area
+        push!(left_scheme_boundaries, :bottom)
+        total_area_scheme_boundaries += boundary_fluxes.bottom_area
     end
 
     # Check top boundary (w velocity)
     top_flux_and_area = initialize_boundary_mass_flux(w, w_bcs.top, Val(:top))
     boundary_fluxes = merge(boundary_fluxes, top_flux_and_area)
     if needs_mass_flux_correction(w_bcs.top)
-        push!(right_matching_scheme_boundaries, :top)
-        total_area_matching_scheme_boundaries += boundary_fluxes.top_area
+        push!(right_scheme_boundaries, :top)
+        total_area_scheme_boundaries += boundary_fluxes.top_area
     end
 
-    boundary_fluxes = merge(boundary_fluxes, (; left_matching_scheme_boundaries,
-                                                right_matching_scheme_boundaries,
-                                                total_area_matching_scheme_boundaries))
+    boundary_fluxes = merge(boundary_fluxes, (; left_scheme_boundaries = Tuple(left_scheme_boundaries),
+                                                right_scheme_boundaries = Tuple(right_scheme_boundaries),
+                                                total_area_scheme_boundaries))
 
-    if length(left_matching_scheme_boundaries) == 0 && length(right_matching_scheme_boundaries) == 0
+    if length(left_scheme_boundaries) == 0 && length(right_scheme_boundaries) == 0
         return nothing
     else
         return boundary_fluxes
@@ -225,7 +225,7 @@ function enforce_open_boundary_mass_conservation!(model, boundary_mass_fluxes)
     u, v, w = model.velocities
 
     ∮udA = open_boundary_mass_inflow(model)
-    A = boundary_mass_fluxes.total_area_matching_scheme_boundaries
+    A = boundary_mass_fluxes.total_area_scheme_boundaries
 
     A⁻¹_∮udA = ∮udA / A
 

src/Oceananigans.jl

@@ -35,6 +35,7 @@ export
     # Boundary conditions
     BoundaryCondition,
     FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition,
+    FlatExtrapolation, PerturbationAdvection,
     FieldBoundaryConditions,
 
     # Fields and field manipulation

test/test_boundary_conditions_integration.jl

@@ -1,8 +1,6 @@
 include("dependencies_for_runtests.jl")
 
 using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction,
-                                       FlatExtrapolationOpenBoundaryCondition,
-                                       PerturbationAdvectionOpenBoundaryCondition
                                        fill_halo_regions!
 
 using Oceananigans: prognostic_fields
@@ -108,22 +106,22 @@ function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
 end
 
 left_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, east = OpenBoundaryCondition(1),
-                                                                    west = FlatExtrapolationOpenBoundaryCondition())
+                                                                    west = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 right_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, west = OpenBoundaryCondition(1),
-                                                                     east = FlatExtrapolationOpenBoundaryCondition())
+                                                                     east = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 left_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, north = OpenBoundaryCondition(1),
-                                                                    south = FlatExtrapolationOpenBoundaryCondition())
+                                                                    south = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 right_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, south = OpenBoundaryCondition(1),
-                                                                     north = FlatExtrapolationOpenBoundaryCondition())
+                                                                     north = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 left_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, top = OpenBoundaryCondition(1),
-                                                                    bottom = FlatExtrapolationOpenBoundaryCondition())
+                                                                    bottom = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 right_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, bottom = OpenBoundaryCondition(1),
-                                                                     top = FlatExtrapolationOpenBoundaryCondition())
+                                                                     top = OpenBoundaryCondition(1, scheme = FlatExtrapolation()))
 
 end_position(::Val{1}, grid) = (grid.Nx+1, 1, 1)
 end_position(::Val{2}, grid) = (1, grid.Ny+1, 1)
@@ -182,7 +180,7 @@ function test_perturbation_advection_open_boundary_conditions(arch, FT)
 
         grid = RectilinearGrid(arch, FT; topology, size = (4, ), x = (0, 4), y = (0, 4), z = (0, 4), halo = (1, ))
 
-        obc = PerturbationAdvectionOpenBoundaryCondition(-1, inflow_timescale = 10.0)
+        obc = OpenBoundaryCondition(-1, scheme = PerturbationAdvection(inflow_timescale = 10.0))
         boundary_conditions = wall_normal_boundary_condition(Val(orientation), obc)
 
         model = NonhydrostaticModel(; grid, boundary_conditions, timestepper = :QuasiAdamsBashforth2)
@@ -195,7 +193,7 @@ function test_perturbation_advection_open_boundary_conditions(arch, FT)
         @test all(view(parent(u), :, :, :) .== -1)
         @test all(interior(u) .== -1)
 
-        obc = PerturbationAdvectionOpenBoundaryCondition((t) -> 0.1*t, inflow_timescale = 0.01, outflow_timescale = 0.5)
+        obc = OpenBoundaryCondition((t) -> 0.1*t, scheme = PerturbationAdvection(inflow_timescale = 0.01, outflow_timescale = 0.5))
         forcing = velocity_forcing(Val(orientation), Forcing((x, t) -> 0.1))
         boundary_conditions = wall_normal_boundary_condition(Val(orientation), obc)
 
@@ -427,27 +425,27 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
             test_flat_extrapolation_open_boundary_conditions(arch, FT)
             test_perturbation_advection_open_boundary_conditions(arch, FT)
 
-            # Only PerturbationAdvectionOpenBoundaryCondition
+            # Only PerturbationAdvection OpenBoundaryCondition
             U₀ = 1
             inflow_timescale = 1e-1
             outflow_timescale = Inf
 
-            u_bcs = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale),
-                                            east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale))
+            u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)),
+                                            east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)))
             boundary_conditions = (; u = u_bcs)
             test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions)
 
-            # Only FlatExtrapolationOpenBoundaryCondition
-            u_bcs = FieldBoundaryConditions(west = FlatExtrapolationOpenBoundaryCondition(),
-                                            east = FlatExtrapolationOpenBoundaryCondition())
+            # Only FlatExtrapolation OpenBoundaryCondition
+            u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation()),
+                                            east = OpenBoundaryCondition(U₀; scheme = FlatExtrapolation()))
             boundary_conditions = (; u = u_bcs)
             test_open_boundary_condition_mass_conservation(arch, FT, boundary_conditions)
 
             # Mixed open boundary conditions
-            u_bcs = FieldBoundaryConditions(west = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale),
-                                            east = PerturbationAdvectionOpenBoundaryCondition(U₀; inflow_timescale, outflow_timescale))
-            v_bcs = FieldBoundaryConditions(south = FlatExtrapolationOpenBoundaryCondition(),
-                                            north = FlatExtrapolationOpenBoundaryCondition())
+            u_bcs = FieldBoundaryConditions(west = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)),
+                                            east = OpenBoundaryCondition(U₀; scheme = PerturbationAdvection(inflow_timescale, outflow_timescale)))
+            v_bcs = FieldBoundaryConditions(south = OpenBoundaryCondition(0; scheme = FlatExtrapolation()),
+                                            north = OpenBoundaryCondition(0; scheme = FlatExtrapolation()))
             w_bcs = FieldBoundaryConditions(bottom = OpenBoundaryCondition(0),
                                             top = OpenBoundaryCondition(0))