(0.95.3) Restores open boundary condition functionality

src/BoundaryConditions/fill_halo_regions.jl

@@ -20,12 +20,14 @@ conditions, possibly recursing into `fields` if it is a nested tuple-of-tuples.
 # Some fields have `nothing` boundary conditions, such as `FunctionField` and `ZeroField`.
 fill_halo_regions!(c::OffsetArray, ::Nothing, args...; kwargs...) = nothing
 
+retrieve_bc(bc) = bc
+
 # Returns the boundary conditions a specific side for `FieldBoundaryConditions` inputs and
 # a tuple of boundary conditions for `NTuple{N, <:FieldBoundaryConditions}` inputs
 for dir in (:west, :east, :south, :north, :bottom, :top)
     extract_side_bc = Symbol(:extract_, dir, :_bc)
     @eval begin
-        @inline $extract_side_bc(bc) = bc.$dir
+        @inline $extract_side_bc(bc) = retrieve_bc(bc.$dir)
         @inline $extract_side_bc(bc::Tuple) = map($extract_side_bc, bc)
     end
 end

src/BoundaryConditions/fill_halo_regions_open.jl

@@ -48,6 +48,8 @@ end
 @inline retrieve_open_bc(bc::OBC) = bc
 @inline retrieve_open_bc(bc) = nothing
 
+@inline retrieve_bc(bc::OBC) = nothing
+
 # for regular halo fills, return nothing if the BC is not an OBC
 @inline left_open_boundary_condition(boundary_condition, loc) = nothing
 @inline left_open_boundary_condition(boundary_conditions, ::Tuple{Face, Center, Center}) = retrieve_open_bc(boundary_conditions.west)
@@ -59,7 +61,7 @@ end
 @inline right_open_boundary_condition(boundary_conditions, ::Tuple{Center, Face, Center}) = retrieve_open_bc(boundary_conditions.north)
 @inline right_open_boundary_condition(boundary_conditions, ::Tuple{Center, Center, Face}) = retrieve_open_bc(boundary_conditions.top)
 
-# Opern boundary fill 
+# Open boundary fill 
 @inline   _fill_west_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[1, j, k]           = getbc(bc, j, k, grid, args...)
 @inline   _fill_east_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[grid.Nx + 1, j, k] = getbc(bc, j, k, grid, args...)
 @inline  _fill_south_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, 1, k]           = getbc(bc, i, k, grid, args...)

src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Grids: xspacing, yspacing, zspacing
+using Oceananigans.Operators: Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ
 
 """
     FlatExtrapolation
@@ -59,12 +59,10 @@ end
     return ϕ + Δϕ
 end
 
-const c = Center()
-
-@inline function _fill_west_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
-    Δx₁ = xspacing(1, j, k, grid, c, c, c)
-    Δx₂ = xspacing(2, j, k, grid, c, c, c)
-    Δx₃ = xspacing(3, j, k, grid, c, c, c)
+@inline function _fill_west_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+    Δx₁ = Δxᶜᶜᶜ(1, j, k, grid)
+    Δx₂ = Δxᶜᶜᶜ(2, j, k, grid)
+    Δx₃ = Δxᶜᶜᶜ(3, j, k, grid)
 
     spacing_factor = Δx₁ / (Δx₂ + Δx₃)
 
@@ -75,12 +73,12 @@ const c = Center()
     return nothing
 end
 
-@inline function _fill_east_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+@inline function _fill_east_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
     i = grid.Nx + 1
 
-    Δx₁ = xspacing(i-1, j, k, grid, c, c, c)
-    Δx₂ = xspacing(i-2, j, k, grid, c, c, c)
-    Δx₃ = xspacing(i-3, j, k, grid, c, c, c)
+    Δx₁ = Δxᶜᶜᶜ(i-1, j, k, grid)
+    Δx₂ = Δxᶜᶜᶜ(i-2, j, k, grid)
+    Δx₃ = Δxᶜᶜᶜ(i-3, j, k, grid)
 
     spacing_factor = Δx₁ / (Δx₂ + Δx₃)
 
@@ -91,10 +89,10 @@ end
     return nothing
 end
 
-@inline function _fill_south_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
-    Δy₁ = yspacing(i, 1, k, grid, c, c, c)
-    Δy₂ = yspacing(i, 2, k, grid, c, c, c)
-    Δy₃ = yspacing(i, 3, k, grid, c, c, c)
+@inline function _fill_south_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+    Δy₁ = Δyᶜᶜᶜ(i, 1, k, grid)
+    Δy₂ = Δyᶜᶜᶜ(i, 2, k, grid)
+    Δy₃ = Δyᶜᶜᶜ(i, 3, k, grid)
 
     spacing_factor = Δy₁ / (Δy₂ + Δy₃)
 
@@ -105,12 +103,12 @@ end
     return nothing
 end
 
-@inline function _fill_north_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+@inline function _fill_north_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
     j = grid.Ny + 1
 
-    Δy₁ = yspacing(i, j-1, k, grid, c, c, c)
-    Δy₂ = yspacing(i, j-2, k, grid, c, c, c)
-    Δy₃ = yspacing(i, j-3, k, grid, c, c, c)
+    Δy₁ = Δyᶜᶜᶜ(i, j-1, k, grid)
+    Δy₂ = Δyᶜᶜᶜ(i, j-2, k, grid)
+    Δy₃ = Δyᶜᶜᶜ(i, j-3, k, grid)
 
     spacing_factor = Δy₁ / (Δy₂ + Δy₃)
 
@@ -121,10 +119,10 @@ end
     return nothing
 end
 
-@inline function _fill_bottom_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
-    Δz₁ = zspacing(i, j, 1, grid, c, c, c)
-    Δz₂ = zspacing(i, j, 2, grid, c, c, c)
-    Δz₃ = zspacing(i, j, 3, grid, c, c, c)
+@inline function _fill_bottom_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+    Δz₁ = Δzᶜᶜᶜ(i, j, 1, grid)
+    Δz₂ = Δzᶜᶜᶜ(i, j, 2, grid)
+    Δz₃ = Δzᶜᶜᶜ(i, j, 3, grid)
 
     spacing_factor = Δz₁ / (Δz₂ + Δz₃)
 
@@ -135,12 +133,12 @@ end
     return nothing
 end
 
-@inline function _fill_top_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
+@inline function _fill_top_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields)
     k = grid.Nz + 1
 
-    Δz₁ = zspacing(i, j, k-1, grid, c, c, c)
-    Δz₂ = zspacing(i, j, k-2, grid, c, c, c)
-    Δz₃ = zspacing(i, j, k-3, grid, c, c, c)
+    Δz₁ = Δzᶜᶜᶜ(i, j, k-1, grid)
+    Δz₂ = Δzᶜᶜᶜ(i, j, k-2, grid)
+    Δz₃ = Δzᶜᶜᶜ(i, j, k-3, grid)
 
     spacing_factor = Δz₁ / (Δz₂ + Δz₃)
 

test/test_boundary_conditions_integration.jl

@@ -1,6 +1,6 @@
 include("dependencies_for_runtests.jl")
 
-using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction
+using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction, FlatExtrapolationOpenBoundaryCondition, fill_halo_regions!
 using Oceananigans: prognostic_fields
 
 function test_boundary_condition(arch, FT, topo, side, field_name, boundary_condition)
@@ -103,6 +103,63 @@ function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
     return isapprox(mean(b) - mean_b₀, flux * model.clock.time / Lz, atol=1e-6)
 end
 
+left_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, east = OpenBoundaryCondition(1),
+                                                                    west = FlatExtrapolationOpenBoundaryCondition())
+
+right_febc(::Val{1}, grid, loc) = FieldBoundaryConditions(grid, loc, west = OpenBoundaryCondition(1),
+                                                                     east = FlatExtrapolationOpenBoundaryCondition())
+
+left_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, north = OpenBoundaryCondition(1),
+                                                                    south = FlatExtrapolationOpenBoundaryCondition())
+
+right_febc(::Val{2}, grid, loc) = FieldBoundaryConditions(grid, loc, south = OpenBoundaryCondition(1),
+                                                                     north = FlatExtrapolationOpenBoundaryCondition())
+
+left_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, top = OpenBoundaryCondition(1),
+                                                                    bottom = FlatExtrapolationOpenBoundaryCondition())
+
+right_febc(::Val{3}, grid, loc) = FieldBoundaryConditions(grid, loc, bottom = OpenBoundaryCondition(1),
+                                                                     top = FlatExtrapolationOpenBoundaryCondition())
+
+end_position(::Val{1}, grid) = (grid.Nx+1, 1, 1)
+end_position(::Val{2}, grid) = (1, grid.Ny+1, 1)
+end_position(::Val{3}, grid) = (1, 1, grid.Nz+1)
+
+function flat_extrapolation_open_boundary_conditions_are_correct(arch, FT)
+    clock = Clock(; time = zero(FT))
+
+    for orientation in 1:3
+        topology = tuple(map(n -> ifelse(n == orientation, Bounded, Flat), 1:3)...)
+
+        normal_location = tuple(map(n -> ifelse(n == orientation, Face, Center), 1:3)...)
+
+        grid = RectilinearGrid(arch, FT; topology, size = (16, ), x = (0, 1), y = (0, 1), z = (0, 1))
+
+        bcs1 = left_febc(Val(orientation), grid, normal_location)
+        bcs2 = right_febc(Val(orientation), grid, normal_location)
+
+        u1 = Field{normal_location...}(grid; boundary_conditions = bcs1)
+        u2 = Field{normal_location...}(grid; boundary_conditions = bcs2)
+
+        set!(u1, (X, ) -> 2-X)
+        set!(u2, (X, ) -> 1+X)
+
+        fill_halo_regions!(u1, clock, (); fill_boundary_normal_velocities = false)
+        fill_halo_regions!(u2, clock, (); fill_boundary_normal_velocities = false)
+
+        # we can stop the wall normal halos being filled after the pressure solve
+        @test interior(u1, 1, 1, 1) .== 2
+        @test interior(u2, end_position(Val(orientation), grid)...) .== 2
+
+        fill_halo_regions!(u1, clock, ())
+        fill_halo_regions!(u2, clock, ())
+
+        # now they should be filled
+        @test interior(u1, 1, 1, 1) .≈ 1.8125
+        @test interior(u2, end_position(Val(orientation), grid)...) .≈ 1.8125
+    end
+end
+
 test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
                                               float_bc(C, FT, ArrayType),
                                               irrational_bc(C, FT, ArrayType),
@@ -273,4 +330,12 @@ test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
             @test fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
         end
     end
+
+    @testset "Open boundary conditions" begin
+        for arch in archs, FT in (Float64,) #float_types
+            A = typeof(arch)
+            @info "  Testing open boundary conditions [$A, $FT]..."
+            @test flat_extrapolation_open_boundary_conditions_are_correct(arch, FT)
+        end
+    end
 end