-
Notifications
You must be signed in to change notification settings - Fork 190
/
apply_flux_bcs.jl
160 lines (128 loc) · 6.73 KB
/
apply_flux_bcs.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
using Oceananigans: instantiated_location
using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Grids
using Oceananigans.Grids: AbstractGrid
#####
##### Algorithm for adding fluxes associated with non-trivial flux boundary conditions.
##### Value and Gradient boundary conditions are handled by filling halos.
#####
# Unpack
apply_x_bcs!(Gc, c, args...) = apply_x_bcs!(Gc, Gc.grid, c, c.boundary_conditions.west, c.boundary_conditions.east, args...)
apply_y_bcs!(Gc, c, args...) = apply_y_bcs!(Gc, Gc.grid, c, c.boundary_conditions.south, c.boundary_conditions.north, args...)
apply_z_bcs!(Gc, c, args...) = apply_z_bcs!(Gc, Gc.grid, c, c.boundary_conditions.bottom, c.boundary_conditions.top, args...)
# Shortcuts for...
#
# Nothing tendencies.
apply_x_bcs!(::Nothing, args...) = NoneEvent()
apply_y_bcs!(::Nothing, args...) = NoneEvent()
apply_z_bcs!(::Nothing, args...) = NoneEvent()
# Not-flux boundary conditions
const NotFluxBC = Union{PBC, VBC, GBC, OBC, ZFBC, Nothing}
apply_x_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = NoneEvent()
apply_y_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = NoneEvent()
apply_z_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = NoneEvent()
# The real deal
"""
Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
apply_x_bcs!(Gc, grid::AbstractGrid, c, west_bc, east_bc, arch::AbstractArchitecture, dep, args...) =
launch!(arch, grid, :yz, _apply_x_bcs!, Gc, instantiated_location(Gc), grid, west_bc, east_bc, args..., dependencies=dep)
"""
Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
apply_y_bcs!(Gc, grid::AbstractGrid, c, south_bc, north_bc, arch::AbstractArchitecture, dep, args...) =
launch!(arch, grid, :xz, _apply_y_bcs!, Gc, instantiated_location(Gc), grid, south_bc, north_bc, args..., dependencies=dep)
"""
Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the top and bottom.
"""
apply_z_bcs!(Gc, grid::AbstractGrid, c, bottom_bc, top_bc, arch::AbstractArchitecture, dep, args...) =
launch!(arch, grid, :xy, _apply_z_bcs!, Gc, instantiated_location(Gc), grid, bottom_bc, top_bc, args..., dependencies=dep)
"""
_apply_x_bcs!(Gc, grid, west_bc, east_bc, args...)
Apply a west and/or east boundary condition to variable `c`.
"""
@kernel function _apply_x_bcs!(Gc, loc, grid, west_bc, east_bc, args...)
j, k = @index(Global, NTuple)
apply_x_west_bc!(Gc, loc, west_bc, j, k, grid, args...)
apply_x_east_bc!(Gc, loc, east_bc, j, k, grid, args...)
end
"""
_apply_y_bcs!(Gc, grid, south_bc, north_bc, args...)
Apply a south and/or north boundary condition to variable `c`.
"""
@kernel function _apply_y_bcs!(Gc, loc, grid, south_bc, north_bc, args...)
i, k = @index(Global, NTuple)
apply_y_south_bc!(Gc, loc, south_bc, i, k, grid, args...)
apply_y_north_bc!(Gc, loc, north_bc, i, k, grid, args...)
end
"""
_apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
Apply a top and/or bottom boundary condition to variable `c`.
"""
@kernel function _apply_z_bcs!(Gc, loc, grid, bottom_bc, top_bc, args...)
i, j = @index(Global, NTuple)
apply_z_bottom_bc!(Gc, loc, bottom_bc, i, j, grid, args...)
apply_z_top_bc!(Gc, loc, top_bc, i, j, grid, args...)
end
# Shortcuts for zero flux or non-flux boundary conditions
@inline apply_x_east_bc!( Gc, loc, ::NotFluxBC, args...) = nothing
@inline apply_x_west_bc!( Gc, loc, ::NotFluxBC, args...) = nothing
@inline apply_y_north_bc!( Gc, loc, ::NotFluxBC, args...) = nothing
@inline apply_y_south_bc!( Gc, loc, ::NotFluxBC, args...) = nothing
@inline apply_z_top_bc!( Gc, loc, ::NotFluxBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
@inline flip(::Center) = Face()
@inline flip(::Face) = Center()
"""
apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
Add the flux divergence associated with a west flux boundary condition on `c`.
Note that because
`tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`
a positive west flux is associated with an *increase* in `Gc` near the west boundary.
If `west_bc.condition` is a function, the function must have the signature
`west_bc.condition(j, k, grid, boundary_condition_args...)`
The same logic holds for south and bottom boundary conditions in `y`, and `z`, respectively.
"""
@inline function apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[1, j, k] = getbc(west_flux, j, k, grid, args...) * Ax(1, j, k, grid, flip(LX), LY, LZ) / volume(1, j, k, grid, LX, LY, LZ)
return nothing
end
@inline function apply_y_south_bc!(Gc, loc, south_flux::BC{<:Flux}, i, k, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[i, 1, k] = getbc(south_flux, i, k, grid, args...) * Ay(i, 1, k, grid, LX, flip(LY), LZ) / volume(i, 1, k, grid, LX, LY, LZ)
return nothing
end
@inline function apply_z_bottom_bc!(Gc, loc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[i, j, 1] = getbc(bottom_flux, i, j, grid, args...) * Az(i, j, 1, grid, LX, LY, flip(LZ)) / volume(i, j, 1, grid, LX, LY, LZ)
return nothing
end
"""
apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
Add the part of flux divergence associated with a east boundary condition on `c`.
Note that because
`tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`
a positive east flux is associated with a *decrease* in `Gc` near the east boundary.
If `east_bc.condition` is a function, the function must have the signature
`east_bc.condition(i, j, grid, boundary_condition_args...)`
The same logic holds for north and top boundary conditions in `y`, and `z`, respectively.
"""
@inline function apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[grid.Nx, j, k] -= getbc(east_flux, j, k, grid, args...) * Ax(grid.Nx 1, j, k, grid, flip(LX), LY, LZ) / volume(grid.Nx, j, k, grid, LX, LY, LZ)
return nothing
end
@inline function apply_y_north_bc!(Gc, loc, north_flux::BC{<:Flux}, i, k, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[i, grid.Ny, k] -= getbc(north_flux, i, k, grid, args...) * Ay(i, grid.Ny 1, k, grid, LX, flip(LY), LZ) / volume(i, grid.Ny, k, grid, LX, LY, LZ)
return nothing
end
@inline function apply_z_top_bc!(Gc, loc, top_flux::BC{<:Flux}, i, j, grid, args...)
LX, LY, LZ = loc
@inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) * Az(i, j, grid.Nz 1, grid, LX, LY, flip(LZ)) / volume(i, j, grid.Nz, grid, LX, LY, LZ)
return nothing
end