Skip to content

Commit

Permalink
Merge pull request #699 from climate-machine/glw/forced-flow-verifica…
Browse files Browse the repository at this point in the history
…tion

More convenient BoundaryFunction functionality
  • Loading branch information
glwagner committed Mar 16, 2020
2 parents bff7abc + b83618c commit 8e3c275
Show file tree
Hide file tree
Showing 9 changed files with 228 additions and 97 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "Oceananigans"
uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
version = "0.26.0"
version = "0.27.0"

[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
Expand Down
2 changes: 1 addition & 1 deletion examples/two_dimensional_turbulence.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ nothing # hide
# `vorticity_operation` we create a field `ω` to store the result of the operation, and a
# `Computation` object for coordinate the computation of vorticity and storage in `ω`:

ω = Field(Face, Face, Cell, model.architecture, model.grid)
ω = Field(Face, Face, Cell, model.architecture, model.grid, TracerBoundaryConditions(model.grid))

vorticity_computation = Computation(vorticity_operation, ω)
nothing # hide
Expand Down
12 changes: 11 additions & 1 deletion src/BoundaryConditions/BoundaryConditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,26 @@ module BoundaryConditions

export
BCType, Flux, Gradient, Value, NoPenetration,

BoundaryCondition, bctype, getbc, setbc!,

PeriodicBoundaryCondition, NoPenetrationBoundaryCondition, NoFluxBoundaryCondition,
FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition,

CoordinateBoundaryConditions,

FieldBoundaryConditions, UVelocityBoundaryConditions, VVelocityBoundaryConditions,
WVelocityBoundaryConditions, TracerBoundaryConditions, PressureBoundaryConditions,

DiffusivityBoundaryConditions,
BoundaryFunction,

BoundaryFunction, TracerBoundaryCondition,
UVelocityBoundaryCondition, VVelocityBoundaryCondition, WVelocityBoundaryCondition,

ParameterizedBoundaryCondition, ParameterizedBoundaryConditionFunction,

apply_z_bcs!, apply_y_bcs!,

fill_halo_regions!, zero_halo_regions!

using CUDAnative
Expand Down
140 changes: 125 additions & 15 deletions src/BoundaryConditions/boundary_function.jl
Original file line number Diff line number Diff line change
@@ -1,33 +1,143 @@
"""
BoundaryFunction{B, X1, X2}(func)
BoundaryFunction{B, X1, X2}(func, parameters=nothing)
A wrapper for user-defined boundary condition functions on the
boundary specified by symbol `B` and at location `(X1, X2)`.
A wrapper for the user-defined boundary condition function `func`, on the
boundary specified by symbol `B` and at location `(X1, X2)`, and with `parameters`.
Example
=======
julia> using Oceananigans: BoundaryCondition, BoundaryFunction, Flux, Cell
Example
=======
julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
(::BoundaryFunction{:z,Cell,Cell,getfield(Main, Symbol("##7#8"))}) (generic function with 1 method)
julia> using Oceananigans, Oceananigans.BoundaryConditions, Oceananigans.Fields
julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);
julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
(::BoundaryFunction{:z,Cell,Cell,var"#7#8",Nothing}) (generic function with 1 method)
julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);
julia> flux_func(x, y, t, p) = cos(p.k * x) * cos(p.ω * t); # function with parameters
julia> parameterized_u_velocity_flux = BoundaryFunction{:z, Face, Cell}(flux_func, (k=4π, ω=3.0))
(::BoundaryFunction{:z,Face,Cell,typeof(flux_func),NamedTuple{(:k, :ω),Tuple{Float64,Float64}}}) (generic function with 1 method)
julia> top_u_bc = BoundaryCondition(Flux, parameterized_u_velocity_flux);
"""
struct BoundaryFunction{B, X1, X2, F} <: Function
struct BoundaryFunction{B, X1, X2, F, P} <: Function
func :: F
parameters :: P

function BoundaryFunction{B, X1, X2}(func) where {B, X1, X2}
function BoundaryFunction{B, X1, X2}(func, parameters=nothing) where {B, X1, X2}
B (:x, :y, :z) || throw(ArgumentError("The boundary B at which the BoundaryFunction is
to be applied must be either :x, :y, or :z."))
return new{B, X1, X2, typeof(func)}(func)
return new{B, X1, X2, typeof(func), typeof(parameters)}(func, parameters)
end
end

@inline call_boundary_function(func, ξ, η, t, ::Nothing) = func(ξ, η, t)
@inline call_boundary_function(func, ξ, η, t, parameters) = func(ξ, η, t, parameters)

@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, clock, state) where {Y, Z} =
bc.func(ynode(Y, j, grid), znode(Z, k, grid), clock.time)
call_boundary_function(bc.func, ynode(Y, j, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, clock, state) where {X, Z} =
bc.func(xnode(X, i, grid), znode(Z, k, grid), clock.time)
call_boundary_function(bc.func, xnode(X, i, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, clock, state) where {X, Y} =
bc.func(xnode(X, i, grid), ynode(Y, j, grid), clock.time)
call_boundary_function(bc.func, xnode(X, i, grid), ynode(Y, j, grid), clock.time, bc.parameters)

#####
##### Convenience constructors
#####

""" Returns the location of `loc` on the boundary `:x`, `:y`, or `:z`. """
boundary_location(::Val{:x}, loc) = (loc[2], loc[3])
boundary_location(::Val{:y}, loc) = (loc[1], loc[3])
boundary_location(::Val{:z}, loc) = (loc[1], loc[2])

""" Returns a boundary function at on the boundary `B` at the appropriate
location for a tracer field. """
TracerBoundaryFunction(B, args...) = BoundaryFunction{B, Cell, Cell}(args...)

""" Returns a boundary function at on the boundary `B` at the appropriate
location for u, the x-velocity field. """
function UVelocityBoundaryFunction(B, args...)
loc = boundary_location(Val(B), (Face, Cell, Cell))
return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

""" Returns a boundary function at on the boundary `B` at the appropriate
location for v, the y-velocity field. """
function VVelocityBoundaryFunction(B, args...)
loc = boundary_location(Val(B), (Cell, Face, Cell))
return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

""" Returns a boundary function at on the boundary `B` at the appropriate
location for w, the z-velocity field. """
function WVelocityBoundaryFunction(B, args...)
loc = boundary_location(Val(B), (Cell, Cell, Face))
return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

"""
TracerBoundaryCondition(bctype, B, args...)
Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to a tracer on the boundary `B`, which is one of `:x, :y, :z`.
The boundary function has the signature
`func(ξ, η, t)`
where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
TracerBoundaryCondition(bctype, B, args...) =
BoundaryCondition(bctype, TracerBoundaryFunction(B, args...))

"""
UVelocityBoundaryCondition(bctype, B, args...)
Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `u`, the `x`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature
`func(ξ, η, t)`
where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
UVelocityBoundaryCondition(bctype, B, args...) =
BoundaryCondition(bctype, UVelocityBoundaryFunction(B, args...))

"""
VVelocityBoundaryCondition(bctype, B, args...)
Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `v`, the `y`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature
`func(ξ, η, t)`
where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
VVelocityBoundaryCondition(bctype, B, args...) =
BoundaryCondition(bctype, VVelocityBoundaryFunction(B, args...))

"""
VVelocityBoundaryCondition(bctype, B, args...)
Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `w`, the `z`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature
`func(ξ, η, t)`
where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
WVelocityBoundaryCondition(bctype, B, args...) =
BoundaryCondition(bctype, WVelocityBoundaryFunction(B, args...))
16 changes: 8 additions & 8 deletions src/BoundaryConditions/field_boundary_conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@ FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

DefaultBoundaryCondition(::Union{Grids.Periodic, Flat}, loc) = PeriodicBoundaryCondition()

DefaultBoundaryCondition(::Bounded, ::Type{Cell}) = NoFluxBoundaryCondition()
DefaultBoundaryCondition(::Bounded, ::Type{Face}) = NoPenetrationBoundaryCondition()
DefaultBoundaryCondition(::Bounded, ::Cell) = NoFluxBoundaryCondition()
DefaultBoundaryCondition(::Bounded, ::Face) = NoPenetrationBoundaryCondition()

function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir)
if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc)
Expand Down Expand Up @@ -69,9 +69,9 @@ function FieldBoundaryConditions(grid, loc;
return FieldBoundaryConditions(x, y, z)
end

UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face, Cell, Cell); user_defined_bcs...)
VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Face, Cell); user_defined_bcs...)
WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Face); user_defined_bcs...)
TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)
PressureBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)
DiffusivityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)
UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face(), Cell(), Cell()); user_defined_bcs...)
VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Face(), Cell()); user_defined_bcs...)
WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Face()); user_defined_bcs...)
TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)
PressureBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)
DiffusivityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell(), Cell(), Cell()); user_defined_bcs...)
115 changes: 53 additions & 62 deletions src/Fields/field.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,110 +49,101 @@ struct Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}
end

"""
Field(X, Y, Z, arch, grid, [ bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
data = zeros(arch, grid, (X, Y, Z)) ] )
Construct a `Field` on `grid` with `data` on architecture `arch` with
boundary conditions `bcs`. Each of `(X, Y, Z)` is either `Cell` or `Face` and determines
the field's location in `(x, y, z)`.
Example
=======
julia> ω = Field(Face, Face, Cell, CPU(), RegularCartesianmodel.grid)
Field(L::Tuple, arch, grid, bcs, [data=zeros(arch, grid)])
Construct a `Field` on some architecture `arch` and a `grid` with some `data`.
The field's location is defined by a tuple `L` of length 3 whose elements are
`Cell` or `Face` and has field boundary conditions `bcs`.
"""
function Field(X, Y, Z, arch, grid,
bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
data = zeros(eltype(grid), arch, grid, (X, Y, Z)))

return Field{X, Y, Z}(data, grid, bcs)
function Field(L::Tuple, arch, grid, bcs,
data=zeros(eltype(grid), arch, grid, (typeof(L[1]), typeof(L[2]), typeof(L[3]))))

return Field{typeof(L[1]), typeof(L[2]), typeof(L[3])}(data, grid, bcs)
end

#####
##### Convenience constructor for Field that uses a 3-tuple of locations rather than a list of locations:
#####
function Field(L::NTuple{3, DataType}, arch, grid, bcs,
data=zeros(eltype(grid), arch, grid, (L[1], L[2], L[3])))

# Type "destantiation": convert Face() to Face and Cell() to Cell if needed.
destantiate(X) = typeof(X)
destantiate(X::DataType) = X
return Field{L[1], L[2], L[3]}(data, grid, bcs)
end

"""
Field(L::Tuple, arch, grid, data, bcs)
Field(X, Y, Z, arch, grid, [data=zeros(arch, grid)], bcs)
Construct a `Field` at the location defined by the 3-tuple `L`,
whose elements are `Cell` or `Face`.
Construct a `Field` on some architecture `arch` and a `grid` with some `data`.
The field's location is defined by `X`, `Y`, `Z` where each is either `Cell` or `Face`
and has field boundary conditions `bcs`.
"""
Field(L::Tuple, args...) = Field(destantiate.(L)..., args...)

#####
##### Special constructors for tracers and velocity fields
#####
Field(X, Y, Z, arch, grid, bcs, data=zeros(eltype(grid), arch, grid, (X, Y, Z))) =
Field((X, Y, Z), arch, grid, bcs, data)

"""
CellField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
[ bcs = TracerBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Cell, Cell) ] )
CellField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=TracerBoundaryConditions(grid),
data=zeros(FT, arch, grid, (Cell, Cell, Cell)))
Return a `Field{Cell, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function CellField(FT::DataType, arch, grid,
bcs = TracerBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Cell, Cell)))
function CellField(arch::AbstractArchitecture, grid,
bcs=TracerBoundaryConditions(grid),
data=zeros(eltype(grid), arch, grid, (Cell, Cell, Cell)))

return Field{Cell, Cell, Cell}(data, grid, bcs)
return Field(Cell, Cell, Cell, arch, grid, bcs, data)
end

"""
XFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
[ bcs = UVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Face, Cell, Cell) ] )
XFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=UVelocityBoundaryConditions(grid),
data=zeros(FT, arch, grid, (Face, Cell, Cell)))
Return a `Field{Face, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function XFaceField(FT::DataType, arch, grid,
bcs = UVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Face, Cell, Cell)))
function XFaceField(arch::AbstractArchitecture, grid,
bcs=UVelocityBoundaryConditions(grid),
data=zeros(eltype(grid), arch, grid, (Face, Cell, Cell)))

return Field{Face, Cell, Cell}(data, grid, bcs)
return Field(Face, Cell, Cell, arch, grid, bcs, data)
end

"""
YFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
[ bcs = VVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Face, Cell)) ] )
YFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=VVelocityBoundaryConditions(grid),
data=zeros(FT, arch, grid, (Cell, Face, Cell)))
Return a `Field{Cell, Face, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function YFaceField(FT::DataType, arch, grid,
bcs = VVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Face, Cell)))
function YFaceField(arch::AbstractArchitecture, grid,
bcs=VVelocityBoundaryConditions(grid),
data=zeros(eltype(grid), arch, grid, (Cell, Face, Cell)))

return Field{Cell, Face, Cell}(data, grid, bcs)
return Field(Cell, Face, Cell, arch, grid, bcs, data)
end

"""
ZFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
[ bcs = WVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Cell, Face)) ] )
ZFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=WVelocityBoundaryConditions(grid),
data=zeros(FT, arch, grid, (Cell, Cell, Face))
Return a `Field{Cell, Cell, Face}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function ZFaceField(FT::DataType, arch, grid,
bcs = WVelocityBoundaryConditions(grid),
data = zeros(FT, arch, grid, (Cell, Cell, Face)))
function ZFaceField(arch::AbstractArchitecture, grid,
bcs=WVelocityBoundaryConditions(grid),
data=zeros(eltype(grid), arch, grid, (Cell, Cell, Face)))

return Field{Cell, Cell, Face}(data, grid, bcs)
return Field(Cell, Cell, Face, arch, grid, bcs, data)
end

CellField(arch::AbstractArchitecture, grid, args...) = CellField(eltype(grid), arch, grid, args...)
XFaceField(arch::AbstractArchitecture, grid, args...) = XFaceField(eltype(grid), arch, grid, args...)
YFaceField(arch::AbstractArchitecture, grid, args...) = YFaceField(eltype(grid), arch, grid, args...)
ZFaceField(arch::AbstractArchitecture, grid, args...) = ZFaceField(eltype(grid), arch, grid, args...)
CellField(FT::DataType, arch, grid, bcs=TracerBoundaryConditions(grid)) =
CellField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Cell, Cell)))

XFaceField(FT::DataType, arch, grid, bcs=UVelocityBoundaryConditions(grid)) =
XFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Face, Cell, Cell)))

YFaceField(FT::DataType, arch, grid, bcs=VVelocityBoundaryConditions(grid)) =
YFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Face, Cell)))

ZFaceField(FT::DataType, arch, grid, bcs=WVelocityBoundaryConditions(grid)) =
ZFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Cell, Face)))

#####
##### Functions for querying fields
Expand Down
Loading

0 comments on commit 8e3c275

Please sign in to comment.