diff --git a/docs/src/model_setup/buoyancy_and_equation_of_state.md b/docs/src/model_setup/buoyancy_and_equation_of_state.md index c5141e5ec5..d746686749 100644 --- a/docs/src/model_setup/buoyancy_and_equation_of_state.md +++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md @@ -23,11 +23,11 @@ end ```jldoctest buoyancy -julia> grid = RectilinearGrid(size=(64, 64, 64), extent=(1, 1, 1)); +julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)); julia> model = NonhydrostaticModel(; grid, buoyancy=nothing) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: () ├── closure: Nothing @@ -35,13 +35,13 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) └── coriolis: Nothing ``` -`buoyancy=nothing` is the default option for`NonhydrostaticModel`, so omitting `buoyancy` -from the `NonhydrostaticModel` constructor yields an identical result: +The option `buoyancy = nothing` is the default for [`NonhydrostaticModel`](@ref), so omitting the +`buoyancy` keyword argument from the `NonhydrostaticModel` constructor yields the same: ```jldoctest buoyancy julia> model = NonhydrostaticModel(; grid) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: () ├── closure: Nothing @@ -50,13 +50,13 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ``` To create a `HydrostaticFreeSurfaceModel` without a buoyancy term we explicitly -specify `buoyancy=nothing` flag. The default tracers `T` and `S` for `HydrostaticFreeSurfaceModel` -may be eliminated when `buoyancy=nothing` by specifying `tracers=()`: +specify `buoyancy = nothing`. The default tracers `T` and `S` for `HydrostaticFreeSurfaceModel` +may be eliminated when `buoyancy = nothing` by specifying `tracers = ()`: ```jldoctest buoyancy julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=nothing, tracers=()) HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: () ├── closure: Nothing @@ -74,24 +74,24 @@ a buoyancy tracer by including `:b` in `tracers` and specifying `buoyancy = Buo ```jldoctest buoyancy julia> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: b ├── closure: Nothing -├── buoyancy: BuoyancyTracer with -ĝ = ZDirection +├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() └── coriolis: Nothing ``` -We follow the same pattern to create a `HydrostaticFreeSurfaceModel` with buoyancy as a tracer: +Similarly for a `HydrostaticFreeSurfaceModel` with buoyancy as a tracer: ```jldoctest buoyancy julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: b ├── closure: Nothing -├── buoyancy: BuoyancyTracer with -ĝ = ZDirection +├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver └── coriolis: Nothing @@ -100,34 +100,35 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = ## Seawater buoyancy `NonhydrostaticModel` and `HydrostaticFreeSurfaceModel` support modeling the buoyancy of seawater -as a function of gravitational acceleration, conservative temperature ``T`` and absolute salinity ``S``. -The relationship between ``T``, ``S``, the geopotential height, and the density perturbation from -a reference value is called the `equation_of_state`. +as a function of the gravitational acceleration, the conservative temperature ``T``, and the absolute +salinity ``S``. The relationship between ``T``, ``S``, the geopotential height, and the density +perturbation from a reference value is called the `equation_of_state`. Specifying `buoyancy = SeawaterBuoyancy()` returns a buoyancy model with a linear equation of state, -[Earth standard](https://en.wikipedia.org/wiki/Standard_gravity) `gravitational_acceleration = 9.80665` (``\text{m}\,\text{s}^{-2}``) and requires the tracers `:T` and `:S`: +[Earth standard](https://en.wikipedia.org/wiki/Standard_gravity) `gravitational_acceleration = 9.80665` (in +S.I. units ``\text{m}\,\text{s}^{-2}``) and requires to add `:T` and `:S` as tracers: ```jldoctest buoyancy julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: (T, S) ├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection +├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() └── coriolis: Nothing ``` - -With `HydrostaticFreeSurfaceModel`, +With `HydrostaticFreeSurfaceModel`, these are the default choices for `buoyancy` and `tracers` so, +either including them or not we get: ```jldoctest buoyancy julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: (T, S) ├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection +├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver └── coriolis: Nothing @@ -138,11 +139,11 @@ is identical to the default, ```jldoctest buoyancy julia> model = HydrostaticFreeSurfaceModel(; grid) HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: (T, S) ├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection +├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() ├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² │ └── solver: FFTImplicitFreeSurfaceSolver └── coriolis: Nothing @@ -159,11 +160,11 @@ SeawaterBuoyancy{Float64}: julia> model = NonhydrostaticModel(; grid, buoyancy=buoyancy, tracers=(:T, :S)) NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: (T, S) ├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection +├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() └── coriolis: Nothing ``` @@ -183,9 +184,9 @@ SeawaterBuoyancy{Float64}: ### Idealized nonlinear equations of state -Instead of a linear equation of state, five idealized (second-order) nonlinear equation of state as described by -[Roquet15Idealized](@cite) may be used. These equations of state are provided via the -[SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package. +Instead of a linear equation of state, five idealized (second-order) nonlinear equation of state +as described by [Roquet15Idealized](@cite) may be used. These equations of state are provided +via the [SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package. ```jldoctest buoyancy julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials @@ -220,7 +221,7 @@ To simulate gravitational accelerations that don't align with the vertical (`z`) we wrap the buoyancy model in `Buoyancy()` function call, which takes the keyword arguments `model` and `gravity_unit_vector`, -```jldoctest buoyancy +```jldoctest buoyancy; filter = r".*@ Oceananigans.BuoyancyModels.*" julia> θ = 45; # degrees julia> g̃ = (0, sind(θ), cosd(θ)); @@ -228,11 +229,15 @@ julia> g̃ = (0, sind(θ), cosd(θ)); julia> model = NonhydrostaticModel(; grid, buoyancy=Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃), tracers=:b) +┌ Warning: The meaning of `gravity_unit_vector` changed in version 0.80.0. +│ In versions 0.79 and earlier, `gravity_unit_vector` indicated the direction _opposite_ to gravity. +│ In versions 0.80.0 and later, `gravity_unit_vector` indicates the direction of gravitational acceleration. +└ @ Oceananigans.BuoyancyModels ~/builds/tartarus-16/clima/oceananigans/src/BuoyancyModels/buoyancy.jl:48 NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: b ├── closure: Nothing -├── buoyancy: BuoyancyTracer with -ĝ = Tuple{Int64, Float64, Float64} +├── buoyancy: BuoyancyTracer with ĝ = Tuple{Float64, Float64, Float64} └── coriolis: Nothing ``` diff --git a/examples/tilted_bottom_boundary_layer.jl b/examples/tilted_bottom_boundary_layer.jl index 2811f574fc..487047c976 100644 --- a/examples/tilted_bottom_boundary_layer.jl +++ b/examples/tilted_bottom_boundary_layer.jl @@ -27,9 +27,9 @@ using Oceananigans -Lx = 400 # m +Lx = 200 # m Lz = 100 # m -Nx = 128 +Nx = 64 Nz = 64 ## Creates a grid with near-constant spacing `refinement * Lz / Nz` @@ -76,12 +76,12 @@ current_figure() # hide # so that ``x`` is the along-slope direction, ``z`` is the across-sloce direction that # is perpendicular to the bottom, and the unit vector anti-aligned with gravity is -ĝ = (sind(θ), 0, cosd(θ)) +ĝ = [sind(θ), 0, cosd(θ)] # Changing the vertical direction impacts both the `gravity_unit_vector` # for `Buoyancy` as well as the `rotation_axis` for Coriolis forces, -buoyancy = Buoyancy(model = BuoyancyTracer(), gravity_unit_vector = ĝ) +buoyancy = Buoyancy(model = BuoyancyTracer(), gravity_unit_vector = -ĝ) coriolis = ConstantCartesianCoriolis(f = 1e-4, rotation_axis = ĝ) # where we have used a constant Coriolis parameter ``f = 10⁻⁴ \rm{s}⁻¹``. @@ -144,7 +144,7 @@ model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure, using Oceananigans.Units using Oceananigans.Grids: min_Δz -simulation = Simulation(model, Δt = 0.5 * min_Δz(grid) / V∞, stop_time = 2days) +simulation = Simulation(model, Δt = 0.5 * min_Δz(grid) / V∞, stop_time = 1days) # We use `TimeStepWizard` to adapt our time-step and print a progress message, diff --git a/src/BuoyancyModels/buoyancy.jl b/src/BuoyancyModels/buoyancy.jl index 82a0f073b8..d59b94bb78 100644 --- a/src/BuoyancyModels/buoyancy.jl +++ b/src/BuoyancyModels/buoyancy.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: ZDirection, validate_unit_vector +using Oceananigans.Grids: NegativeZDirection, validate_unit_vector struct Buoyancy{M, G} model :: M @@ -6,23 +6,23 @@ struct Buoyancy{M, G} end """ - Buoyancy(; model, gravity_unit_vector=ZDirection()) + Buoyancy(; model, gravity_unit_vector=NegativeZDirection()) -Uses a given buoyancy `model` to create buoyancy in a model. The optional keyword argument -`gravity_unit_vector` can be used to specify the direction opposite to the gravitational -acceleration (which we take here to mean the "vertical" direction). +Construct a `buoyancy` given a buoyancy `model`. Optional keyword argument `gravity_unit_vector` +can be used to specify the direction of gravity (default `NegativeZDirection()`). +The buoyancy acceleration acts in the direction opposite to gravity. Example ======= -```jldoctest +```jldoctest; filter = r"└ @ Oceananigans.BuoyancyModels.*" using Oceananigans -grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1000, 100)) +grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1, 1)) θ = 45 # degrees -g̃ = (0, sind(θ), cosd(θ)) +g̃ = (0, sind(θ), cosd(θ)); buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃) @@ -30,28 +30,36 @@ model = NonhydrostaticModel(grid=grid, buoyancy=buoyancy, tracers=:b) # output +┌ Warning: The meaning of `gravity_unit_vector` changed in version 0.80.0. +│ In versions 0.79 and earlier, `gravity_unit_vector` indicated the direction _opposite_ to gravity. +│ In versions 0.80.0 and later, `gravity_unit_vector` indicates the direction of gravitational acceleration. +└ @ Oceananigans.BuoyancyModels ~/builds/tartarus-16/clima/oceananigans/src/BuoyancyModels/buoyancy.jl:48 NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) ├── grid: 1×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo ├── timestepper: QuasiAdamsBashforth2TimeStepper ├── tracers: b ├── closure: Nothing -├── buoyancy: BuoyancyTracer with -ĝ = Tuple{Int64, Float64, Float64} +├── buoyancy: BuoyancyTracer with ĝ = Tuple{Float64, Float64, Float64} └── coriolis: Nothing ``` """ -function Buoyancy(; model, gravity_unit_vector=ZDirection()) +function Buoyancy(; model, gravity_unit_vector=NegativeZDirection()) + gravity_unit_vector != NegativeZDirection() && + @warn """The meaning of `gravity_unit_vector` changed in version 0.80.0. + In versions 0.79 and earlier, `gravity_unit_vector` indicated the direction _opposite_ to gravity. + In versions 0.80.0 and later, `gravity_unit_vector` indicates the direction of gravitational acceleration.""" gravity_unit_vector = validate_unit_vector(gravity_unit_vector) return Buoyancy(model, gravity_unit_vector) end -@inline ĝ_x(buoyancy) = @inbounds buoyancy.gravity_unit_vector[1] -@inline ĝ_y(buoyancy) = @inbounds buoyancy.gravity_unit_vector[2] -@inline ĝ_z(buoyancy) = @inbounds buoyancy.gravity_unit_vector[3] +@inline ĝ_x(buoyancy) = @inbounds -buoyancy.gravity_unit_vector[1] +@inline ĝ_y(buoyancy) = @inbounds -buoyancy.gravity_unit_vector[2] +@inline ĝ_z(buoyancy) = @inbounds -buoyancy.gravity_unit_vector[3] -@inline ĝ_x(::Buoyancy{M, ZDirection}) where M = 0 -@inline ĝ_y(::Buoyancy{M, ZDirection}) where M = 0 -@inline ĝ_z(::Buoyancy{M, ZDirection}) where M = 1 +@inline ĝ_x(::Buoyancy{M, NegativeZDirection}) where M = 0 +@inline ĝ_y(::Buoyancy{M, NegativeZDirection}) where M = 0 +@inline ĝ_z(::Buoyancy{M, NegativeZDirection}) where M = 1 ##### ##### For convenience @@ -70,7 +78,6 @@ end regularize_buoyancy(b) = b regularize_buoyancy(b::AbstractBuoyancyModel) = Buoyancy(model=b) -Base.summary(buoyancy::Buoyancy) = string(summary(buoyancy.model), " with -ĝ = ", summary(buoyancy.gravity_unit_vector)) - -Base.show(io::IO, buoyancy::Buoyancy) = summary(buoyancy) +Base.summary(buoyancy::Buoyancy) = string(summary(buoyancy.model), " with ĝ = ", summary(buoyancy.gravity_unit_vector)) +Base.show(io::IO, buoyancy::Buoyancy) = print(io, sprint(show, buoyancy.model), "\nwith `gravity_unit_vector` = ", summary(buoyancy.gravity_unit_vector)) diff --git a/src/BuoyancyModels/g_dot_b.jl b/src/BuoyancyModels/g_dot_b.jl index 2b901cabfe..9487fc199c 100644 --- a/src/BuoyancyModels/g_dot_b.jl +++ b/src/BuoyancyModels/g_dot_b.jl @@ -2,6 +2,6 @@ @inline y_dot_g_b(i, j, k, grid, buoyancy, C) = ĝ_y(buoyancy) * buoyancy_perturbation(i, j, k, grid, buoyancy.model, C) @inline z_dot_g_b(i, j, k, grid, buoyancy, C) = ĝ_z(buoyancy) * buoyancy_perturbation(i, j, k, grid, buoyancy.model, C) -@inline x_dot_g_b(i, j, k, grid, ::Buoyancy{M, ZDirection}, C) where M = 0 -@inline y_dot_g_b(i, j, k, grid, ::Buoyancy{M, ZDirection}, C) where M = 0 -@inline z_dot_g_b(i, j, k, grid, b::Buoyancy{M, ZDirection}, C) where M = buoyancy_perturbation(i, j, k, grid, b.model, C) +@inline x_dot_g_b(i, j, k, grid, ::Buoyancy{M, NegativeZDirection}, C) where M = 0 +@inline y_dot_g_b(i, j, k, grid, ::Buoyancy{M, NegativeZDirection}, C) where M = 0 +@inline z_dot_g_b(i, j, k, grid, b::Buoyancy{M, NegativeZDirection}, C) where M = buoyancy_perturbation(i, j, k, grid, b.model, C) diff --git a/src/Coriolis/constant_cartesian_coriolis.jl b/src/Coriolis/constant_cartesian_coriolis.jl index ccf01e84be..517977a408 100644 --- a/src/Coriolis/constant_cartesian_coriolis.jl +++ b/src/Coriolis/constant_cartesian_coriolis.jl @@ -41,7 +41,7 @@ function ConstantCartesianCoriolis(FT=Float64; fx=nothing, fy=nothing, fz=nothin elseif !isnothing(f) all(isnothing.((fx, fy, fz, latitude))) || throw(ArgumentError("Only `rotation_axis` can be specified when using `f`.")) - rotation_axis = validate_unit_vector(rotation_axis) + rotation_axis = validate_unit_vector(rotation_axis, FT) if rotation_axis isa ZDirection fx = fy = 0 fz = f diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl index 0c0c0892ec..7b2fc868ab 100644 --- a/src/Grids/Grids.jl +++ b/src/Grids/Grids.jl @@ -22,7 +22,7 @@ using OffsetArrays using Oceananigans using Oceananigans.Architectures -import Base: size, length, eltype, show +import Base: size, length, eltype, show, - import Oceananigans.Architectures: architecture ##### diff --git a/src/Grids/grid_utils.jl b/src/Grids/grid_utils.jl index b156edaa72..a7a6b16250 100644 --- a/src/Grids/grid_utils.jl +++ b/src/Grids/grid_utils.jl @@ -416,7 +416,17 @@ end ##### struct ZDirection end -Base.summary(::ZDirection) = "ZDirection" + +Base.summary(::ZDirection) = "ZDirection()" +Base.show(io::IO, zdir::ZDirection) = print(io, summary(zdir)) + +struct NegativeZDirection end + +Base.summary(::NegativeZDirection) = "NegativeZDirection()" +Base.show(io::IO, zdir::NegativeZDirection) = print(io, summary(zdir)) + +-(::NegativeZDirection) = ZDirection() +-(::ZDirection) = NegativeZDirection() ##### ##### Show utils diff --git a/src/Grids/input_validation.jl b/src/Grids/input_validation.jl index 8fdc184074..23ee67bc40 100644 --- a/src/Grids/input_validation.jl +++ b/src/Grids/input_validation.jl @@ -173,17 +173,18 @@ function validate_vertically_stretched_grid_xy(TX, TY, FT, x, y) return FT(Lx), FT(Ly), FT.(x), FT.(y) end -validate_unit_vector(ê::ZDirection) = ê +validate_unit_vector(ê::ZDirection, FT::DataType=Float64) = ê +validate_unit_vector(ê::NegativeZDirection, FT::DataType=Float64) = ê -function validate_unit_vector(ê) +function validate_unit_vector(ê, FT::DataType=Float64) length(ê) == 3 || throw(ArgumentError("unit vector must have length 3")) ex, ey, ez = ê ex^2 + ey^2 + ez^2 ≈ 1 || - throw(ArgumentError("unit vector `ê` must have ê[1]² + ê[2]² + ê[3]² ≈ 1")) + throw(ArgumentError("unit vector `ê` must satisfy ê[1]² + ê[2]² + ê[3]² ≈ 1")) - return tuple(ê...) + return tuple(FT(ex), FT(ey), FT(ez)) end function validate_index(idx, loc, topo, N, H) diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index de9346544a..401b87dfb1 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -1,56 +1,58 @@ +include("dependencies_for_runtests.jl") + using Oceananigans.Coriolis: Ω_Earth using Oceananigans.Advection: EnergyConservingScheme, EnstrophyConservingScheme function instantiate_fplane_1(FT) coriolis = FPlane(FT, f=π) - return coriolis.f == FT(π) + return coriolis.f ≈ FT(π) end function instantiate_fplane_2(FT) coriolis = FPlane(FT, rotation_rate=2, latitude=30) - return coriolis.f == FT(2) + return coriolis.f ≈ FT(2) end function instantiate_constant_coriolis_1(FT) coriolis = ConstantCartesianCoriolis(FT, f=1, rotation_axis=[0, cosd(45), sind(45)]) - @test coriolis.fy == FT(cosd(45)) - @test coriolis.fz == FT(sind(45)) + @test coriolis.fy ≈ FT(cosd(45)) + @test coriolis.fz ≈ FT(sind(45)) end function instantiate_constant_coriolis_2(FT) coriolis = ConstantCartesianCoriolis(FT, f=10, rotation_axis=[√(1/3),√(1/3),√(1/3)]) - @test coriolis.fx == FT(10*√(1/3)) - @test coriolis.fy == FT(10*√(1/3)) - @test coriolis.fz == FT(10*√(1/3)) + @test coriolis.fx ≈ FT(10*√(1/3)) + @test coriolis.fy ≈ FT(10*√(1/3)) + @test coriolis.fz ≈ FT(10*√(1/3)) end function instantiate_betaplane_1(FT) coriolis = BetaPlane(FT, f₀=π, β=2π) - @test coriolis.f₀ == FT(π) - @test coriolis.β == FT(2π) + @test coriolis.f₀ ≈ FT(π) + @test coriolis.β ≈ FT(2π) end function instantiate_betaplane_2(FT) coriolis = BetaPlane(FT, latitude=70, radius=2π, rotation_rate=3π) - @test coriolis.f₀ == FT(6π * sind(70)) - @test coriolis.β == FT(6π * cosd(70) / 2π) + @test coriolis.f₀ ≈ FT(6π * sind(70)) + @test coriolis.β ≈ FT(6π * cosd(70) / 2π) end function instantiate_ntbetaplane_1(FT) coriolis = NonTraditionalBetaPlane(FT, fz=π, fy=ℯ, β=1//7, γ=5) - @test coriolis.fz == FT(π) - @test coriolis.fy == FT(ℯ) - @test coriolis.β == FT(1//7) - @test coriolis.γ == FT(5) + @test coriolis.fz ≈ FT(π) + @test coriolis.fy ≈ FT(ℯ) + @test coriolis.β ≈ FT(1//7) + @test coriolis.γ ≈ FT(5) end function instantiate_ntbetaplane_2(FT) Ω, φ, R = π, 17, ℯ coriolis = NonTraditionalBetaPlane(FT, rotation_rate=Ω, latitude=φ, radius=R) - @test coriolis.fz == FT(+2Ω*sind(φ)) - @test coriolis.fy == FT(+2Ω*cosd(φ)) - @test coriolis.β == FT(+2Ω*cosd(φ)/R) - @test coriolis.γ == FT(-4Ω*sind(φ)/R) + @test coriolis.fz ≈ FT(+2Ω*sind(φ)) + @test coriolis.fy ≈ FT(+2Ω*cosd(φ)) + @test coriolis.β ≈ FT(+2Ω*cosd(φ)/R) + @test coriolis.γ ≈ FT(-4Ω*sind(φ)/R) end function instantiate_hydrostatic_spherical_coriolis1(FT) diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl index 21993c3213..84cca293e4 100644 --- a/test/test_dynamics.jl +++ b/test/test_dynamics.jl @@ -260,8 +260,8 @@ function stratified_fluid_remains_at_rest_with_tilted_gravity_buoyancy_tracer(ar topo = (Periodic, Bounded, Bounded) grid = RectilinearGrid(arch, FT, topology=topo, size=(1, N, N), extent=(L, L, L)) - g̃ = (0, sind(θ), cosd(θ)) - buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃) + g̃ = [0, sind(θ), cosd(θ)] + buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=-g̃) y_bc = GradientBoundaryCondition(N² * g̃[2]) z_bc = GradientBoundaryCondition(N² * g̃[3])