Flips sign for gravity_unit_vector to match its description (attempt …

…#2) (#2990)

* ZDirection -> NegativeZDirection (doctests pass locally atm)

* @info -> @warn

* comment parts of the documentation

* uncomment tilted bbl example (docs still build locally w this change)

* bring back other doc pages (still builds)

* full docs

* adds review suggestions

* bugfix

* better unit-vector validation + fix coriolis tests

---------

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>

docs/src/model_setup/buoyancy_and_equation_of_state.md

@@ -23,25 +23,25 @@ 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
 ├── buoyancy: Nothing
 └── 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,19 +221,23 @@ 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(θ));
 
 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
 ```

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,
 

src/BuoyancyModels/buoyancy.jl

@@ -1,57 +1,65 @@
-using Oceananigans.Grids: ZDirection, validate_unit_vector
+using Oceananigans.Grids: NegativeZDirection, validate_unit_vector
 
 struct Buoyancy{M, G}
     model :: M
     gravity_unit_vector :: 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̃)
 
 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))

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)

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

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
 
 #####

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

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)