[WIP] Adds a pertubation advection open boundary matching scheme

src/BoundaryConditions/BoundaryConditions.jl

@@ -40,4 +40,5 @@ 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")
 end # module

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl

@@ -0,0 +1,112 @@
+using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ, Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ, Az_qᶜᶜᶠ
+
+"""
+    PerturbationAdvection
+
+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:
+...
+see latex document for now
+
+TODO: check what the coriolis is doing, and check what happens if U is the mean velocity
+"""
+struct PerturbationAdvection{FT}
+    inflow_timescale :: FT
+   outflow_timescale :: FT
+end
+
+Adapt.adapt_structure(to, pe::PerturbationAdvection) = 
+    PerturbationAdvection(adapt(to, pe.outflow_timescale),
+                          adapt(to, pe.inflow_timescale))
+
+function PerturbationAdvectionOpenBoundaryCondition(val, FT = Float64; 
+                                                    outflow_timescale = Inf, 
+                                                    inflow_timescale = 300.0, kwargs...)
+
+    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, l, m, boundary_indices, boundary_adjacent_indices, 
+                                      grid, u, clock, model_fields, ΔX)
+    Δt = clock.last_stage_Δt
+
+    Δt = ifelse(isinf(Δt), 0, Δt)
+
+    ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields)
+
+    uᵢⁿ     = @inbounds getindex(u, boundary_indices...)
+    uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...)
+
+    U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹))
+
+    pa = bc.classification.matching_scheme
+
+    τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale)
+
+    τ̃ = Δt / τ
+
+    uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U)
+
+    @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...)
+
+    return nothing
+end
+
+@inline function step_left_boundary!(bc, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, 
+                                     grid, u, clock, model_fields, ΔX)
+    Δt = clock.last_stage_Δt
+
+    Δt = ifelse(isinf(Δt), 0, Δt)
+
+    ūⁿ⁺¹ = getbc(bc, l, m, grid, clock, model_fields)
+
+    uᵢⁿ     = @inbounds getindex(u, boundary_secret_storage_indices...)
+    uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, boundary_adjacent_indices...)
+
+    U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹))
+
+    pa = bc.classification.matching_scheme
+
+    τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale)
+
+    τ̃ = Δt / τ
+
+    u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U)
+
+    @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...)
+    @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...)
+
+    return nothing
+end
+
+@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, ::Tuple{Face, Any, Any}, clock, model_fields)
+    i = grid.Nx + 1
+
+    boundary_indices = (i, j, k)
+    boundary_adjacent_indices = (i-1, j, k)
+
+    Δx = Δxᶠᶜᶜ(i, j, k, grid)
+
+    step_right_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx)
+
+    return nothing
+end
+
+@inline function _fill_west_halo!(j, k, grid, u, bc::PAOBC, ::Tuple{Face, Any, Any}, clock, model_fields)
+    boundary_indices = (1, j, k)
+    boundary_adjacent_indices = (2, j, k)
+    boundary_secret_storage_indices = (0, j, k)
+
+    Δx = Δxᶠᶜᶜ(1, j, k, grid)
+
+    step_left_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δx)
+
+    return nothing
+end

src/Oceananigans.jl

@@ -221,6 +221,8 @@ include("Simulations/Simulations.jl")
 # Abstractions for distributed and multi-region models
 include("MultiRegion/MultiRegion.jl")
 
+include("boundary_normal_mean_velocity.jl")
+
 #####
 ##### Needed so we can export names from sub-modules at the top-level
 #####

src/boundary_normal_mean_velocity.jl

@@ -0,0 +1,46 @@
+using Adapt
+using Oceananigans.AbstractOperations: GridMetricOperation, Ax, Ay, Az
+using Oceananigans.BoundaryConditions: BoundaryCondition, Open, PerturbationAdvection
+
+import Adapt: adapt_structure
+import Oceananigans.BoundaryConditions: update_boundary_condition!
+
+struct BoundaryNormalMeanVelocity{BV}
+   boundary_velocity :: BV
+
+   BoundaryNormalMeanVelocity(; value::BV = Ref(0.0)) where BV = new{BV}(value)
+end
+
+@inline (value::BoundaryNormalMeanVelocity)(args...) = value.boundary_velocity[]
+
+Adapt.adapt_structure(to, mo::BoundaryNormalMeanVelocity) = 
+    BoundaryNormalMeanVelocity(adapt(to, mo.mean_outflow_velocity[]))
+
+const MOPABC = BoundaryCondition{<:Open{<:PerturbationAdvection}, <:BoundaryNormalMeanVelocity}
+
+@inline boundary_normal_area(::Union{Val{:west}, Val{:east}}, grid)   = GridMetricOperation((Face, Center, Center), Ax, grid)
+@inline boundary_normal_area(::Union{Val{:south}, Val{:north}}, grid) = GridMetricOperation((Center, Face, Center), Ay, grid)
+@inline boundary_normal_area(::Union{Val{:bottom}, Val{:top}}, grid)  = GridMetricOperation((Center, Center, Face), Az, grid)
+
+@inline boundary_adjacent_index(::Val{:east}, grid) = (size(grid, 1), 1, 1), (2, 3)
+@inline boundary_adjacent_index(::Val{:west}, grid) = (2, 1, 1), (2, 3)
+
+@inline boundary_adjacent_index(::Val{:north}, grid) = (1, size(grid, 2), 1), (1, 3)
+@inline boundary_adjacent_index(::Val{:south}, grid) = (1, 2, 1), (1, 3)
+
+@inline boundary_adjacent_index(::Val{:top}, grid)    = (1, 1, size(grid, 3)), (2, 3)
+@inline boundary_adjacent_index(::Val{:bottom}, grid) = (1, 1, 2), (2, 3)
+
+function update_boundary_condition!(bc::MOPABC, val_side, u, model)
+    grid = model.grid
+
+    An = boundary_normal_area(val_side, grid)
+
+    (i, j, k), dims = boundary_adjacent_index(val_side, grid)
+
+    total_area = sum(An; dims)[i, j, k]
+
+    Ū = sum(u * An; dims) / total_area
+
+    bc.condition.boundary_velocity[] = Ū[i, j, k]
+end