From 19a6887e3375f87a230df00dc5ea1edee71c2f15 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 16:39:21 +0100 Subject: [PATCH 1/3] polar boundary conditions --- src/BoundaryConditions/BoundaryConditions.jl | 1 + .../field_boundary_conditions.jl | 20 ++- .../polar_boundary_condition.jl | 121 ++++++++++++++++++ test/test_boundary_conditions.jl | 14 +- 4 files changed, 148 insertions(+), 8 deletions(-) create mode 100644 src/BoundaryConditions/polar_boundary_condition.jl diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl index 2011bd3df4..95c098f17a 100644 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ b/src/BoundaryConditions/BoundaryConditions.jl @@ -37,6 +37,7 @@ include("fill_halo_regions_nothing.jl") include("apply_flux_bcs.jl") include("update_boundary_conditions.jl") +include("polar_boundary_condition.jl") include("flat_extrapolation_open_boundary_matching_scheme.jl") end # module diff --git a/src/BoundaryConditions/field_boundary_conditions.jl b/src/BoundaryConditions/field_boundary_conditions.jl index 92c749af6c..38202982ef 100644 --- a/src/BoundaryConditions/field_boundary_conditions.jl +++ b/src/BoundaryConditions/field_boundary_conditions.jl @@ -59,7 +59,6 @@ end FieldBoundaryConditions(indices::Tuple, bcs::FieldBoundaryConditions) = FieldBoundaryConditions(indices, (getproperty(bcs, side) for side in propertynames(bcs))...) - FieldBoundaryConditions(indices::Tuple, ::Nothing) = nothing window_boundary_conditions(::Colon, left, right) = left, right @@ -178,6 +177,13 @@ function regularize_immersed_boundary_condition(ibc, grid, loc, field_name, args return NoFluxBoundaryCondition() end + regularize_west_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) + regularize_east_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) + regularize_south_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) + regularize_north_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) +regularize_bottom_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) + regularize_top_boundary_condition(bc, args...) = regularize_boundary_condition(bc, args...) + # regularize default boundary conditions function regularize_boundary_condition(default::DefaultBoundaryCondition, grid, loc, dim, args...) default_bc = default_prognostic_bc(topology(grid, dim)(), loc[dim](), default) @@ -212,12 +218,12 @@ function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions, loc = assumed_field_location(field_name) - west = regularize_boundary_condition(bcs.west, grid, loc, 1, LeftBoundary, prognostic_names) - east = regularize_boundary_condition(bcs.east, grid, loc, 1, RightBoundary, prognostic_names) - south = regularize_boundary_condition(bcs.south, grid, loc, 2, LeftBoundary, prognostic_names) - north = regularize_boundary_condition(bcs.north, grid, loc, 2, RightBoundary, prognostic_names) - bottom = regularize_boundary_condition(bcs.bottom, grid, loc, 3, LeftBoundary, prognostic_names) - top = regularize_boundary_condition(bcs.top, grid, loc, 3, RightBoundary, prognostic_names) + west = regularize_west_boundary_condition(bcs.west, grid, loc, 1, LeftBoundary, prognostic_names) + east = regularize_east_boundary_condition(bcs.east, grid, loc, 1, RightBoundary, prognostic_names) + south = regularize_south_boundary_condition(bcs.south, grid, loc, 2, LeftBoundary, prognostic_names) + north = regularize_north_boundary_condition(bcs.north, grid, loc, 2, RightBoundary, prognostic_names) + bottom = regularize_bottom_boundary_condition(bcs.bottom, grid, loc, 3, LeftBoundary, prognostic_names) + top = regularize_top_boundary_condition(bcs.top, grid, loc, 3, RightBoundary, prognostic_names) immersed = regularize_immersed_boundary_condition(bcs.immersed, grid, loc, field_name, prognostic_names) diff --git a/src/BoundaryConditions/polar_boundary_condition.jl b/src/BoundaryConditions/polar_boundary_condition.jl new file mode 100644 index 0000000000..160695d8f7 --- /dev/null +++ b/src/BoundaryConditions/polar_boundary_condition.jl @@ -0,0 +1,121 @@ +using Oceananigans.Architectures: device +using Oceananigans.Grids: inactive_node +using CUDA: @allowscalar + +struct PolarValue{D, S} + data :: D + side :: S +end + +Adapt.adapt_structure(to, pv::PolarValue) = PolarValue(Adapt.adapt(to, pv.data), nothing) + +const PolarBoundaryCondition{V} = BoundaryCondition{<:Value, <:PolarValue} + +function PolarBoundaryCondition(grid, side) + FT = eltype(grid) + arch = architecture(grid) + data = on_architecture(arch, zeros(FT, grid.Nz)) + return ValueBoundaryCondition(PolarValue(data, side)) +end + +# Just a column +@inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[k] + +# TODO: vectors should have a different treatment since vector components should account for the frame of reference +# North - South flux boundary conditions are not valid on a Latitude-Longitude grid if the last / first rows represent the poles +function latitude_north_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondition()) + # Check if the halo lies beyond the north pole + φmax = @allowscalar φnode(grid.Ny+1, grid, Center()) + + # No problem! + if φmax < 90 || loc[1] == Nothing + return default_bc + end + + return PolarBoundaryCondition(grid, :north) +end + +# North - South flux boundary conditions are not valid on a Latitude-Longitude grid if the last / first rows represent the poles +function latitude_south_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondition()) + # Check if the halo lies beyond the south pole + φmin = @allowscalar φnode(0, grid, Face()) + + # No problem! + if φmin > -90 || loc[1] == Nothing + return default_bc + end + + return PolarBoundaryCondition(grid, :south) +end + +regularize_north_boundary_condition(bc::DefaultBoundaryCondition, grid::LatitudeLongitudeGrid, loc, args...) = + regularize_boundary_condition(latitude_north_auxiliary_bc(grid, loc, bc), grid, loc, args...) + +regularize_south_boundary_condition(bc::DefaultBoundaryCondition, grid::LatitudeLongitudeGrid, loc, args...) = + regularize_boundary_condition(latitude_south_auxiliary_bc(grid, loc, bc), grid, loc, args...) + +@kernel function _average_pole_value!(data, c, j, grid, loc) + k = @index(Global, Linear) + c̄ = zero(grid) + n = 0 + @inbounds for i in 1:grid.Nx + inactive = inactive_node(i, j, k, grid, loc...) + c̄ += ifelse(inactive, 0, c[i, j, k]) + n += ifelse(inactive, 0, 1) + end + @inbounds data[k] = ifelse(n == 0, 0, c̄ / n) +end + +function update_pole_value!(bc::PolarValue, c, grid, loc) + j = bc.side == :north ? grid.Ny : 1 + dev = device(architecture(grid)) + _average_pole_value!(dev, 16, grid.Nz)(bc.data, c, j, grid, loc) + return nothing +end + +function fill_south_halo!(c, bc::PolarBoundaryCondition, size, offset, loc, arch, grid, args...; only_local_halos = false, kwargs...) + update_pole_value!(bc.condition, c, grid, loc) + return launch!(arch, grid, KernelParameters(size, offset), + _fill_only_south_halo!, c, bc, loc, grid, Tuple(args); kwargs...) +end + +function fill_north_halo!(c, bc::PolarBoundaryCondition, size, offset, loc, arch, grid, args...; only_local_halos = false, kwargs...) + update_pole_value!(bc.condition, c, grid, loc) + return launch!(arch, grid, KernelParameters(size, offset), + _fill_only_north_halo!, c, bc, loc, grid, Tuple(args); kwargs...) +end + +function fill_south_and_north_halo!(c, south_bc::PolarBoundaryCondition, north_bc, size, offset, loc, arch, grid, args...; only_local_halos = false, kwargs...) + update_pole_value!(south_bc.condition, c, grid, loc) + return launch!(arch, grid, KernelParameters(size, offset), + _fill_south_and_north_halo!, c, south_bc, north_bc, loc, grid, Tuple(args); kwargs...) +end + +function fill_south_and_north_halo!(c, south_bc, north_bc::PolarBoundaryCondition, size, offset, loc, arch, grid, args...; only_local_halos = false, kwargs...) + update_pole_value!(north_bc.condition, c, grid, loc) + return launch!(arch, grid, KernelParameters(size, offset), + _fill_south_and_north_halo!, c, south_bc, north_bc, loc, grid, Tuple(args); kwargs...) +end + +function fill_south_and_north_halo!(c, south_bc::PolarBoundaryCondition, north_bc::PolarBoundaryCondition, size, offset, loc, arch, grid, args...; only_local_halos = false, kwargs...) + update_pole_value!(south_bc.condition, c, grid, loc) + update_pole_value!(north_bc.condition, c, grid, loc) + return launch!(arch, grid, KernelParameters(size, offset), + _fill_south_and_north_halo!, c, south_bc, north_bc, loc, grid, Tuple(args); kwargs...) +end + +# If it is a LatitudeLongitudeGrid, we include the PolarBoundaryConditions +function FieldBoundaryConditions(grid::LatitudeLongitudeGrid, location, indices=(:, :, :); + west = default_auxiliary_bc(topology(grid, 1)(), location[1]()), + east = default_auxiliary_bc(topology(grid, 1)(), location[1]()), + south = default_auxiliary_bc(topology(grid, 2)(), location[2]()), + north = default_auxiliary_bc(topology(grid, 2)(), location[2]()), + bottom = default_auxiliary_bc(topology(grid, 3)(), location[3]()), + top = default_auxiliary_bc(topology(grid, 3)(), location[3]()), + immersed = NoFluxBoundaryCondition()) + + north = latitude_north_auxiliary_bc(grid, location, north) + south = latitude_south_auxiliary_bc(grid, location, south) + + return FieldBoundaryConditions(indices, west, east, south, north, bottom, top, immersed) +end \ No newline at end of file diff --git a/test/test_boundary_conditions.jl b/test/test_boundary_conditions.jl index 0cb88f4330..b9899cb25b 100644 --- a/test/test_boundary_conditions.jl +++ b/test/test_boundary_conditions.jl @@ -1,6 +1,6 @@ include("dependencies_for_runtests.jl") -using Oceananigans.BoundaryConditions: PBC, ZFBC, OBC, ContinuousBoundaryFunction, DiscreteBoundaryFunction, regularize_field_boundary_conditions +using Oceananigans.BoundaryConditions: PBC, ZFBC, VBC, OBC, ContinuousBoundaryFunction, DiscreteBoundaryFunction, regularize_field_boundary_conditions using Oceananigans.Fields: Face, Center simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t) @@ -252,5 +252,17 @@ end @test T_bcs.south === one_bc @test T_bcs.top === one_bc @test T_bcs.bottom === one_bc + + grid = LatitudeLongitudeGrid(size=(10, 10, 10), latitude=(-90, 90), longitude=(0, 360), z = (0, 1)) + f = CenterField(grid) + + @test f.boundary_conditions.north isa VBC + @test f.boundary_conditions.south isa VBC + + set!(f, (x, y, z) -> x) + fill_halo_regions!(f) + + @test all(f.data[1:10, 0, 1:10] .== 2 * mean(f.data[1:10, 1, 1:10]) .- f.data[1:10, 1, 1:10]) + @test all(f.data[1:10, 11, 1:10] .== 2 * mean(f.data[1:10, 10, 1:10]) .- f.data[1:10, 10, 1:10]) end end From 9a8e90cc023f170a30053fa4c5d3f5da2d9a74ad Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 19:04:38 +0100 Subject: [PATCH 2/3] do not touch face fields --- src/BoundaryConditions/polar_boundary_condition.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/BoundaryConditions/polar_boundary_condition.jl b/src/BoundaryConditions/polar_boundary_condition.jl index 160695d8f7..9e7bfc49d1 100644 --- a/src/BoundaryConditions/polar_boundary_condition.jl +++ b/src/BoundaryConditions/polar_boundary_condition.jl @@ -21,14 +21,15 @@ end # Just a column @inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[k] -# TODO: vectors should have a different treatment since vector components should account for the frame of reference +# TODO: vectors should have a different treatment since vector components should account for the frame of reference. +# For the moment, the `PolarBoundaryConditions` is implemented only for fields that have `loc[2] == Center()` # North - South flux boundary conditions are not valid on a Latitude-Longitude grid if the last / first rows represent the poles function latitude_north_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondition()) # Check if the halo lies beyond the north pole φmax = @allowscalar φnode(grid.Ny+1, grid, Center()) # No problem! - if φmax < 90 || loc[1] == Nothing + if φmax < 90 || loc[1] == Nothing || !(loc[2] == Center) return default_bc end @@ -41,7 +42,7 @@ function latitude_south_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondit φmin = @allowscalar φnode(0, grid, Face()) # No problem! - if φmin > -90 || loc[1] == Nothing + if φmin > -90 || loc[1] == Nothing || !(loc[2] == Center) return default_bc end From ce0d9f0006a497a534b065fbb7bcb47f8d7a7ab7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 3 Dec 2024 10:51:10 +0100 Subject: [PATCH 3/3] correct tests --- .../polar_boundary_condition.jl | 41 +++++++++++-------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/src/BoundaryConditions/polar_boundary_condition.jl b/src/BoundaryConditions/polar_boundary_condition.jl index 9e7bfc49d1..2d3cef673d 100644 --- a/src/BoundaryConditions/polar_boundary_condition.jl +++ b/src/BoundaryConditions/polar_boundary_condition.jl @@ -1,5 +1,4 @@ -using Oceananigans.Architectures: device -using Oceananigans.Grids: inactive_node +using Oceananigans.Grids: inactive_node, new_data using CUDA: @allowscalar struct PolarValue{D, S} @@ -11,29 +10,34 @@ Adapt.adapt_structure(to, pv::PolarValue) = PolarValue(Adapt.adapt(to, pv.data), const PolarBoundaryCondition{V} = BoundaryCondition{<:Value, <:PolarValue} -function PolarBoundaryCondition(grid, side) - FT = eltype(grid) - arch = architecture(grid) - data = on_architecture(arch, zeros(FT, grid.Nz)) +function PolarBoundaryCondition(grid, side, zloc) + FT = eltype(grid) + loc = (Nothing, Nothing, zloc) + data = new_data(FT, grid, loc) return ValueBoundaryCondition(PolarValue(data, side)) end # Just a column -@inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[k] +@inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[1, 1, k] # TODO: vectors should have a different treatment since vector components should account for the frame of reference. -# For the moment, the `PolarBoundaryConditions` is implemented only for fields that have `loc[2] == Center()` +# For the moment, the `PolarBoundaryConditions` is implemented only for fields that have `loc[1] == loc[2] == Center()`, which +# we assume are not components of horizontal vectors that would require rotation. (The `w` velocity if not a tracer, but it does +# not require rotation since it is a scalar field.) # North - South flux boundary conditions are not valid on a Latitude-Longitude grid if the last / first rows represent the poles function latitude_north_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondition()) # Check if the halo lies beyond the north pole φmax = @allowscalar φnode(grid.Ny+1, grid, Center()) + # Assumption: fields at `Center`s in x and y are not vector components + rotated_field = loc[1] != Center || loc[2] != Center + # No problem! - if φmax < 90 || loc[1] == Nothing || !(loc[2] == Center) + if φmax < 90 || rotated_field return default_bc end - return PolarBoundaryCondition(grid, :north) + return PolarBoundaryCondition(grid, :north, loc[3]) end # North - South flux boundary conditions are not valid on a Latitude-Longitude grid if the last / first rows represent the poles @@ -41,12 +45,15 @@ function latitude_south_auxiliary_bc(grid, loc, default_bc=DefaultBoundaryCondit # Check if the halo lies beyond the south pole φmin = @allowscalar φnode(0, grid, Face()) + # Assumption: fields at `Center`s in x and y are not vector components + rotated_field = loc[1] != Center || loc[2] != Center + # No problem! - if φmin > -90 || loc[1] == Nothing || !(loc[2] == Center) + if φmin > -90 || rotated_field return default_bc end - return PolarBoundaryCondition(grid, :south) + return PolarBoundaryCondition(grid, :south, loc[3]) end regularize_north_boundary_condition(bc::DefaultBoundaryCondition, grid::LatitudeLongitudeGrid, loc, args...) = @@ -56,7 +63,7 @@ regularize_south_boundary_condition(bc::DefaultBoundaryCondition, grid::Latitude regularize_boundary_condition(latitude_south_auxiliary_bc(grid, loc, bc), grid, loc, args...) @kernel function _average_pole_value!(data, c, j, grid, loc) - k = @index(Global, Linear) + i′, j′, k = @index(Global, NTuple) c̄ = zero(grid) n = 0 @inbounds for i in 1:grid.Nx @@ -64,13 +71,15 @@ regularize_south_boundary_condition(bc::DefaultBoundaryCondition, grid::Latitude c̄ += ifelse(inactive, 0, c[i, j, k]) n += ifelse(inactive, 0, 1) end - @inbounds data[k] = ifelse(n == 0, 0, c̄ / n) + @inbounds data[i′, j′, k] = ifelse(n == 0, 0, c̄ / n) end function update_pole_value!(bc::PolarValue, c, grid, loc) j = bc.side == :north ? grid.Ny : 1 - dev = device(architecture(grid)) - _average_pole_value!(dev, 16, grid.Nz)(bc.data, c, j, grid, loc) + Nz = size(c, 3) + Oz = c.offsets[3] + params = KernelParameters(1:1, 1:1, 1+Oz:Nz+Oz) + launch!(architecture(grid), grid, params, _average_pole_value!, bc.data, c, j, grid, loc) return nothing end