Merge remote-tracking branch 'origin/main' into ss/mask-inactive-points

.buildkite/pipeline.yml

@@ -319,13 +319,14 @@ steps:
 ##### Turbulence Closures
 #####
 
+
   - label: "🎣 gpu turbulence closures"
     env:
       JULIA_DEPOT_PATH: "$SVERDRUP_HOME/.julia-$BUILDKITE_BUILD_NUMBER"
       TEST_GROUP: "turbulence_closures"
       GPU_TEST: "true"
     commands:
-      - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project -e 'using Pkg; Pkg.test()'"
+      - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia --color=yes --project -e 'using Pkg; Pkg.test()'"
     agents:
       queue: Oceananigans
       architecture: GPU
@@ -335,20 +336,7 @@ steps:
           limit: 1
     depends_on: "init_gpu"
 
-  - label: "🎏 cpu turbulence closures"
-    env:
-      JULIA_DEPOT_PATH: "$TARTARUS_HOME/.julia-$BUILDKITE_BUILD_NUMBER"
-      TEST_GROUP: "turbulence_closures"
-    commands:
-      - "$TARTARUS_HOME/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project -e 'using Pkg; Pkg.test()'"
-    agents:
-      queue: Oceananigans
-      architecture: CPU
-    retry:
-      automatic:
-        - exit_status: 1
-          limit: 1
-    depends_on: "init_cpu"
+#  The CPU turbulence closures test used to be here, but was moved to Github Actions. See https://github.com/CliMA/Oceananigans.jl/pull/3926
 
 #####
 ##### HydrostaticFreeSurfaceModel

.github/workflows/ci.yml

@@ -0,0 +1,36 @@
+name: CI
+on:
+  - push
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.10'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          TEST_GROUP: "turbulence_closures"

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.95.0"
+version = "0.95.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

benchmark/benchmark_equations_of_state.jl

@@ -3,7 +3,7 @@ push!(LOAD_PATH, joinpath(@__DIR__, ".."))
 using BenchmarkTools
 using CUDA
 using Oceananigans
-using Oceananigans.BuoyancyModels
+using Oceananigans.BuoyancyFormulations
 using SeawaterPolynomials
 using Benchmarks
 

docs/src/appendix/library.md

@@ -40,7 +40,7 @@ Private = false
 ## Buoyancy models
 
 ```@autodocs
-Modules = [Oceananigans.BuoyancyModels]
+Modules = [Oceananigans.BuoyancyFormulations]
 Private = false
 ```
 

docs/src/model_setup/buoyancy_and_equation_of_state.md

@@ -220,8 +220,9 @@ BoussinesqEquationOfState{Float64}:
 ## The direction of gravitational acceleration
 
 To simulate gravitational accelerations that don't align with the vertical (`z`) coordinate,
-we wrap the buoyancy model in
-`Buoyancy()` function call, which takes the keyword arguments `model` and `gravity_unit_vector`,
+we use `BuoyancyForce(formulation; gravity_unit_vector)`, wherein the buoyancy `formulation`
+can be `BuoyancyTracer`, `SeawaterBuoyancy`, etc, in addition to the `gravity_unit_vector`.
+For example,
 
 ```jldoctest buoyancy
 julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1));
@@ -230,9 +231,9 @@ julia> θ = 45; # degrees
 
 julia> g̃ = (0, sind(θ), cosd(θ));
 
-julia> buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃)
-Buoyancy:
-├── model: BuoyancyTracer
+julia> buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃)
+BuoyancyForce:
+├── formulation: BuoyancyTracer
 └── gravity_unit_vector: (0.0, 0.707107, 0.707107)
 
 julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=:b)

examples/langmuir_turbulence.jl

@@ -44,7 +44,7 @@ grid = RectilinearGrid(size=(32, 32, 32), extent=(128, 128, 64))
 # (half the distance from wave crest to wave trough), which determine the wave
 # frequency and the vertical scale of the Stokes drift profile.
 
-using Oceananigans.BuoyancyModels: g_Earth
+using Oceananigans.BuoyancyFormulations: g_Earth
 
  amplitude = 0.8 # m
 wavelength = 60  # m

examples/tilted_bottom_boundary_layer.jl

@@ -74,9 +74,9 @@ current_figure() #hide
 ĝ = [sind(θ), 0, cosd(θ)]
 
 # Changing the vertical direction impacts both the `gravity_unit_vector`
-# for `Buoyancy` as well as the `rotation_axis` for Coriolis forces,
+# for `BuoyancyForce` as well as the `rotation_axis` for Coriolis forces,
 
-buoyancy = Buoyancy(model = BuoyancyTracer(), gravity_unit_vector = -ĝ)
+buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector = -ĝ)
 coriolis = ConstantCartesianCoriolis(f = 1e-4, rotation_axis = ĝ)
 
 # where above we used a constant Coriolis parameter ``f = 10^{-4} \, \rm{s}^{-1}``.

src/Advection/flat_advective_fluxes.jl

@@ -23,16 +23,22 @@ end
 
 Grids = [:XFlatGrid, :YFlatGrid, :ZFlatGrid, :XFlatGrid, :YFlatGrid, :ZFlatGrid]
 
-for side in (:left_biased, :right_biased, :symmetric)
-    for (dir, Grid) in zip([:xᶠᵃᵃ, :yᵃᶠᵃ, :zᵃᵃᶠ, :xᶜᵃᵃ, :yᵃᶜᵃ, :zᵃᵃᶜ], Grids)
-        interp_function = Symbol(side, :_interpolate_, dir)
-        @eval begin
-            @inline $interp_function(i, j, k, grid::$Grid, scheme, ψ, args...)           = @inbounds ψ[i, j, k]
-            @inline $interp_function(i, j, k, grid::$Grid, scheme, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
-
-            @inline $interp_function(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ, args...)           = @inbounds ψ[i, j, k]
-            @inline $interp_function(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
-            @inline $interp_function(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ::Function, S::AbstractSmoothnessStencil, args...) = @inbounds ψ(i, j, k, grid, args...)
-        end
+for (dir, Grid) in zip([:xᶠᵃᵃ, :yᵃᶠᵃ, :zᵃᵃᶠ, :xᶜᵃᵃ, :yᵃᶜᵃ, :zᵃᵃᶜ], Grids)
+    bias_interp = Symbol(:biased_interpolate_, dir)
+    symm_interp = Symbol(:symmetric_interpolate_, dir)
+    @eval begin
+        @inline $symm_interp(i, j, k, grid::$Grid, scheme, ψ, args...)           = @inbounds ψ[i, j, k]
+        @inline $symm_interp(i, j, k, grid::$Grid, scheme, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
+
+        @inline $symm_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ, args...)           = @inbounds ψ[i, j, k]
+        @inline $symm_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
+        @inline $symm_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, ψ::Function, S::AbstractSmoothnessStencil, args...) = @inbounds ψ(i, j, k, grid, args...)
+
+        @inline $bias_interp(i, j, k, grid::$Grid, scheme, bias, ψ, args...)           = @inbounds ψ[i, j, k]
+        @inline $bias_interp(i, j, k, grid::$Grid, scheme, bias, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
+
+        @inline $bias_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, bias, ψ, args...)           = @inbounds ψ[i, j, k]
+        @inline $bias_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, bias, ψ::Function, args...) = @inbounds ψ(i, j, k, grid, args...)
+        @inline $bias_interp(i, j, k, grid::$Grid, scheme::AbstractUpwindBiasedAdvectionScheme, bias, ψ::Function, S::AbstractSmoothnessStencil, args...) = @inbounds ψ(i, j, k, grid, args...)
     end
 end

src/Advection/weno_interpolants.jl

@@ -327,7 +327,7 @@ where
 
 The ``α`` values are normalized before returning
 """
-@inline function biased_weno_weights(ψ, scheme::WENO{N, FT}, args...) where {N, FT}
+@inline function biased_weno_weights(ψ, grid, scheme::WENO{N, FT}, args...) where {N, FT}
     β = beta_loop(scheme, ψ)
 
     τ = global_smoothness_indicator(Val(N), β)
@@ -336,11 +336,11 @@ The ``α`` values are normalized before returning
     return α ./ sum(α)
 end
 
-@inline function biased_weno_weights(ijk, scheme::WENO{N, FT}, bias, dir, ::VelocityStencil, u, v) where {N, FT}
+@inline function biased_weno_weights(ijk, grid, scheme::WENO{N, FT}, bias, dir, ::VelocityStencil, u, v) where {N, FT}
     i, j, k = ijk
 
-    uₛ = tangential_stencil_u(i, j, k, scheme, bias, dir, u)
-    vₛ = tangential_stencil_v(i, j, k, scheme, bias, dir, v)
+    uₛ = tangential_stencil_u(i, j, k, grid, scheme, bias, dir, u)
+    vₛ = tangential_stencil_v(i, j, k, grid, scheme, bias, dir, v)
     βᵤ = beta_loop(scheme, uₛ)
     βᵥ = beta_loop(scheme, vₛ)
     β  =  beta_sum(scheme, βᵤ, βᵥ)
@@ -379,10 +379,10 @@ julia> load_weno_stencil(2, :x)
     for (idx, c) in enumerate(-buffer:buffer-1)
         if func 
             stencil[idx] =  dir == :x ? 
-                            :(ψ(i + $c, j, k, args...)) :
+                            :(ψ(i + $c, j, k, grid, args...)) :
                             dir == :y ?
-                            :(ψ(i, j + $c, k, args...)) :
-                            :(ψ(i, j, k + $c, args...))
+                            :(ψ(i, j + $c, k, grid, args...)) :
+                            :(ψ(i, j, k + $c, grid, args...))
         else    
             stencil[idx] =  dir == :x ? 
                             :(ψ[i + $c, j, k]) :
@@ -400,27 +400,27 @@ end
 for dir in (:x, :y, :z), (T, f) in zip((:Any, :Function), (false, true))
     stencil = Symbol(:weno_stencil_, dir)
     @eval begin
-        @inline function $stencil(i, j, k, ::WENO{2}, bias, ψ::$T, args...) 
+        @inline function $stencil(i, j, k, grid, ::WENO{2}, bias, ψ::$T, args...) 
             S = @inbounds $(load_weno_stencil(2, dir, f))
             return S₀₂(S, bias), S₁₂(S, bias)
         end
 
-        @inline function $stencil(i, j, k, ::WENO{3}, bias, ψ::$T, args...) 
+        @inline function $stencil(i, j, k, grid, ::WENO{3}, bias, ψ::$T, args...) 
             S = @inbounds $(load_weno_stencil(3, dir, f))
             return S₀₃(S, bias), S₁₃(S, bias), S₂₃(S, bias)
         end
 
-        @inline function $stencil(i, j, k, ::WENO{4}, bias, ψ::$T, args...) 
+        @inline function $stencil(i, j, k, grid, ::WENO{4}, bias, ψ::$T, args...) 
             S = @inbounds $(load_weno_stencil(4, dir, f))
             return S₀₄(S, bias), S₁₄(S, bias), S₂₄(S, bias), S₃₄(S, bias)
         end
 
-        @inline function $stencil(i, j, k, ::WENO{5}, bias, ψ::$T, args...) 
+        @inline function $stencil(i, j, k, grid, ::WENO{5}, bias, ψ::$T, args...) 
             S = @inbounds $(load_weno_stencil(5, dir, f))
             return S₀₅(S, bias), S₁₅(S, bias), S₂₅(S, bias), S₃₅(S, bias), S₄₅(S, bias)
         end
 
-        @inline function $stencil(i, j, k, ::WENO{6}, bias, ψ::$T, args...) 
+        @inline function $stencil(i, j, k, grid, ::WENO{6}, bias, ψ::$T, args...) 
             S = @inbounds $(load_weno_stencil(6, dir, f))
             return S₀₆(S, bias), S₁₆(S, bias), S₂₆(S, bias), S₃₆(S, bias), S₄₆(S, bias), S₅₆(S, bias)
         end
@@ -455,10 +455,10 @@ end
 
 # Stencil for vector invariant calculation of smoothness indicators in the horizontal direction
 # Parallel to the interpolation direction! (same as left/right stencil)
-@inline tangential_stencil_u(i, j, k, scheme, bias, ::Val{1}, u) = @inbounds weno_stencil_x(i, j, k, scheme, bias, ℑyᵃᶠᵃ, u)
-@inline tangential_stencil_u(i, j, k, scheme, bias, ::Val{2}, u) = @inbounds weno_stencil_y(i, j, k, scheme, bias, ℑyᵃᶠᵃ, u)
-@inline tangential_stencil_v(i, j, k, scheme, bias, ::Val{1}, v) = @inbounds weno_stencil_x(i, j, k, scheme, bias, ℑxᶠᵃᵃ, v)
-@inline tangential_stencil_v(i, j, k, scheme, bias, ::Val{2}, v) = @inbounds weno_stencil_y(i, j, k, scheme, bias, ℑxᶠᵃᵃ, v)
+@inline tangential_stencil_u(i, j, k, grid, scheme, bias, ::Val{1}, u) = @inbounds weno_stencil_x(i, j, k, grid, scheme, bias, ℑyᵃᶠᵃ, u)
+@inline tangential_stencil_u(i, j, k, grid, scheme, bias, ::Val{2}, u) = @inbounds weno_stencil_y(i, j, k, grid, scheme, bias, ℑyᵃᶠᵃ, u)
+@inline tangential_stencil_v(i, j, k, grid, scheme, bias, ::Val{1}, v) = @inbounds weno_stencil_x(i, j, k, grid, scheme, bias, ℑxᶠᵃᵃ, v)
+@inline tangential_stencil_v(i, j, k, grid, scheme, bias, ::Val{2}, v) = @inbounds weno_stencil_y(i, j, k, grid, scheme, bias, ℑxᶠᵃᵃ, v)
 
 # Trick to force compilation of Val(stencil-1) and avoid loops on the GPU
 @inline function metaprogrammed_weno_reconstruction(buffer)
@@ -509,36 +509,36 @@ for (interp, dir, val, cT) in zip([:xᶠᵃᵃ, :yᵃᶠᵃ, :zᵃᵃᶠ], [:x,
                                             scheme::WENO{N, FT, XT, YT, ZT}, bias,
                                             ψ, idx, loc, args...) where {N, FT, XT, YT, ZT}
 
-            ψₜ = $stencil(i, j, k, scheme, bias, ψ, grid, args...)
-            ω = biased_weno_weights(ψₜ, scheme, bias, Val($val), Nothing, args...)
+            ψₜ = $stencil(i, j, k, grid, scheme, bias, ψ, args...)
+            ω = biased_weno_weights(ψₜ, grid, scheme, bias, Val($val), Nothing, args...)
             return weno_reconstruction(scheme, bias, ψₜ, ω, $cT, $val, idx, loc)
         end
 
         @inline function $interpolate_func(i, j, k, grid, 
                                             scheme::WENO{N, FT, XT, YT, ZT}, bias, 
                                             ψ, idx, loc, VI::AbstractSmoothnessStencil, args...) where {N, FT, XT, YT, ZT}
 
-            ψₜ = $stencil(i, j, k, scheme, bias, ψ, grid, args...)
-            ω = biased_weno_weights(ψₜ, scheme, bias, Val($val), VI, args...)
+            ψₜ = $stencil(i, j, k, grid, scheme, bias, ψ, args...)
+            ω = biased_weno_weights(ψₜ, grid, scheme, bias, Val($val), VI, args...)
             return weno_reconstruction(scheme, bias, ψₜ, ω, $cT, $val, idx, loc)
         end
 
         @inline function $interpolate_func(i, j, k, grid, 
                                             scheme::WENO{N, FT, XT, YT, ZT}, bias, 
                                             ψ, idx, loc, VI::VelocityStencil, u, v, args...) where {N, FT, XT, YT, ZT}
 
-            ψₜ = $stencil(i, j, k, scheme, bias, ψ, grid, u, v, args...)
-            ω = biased_weno_weights((i, j, k), scheme, bias, Val($val), VI, u, v)
+            ψₜ = $stencil(i, j, k, grid, scheme, bias, ψ, u, v, args...)
+            ω = biased_weno_weights((i, j, k), grid, scheme, bias, Val($val), VI, u, v)
             return weno_reconstruction(scheme, bias, ψₜ, ω, $cT, $val, idx, loc)
         end
 
         @inline function $interpolate_func(i, j, k, grid, 
                                             scheme::WENO{N, FT, XT, YT, ZT}, bias, 
                                             ψ, idx, loc, VI::FunctionStencil, args...) where {N, FT, XT, YT, ZT}
 
-            ψₜ = $stencil(i, j, k, scheme, bias, ψ,       grid, args...)
-            ψₛ = $stencil(i, j, k, scheme, bias, VI.func, grid, args...)
-            ω = biased_weno_weights(ψₛ, scheme, bias, Val($val), VI, args...)
+            ψₜ = $stencil(i, j, k, grid, scheme, bias, ψ,       args...)
+            ψₛ = $stencil(i, j, k, grid, scheme, bias, VI.func, args...)
+            ω = biased_weno_weights(ψₛ, grid, scheme, bias, Val($val), VI, args...)
             return weno_reconstruction(scheme, bias, ψₜ, ω, $cT, $val, idx, loc)
         end
     end

src/BoundaryConditions/polar_boundary_condition.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Grids: inactive_node, new_data
+using Oceananigans.Grids: inactive_node, new_data, YFlatGrid
 using CUDA: @allowscalar
 
 struct PolarValue{D, S}
@@ -20,6 +20,10 @@ end
 # Just a column
 @inline getbc(pv::BC{<:Value, <:PolarValue}, i, k, args...) = @inbounds pv.condition.data[1, 1, k]
 
+# YFlat grids do not have boundary conditions!
+latitude_north_auxiliary_bc(::YFlatGrid, args...) = nothing
+latitude_south_auxiliary_bc(::YFlatGrid, args...) = nothing
+
 # TODO: vectors should have a different treatment since vector components should account for the frame of reference.
 # For the moment, the `PolarBoundaryConditions` is implemented only for fields that have `loc[1] == loc[2] == Center()`, which
 # we assume are not components of horizontal vectors that would require rotation. (The `w` velocity if not a tracer, but it does

src/BuoyancyModels/BuoyancyModels.jl → src/BuoyancyFormulations/BuoyancyFormulations.jl

@@ -1,12 +1,9 @@
-module BuoyancyModels
+module BuoyancyFormulations
 
 export
-    Buoyancy, BuoyancyTracer, SeawaterBuoyancy, buoyancy_perturbationᶜᶜᶜ,
-    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState, TEOS10,
+    BuoyancyForce, BuoyancyTracer, SeawaterBuoyancy, LinearEquationOfState,
     ∂x_b, ∂y_b, ∂z_b, buoyancy_perturbationᶜᶜᶜ, x_dot_g_bᶠᶜᶜ, y_dot_g_bᶜᶠᶜ, z_dot_g_bᶜᶜᶠ,
-    top_buoyancy_flux,
-    buoyancy_frequency_squared,
-    BuoyancyField
+    top_buoyancy_flux, buoyancy, buoyancy_frequency, BuoyancyField
 
 using Printf
 using Oceananigans.Grids
@@ -19,11 +16,11 @@ import SeawaterPolynomials: ρ′, thermal_expansion, haline_contraction
 const g_Earth = 9.80665    # [m s⁻²] conventional standard value for Earth's gravity https://en.wikipedia.org/wiki/Gravitational_acceleration#Gravity_model_for_Earth
 
 """
-    AbstractBuoyancyModel{EOS}
+    AbstractBuoyancyFormulation{EOS}
 
 Abstract supertype for buoyancy models.
 """
-abstract type AbstractBuoyancyModel{EOS} end
+abstract type AbstractBuoyancyFormulation{EOS} end
 
 """
     AbstractEquationOfState
@@ -41,7 +38,7 @@ function validate_buoyancy(buoyancy, tracers)
     return nothing
 end
 
-include("buoyancy.jl")
+include("buoyancy_force.jl")
 include("no_buoyancy.jl")
 include("buoyancy_tracer.jl")
 include("seawater_buoyancy.jl")

src/BuoyancyFormulations/buoyancy_field.jl

@@ -0,0 +1,20 @@
+using Oceananigans.AbstractOperations: KernelFunctionOperation
+using Oceananigans.Fields: Field, ZeroField
+
+buoyancy(::Nothing, args...) = ZeroField()
+buoyancy(::BuoyancyTracer, grid, tracers) = tracers.b
+
+# TODO: move to Models
+buoyancy(model) = buoyancy(model.buoyancy, model.grid, model.tracers)
+buoyancy(bf::BuoyancyForce, grid, tracers) = buoyancy(bf.formulation, grid, tracers)
+
+buoyancy(bm::AbstractBuoyancyFormulation, grid, tracers) =
+    KernelFunctionOperation{Center, Center, Center}(buoyancy_perturbationᶜᶜᶜ, grid, bm, tracers)
+
+BuoyancyField(model) = Field(buoyancy(model))
+
+buoyancy_frequency(model) = buoyancy_frequency(model.buoyancy, model.grid, model.tracers)
+buoyancy_frequency(b::BuoyancyForce, grid, tracers) = buoyancy_frequency(b.formulation, grid, tracers)
+buoyancy_frequency(bm::AbstractBuoyancyFormulation, grid, tracers) =
+    KernelFunctionOperation{Center, Center, Face}(∂z_b, grid, bm, tracers)
+