Skip to content

Commit

Permalink
Refactoring for ImmersedBoundaries to take precedence over Solvers (
Browse files Browse the repository at this point in the history
#3847)

* try it out

* precompiles

* couple of bugfixes

* using inactive_node

* import regularize_boundary_condition

* flip the location

* flux form advection

* import conditional_flux

* disambiguation

* Update src/Coriolis/hydrostatic_spherical_coriolis.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* update importing

---------

Co-authored-by: Gregory L. Wagner <[email protected]>
  • Loading branch information
simone-silvestri and glwagner authored Oct 21, 2024
1 parent e8eccce commit 7cbf013
Show file tree
Hide file tree
Showing 18 changed files with 498 additions and 577 deletions.
1 change: 1 addition & 0 deletions src/Advection/Advection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ include("flux_form_advection.jl")

include("flat_advective_fluxes.jl")
include("topologically_conditional_interpolation.jl")
include("immersed_advective_fluxes.jl")
include("momentum_advection_operators.jl")
include("tracer_advection_operators.jl")
include("positivity_preserving_tracer_advection_operators.jl")
Expand Down
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
using Oceananigans.Advection: AbstractAdvectionScheme, advection_buffers
using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑxᶜᵃᵃ, ℑyᵃᶠᵃ, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, ℑzᵃᵃᶜ
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractTimeDiscretization
using Oceananigans.Advection: LOADV, HOADV, WENO, FluxFormAdvection
using Oceananigans.ImmersedBoundaries
using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node
using Oceananigans.Fields: ZeroField

const ATC = AbstractTurbulenceClosure
const ATD = AbstractTimeDiscretization
const IBG = ImmersedBoundaryGrid

const c = Center()
const f = Face()

"""
conditional_flux(i, j, k, ibg::IBG, ℓx, ℓy, ℓz, qᴮ, qᴵ, nc)
Expand Down Expand Up @@ -33,31 +33,7 @@ end
@inline conditional_flux_ccf(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, c, c, f, qᴮ, qᴵ)

#####
##### Diffusive fluxes
#####

# ccc, ffc, fcf
@inline _viscous_flux_ux(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), viscous_flux_ux(i, j, k, ibg, args...))
@inline _viscous_flux_uy(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), viscous_flux_uy(i, j, k, ibg, args...))
@inline _viscous_flux_uz(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), viscous_flux_uz(i, j, k, ibg, args...))

# ffc, ccc, cff
@inline _viscous_flux_vx(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), viscous_flux_vx(i, j, k, ibg, args...))
@inline _viscous_flux_vy(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), viscous_flux_vy(i, j, k, ibg, args...))
@inline _viscous_flux_vz(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), viscous_flux_vz(i, j, k, ibg, args...))

# fcf, cff, ccc
@inline _viscous_flux_wx(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), viscous_flux_wx(i, j, k, ibg, args...))
@inline _viscous_flux_wy(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), viscous_flux_wy(i, j, k, ibg, args...))
@inline _viscous_flux_wz(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), viscous_flux_wz(i, j, k, ibg, args...))

# fcc, cfc, ccf
@inline _diffusive_flux_x(i, j, k, ibg::IBG, args...) = conditional_flux_fcc(i, j, k, ibg, zero(ibg), diffusive_flux_x(i, j, k, ibg, args...))
@inline _diffusive_flux_y(i, j, k, ibg::IBG, args...) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), diffusive_flux_y(i, j, k, ibg, args...))
@inline _diffusive_flux_z(i, j, k, ibg::IBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), diffusive_flux_z(i, j, k, ibg, args...))

#####
##### Advective fluxes
##### Immersed Advective fluxes
#####

# dx(uu), dy(vu), dz(wu)
Expand All @@ -82,69 +58,27 @@ end
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, args...) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_tracer_flux_y(i, j, k, ibg, args...))
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), advective_tracer_flux_z(i, j, k, ibg, args...))

# Disambiguation for tracer fluxes....
@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_x(i, j, k, ibg, advection.x, args...)

@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_y(i, j, k, ibg, advection.y, args...)

@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_z(i, j, k, ibg, advection.z, args...)

# Disambiguation for momentum fluxes in x....
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Uu(i, j, k, ibg, advection.x, args...)

@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Uv(i, j, k, ibg, advection.x, args...)

@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Uw(i, j, k, ibg, advection.x, args...)

# Disambiguation for momentum fluxes in y....
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Vu(i, j, k, ibg, advection.y, args...)

@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Vv(i, j, k, ibg, advection.y, args...)

@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Vw(i, j, k, ibg, advection.y, args...)

# Disambiguation for momentum fluxes in z....
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Wu(i, j, k, ibg, advection.z, args...)

@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Wv(i, j, k, ibg, advection.z, args...)

@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_momentum_flux_Ww(i, j, k, ibg, advection.z, args...)

# Fallback for `nothing` Advection
# Fallback for `nothing` advection
@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)

# dx(uu), dy(vu), dz(wu)
# ccc, ffc, fcf
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
# Disambiguation for `FluxForm` advection
@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uu(i, j, k, ibg, advection.x, args...))
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vu(i, j, k, ibg, advection.y, args...))
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wu(i, j, k, ibg, advection.z, args...))

# dx(uv), dy(vv), dz(wv)
# ffc, ccc, cff
@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uv(i, j, k, ibg, advection.x, args...))
@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vv(i, j, k, ibg, advection.y, args...))
@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wv(i, j, k, ibg, advection.z, args...))

# dx(uw), dy(vw), dz(ww)
# fcf, cff, ccc
@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uw(i, j, k, ibg, advection.x, args...))
@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vw(i, j, k, ibg, advection.y, args...))
@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Ww(i, j, k, ibg, advection.z, args...))

@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg)
@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_fcc(i, j, k, ibg, zero(ibg), advective_tracer_flux_x(i, j, k, ibg, advection.x, args...))
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_tracer_flux_y(i, j, k, ibg, advection.y, args...))
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), advective_tracer_flux_z(i, j, k, ibg, advection.z, args...))

#####
##### "Boundary-aware" reconstruct
Expand Down
10 changes: 10 additions & 0 deletions src/Advection/tracer_advection_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,16 @@
##### Fallback tracer fluxes!
#####

# Disambiguation for tracer fluxes....
@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_x(i, j, k, ibg, advection.x, args...)

@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_y(i, j, k, ibg, advection.y, args...)

@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) =
_advective_tracer_flux_z(i, j, k, ibg, advection.z, args...)

# Fallback for `nothing` advection
@inline _advective_tracer_flux_x(i, j, k, grid, ::Nothing, args...) = zero(grid)
@inline _advective_tracer_flux_y(i, j, k, grid, ::Nothing, args...) = zero(grid)
Expand Down
2 changes: 2 additions & 0 deletions src/Coriolis/hydrostatic_spherical_coriolis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ using Oceananigans.Advection: EnergyConserving, EnstrophyConserving
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.ImmersedBoundaries

"""
struct ActiveCellEnstrophyConserving
Expand Down Expand Up @@ -49,6 +50,7 @@ Adapt.adapt_structure(to, coriolis::HydrostaticSphericalCoriolis) =

@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = φnode(j, grid, Face())
@inline φᶠᶠᵃ(i, j, k, grid::OrthogonalSphericalShellGrid) = φnode(i, j, grid, Face(), Face())
@inline φᶠᶠᵃ(i, j, k, grid::ImmersedBoundaryGrid) = φᶠᶠᵃ(i, j, k, grid.underlying_grid)

@inline fᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis) =
2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))
Expand Down
1 change: 1 addition & 0 deletions src/DistributedComputations/DistributedComputations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ using Oceananigans.Grids
include("distributed_architectures.jl")
include("partition_assemble.jl")
include("distributed_grids.jl")
include("distributed_immersed_boundaries.jl")
include("distributed_on_architecture.jl")
include("distributed_kernel_launching.jl")
include("halo_communication_bcs.jl")
Expand Down
154 changes: 154 additions & 0 deletions src/DistributedComputations/distributed_immersed_boundaries.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
using Oceananigans.Utils: getnamewrapper
using Oceananigans.ImmersedBoundaries
using Oceananigans.ImmersedBoundaries: AbstractGridFittedBottom,
GridFittedBottom,
GridFittedBoundary,
compute_mask,
interior_active_indices

import Oceananigans.ImmersedBoundaries: map_interior_active_cells

# For the moment we extend distributed in the `ImmersedBoundaryGrids` module.
# When we fix the immersed boundary module to remove all the `TurbulenceClosure` stuff
# we can move this file back to `DistributedComputations` if we want `ImmersedBoundaries`
# to take precedence
const DistributedImmersedBoundaryGrid = ImmersedBoundaryGrid{FT, TX, TY, TZ, <:DistributedGrid, I, M, <:Distributed} where {FT, TX, TY, TZ, I, M}

function reconstruct_global_grid(grid::ImmersedBoundaryGrid)
arch = grid.architecture
local_ib = grid.immersed_boundary
global_ug = reconstruct_global_grid(grid.underlying_grid)
global_ib = getnamewrapper(local_ib)(construct_global_array(arch, local_ib.bottom_height, size(grid)))
return ImmersedBoundaryGrid(global_ug, global_ib)
end

function with_halo(new_halo, grid::DistributedImmersedBoundaryGrid)
immersed_boundary = grid.immersed_boundary
underlying_grid = grid.underlying_grid
new_underlying_grid = with_halo(new_halo, underlying_grid)
new_immersed_boundary = resize_immersed_boundary(immersed_boundary, new_underlying_grid)
return ImmersedBoundaryGrid(new_underlying_grid, new_immersed_boundary)
end

function scatter_local_grids(global_grid::ImmersedBoundaryGrid, arch::Distributed, local_size)
ib = global_grid.immersed_boundary
ug = global_grid.underlying_grid

local_ug = scatter_local_grids(ug, arch, local_size)

# Kinda hacky
local_bottom_height = partition(ib.bottom_height, arch, local_size)
ImmersedBoundaryConstructor = getnamewrapper(ib)
local_ib = ImmersedBoundaryConstructor(local_bottom_height)

return ImmersedBoundaryGrid(local_ug, local_ib)
end

"""
function resize_immersed_boundary!(ib, grid)
If the immersed condition is an `OffsetArray`, resize it to match
the total size of `grid`
"""
resize_immersed_boundary(ib::AbstractGridFittedBottom, grid) = ib
resize_immersed_boundary(ib::GridFittedBoundary, grid) = ib

function resize_immersed_boundary(ib::GridFittedBoundary{<:OffsetArray}, grid)

Nx, Ny, Nz = size(grid)
Hx, Hy, Nz = halo_size(grid)

mask_size = (Nx, Ny, Nz) .+ 2 .* (Hx, Hy, Hz)

# Check that the size of a bottom field are
# consistent with the size of the grid
if any(size(ib.mask) .!= mask_size)
@warn "Resizing the mask to match the grids' halos"
mask = compute_mask(grid, ib)
return getnamewrapper(ib)(mask)
end

return ib
end

function resize_immersed_boundary(ib::AbstractGridFittedBottom{<:OffsetArray}, grid)

Nx, Ny, _ = size(grid)
Hx, Hy, _ = halo_size(grid)

bottom_heigth_size = (Nx, Ny) .+ 2 .* (Hx, Hy)

# Check that the size of a bottom field are
# consistent with the size of the grid
if any(size(ib.bottom_height) .!= bottom_heigth_size)
@warn "Resizing the bottom field to match the grids' halos"
bottom_field = Field((Center, Center, Nothing), grid)
cpu_bottom = on_architecture(CPU(), ib.bottom_height)[1:Nx, 1:Ny]
set!(bottom_field, cpu_bottom)
fill_halo_regions!(bottom_field)
offset_bottom_array = dropdims(bottom_field.data, dims=3)

return getnamewrapper(ib)(offset_bottom_array)
end

return ib
end

# A distributed grid with split interior map
const DistributedActiveCellsIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:DistributedGrid, <:Any, <:NamedTuple}

# In case of a `DistributedGrid` we want to have different maps depending on the partitioning of the domain:
#
# If we partition the domain in the x-direction, we typically want to have the option to split three-dimensional
# kernels in a `halo-independent` part in the range Hx+1:Nx-Hx, 1:Ny, 1:Nz and two `halo-dependent` computations:
# a west one spanning 1:Hx, 1:Ny, 1:Nz and an east one spanning Nx-Hx+1:Nx, 1:Ny, 1:Nz.
# For this reason we need three different maps, one containing the `halo_independent` active region, a `west` map and an `east` map.
# For the same reason we need to construct `south` and `north` maps if we partition the domain in the y-direction.
# Therefore, the `interior_active_cells` in this case is a `NamedTuple` containing 5 elements.
# Note that boundary-adjacent maps corresponding to non-partitioned directions are set to `nothing`
function map_interior_active_cells(ibg::ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:DistributedGrid})

arch = architecture(ibg)

# If we using a synchronized architecture, nothing
# changes with serial execution.
if arch isa SynchronizedDistributed
return interior_active_indices(ibg; parameters = :xyz)
end

Rx, Ry, _ = arch.ranks
Tx, Ty, _ = topology(ibg)
Nx, Ny, Nz = size(ibg)
Hx, Hy, _ = halo_size(ibg)

x_boundary = (Hx, Ny, Nz)
y_boundary = (Nx, Hy, Nz)

left_offsets = (0, 0, 0)
right_x_offsets = (Nx-Hx, 0, 0)
right_y_offsets = (0, Ny-Hy, 0)

include_west = !isa(ibg, XFlatGrid) && (Rx != 1) && !(Tx == RightConnected)
include_east = !isa(ibg, XFlatGrid) && (Rx != 1) && !(Tx == LeftConnected)
include_south = !isa(ibg, YFlatGrid) && (Ry != 1) && !(Ty == RightConnected)
include_north = !isa(ibg, YFlatGrid) && (Ry != 1) && !(Ty == LeftConnected)

west_halo_dependent_cells = include_west ? interior_active_indices(ibg; parameters = KernelParameters(x_boundary, left_offsets)) : nothing
east_halo_dependent_cells = include_east ? interior_active_indices(ibg; parameters = KernelParameters(x_boundary, right_x_offsets)) : nothing
south_halo_dependent_cells = include_south ? interior_active_indices(ibg; parameters = KernelParameters(y_boundary, left_offsets)) : nothing
north_halo_dependent_cells = include_north ? interior_active_indices(ibg; parameters = KernelParameters(y_boundary, right_y_offsets)) : nothing

nx = Rx == 1 ? Nx : (Tx == RightConnected || Tx == LeftConnected ? Nx - Hx : Nx - 2Hx)
ny = Ry == 1 ? Ny : (Ty == RightConnected || Ty == LeftConnected ? Ny - Hy : Ny - 2Hy)

ox = Rx == 1 || Tx == RightConnected ? 0 : Hx
oy = Ry == 1 || Ty == RightConnected ? 0 : Hy

halo_independent_cells = interior_active_indices(ibg; parameters = KernelParameters((nx, ny, Nz), (ox, oy, 0)))

return (; halo_independent_cells,
west_halo_dependent_cells,
east_halo_dependent_cells,
south_halo_dependent_cells,
north_halo_dependent_cells)
end
Loading

0 comments on commit 7cbf013

Please sign in to comment.