Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FINALLY reverse the k index #462

Merged
merged 36 commits into from
Oct 22, 2019
Merged
Show file tree
Hide file tree
Changes from 32 commits
Commits
Show all changes
36 commits
Select commit Hold shift + click to select a range
ed46bb2
Reverse index for operators
ali-ramadhan Sep 20, 2019
0a929bd
Reverse index for turbulence operators
ali-ramadhan Sep 20, 2019
826d302
Reverse index for halos
ali-ramadhan Sep 21, 2019
9ec27d3
Reverse index for compute_w_from_continuity
ali-ramadhan Sep 21, 2019
534c4c1
Reverse zC and zF
ali-ramadhan Sep 21, 2019
d362be6
Fix grid tests
ali-ramadhan Oct 15, 2019
6f47019
Model is incompressible again
ali-ramadhan Oct 15, 2019
8e26cf8
Update k indexing on turbulence closure tests
ali-ramadhan Oct 15, 2019
2343bdf
Add test for grid range lengths
ali-ramadhan Oct 15, 2019
79d055d
Update variable names in ancient grid tests
ali-ramadhan Oct 15, 2019
43bd7e9
Add test for roundoff error in ranges that fails
ali-ramadhan Oct 15, 2019
b08f6fa
Fix grid ranges.
ali-ramadhan Oct 15, 2019
52d85ea
Rename Pearson vortex to Taylor-Green vortex
ali-ramadhan Oct 15, 2019
05fb3ac
Update halo filling for no-penetration boundary conditions.
ali-ramadhan Oct 16, 2019
1377c34
Integrate from bottom upwards to recompute w from continuity.
ali-ramadhan Oct 16, 2019
d1c906f
Fix indexing on turbulence closure tests.
ali-ramadhan Oct 16, 2019
7c1fa1c
Perform buoyancy integral from the surface downwards.
ali-ramadhan Oct 16, 2019
b1542c9
Only have to zero out `div_u` at `k=Nz`.
ali-ramadhan Oct 16, 2019
f462471
Generate all ranges with function
ali-ramadhan Oct 17, 2019
9a8aa2f
Correct buoyancy integral
ali-ramadhan Oct 17, 2019
bc049e1
Reverse output for rising thermal bubble regression test.
ali-ramadhan Oct 17, 2019
20be6f4
Remapping `:bottom -> :left` and `:top -> :right`.
ali-ramadhan Oct 18, 2019
0445cd4
Reverse k index and remap top/bottom for filling VBC and GBC halos
ali-ramadhan Oct 18, 2019
5ea35a3
Reverse output for Rayleigh-Benard convection regression test.
ali-ramadhan Oct 18, 2019
0de2765
Only compare `k = 2:Nz` as `w[:, :, 1] = 0`.
ali-ramadhan Oct 18, 2019
fb8d906
Fix bottom/top alias when using `CoordinateBoundaryConditions`.
ali-ramadhan Oct 18, 2019
8f09f58
Fix flux boundary conditions after reversing k index.
ali-ramadhan Oct 18, 2019
f0a91d6
Merge branch 'master' into reverse-k-index
ali-ramadhan Oct 18, 2019
a8b5a9a
Update old regression tests and fix post-merge typo.
ali-ramadhan Oct 18, 2019
5dff125
Reverse output for LES regression tests.
ali-ramadhan Oct 18, 2019
6ee988b
Fix bottom/top mixup in applying flux boundary conditions.
ali-ramadhan Oct 21, 2019
ed452ff
`return` -> `return nothing`
ali-ramadhan Oct 21, 2019
861a0ea
Switch `_fill_bottom_halo!` and `_fill_top_halo!` where appropriate.
ali-ramadhan Oct 21, 2019
4170bbd
Merge branch 'master' into reverse-k-index
ali-ramadhan Oct 21, 2019
41575ad
Need to export NetCDFOutputWriter.
ali-ramadhan Oct 22, 2019
419f176
Fix thermal bubble regression test by directly using NCDatasets.jl
ali-ramadhan Oct 22, 2019
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 19 additions & 20 deletions src/Operators/ops_regular_cartesian_grid.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,10 @@ using Oceananigans: AbstractGrid, RegularCartesianGrid
@inline δy_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j+1, k] - f[i, j, k]
@inline δy_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, j-1, k]

@inline δz_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k-1] - f[i, j, k]
@inline δz_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, j, k+1]
@inline δz_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, j, k+1]
@inline δz_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k-1] - f[i, j, k]
@inline δz_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, j, k-1]
@inline δz_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j, k]
@inline δz_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j, k]
@inline δz_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, j, k-1]

@inline avgx_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i-1, j, k])
@inline avgx_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i+1, j, k] + f[i, j, k])
Expand All @@ -23,9 +23,9 @@ using Oceananigans: AbstractGrid, RegularCartesianGrid
@inline avgy_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j+1, k] + f[i, j, k])
@inline avgy_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j-1, k])

@inline avgz_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])
@inline avgz_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k+1] + f[i, j, k])
@inline avgz_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])
@inline avgz_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])
@inline avgz_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k+1] + f[i, j, k])
@inline avgz_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])

@inline ∂x_p(i, j, k, grid::RegularCartesianGrid, p) = δx_c2f(grid, p, i, j, k) / grid.Δx
@inline ∂y_p(i, j, k, grid::RegularCartesianGrid, p) = δy_c2f(grid, p, i, j, k) / grid.Δy
Expand All @@ -48,12 +48,12 @@ end
end

@inline function δz_f2c_ab̄ᶻ(g::RegularCartesianGrid, a, b, i, j, k)
@inbounds (a[i, j, k] * avgz_c2f(g, b, i, j, k) -
a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1))
@inbounds (a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1) -
a[i, j, k] * avgz_c2f(g, b, i, j, k))
end

@inline div_flux(g::RegularCartesianGrid, u, v, w, Q, i, j, k) =
(δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)
(δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)

@inline δx_c2f_ūˣūˣ(g::RegularCartesianGrid, u, i, j, k) =
avgx_f2c(g, u, i, j, k)^2 - avgx_f2c(g, u, i-1, j, k)^2
Expand All @@ -64,8 +64,8 @@ end
end

@inline function δz_e2f_w̄ˣūᶻ(g::RegularCartesianGrid, u, w, i, j, k)
@inbounds avgx_f2e(g, w, i, j, k) * avgz_f2e(g, u, i, j, k) -
avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1)
@inbounds avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1) -
avgx_f2e(g, w, i, j, k) * avgz_f2e(g, u, i, j, k)
end

@inline u∇u(g::RegularCartesianGrid, u, v, w, i, j, k) =
Expand All @@ -80,8 +80,8 @@ end
avgy_f2c(g, v, i, j, k)^2 - avgy_f2c(g, v, i, j-1, k)^2

@inline function δz_e2f_w̄ʸv̄ᶻ(g::RegularCartesianGrid, v, w, i, j, k)
@inbounds avgy_f2e(g, w, i, j, k) * avgz_f2e(g, v, i, j, k) -
avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1)
@inbounds avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1) -
avgy_f2e(g, w, i, j, k) * avgz_f2e(g, v, i, j, k)
end

@inline u∇v(g::RegularCartesianGrid, u, v, w, i, j, k) =
Expand All @@ -98,26 +98,25 @@ end
end

@inline δz_c2f_w̄ᶻw̄ᶻ(g::RegularCartesianGrid, w, i, j, k) =
avgz_f2c(g, w, i, j, k-1)^2 - avgz_f2c(g, w, i, j, k)^2
avgz_f2c(g, w, i, j, k)^2 - avgz_f2c(g, w, i, j, k-1)^2

@inline u∇w(g::RegularCartesianGrid, u, v, w, i, j, k) =
(δx_e2f_ūᶻw̄ˣ(g, u, w, i, j, k) / g.Δx) + (δy_e2f_v̄ᶻw̄ʸ(g, v, w, i, j, k) / g.Δy) + (δz_c2f_w̄ᶻw̄ᶻ(g, w, i, j, k) / g.Δz)

@inline δx²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δx_c2f(g, f, i+1, j, k) - δx_c2f(g, f, i, j, k)
@inline δy²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δy_c2f(g, f, i, j+1, k) - δy_c2f(g, f, i, j, k)
@inline δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δz_c2f(g, f, i, j, k) - δz_c2f(g, f, i, j, k+1)
@inline δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δz_c2f(g, f, i, j, k+1) - δz_c2f(g, f, i, j, k)

@inline κ∇²(g::RegularCartesianGrid, Q, κh, κv, i, j, k) =
((κh/g.Δx^2) * δx²_c2f2c(g, Q, i, j, k)) + ((κh/g.Δy^2) * δy²_c2f2c(g, Q, i, j, k)) + ((κv/g.Δz^2) * δz²_c2f2c(g, Q, i, j, k))

@inline δx²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2c(g, f, i, j, k) - δx_f2c(g, f, i-1, j, k)
@inline δy²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2c(g, f, i, j, k) - δy_f2c(g, f, i, j-1, k)
@inline δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2c(g, f, i, j, k) - δz_f2c(g, f, i, j, k-1)

@inline δx²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2e(g, f, i+1, j, k) - δx_f2e(g, f, i, j, k)
@inline δy²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2e(g, f, i, j+1, k) - δy_f2e(g, f, i, j, k)

@inline δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2e(g, f, i, j, k) - δz_f2e(g, f, i, j, k+1)
@inline δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2c(g, f, i, j, k-1) - δz_f2c(g, f, i, j, k)
@inline δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2e(g, f, i, j, k+1) - δz_f2e(g, f, i, j, k)

@inline 𝜈∇²u(g::RegularCartesianGrid, u, 𝜈h, 𝜈v, i, j, k) =
((𝜈h/g.Δx^2) * δx²_f2c2f(g, u, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, u, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, u, i, j, k))
Expand All @@ -129,6 +128,6 @@ end
((𝜈h/g.Δx^2) * δx²_f2e2f(g, w, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, w, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2c2f(g, w, i, j, k))

@inline ∇²(g::RegularCartesianGrid, f, i, j, k) =
(δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)
(δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)

@inline ∇h_u(i, j, k, grid, u, v) = δx_f2c(grid, u, i, j, k) / grid.Δx + δy_f2c(grid, v, i, j, k) / grid.Δy
8 changes: 4 additions & 4 deletions src/TurbulenceClosures/closure_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@
@inline ∂y_aca(i, j, k, grid, v::AbstractArray) = @inbounds (v[i, j+1, k] - v[i, j, k]) / grid.Δy
@inline ∂y_afa(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i, j-1, k]) / grid.Δy

@inline ∂z_aac(i, j, k, grid, w::AbstractArray) = @inbounds (w[i, j, k] - w[i, j, k+1]) / grid.Δz
@inline ∂z_aaf(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k-1] - c[i, j, k]) / grid.Δz
@inline ∂z_aac(i, j, k, grid, w::AbstractArray) = @inbounds (w[i, j, k+1] - w[i, j, k]) / grid.Δz
@inline ∂z_aaf(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i, j, k-1]) / grid.Δz

#
# Differentiation and interpolation operators for functions
Expand Down Expand Up @@ -108,7 +108,7 @@ Differentiate the function or callable object
located at `aac` in `z`, across `aaf`.
"""
@inline ∂z_aaf(i, j, k, grid::AbstractGrid, F::TF, args...) where TF<:Function =
(F(i, j, k-1, grid, args...) - F(i, j, k, grid, args...)) / grid.Δz
(F(i, j, k, grid, args...) - F(i, j, k-1, grid, args...)) / grid.Δz

"""
∂z_aac(i, j, k, grid, F, args...)
Expand All @@ -120,7 +120,7 @@ Differentiate the function or callable object
located at `aaf` in `z`, across `aac`.
"""
@inline ∂z_aac(i, j, k, grid, F::TF, args...) where TF<:Function =
(F(i, j, k, grid, args...) - F(i, j, k+1, grid, args...)) / grid.Δz
(F(i, j, k+1, grid, args...) - F(i, j, k, grid, args...)) / grid.Δz

#####
##### Double differentiation
Expand Down
47 changes: 22 additions & 25 deletions src/boundary_conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -139,20 +139,17 @@ const CBC = CoordinateBoundaryConditions
PeriodicBCs() = CBC(PeriodicBC(), PeriodicBC())

# Here we overload setproperty! and getproperty to permit users to call
# the 'right' and 'left' bcs in the z-direction 'bottom' and 'top'.
#
# Note that 'right' technically corresponds to face point N+1. Thus
# the fact that right == bottom is associated with the reverse z-indexing
# convention. With ordinary indexing, right == top.
# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top'.
# Note that 'right' technically corresponds to face point N+1.
Base.setproperty!(cbc::CBC, side::Symbol, bc) = setbc!(cbc, Val(side), bc)
setbc!(cbc::CBC, ::Val{S}, bc) where S = setfield!(cbc, S, bc)
setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :right, bc)
setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :left, bc)
setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :left, bc)
setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :right, bc)

Base.getproperty(cbc::CBC, side::Symbol) = getbc(cbc, Val(side))
getbc(cbc::CBC, ::Val{S}) where S = getfield(cbc, S)
getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :right)
getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :right)

#####
##### Boundary conditions for Fields
Expand Down Expand Up @@ -189,7 +186,7 @@ function HorizontallyPeriodicBCs(; top = BoundaryCondition(Flux, nothing),

x = PeriodicBCs()
y = PeriodicBCs()
z = CoordinateBoundaryConditions(top, bottom)
z = CoordinateBoundaryConditions(bottom, top)

return FieldBoundaryConditions(x, y, z)
end
Expand All @@ -214,7 +211,7 @@ function ChannelBCs(; north = BoundaryCondition(Flux, nothing),

x = PeriodicBCs()
y = CoordinateBoundaryConditions(south, north)
z = CoordinateBoundaryConditions(top, bottom)
z = CoordinateBoundaryConditions(bottom, top)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very satisfying.


return FieldBoundaryConditions(x, y, z)
end
Expand All @@ -228,12 +225,12 @@ default_tracer_bcs(tracers, solution_bcs) = DefaultTracerBoundaryConditions(solu
"""
SolutionBoundaryConditions(tracers, proposal_bcs)

Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with
fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions
`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w`
Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with
fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions
`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w`
and may contain some or all of the boundary conditions on `tracers`.
"""
SolutionBoundaryConditions(tracers, proposal_bcs) =
SolutionBoundaryConditions(tracers, proposal_bcs) =
with_tracers(tracers, proposal_bcs, default_tracer_bcs, with_velocities=true)

"""
Expand Down Expand Up @@ -282,7 +279,7 @@ const BoundaryConditions = HorizontallyPeriodicSolutionBCs
#####
##### Tracer, tendency and pressure boundary condition "translators":
#####
##### * Default boundary conditions on tracers are periodic or no flux and
##### * Default boundary conditions on tracers are periodic or no flux and
##### can be derived from boundary conditions on any field
#####
##### * Boundary conditions on tendency terms are
Expand All @@ -296,7 +293,7 @@ const BoundaryConditions = HorizontallyPeriodicSolutionBCs
DefaultTracerBC(::BC) = BoundaryCondition(Flux, nothing)
DefaultTracerBC(::PBC) = PeriodicBC()

DefaultTracerCoordinateBCs(bcs) =
DefaultTracerCoordinateBCs(bcs) =
CoordinateBoundaryConditions(DefaultTracerBC(bcs.left), DefaultTracerBC(bcs.right))

DefaultTracerBoundaryConditions(field_bcs) =
Expand All @@ -310,7 +307,7 @@ TendencyBC(::NPBC) = NoPenetrationBC()
TendencyCoordinateBCs(bcs) =
CoordinateBoundaryConditions(TendencyBC(bcs.left), TendencyBC(bcs.right))

TendencyFieldBCs(field_bcs) =
TendencyFieldBCs(field_bcs) =
FieldBoundaryConditions(Tuple(TendencyCoordinateBCs(bcs) for bcs in field_bcs))

TendenciesBoundaryConditions(solution_bcs) =
Expand Down Expand Up @@ -365,7 +362,7 @@ the source term `Gc` at the top and bottom.
"""
function apply_z_bcs!(Gc, arch, grid, top_bc, bottom_bc, args...)
@launch device(arch) config=launch_config(grid, 2) _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
return
return nothing
end

# Fall back functions for boundary conditions that are not of type Flux.
Expand All @@ -384,13 +381,13 @@ Note that because

`tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive top flux is associated with a *decrease* in `Gc` near the top boundary.
a positive top flux is associated with an *increase* in `Gc` near the top boundary.
If `top_bc.condition` is a function, the function must have the signature

`top_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...) =
@inbounds Gc[i, j, 1] -= getbc(top_flux, i, j, grid, args...) / grid.Δz
@inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) / grid.Δz

"""
apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
Expand All @@ -400,24 +397,24 @@ Note that because

`tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive bottom flux is associated with an *increase* in `Gc` near the bottom boundary.
a positive bottom flux is associated with a *decrease* in `Gc` near the bottom boundary.
If `bottom_bc.condition` is a function, the function must have the signature

`bottom_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) =
@inbounds Gc[i, j, grid.Nz] += getbc(bottom_flux, i, j, grid, args...) / grid.Δz
@inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / grid.Δz

"""
_apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)

Apply a top and/or bottom boundary condition to variable `c`.
"""
function _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
function _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
@loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
@loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
apply_z_top_bc!(Gc, top_bc, i, j, grid, args...)
apply_z_bottom_bc!(Gc, bottom_bc, i, j, grid, args...)
apply_z_top_bc!(Gc, top_bc, i, j, grid, args...)
end
end
end
4 changes: 2 additions & 2 deletions src/diagnostics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ function run_diagnostic(model, havg::HorizontalAverage{NTuple{1}})
zero_halo_regions!(havg.fields[1], model.grid)
sum!(havg.profile, havg.fields[1])
normalize_horizontal_sum!(havg, model.grid)
return
return nothing
end

function run_diagnostic(model::Model, havg::HorizontalAverage)
Expand All @@ -93,7 +93,7 @@ function run_diagnostic(model::Model, havg::HorizontalAverage)
sum!(havg.profile, tmp)
normalize_horizontal_sum!(havg, model.grid)

return
return nothing
end

function (havg::HorizontalAverage{F, Nothing})(model) where F
Expand Down
12 changes: 6 additions & 6 deletions src/grids.jl
Original file line number Diff line number Diff line change
Expand Up @@ -106,13 +106,13 @@ function RegularCartesianGrid(T, N, L)

V = Δx*Δy*Δz

xC = Δx/2:Δx:Lx
yC = Δy/2:Δy:Ly
zC = -Δz/2:-Δz:-Lz
xC = range(Δx/2, Lx-Δx/2; length=Nx)
yC = range(Δy/2, Ly-Δy/2; length=Ny)
zC = range(-Lz+Δz/2, -Δz/2; length=Nz)

xF = 0:Δx:Lx
yF = 0:Δy:Ly
zF = 0:-Δz:-Lz
xF = range(0, Lx; length=Nx+1)
yF = range(0, Ly; length=Ny+1)
zF = range(-Lz, 0; length=Nz+1)

RegularCartesianGrid{T, typeof(xC)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz,
Lx, Ly, Lz, Δx, Δy, Δz, Ax, Ay, Az, V,
Expand Down
Loading