Remove the use of indices for storage in PerturbationAdvection (#4578)

* remove secret_indices

* different index notation

* cleanup and improvement of validation scripts

* uncomment flat_extrapolation

* change notation for setp_left_boundary!

* follow nomenclature already in fill_halo_regions_value_gradient.jl

* slightly different notation

Author: tomchor

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl

@@ -26,13 +26,13 @@ relaxation to the mean velocity (i.e. u′→0) then Fᵤ = -u′ / τ then we c
     u′ⁿ⁺¹ = (uⁿ + Ũu′ⁿ⁺¹ᵢ₋₁ - Uⁿ⁺¹) / (1 + Ũ + Δt/τ),
 
 where Ũ = U Δt / Δx, then uⁿ⁺¹ is:
-    uⁿ⁺¹ = (uᵢⁿ + Ũuᵢ₋₁ⁿ⁺¹ + Uⁿ⁺¹τ̃) / (1 + τ̃ + U)
+    uⁿ⁺¹ = (uᵢⁿ + Ũuᵢ₋₁ⁿ⁺¹ + Uⁿ⁺¹τ̃) / (1 + τ̃ + Ũ)
 
 where τ̃ = Δt/τ.
 
 The same operation can be repeated for left boundaries.
 
-The relaxation timescale ``τ`` can be set to different values depending on whether the 
+The relaxation timescale ``τ`` can be set to different values depending on whether
 ``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
@@ -64,6 +64,8 @@ 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))
 
@@ -76,45 +78,48 @@ const PAOBC = BoundaryCondition{<:Open{<:PerturbationAdvection}}
 
 @inline function step_right_boundary!(bc::PAOBC, l, m, boundary_indices, boundary_adjacent_indices,
                                       grid, u, clock, model_fields, ΔX)
+    iᴮ, jᴮ, kᴮ = boundary_indices
+    iᴬ, jᴬ, kᴬ = boundary_adjacent_indices
     Δ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...)
+    ūⁿ⁺¹    = getbc(bc, l, m, grid, clock, model_fields)
+    uᵢⁿ     = @inbounds getindex(u, iᴮ, jᴮ, kᴮ)
+    uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ)
     U = max(0, min(1, Δt / ΔX * ūⁿ⁺¹))
 
     pa = bc.classification.matching_scheme
     τ = ifelse(ūⁿ⁺¹ >= 0, pa.outflow_timescale, pa.inflow_timescale)
-    τ̃ = Δt / τ
+    τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale
 
     relaxed_uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U)
-    uᵢⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_uᵢⁿ⁺¹)
+    uᵢⁿ⁺¹         = ifelse(τ == 0, ūⁿ⁺¹, relaxed_uᵢⁿ⁺¹)
 
-    @inbounds setindex!(u, uᵢⁿ⁺¹, boundary_indices...)
+    @inbounds setindex!(u, uᵢⁿ⁺¹, iᴮ, jᴮ, kᴮ)
 
     return nothing
 end
 
-@inline function step_left_boundary!(bc::PAOBC, 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,
                                      grid, u, clock, model_fields, ΔX)
+    iᴮ, jᴮ, kᴮ = boundary_indices
+    iᴬ, jᴬ, kᴬ = boundary_adjacent_indices
     Δ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...)
+    ūⁿ⁺¹    = getbc(bc, l, m, grid, clock, model_fields)
+    uᵢⁿ     = @inbounds getindex(u, iᴮ, jᴮ, kᴮ)
+    uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, iᴬ, jᴬ, kᴬ)
     U = min(0, max(-1, Δt / ΔX * ūⁿ⁺¹))
 
     pa = bc.classification.matching_scheme
     τ = ifelse(ūⁿ⁺¹ <= 0, pa.outflow_timescale, pa.inflow_timescale)
-    τ̃ = Δt / τ
+    τ̃ = Δt / τ # last stage Δt normalized by the inflow/output timescale
 
     relaxed_u₁ⁿ⁺¹ = (uᵢⁿ - U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ - U)
-    u₁ⁿ⁺¹ = ifelse(τ == 0, ūⁿ⁺¹, relaxed_u₁ⁿ⁺¹)
+    u₁ⁿ⁺¹         = ifelse(τ == 0, ūⁿ⁺¹, relaxed_u₁ⁿ⁺¹)
 
-    @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_indices...)
-    @inbounds setindex!(u, u₁ⁿ⁺¹, boundary_secret_storage_indices...)
+    @inbounds setindex!(u, u₁ⁿ⁺¹, iᴮ, jᴮ, kᴮ)
 
     return nothing
 end
@@ -134,11 +139,9 @@ 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)
+    step_left_boundary!(bc, j, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δx)
 
     return nothing
 end
@@ -149,7 +152,6 @@ end
     boundary_adjacent_indices = (i, j-1, k)
 
     Δy = Δyᶜᶠᶜ(i, j, k, grid)
-
     step_right_boundary!(bc, i, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δy)
 
     return nothing
@@ -158,11 +160,9 @@ end
 @inline function _fill_south_halo!(i, k, grid, u, bc::PAOBC, ::Tuple{Any, Face, Any}, clock, model_fields)
     boundary_indices = (i, 1, k)
     boundary_adjacent_indices = (i, 2, k)
-    boundary_secret_storage_indices = (i, 0, k)
 
     Δy = Δyᶜᶠᶜ(i, 1, k, grid)
-
-    step_left_boundary!(bc, i, k, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δy)
+    step_left_boundary!(bc, i, k, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δy)
 
     return nothing
 end
@@ -173,7 +173,6 @@ end
     boundary_adjacent_indices = (i, j, k-1)
 
     Δz = Δzᶜᶜᶠ(i, j, k, grid)
-
     step_right_boundary!(bc, i, j, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δz)
 
     return nothing
@@ -182,11 +181,9 @@ end
 @inline function _fill_bottom_halo!(i, j, grid, u, bc::PAOBC, ::Tuple{Any, Any, Face}, clock, model_fields)
     boundary_indices = (i, j, 1)
     boundary_adjacent_indices = (i, j, 2)
-    boundary_secret_storage_indices = (i, j, 0)
 
     Δz = Δzᶜᶜᶠ(i, j, 1, grid)
-
-    step_left_boundary!(bc, i, j, boundary_indices, boundary_adjacent_indices, boundary_secret_storage_indices, grid, u, clock, model_fields, Δz)
+    step_left_boundary!(bc, i, j, boundary_indices, boundary_adjacent_indices, grid, u, clock, model_fields, Δz)
 
     return nothing
 end

validation/open_boundaries/cylinder.jl

@@ -54,49 +54,34 @@ function drag(model;
     ∂₁uᵣ = Field(∂x(uᶜ), indices = (i₂, j₁:j₂, 1))
 
     ∂ₜuᶜ = Field(@at (Center, Center, Center) model.timestepper.Gⁿ.u)
-
     ∂ₜu = Field(∂ₜuᶜ, indices = (i₁:i₂, j₁:j₂, 1))
 
     p = model.pressures.pNHS
-
     ∫∂ₓp = Field(∂x(p), indices = (i₁:i₂, j₁:j₂, 1))
 
     a_local = Field(Integral(∂ₜu))
-
     a_flux = Field(Integral(uᵣ²)) - Field(Integral(uₗ²)) + Field(Integral(uvᵣ)) - Field(Integral(uvₗ))
-
     a_viscous_stress = 2ν * (Field(Integral(∂₁uᵣ)) - Field(Integral(∂₁uₗ)))
-
     a_pressure = Field(Integral(∫∂ₓp))
 
     return a_local + a_flux + a_pressure - a_viscous_stress
 end
 
 function cylinder_model(open_boundaries;
-
                         obc_name = "",
-
                         u∞ = 1,
                         r = 1/2,
-
                         cylinder = (x, y) -> ((x^2 + y^2) ≤ r^2),
-
                         stop_time = 100,
-
                         arch = GPU(),
-
                         Re = 100,
-                        Ny = 512,
+                        Ny = 128,
                         Nx = Ny,
-
                         ϵ = 0, # break up-down symmetry
                         x = (-6, 12), # 18
                         y = (-6 + ϵ, 6 + ϵ),  # 12
-
                         grid_kwargs = (; size=(Nx, Ny), x, y, halo=(6, 6), topology=(Bounded, Bounded, Flat)),
-
                         prefix = "flow_around_cylinder_Re$(Re)_Ny$(Ny)_$(obc_name)",
-
                         drag_averaging_window = 29)
 
     grid = RectilinearGrid(arch; grid_kwargs...)
@@ -108,9 +93,7 @@ function cylinder_model(open_boundaries;
     closure = ScalarDiffusivity(ν=1/Re)
 
     no_slip = ValueBoundaryCondition(0)
-
     u_bcs = FieldBoundaryConditions(immersed=no_slip, east=open_boundaries.east, west=open_boundaries.west)
-
     v_bcs = FieldBoundaryConditions(immersed=no_slip,
                                     east=GradientBoundaryCondition(0),
                                     west=ValueBoundaryCondition(0))
@@ -221,8 +204,8 @@ function cylinder_model(open_boundaries;
 
     supertitle = Label(fig[0, :], title)
 
-    CairoMakie.record(fig, prefix"*.mp4", 1:length(u.times), framerate=20) do i
-        @info "$i"
+    CairoMakie.record(fig, prefix * ".mp4", 1:length(u.times), framerate=20) do i
+        i % 50 == 0 && @info "$i"
         n[] = i
     end
 
@@ -231,15 +214,16 @@ end
 
 u∞ = 1
 
-feobc = (east = FlatExtrapolationOpenBoundaryCondition(), west = OpenBoundaryCondition(u∞))
+feobcs = (west = OpenBoundaryCondition(u∞), east = FlatExtrapolationOpenBoundaryCondition(),)
 
-paobcs = (east = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 1/4, outflow_timescale = Inf),
-          west = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 0.1, outflow_timescale = 0.1))
+paobcs = (west = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 0.1, outflow_timescale = 0.1),
+          east = PerturbationAdvectionOpenBoundaryCondition(u∞; inflow_timescale = 1/4, outflow_timescale = Inf),
+          )
 
-obcs = (; flat_extrapolation=feobc,
+obcs = (; flat_extrapolation=feobcs,
           perturbation_advection=paobcs)
 
 for (obc_name, obc) in pairs(obcs)
     @info "Running $(obc_name)"
-    cylinder_model(obc; obc_name, u∞)
+    cylinder_model(obc; obc_name, u∞, arch=CPU(), Ny=16)
 end

validation/open_boundaries/oscillating_flow.jl

@@ -8,7 +8,7 @@
 # This case also has a stretched grid to validate the matching scheme on a stretched grid.
 
 using Oceananigans, CairoMakie
-using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition
+using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition, PerturbationAdvectionOpenBoundaryCondition
 
 @kwdef struct Cylinder{FT}
     D :: FT = 1.0
@@ -122,6 +122,7 @@ function run_cylinder(grid, boundary_conditions; plot=true, stop_time = 50, simn
     end
 end
 
+inflow_timescale = outflow_timescale = 1/4
 matching_scheme_name(obc) = string(nameof(typeof(obc.classification.matching_scheme)))
 for grid in (xygrid, xzgrid)
 
@@ -132,14 +133,26 @@ for grid in (xygrid, xzgrid)
     u_boundaries_fe = FieldBoundaryConditions(west = u_fe, east = u_fe)
     v_boundaries_fe = FieldBoundaryConditions(south = v_fe, north = v_fe)
     w_boundaries_fe = FieldBoundaryConditions(bottom = w_fe, top = w_fe)
+    feobcs = (u = u_boundaries_fe, v = v_boundaries_fe, w = w_boundaries_fe)
 
-    if grid isa Oceananigans.Grids.ZFlatGrid
-        boundary_conditions = (u = u_boundaries_fe, v = v_boundaries_fe)
-        simname = "xy_" * matching_scheme_name(u_boundaries_fe.east)
-    elseif grid isa Oceananigans.Grids.YFlatGrid
-        boundary_conditions = (u = u_boundaries_fe, w = w_boundaries_fe)
-        simname = "xz_" * matching_scheme_name(u_boundaries_fe.east)
+    u_boundaries_pa = FieldBoundaryConditions(west   = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale),
+                                              east   = PerturbationAdvectionOpenBoundaryCondition(u∞; parameters = (; U, T), inflow_timescale, outflow_timescale))
+    v_boundaries_pa = FieldBoundaryConditions(south  = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale),
+                                              north  = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale))
+    w_boundaries_pa = FieldBoundaryConditions(bottom = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale),
+                                              top    = PerturbationAdvectionOpenBoundaryCondition(v∞; parameters = (; U, T), inflow_timescale, outflow_timescale))
+    paobcs = (u = u_boundaries_pa, v = v_boundaries_pa, w = w_boundaries_pa)
+
+    for obcs in (feobcs, paobcs,)
+        if grid isa Oceananigans.Grids.ZFlatGrid
+            boundary_conditions = (u = obcs.u, v = obcs.v)
+            simname = "xy_" * matching_scheme_name(boundary_conditions.u.east)
+        elseif grid isa Oceananigans.Grids.YFlatGrid
+            boundary_conditions = (u = obcs.u, w = obcs.w)
+            simname = "xz_" * matching_scheme_name(boundary_conditions.u.east)
+        end
+        @info "Running $simname"
+        run_cylinder(grid, boundary_conditions, simname = simname, stop_time = T)
     end
-    run_cylinder(grid, boundary_conditions, simname = simname, stop_time = T)
 end