Small improvements to PerturbationAdvectionOpenBoundaryCondition interface (#4394)

* move docstring ti user-faacing method

* add docstring

* guard agaisnt infinity for backward step

* remove extra lines

* use ifelse function

* Update perturbation_advection_open_boundary_matching_scheme.jl

* Update perturbation_advection_open_boundary_matching_scheme.jl

* remove forward step

* don't hardcode float type

* move intuitive testing

* avoid calculations inside ifelse statements

* typo andf clarification

* import defaults first

* bugfix

---------

Co-authored-by: Jago Strong-Wright <[email protected]>

Author: tomchor

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl

@@ -1,4 +1,5 @@
 using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ, Ax_qᶠᶜᶜ, Ay_qᶜᶠᶜ, Az_qᶜᶜᶠ
+using Oceananigans: defaults
 
 """
     PerturbationAdvection
@@ -30,24 +31,41 @@ where Ũ = U Δt / Δx, then uⁿ⁺¹ is:
 where τ̃ = Δt/τ.
 
 The same operation can be repeated for left boundaries.
+
+The relaxation timescale ``τ`` can be set to different values depending on whether the 
+``U`` points in or out of the domain (`inflow_timescale`/`outflow_timescale`). Since the
+scheme is only valid when the flow is directed out of the domain the boundary condition
+falls back to relaxation to the prescribed value. By default this happens instantly but
+if the direction varies this may not be preferable. It is beneficial to relax the outflow
+(i.e. non-zero `outflow_timescale`) to reduce the shock when the flow changes direction
+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{VT, FT}
-       backward_step :: VT
+struct PerturbationAdvection{FT}
     inflow_timescale :: FT
    outflow_timescale :: FT
 end
 
 Adapt.adapt_structure(to, pe::PerturbationAdvection) =
-    PerturbationAdvection(adapt(to, pe.backward_step),
-                          adapt(to, pe.inflow_timescale),
+    PerturbationAdvection(adapt(to, pe.inflow_timescale),
                           adapt(to, pe.outflow_timescale))
 
-function PerturbationAdvectionOpenBoundaryCondition(val, FT = Float64;
-                                                    backward_step = true,
+"""
+    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 = 300.0, kwargs...)
+                                                    inflow_timescale = 0, kwargs...)
 
-    classification = Open(PerturbationAdvection(Val(backward_step), inflow_timescale, outflow_timescale))
+    classification = Open(PerturbationAdvection(inflow_timescale, outflow_timescale))
 
     @warn "`PerturbationAdvection` open boundaries matching scheme is experimental and un-tested/validated"
 
@@ -56,100 +74,45 @@ end
 
 const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}}
 
-const BPAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection{Val{true}}}}
-const FPAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection{Val{false}}}}
-
-@inline function step_right_boundary!(bc::BPAOBC, l, m, boundary_indices, boundary_adjacent_indices,
+@inline function step_right_boundary!(bc::PAOBC, 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)
+    relaxed_uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U)
+    uᵢⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_uᵢⁿ⁺¹)
 
     @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...)
 
     return nothing
 end
 
-@inline function step_left_boundary!(bc::BPAOBC, l, m, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices,
+@inline function step_left_boundary!(bc::PAOBC, 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 step_right_boundary!(bc::FPAOBC, 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ᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + (ūⁿ⁺¹ - uᵢⁿ) * τ̃
-    @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...)
-
-    return nothing
-end
-
-@inline function step_left_boundary!(bc::FPAOBC, 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 / τ
+    relaxed_u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U)
+    u₁ⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_u₁ⁿ⁺¹)
 
-    u₁ⁿ⁺¹ = uᵢⁿ - U * (uᵢ₋₁ⁿ⁺¹ - ūⁿ⁺¹) + (ūⁿ⁺¹ - uᵢⁿ) * τ̃
     @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...)
     @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...)
 

test/test_boundary_conditions_integration.jl

@@ -176,7 +176,7 @@ velocity_forcing(::Val{1}, forcing) = (; u = forcing)
 velocity_forcing(::Val{2}, forcing) = (; v = forcing)
 velocity_forcing(::Val{3}, forcing) = (; w = forcing)
 
-function test_pertubation_advection_open_boundary_conditions(arch, FT)
+function test_perturbation_advection_open_boundary_conditions(arch, FT)
     for orientation in 1:3
         topology = tuple(map(n -> ifelse(n == orientation, Bounded, Flat), 1:3)...)
 
@@ -447,7 +447,7 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
             A = typeof(arch)
             @info "  Testing open boundary conditions [$A, $FT]..."
             test_flat_extrapolation_open_boundary_conditions(arch, FT)
-            test_pertubation_advection_open_boundary_conditions(arch, FT)
+            test_perturbation_advection_open_boundary_conditions(arch, FT)
         end
     end
 end