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

Project.toml

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

src/BoundaryConditions/BoundaryConditions.jl

@@ -5,6 +5,7 @@ export
     BoundaryCondition, getbc, setbc!,
     PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition, MultiRegionCommunicationBoundaryCondition,
     FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, DistributedCommunicationBoundaryCondition,
+    PerturbationAdvection,
     validate_boundary_condition_topology, validate_boundary_condition_architecture,
     FieldBoundaryConditions,
     compute_x_bcs!, compute_y_bcs!, compute_z_bcs!,
@@ -40,6 +41,5 @@ 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("perturbation_advection.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/perturbation_advection_open_boundary_matching_scheme.jl → src/BoundaryConditions/perturbation_advection.jl

@@ -1,8 +1,19 @@
 using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ, Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ, Az_qᶜᶜᶠ
 using Oceananigans: defaults
 
+struct PerturbationAdvection{FT}
+    inflow_timescale :: FT
+   outflow_timescale :: FT
+end
+
 """
-    PerturbationAdvection
+    PerturbationAdvection(FT = defaults.FloatType;
+                          outflow_timescale = Inf,
+                          inflow_timescale = 0)
+
+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
@@ -43,36 +54,18 @@ to point into the domain.
 The ideal value of the timescales probably depend on the grid spacing and details of the
 boundary flow.
 """
-struct PerturbationAdvection{FT}
-    inflow_timescale :: FT
-   outflow_timescale :: FT
+function PerturbationAdvection(FT = defaults.FloatType;
+                               outflow_timescale = Inf,
+                               inflow_timescale = 0)
+    inflow_timescale = convert(FT, inflow_timescale)
+    outflow_timescale = convert(FT, outflow_timescale)
+    return PerturbationAdvection(inflow_timescale, outflow_timescale)
 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...)
-    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...)
-end
-
 const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}}
 
 @inline function step_right_boundary!(bc::PAOBC, l, m, boundary_indices, boundary_adjacent_indices,
@@ -87,7 +80,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 +104,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

@@ -1,10 +1,10 @@
-using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection, FlatExtrapolation
+using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection
 using Oceananigans.AbstractOperations: Integral, Ax, Ay, Az, GridMetricOperation
 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(::OBC) = true
 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(boundary_fluxes.left_scheme_boundaries) == 0 && length(boundary_fluxes.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