From f8ad49607749e42af0c676f94833826030b3bca0 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 3 Oct 2019 09:21:16 -0400 Subject: [PATCH 01/18] Provides Fplane with rotation rate and latitude --- src/coriolis.jl | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 64dad24d4e..5cd9e6cd01 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -24,19 +24,47 @@ end FPlane([T=Float64;] f) Returns a parameter object for constant rotation at the angular frequency -`2f`, and therefore with background vorticity `f`, around a vertical axis. +`f`, and therefore with background vorticity `2f`, around a vertical axis. -Also called `FPlane`, after the "f-plane" approximation for the local effect of +Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ -function FPlane(T=Float64; f) +function FPlane(T=Float64; f) return FPlane{T}(f) end -@inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = +""" + FPlane([T=Float64;] Ω, lat) + +Returns a parameter object for constant rotation at the angular frequency +`f = 2Ωsin(lat)`, and therefore with background vorticity `f = 4Ωsin(lat)`, +around a vertical axis. + +Also called `FPlane`, after the "f-plane" approximation for the local effect of +Earth's rotation in a planar coordinate system tangent to the Earth's surface. +""" +function FPlane(T=Float64; Ω, lat) + return FPlane{T}(2*Ω*sind(lat)) +end + + +# """ +# βPlane{T} <: AbstractRotation +# +# A parameter object for linearly increasing rotation in meridional direction. +# """ +# struct βPlane{T} <: AbstractRotation +# f :: T +# end +# +# function βPlane(T=Float64; f₀, β) +# return βPlane{T}(f₀ + β*y) +# end + +@inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = T(0.5) * f * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) -@inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = +@inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) From 24fd4e35a99e3ae4d25abeaf035629f86742136e Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 3 Oct 2019 09:33:01 -0400 Subject: [PATCH 02/18] Adds beta-plane --- src/coriolis.jl | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 5cd9e6cd01..a4f231d902 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -48,18 +48,19 @@ function FPlane(T=Float64; Ω, lat) end -# """ -# βPlane{T} <: AbstractRotation -# -# A parameter object for linearly increasing rotation in meridional direction. -# """ -# struct βPlane{T} <: AbstractRotation -# f :: T -# end -# -# function βPlane(T=Float64; f₀, β) -# return βPlane{T}(f₀ + β*y) -# end +""" + βPlane{T} <: AbstractRotation + +A parameter object for meridionally increasing rotation (`f = f₀ + βy`). +""" +struct βPlane{T} <: AbstractRotation + f₀ :: T + β :: T +end + +function βPlane(T=Float64; f₀, β) + return βPlane{T}(f₀, β) +end @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = T(0.5) * f * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) @@ -67,7 +68,15 @@ end @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) +@inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = + T(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) + +@inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = + T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) + @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) +@inline x_f_cross_U(i, j, k, grid, rotation::βPlane, U) = -fv(i, j, k, grid, rotation.f, rotation.β, U.v) +@inline y_f_cross_U(i, j, k, grid, rotation::βPlane, U) = fu(i, j, k, grid, rotation.f, rotation.β, U.u) @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) From 691635807fa09f8d03cbf3c87a57c9f4896f7024 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 3 Oct 2019 10:51:45 -0400 Subject: [PATCH 03/18] Fixes docstrings, whitespace, uses inbounds --- src/coriolis.jl | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index a4f231d902..88b9860c97 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -24,7 +24,7 @@ end FPlane([T=Float64;] f) Returns a parameter object for constant rotation at the angular frequency -`f`, and therefore with background vorticity `2f`, around a vertical axis. +`f/2`, and therefore with background vorticity `f`, around a vertical axis. Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. @@ -34,32 +34,32 @@ function FPlane(T=Float64; f) end """ - FPlane([T=Float64;] Ω, lat) + FPlane([T=Float64;] Ω, latitude) Returns a parameter object for constant rotation at the angular frequency -`f = 2Ωsin(lat)`, and therefore with background vorticity `f = 4Ωsin(lat)`, +`Ωsin(latitude), and therefore with background vorticity `f = 2Ωsin(latitude), around a vertical axis. Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ -function FPlane(T=Float64; Ω, lat) - return FPlane{T}(2*Ω*sind(lat)) +function FPlane(T=Float64; Ω, latitude) + return FPlane{T}(2*Ω*sind(latitude)) end """ - βPlane{T} <: AbstractRotation + BetaPlane{T} <: AbstractRotation A parameter object for meridionally increasing rotation (`f = f₀ + βy`). """ -struct βPlane{T} <: AbstractRotation +struct BetaPlane{T} <: AbstractRotation f₀ :: T β :: T end -function βPlane(T=Float64; f₀, β) - return βPlane{T}(f₀, β) +function BetaPlane(T=Float64; f₀, β) + return BetaPlane{T}(f₀, β) end @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = @@ -69,14 +69,16 @@ end T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = - T(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) + T(0.5) * (f₀ + β * @inbounds(grid.yC[j])) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = - T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) + T(0.5) * (f₀ + β * @inbounds(grid.yF[j])) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) -@inline x_f_cross_U(i, j, k, grid, rotation::βPlane, U) = -fv(i, j, k, grid, rotation.f, rotation.β, U.v) -@inline y_f_cross_U(i, j, k, grid, rotation::βPlane, U) = fu(i, j, k, grid, rotation.f, rotation.β, U.u) + +@inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f, rotation.β, U.v) +@inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f, rotation.β, U.u) + @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) From 16c784332a093bcd5e856d35fcbbd88a10eda4b9 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 10 Oct 2019 11:14:41 -0400 Subject: [PATCH 04/18] Additional format for specifying beta plane --- src/Oceananigans.jl | 2 +- src/coriolis.jl | 32 +++++++++++++++++++++++++++----- 2 files changed, 28 insertions(+), 6 deletions(-) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index cbbc88914f..e8a205508e 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -12,7 +12,7 @@ export CPU, GPU, # Constants - FPlane, + FPlane, BetaPlane, second, minute, hour, day, # Grids diff --git a/src/coriolis.jl b/src/coriolis.jl index 88b9860c97..e2f24e499d 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -51,28 +51,50 @@ end """ BetaPlane{T} <: AbstractRotation -A parameter object for meridionally increasing rotation (`f = f₀ + βy`). +A parameter object for meridionally increasing Coriolis +parameter (`f = f₀ + βy`). """ struct BetaPlane{T} <: AbstractRotation f₀ :: T β :: T end +""" + BetaPlane([T=Float64;] f₀, β) + +A parameter object for meridionally increasing Coriolis +parameter (`f = f₀ + βy`). +""" function BetaPlane(T=Float64; f₀, β) return BetaPlane{T}(f₀, β) end +""" + BetaPlane([T=Float64;] Ω, latitude, R) + +Returns a parameter object for meridionally increasing rotation at the +angular frequency `Ωsin(latitude), and therefore with Coriolis parameter + `f₀ = 2Ωsin(latitude), around a vertical axis. +The Coriolis parameter increases meridionally as `f = f₀ + βy`, where +`β = 2Ωcos(latitude)/R`. +""" +function BetaPlane(T=Float64; Ω, latitude, R) + f₀ = 2*Ω*sind(latitude) + β = 2*Ω*cosd(latitude)/R + return BetaPlane{T}(f₀, β) +end + @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = T(0.5) * f * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) -@inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = - T(0.5) * (f₀ + β * @inbounds(grid.yC[j])) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) +@inbounds @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = + T(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) -@inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = - T(0.5) * (f₀ + β * @inbounds(grid.yF[j])) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) +@inbounbds @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = + T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) From 29cc97da6a8067a3c2eb0cee8dfd0eecf618165b Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 10 Oct 2019 11:25:07 -0400 Subject: [PATCH 05/18] Fixes a typo --- src/coriolis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index beb18115af..2392efe1c7 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -93,7 +93,7 @@ end @inbounds @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = T(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) -@inbounbds @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = +@inbounds @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) From a9489c3d815c37ab275a36efbbeefd678b6cbc6e Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 10 Oct 2019 11:34:35 -0400 Subject: [PATCH 06/18] Now rotation calls f0 instead of f for BetaPlane --- src/coriolis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 2392efe1c7..236b10afc7 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -99,8 +99,8 @@ end @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) -@inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f, rotation.β, U.v) -@inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f, rotation.β, U.u) +@inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f₀, rotation.β, U.v) +@inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f₀, rotation.β, U.u) @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) From d56d755a69f9880cc7f9a09dc59c5b1442a9c75c Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 10 Oct 2019 11:57:22 -0400 Subject: [PATCH 07/18] zfcrossu for betaplane --- src/coriolis.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 236b10afc7..4df4e14af2 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -29,7 +29,7 @@ Returns a parameter object for constant rotation at the angular frequency Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ -function FPlane(T::DataType=Float64; f) +function FPlane(T::DataType=Float64; f) return FPlane{T}(f) end @@ -98,9 +98,9 @@ end @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) + @inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f₀, rotation.β, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f₀, rotation.β, U.u) - -@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) - +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::BetaPlane, U) where T = zero(T) From 9cb54b1e07f3f2be72ff0b881508455bd9801820 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Thu, 10 Oct 2019 12:03:58 -0400 Subject: [PATCH 08/18] z_f_cross_U for betaplane --- src/coriolis.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 4df4e14af2..cbad4f6684 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -98,9 +98,8 @@ end @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) -@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) - @inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f₀, rotation.β, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f₀, rotation.β, U.u) -@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::BetaPlane, U) where T = zero(T) + +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::Union{FPlane, BetaPlane}, U) where T = zero(T) From 59ee18a923718841a5b1acb82cafffb66a10f2a9 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:11:59 -0400 Subject: [PATCH 09/18] Adds BetaPlane as possible coriolis option --- src/coriolis.jl | 60 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 58 insertions(+), 2 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index d4e4f1a891..6d946ee250 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -1,3 +1,5 @@ +using .TurbulenceClosures: ▶xy_cfa, ▶xy_fca + ##### ##### Functions for non-rotating models ##### @@ -39,7 +41,61 @@ end @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) -@inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) -@inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) +@inbounds @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = + +@inbounds @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = + T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) + +@inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = -fv(i, j, k, grid, coriolis.f, U.v) +@inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = fu(i, j, k, grid, coriolis.f, U.u) @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::FPlane, U) where T = zero(T) +##### +##### The Beta Plane +##### + +""" + BetaPlane{T} <: AbstractRotation + +A parameter object for meridionally increasing Coriolis +parameter (`f = f₀ + βy`). +""" +struct BetaPlane{T} <: AbstractRotation + f₀ :: T + β :: T +end + +""" + BetaPlane([T=Float64;] f₀=nothing, β=nothing, + rotation_rate=nothing, latitude=nothing, radius=nothing) + +A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`). + +The user may specify both `f₀` and `β`, or the three parameters +`rotation_rate`, `latitude`, and `radius` that specify the rotation rate and radius +of a planet, and the central latitude at which the `β`-plane approximation is to be made. +""" +function BetaPlane(T=Float64; f₀=nothing, β=nothing, + rotation_rate=nothing, latitude=nothing, radius=nothing) + + f_and_β = f₀ != nothing && β != nothing + planet_parameters = rotation_rate != nothing && latitude != nothing && radius != nothing + + (f_and_β || planet_parameters) && !(f_and_β && planet_parameters) || throw(ArgumentError( + "Either keywords f₀ and β must be specified, *or* all of rotation_rate, latitude, and radius.")) + + if planet_parameters + f₀ = 2rotation_rate * sind(latitude) + β = 2rotation_rate * cosd(latitude) / radius + end + + return BetaPlane{T}(f₀, β) +end + +@inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = + @inbounds (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, v) + +@inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = + @inbounds (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, v) + +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::BetaPlane, U) where T = zero(T) From 966ec68409a07563a000212388a48febe3fce6bc Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:12:52 -0400 Subject: [PATCH 10/18] Adds tests for instantiating beta plane and exports from top level --- src/Oceananigans.jl | 2 +- test/test_coriolis.jl | 17 +++++++++++++++-- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index 66b36566ae..437a07d42b 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -12,7 +12,7 @@ export CPU, GPU, # Constants - FPlane, + FPlane, BetaPlane, second, minute, hour, day, # Grids diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index 0317fb2e56..8612ba08a3 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -1,14 +1,27 @@ -function instantiate_coriolis(T) +function instantiate_fplane(T) coriolis = FPlane(T, f=π) return coriolis.f == T(π) end +function instantiate_betaplane_1(T) + coriolis = BetaPlane(T, f₀=π, β=2π) + return coriolis.f₀ == T(π) +end + +function instantiate_betaplane_2(T) + coriolis = BetaPlane(T, latitude=70, radius=2π, rotation_rate=3π) + return coriolis.f₀ == T(6π * sind(70)) +end + + @testset "Coriolis" begin println("Testing Coriolis...") @testset "Coriolis" begin for T in float_types - @test instantiate_coriolis(T) + @test instantiate_fplane(T) + @test instantiate_betaplane_1(T) + @test instantiate_betaplane_2(T) end end end From 5d8b480181e5da0efe05919d6d9d29c52789416a Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:28:04 -0400 Subject: [PATCH 11/18] Adds test for argument error on instatntiation of betaplane --- test/test_coriolis.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index 8612ba08a3..64a0165a92 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -13,7 +13,6 @@ function instantiate_betaplane_2(T) return coriolis.f₀ == T(6π * sind(70)) end - @testset "Coriolis" begin println("Testing Coriolis...") @@ -22,6 +21,11 @@ end @test instantiate_fplane(T) @test instantiate_betaplane_1(T) @test instantiate_betaplane_2(T) + # Non-exhaustively test that BetaPlane throws an ArgumentError + @test_throws ArgumentError BetaPlane(T, f₀=1) + @test_throws ArgumentError BetaPlane(T, β=1) + @test_throws ArgumentError BetaPlane(T, latitude=70) + @test_throws ArgumentError BetaPlane(T, f₀=1e-4, β=1e-11, latitude=70) end end end From 23e44f9fb75bc5edce1529c0b2775ce027700f56 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:28:19 -0400 Subject: [PATCH 12/18] Fixes bugs in coriolis force and beta plane arguments --- src/coriolis.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 6d946ee250..f1bdc437a8 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -81,8 +81,10 @@ function BetaPlane(T=Float64; f₀=nothing, β=nothing, f_and_β = f₀ != nothing && β != nothing planet_parameters = rotation_rate != nothing && latitude != nothing && radius != nothing - (f_and_β || planet_parameters) && !(f_and_β && planet_parameters) || throw(ArgumentError( - "Either keywords f₀ and β must be specified, *or* all of rotation_rate, latitude, and radius.")) + (f_and_β && all(Tuple(p === nothing for p in (rotation_rate, latitude, radius)))) || + (planet_parameters && all(Tuple(p === nothing for p in (f₀, β)))) || + throw(ArgumentError("Either both keywords f₀ and β must be specified, + *or* all of rotation_rate, latitude, and radius.")) if planet_parameters f₀ = 2rotation_rate * sind(latitude) @@ -93,9 +95,9 @@ function BetaPlane(T=Float64; f₀=nothing, β=nothing, end @inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = - @inbounds (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, v) + @inbounds (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, U.v) @inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = - @inbounds (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, v) + @inbounds (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, U.v) @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::BetaPlane, U) where T = zero(T) From 383e103a19e3f25806b8310d2cf271db33158d9c Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:28:36 -0400 Subject: [PATCH 13/18] Adds test for time-stepping with beta plane coriolis --- test/test_time_stepping.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index 05c284a1e2..5cfffb1098 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -5,6 +5,17 @@ function time_stepping_works(arch, FT, Closure) return true # test that no errors/crashes happen when time stepping. end +function time_stepping_works_with_beta_plane(arch, FT) + model = Model( + float_type = FT, + architecture = arch, + grid = RegularCartesianGrid(FT, (16, 16, 16), (1, 2, 3)), + coriolis = BetaPlane(FT, f₀=1e-4, β=1e-11) + ) + time_step!(model, 1, 1) + return true # test that no errors/crashes happen when time stepping. +end + function run_first_AB2_time_step_tests(arch, FT) add_ones(args...) = 1.0 model = BasicModel(N=(16, 16, 16), L=(1, 2, 3), architecture=arch, float_type=FT, @@ -151,6 +162,10 @@ Closures = (ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity, @test time_stepping_works(arch, FT, Closure) end + for arch in archs, FT in float_types + @test time_stepping_works_with_beta_plane(arch, FT) + end + @testset "2nd-order Adams-Bashforth" begin println(" Testing 2nd-order Adams-Bashforth...") for arch in archs, FT in float_types From de5764b9810dfd1a178cb3c4f04286891890092b Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 21 Oct 2019 07:30:28 -0400 Subject: [PATCH 14/18] Another bug fix in beta plane coriolis implementation --- src/coriolis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index f1bdc437a8..cb0471921f 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -95,9 +95,9 @@ function BetaPlane(T=Float64; f₀=nothing, β=nothing, end @inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = - @inbounds (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, U.v) + @inbounds - (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, U.v) @inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = - @inbounds (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, U.v) + @inbounds (coriolis.f₀ + coriolis.β * grid.yF[j]) * ▶xy_cfa(i, j, k, grid, U.u) @inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::BetaPlane, U) where T = zero(T) From efa43fcdaa122435f5203002d7a0e52afc18c48c Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Mon, 21 Oct 2019 13:21:14 -0400 Subject: [PATCH 15/18] Few minor changes --- src/coriolis.jl | 64 ++++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index cbad4f6684..b0ad28df9e 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -2,9 +2,9 @@ ##### Functions for non-rotating models ##### -@inline x_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::Nothing, U) where T = zero(T) -@inline y_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::Nothing, U) where T = zero(T) -@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::Nothing, U) where T = zero(T) +@inline x_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT) +@inline y_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT) +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U) where FT = zero(FT) ##### ##### The 'FPlane' approximation. This is equivalent to a model with a constant @@ -12,16 +12,16 @@ ##### """ - FPlane{T} <: AbstractRotation + FPlane{FT} <: AbstractRotation A parameter object for constant rotation around a vertical axis. """ -struct FPlane{T} <: AbstractRotation - f :: T +struct FPlane{FT} <: AbstractRotation + f :: FT end """ - FPlane([T=Float64;] f) + FPlane([FT=Float64;] f) Returns a parameter object for constant rotation at the angular frequency `f/2`, and therefore with background vorticity `f`, around a vertical axis. @@ -29,12 +29,12 @@ Returns a parameter object for constant rotation at the angular frequency Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ -function FPlane(T::DataType=Float64; f) - return FPlane{T}(f) +function FPlane(FT=Float64; f) + return FPlane{FT}(f) end """ - FPlane([T=Float64;] Ω, latitude) + FPlane([FT=Float64;] Ω, latitude) Returns a parameter object for constant rotation at the angular frequency `Ωsin(latitude), and therefore with background vorticity `f = 2Ωsin(latitude), @@ -43,58 +43,58 @@ around a vertical axis. Also called `FPlane`, after the "f-plane" approximation for the local effect of Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ -function FPlane(T=Float64; Ω, latitude) - return FPlane{T}(2*Ω*sind(latitude)) +function FPlane(FT=Float64; Ω, latitude) + return FPlane{FT}(2*Ω*sind(latitude)) end """ - BetaPlane{T} <: AbstractRotation + BetaPlane{FT} <: AbstractRotation A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`). """ -struct BetaPlane{T} <: AbstractRotation - f₀ :: T - β :: T +struct BetaPlane{FT} <: AbstractRotation + f₀ :: FT + β :: FT end """ - BetaPlane([T=Float64;] f₀, β) + BetaPlane([FT=Float64;] f₀, β) A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`). """ -function BetaPlane(T=Float64; f₀, β) - return BetaPlane{T}(f₀, β) +function BetaPlane(FT=Float64; f₀, β) + return BetaPlane{FT}(f₀, β) end """ - BetaPlane([T=Float64;] Ω, latitude, R) + BetaPlane([FT=Float64;] Ω, latitude, R) Returns a parameter object for meridionally increasing rotation at the angular frequency `Ωsin(latitude), and therefore with Coriolis parameter `f₀ = 2Ωsin(latitude), around a vertical axis. The Coriolis parameter increases meridionally as `f = f₀ + βy`, where -`β = 2Ωcos(latitude)/R`. +`β = 2Ωcos(latitude)/R`, and `R` is the radius of the planet. """ -function BetaPlane(T=Float64; Ω, latitude, R) +function BetaPlane(FT=Float64; Ω, latitude, R) f₀ = 2*Ω*sind(latitude) β = 2*Ω*cosd(latitude)/R - return BetaPlane{T}(f₀, β) + return BetaPlane{FT}(f₀, β) end -@inline fv(i, j, k, grid::RegularCartesianGrid{T}, f, v) where T = - T(0.5) * f * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) +@inline fv(i, j, k, grid::RegularCartesianGrid{FT}, f, v) where FT = + FT(0.5) * f * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) -@inline fu(i, j, k, grid::RegularCartesianGrid{T}, f, u) where T = - T(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) +@inline fu(i, j, k, grid::RegularCartesianGrid{FT}, f, u) where FT = + FT(0.5) * f * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) -@inbounds @inline fv(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, v) where T = - T(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) +@inbounds @inline fv(i, j, k, grid::RegularCartesianGrid{FT}, f₀, β, v) where FT = + FT(0.5) * (f₀ + β * grid.yC[j]) * (avgy_f2c(grid, v, i-1, j, k) + avgy_f2c(grid, v, i, j, k)) -@inbounds @inline fu(i, j, k, grid::RegularCartesianGrid{T}, f₀, β, u) where T = - T(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) +@inbounds @inline fu(i, j, k, grid::RegularCartesianGrid{FT}, f₀, β, u) where FT = + FT(0.5) * (f₀ + β * grid.yF[j]) * (avgx_f2c(grid, u, i, j-1, k) + avgx_f2c(grid, u, i, j, k)) @inline x_f_cross_U(i, j, k, grid, rotation::FPlane, U) = -fv(i, j, k, grid, rotation.f, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::FPlane, U) = fu(i, j, k, grid, rotation.f, U.u) @@ -102,4 +102,4 @@ end @inline x_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = -fv(i, j, k, grid, rotation.f₀, rotation.β, U.v) @inline y_f_cross_U(i, j, k, grid, rotation::BetaPlane, U) = fu(i, j, k, grid, rotation.f₀, rotation.β, U.u) -@inline z_f_cross_U(i, j, k, grid::AbstractGrid{T}, ::Union{FPlane, BetaPlane}, U) where T = zero(T) +@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Union{FPlane, BetaPlane}, U) where FT = zero(FT) From 282db77b1298615b88e9efbc7c694efc354c4f06 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Wed, 23 Oct 2019 12:09:20 -0400 Subject: [PATCH 16/18] Fixes tests --- src/coriolis.jl | 2 +- test/test_time_stepping.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 48621bbdcd..0d20139b23 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -35,7 +35,7 @@ Earth's rotation in a planar coordinate system tangent to the Earth's surface. """ function FPlane(FT=Float64; f=nothing, rotation_rate=nothing, latitude=nothing) - if f==nothing && rotation_rate != nothing && latitude != nothing + if f == nothing && rotation_rate != nothing && latitude != nothing return FPlane{FT}(2rotation_rate*sind(latitude)) elseif f != nothing && rotation_rate == nothing && latitude == nothing return FPlane{FT}(f) diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index e5be62c9c3..906c9031ec 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -9,7 +9,7 @@ function time_stepping_works_with_beta_plane(arch, FT) model = Model( float_type = FT, architecture = arch, - grid = RegularCartesianGrid(FT, (16, 16, 16), (1, 2, 3)), + grid = RegularCartesianGrid(FT, size=(16, 16, 16), length=(1, 2, 3)), coriolis = BetaPlane(FT, f₀=1e-4, β=1e-11) ) time_step!(model, 1, 1) From d84656aa04849ee2d1b87f0c7ab9650ffc88965e Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Wed, 23 Oct 2019 13:48:40 -0400 Subject: [PATCH 17/18] Adds more tests --- test/test_coriolis.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index b4b35a3cc3..9c52138873 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -4,8 +4,8 @@ function instantiate_fplane_1(T) end function instantiate_fplane_2(T) - coriolis = FPlane(T, rotation_rate=7.29e-5, latitude=45) - return coriolis.f == T(2rotation_rate*sind(latitude)) + coriolis = FPlane(T, rotation_rate=2, latitude=30) + return coriolis.f == T(2) end function instantiate_betaplane_1(T) @@ -24,9 +24,12 @@ end @testset "Coriolis" begin for T in float_types @test instantiate_fplane_1(T) + @test instantiate_fplane_2(T) @test instantiate_betaplane_1(T) @test instantiate_betaplane_2(T) # Test that FPlane throws an ArgumentError + @test_throws ArgumentError FPlane(T, rotation_rate=7e-5) + @test_throws ArgumentError FPlane(T, latitude=40) @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5) @test_throws ArgumentError FPlane(T, f=1, latitude=40) @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5, latitude=40) From eecb2976edb2517e8392cb6e8705c0b5f55f0819 Mon Sep 17 00:00:00 2001 From: suyashbire1 Date: Wed, 23 Oct 2019 13:49:47 -0400 Subject: [PATCH 18/18] Eliminates whitespace --- src/coriolis.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/coriolis.jl b/src/coriolis.jl index 0d20139b23..c14aa88b1c 100644 --- a/src/coriolis.jl +++ b/src/coriolis.jl @@ -47,10 +47,9 @@ function FPlane(FT=Float64; f=nothing, rotation_rate=nothing, latitude=nothing) end @inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = - coriolis.f * ▶xy_fca(i, j, k, grid, U.v) - @inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = coriolis.f * ▶xy_cfa(i, j, k, grid, U.u) - @inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::FPlane, U) where FT = zero(FT) + ##### ##### The Beta Plane ##### @@ -95,8 +94,6 @@ function BetaPlane(T=Float64; f₀=nothing, β=nothing, return BetaPlane{T}(f₀, β) end - - @inline x_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) = @inbounds - (coriolis.f₀ + coriolis.β * grid.yC[j]) * ▶xy_fca(i, j, k, grid, U.v) @inline y_f_cross_U(i, j, k, grid, coriolis::BetaPlane, U) =