-
Notifications
You must be signed in to change notification settings - Fork 191
/
grid_fitted_immersed_boundaries.jl
211 lines (163 loc) · 8.26 KB
/
grid_fitted_immersed_boundaries.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
using Adapt
using CUDA: CuArray
using OffsetArrays: OffsetArray
using Oceananigans.Utils: getnamewrapper
using Oceananigans.Grids: total_size
using Oceananigans.Fields: fill_halo_regions!
using Oceananigans.Architectures: arch_array
using Oceananigans.BoundaryConditions: FBC
using Printf
import Oceananigans.TurbulenceClosures: ivd_upper_diagonal,
ivd_lower_diagonal,
z_bottom
import Oceananigans.TurbulenceClosures: immersed_∂ⱼ_τ₁ⱼ,
immersed_∂ⱼ_τ₂ⱼ,
immersed_∂ⱼ_τ₃ⱼ,
immersed_∇_dot_qᶜ
#####
##### Some conveniences for grid fitted boundaries
#####
abstract type AbstractGridFittedBoundary <: AbstractImmersedBoundary end
const GFIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:AbstractGridFittedBoundary}
#####
##### ImmersedBoundaryGrids require one additional halo to check `inactive_node` for
##### Faces on the first halo
#####
#####
##### GridFittedBottom (2.5D immersed boundary with modified bottom height)
#####
abstract type AbstractGridFittedBottom{H} <: AbstractGridFittedBoundary end
struct CenterImmersedCondition end
struct InterfaceImmersedCondition end
"""
GridFittedBottom(bottom_height, [immersed_condition=CenterImmersedCondition()])
Return an immersed boundary with an irregular bottom fit to the underlying grid.
"""
struct GridFittedBottom{H, I} <: AbstractGridFittedBottom{H}
bottom_height :: H
immersed_condition :: I
end
GridFittedBottom(bottom_height) = GridFittedBottom(bottom_height, CenterImmersedCondition())
function Base.summary(ib::GridFittedBottom)
hmax = maximum(parent(ib.bottom_height))
hmin = minimum(parent(ib.bottom_height))
return @sprintf("GridFittedBottom(min(h)=%.2e, max(h)=%.2e)", hmin, hmax)
end
Base.summary(ib::GridFittedBottom{<:Function}) = @sprintf("GridFittedBottom(%s)", ib.bottom_height)
Base.show(io::IO, ib::GridFittedBottom) = print(io, summary(ib))
"""
ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
Return a grid with `GridFittedBottom` immersed boundary.
Computes ib.bottom_height and wraps in an array.
"""
function ImmersedBoundaryGrid(grid, ib::AbstractGridFittedBottom)
bottom_field = Field((Center, Center, Nothing), grid)
set!(bottom_field, ib.bottom_height)
fill_halo_regions!(bottom_field)
offset_bottom_array = dropdims(bottom_field.data, dims=3)
# TODO: maybe clean this up
new_ib = getnamewrapper(ib)(offset_bottom_array)
return ImmersedBoundaryGrid(grid, new_ib)
end
function ImmersedBoundaryGrid(grid, ib::AbstractGridFittedBottom{<:OffsetArray})
TX, TY, TZ = topology(grid)
validate_ib_size(grid, ib)
return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end
function validate_ib_size(grid, ib)
bottom_height_size = total_size(grid, (Center, Center, Nothing))[1:2]
size(ib.bottom_height) != bottom_height_size &&
throw(ArgumentError("The dimensions of the immersed boundary $(size(ib.bottom_height)) do not match the grid size $(bottom_height_size)"))
return nothing
end
@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition})
z = znode(i, j, k+1, underlying_grid, c, c, f)
h = @inbounds ib.bottom_height[i, j]
return z <= h
end
@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:CenterImmersedCondition})
z = znode(i, j, k, underlying_grid, c, c, c)
h = @inbounds ib.bottom_height[i, j]
return z <= h
end
@inline z_bottom(i, j, ibg::GFIBG) = @inbounds ibg.immersed_boundary.bottom_height[i, j]
on_architecture(arch, ib::GridFittedBottom) = GridFittedBottom(arch_array(arch, ib.bottom_height))
Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom_height))
#####
##### Implicit vertical diffusion
#####
#####
##### For a center solver we have to check the interface "solidity" at faces k+1 in both the
##### Upper diagonal and the Lower diagonal
##### (because of tridiagonal convention where lower_diagonal on row k is found at k-1)
##### Same goes for the face solver, where we check at centers k in both Upper and lower diagonal
#####
@inline immersed_ivd_peripheral_node(i, j, k, ibg, LX, LY, ::Center) = immersed_peripheral_node(i, j, k+1, ibg, LX, LY, Face())
@inline immersed_ivd_peripheral_node(i, j, k, ibg, LX, LY, ::Face) = immersed_peripheral_node(i, j, k, ibg, LX, LY, Center())
# Extend the upper and lower diagonal functions of the batched tridiagonal solver
for location in (:upper_, :lower_)
immersed_func = Symbol(:immersed_ivd_, location, :diagonal)
ordinary_func = Symbol(:ivd_ , location, :diagonal)
@eval begin
# Disambiguation
@inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Face, clock, Δt, κz) =
$immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
@inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Center, clock, Δt, κz) =
$immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
@inline function $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
return ifelse(immersed_ivd_peripheral_node(i, j, k, ibg, ℓx, ℓy, ℓz),
zero(eltype(ibg.underlying_grid)),
$ordinary_func(i, j, k, ibg.underlying_grid, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz))
end
end
end
#####
##### GridFittedBoundary (experimental 3D immersed boundary)
#####
struct GridFittedBoundary{M} <: AbstractGridFittedBoundary
mask :: M
end
@inline _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBoundary{<:AbstractArray}) = @inbounds ib.mask[i, j, k]
@inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBoundary)
x, y, z = node(i, j, k, underlying_grid, c, c, c)
return ib.mask(x, y, z)
end
function compute_mask(grid, ib)
mask_field = Field{Center, Center, Center}(grid, Bool)
set!(mask_field, ib.mask)
fill_halo_regions!(mask_field)
return mask_field
end
function ImmersedBoundaryGrid(grid, ib::GridFittedBoundary; precompute_mask=true)
TX, TY, TZ = topology(grid)
if precompute_mask
mask_field = compute_mask(grid, ib)
new_ib = GridFittedBoundary(mask_field)
return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib)
else
return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end
end
function ImmersedBoundaryGrid(grid, ib::GridFittedBoundary{<:OffsetArray}; kw...)
TX, TY, TZ = topology(grid)
return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end
on_architecture(arch, ib::GridFittedBoundary{<:AbstractArray}) = GridFittedBoundary(arch_array(arch, ib.mask))
on_architecture(arch, ib::GridFittedBoundary{<:Field}) = GridFittedBoundary(compute_mask(on_architecture(arch, ib.mask.grid), ib))
on_architecture(arch, ib::GridFittedBoundary) = ib # need a workaround...
Adapt.adapt_structure(to, ib::AbstractGridFittedBoundary) = GridFittedBoundary(adapt(to, ib.mask))
# fallback
immersed_cell(i, j, k, grid, ib) = _immersed_cell(i, j, k, grid, ib)
# support for Flat grids
using Oceananigans.Grids: AbstractGrid
for IB in [:GridFittedBottom, :GridFittedBoundary]
@eval begin
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, <:Any, <:Any}, ib::$IB) = _immersed_cell(1, j, k, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, <:Any, Flat, <:Any}, ib::$IB) = _immersed_cell(i, 1, k, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, <:Any, <:Any, Flat}, ib::$IB) = _immersed_cell(i, j, 1, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, Flat, <:Any}, ib::$IB) = _immersed_cell(1, 1, k, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, <:Any, Flat}, ib::$IB) = _immersed_cell(1, j, 1, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, <:Any, Flat, Flat}, ib::$IB) = _immersed_cell(i, 1, 1, grid, ib)
@inline immersed_cell(i, j, k, grid::AbstractGrid{<:Any, Flat, Flat, Flat}, ib::$IB) = _immersed_cell(1, 1, 1, grid, ib)
end
end