diff --git a/Project.toml b/Project.toml index b7de074bd9..c41dd2bfc0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" authors = ["Climate Modeling Alliance and contributors"] -version = "0.95.2" +version = "0.95.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/BoundaryConditions/fill_halo_regions.jl b/src/BoundaryConditions/fill_halo_regions.jl index 7ac3e79ea5..5eb29f687e 100644 --- a/src/BoundaryConditions/fill_halo_regions.jl +++ b/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 diff --git a/src/BoundaryConditions/fill_halo_regions_open.jl b/src/BoundaryConditions/fill_halo_regions_open.jl index 039de899dd..e19bd3245a 100644 --- a/src/BoundaryConditions/fill_halo_regions_open.jl +++ b/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...) diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index fbbcb748f3..1325b26173 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -1,4 +1,4 @@ -using Oceananigans.Operators: Δx, Δy, Δz +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₁ = Δx(1, j, k, grid, c, c, c) - Δx₂ = Δx(2, j, k, grid, c, c, c) - Δx₃ = Δx(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₁ = Δx(i-1, j, k, grid, c, c, c) - Δx₂ = Δx(i-2, j, k, grid, c, c, c) - Δx₃ = Δx(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₁ = Δy(i, 1, k, grid, c, c, c) - Δy₂ = Δy(i, 2, k, grid, c, c, c) - Δy₃ = Δy(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₁ = Δy(i, j-1, k, grid, c, c, c) - Δy₂ = Δy(i, j-2, k, grid, c, c, c) - Δy₃ = Δy(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₁ = Δz(i, j, 1, grid, c, c, c) - Δz₂ = Δz(i, j, 2, grid, c, c, c) - Δz₃ = Δz(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₁ = Δz(i, j, k-1, grid, c, c, c) - Δz₂ = Δz(i, j, k-2, grid, c, c, c) - Δz₃ = Δz(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₃) diff --git a/test/test_boundary_conditions_integration.jl b/test/test_boundary_conditions_integration.jl index 16efa34729..de169e0a2d 100644 --- a/test/test_boundary_conditions_integration.jl +++ b/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 test_flat_extrapolation_open_boundary_conditions(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(arch, FT) + end + end end