Merge branch 'main' into glw/divergent-lagrangian-mean

.buildkite/pipeline.yml

@@ -1,7 +1,7 @@
 env:
-  JULIA_VERSION: "1.10.6"
+  JULIA_VERSION: "1.10.8"
   JULIA_MINOR_VERSION: "1.10"
-  SVERDRUP_HOME: "/data5/glwagner"
+  SVERDRUP_HOME: "/home/users/glwagner"
   TARTARUS_HOME: "/storage5/buildkite-agent"
   JULIA_PKG_SERVER_REGISTRY_PREFERENCE: eager
   JULIA_NUM_PRECOMPILE_TASKS: 8
@@ -649,6 +649,42 @@ steps:
           limit: 1
     depends_on: "init_cpu"
 
+#####
+##### Vertical Coordinates tests
+#####
+
+  - label: "🥑 gpu vertical coordinate"
+    env:
+      JULIA_DEPOT_PATH: "$SVERDRUP_HOME/.julia-$BUILDKITE_BUILD_NUMBER"
+      TEST_GROUP: "vertical_coordinate"
+      GPU_TEST: "true"
+    commands:
+      - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project -e 'using Pkg; Pkg.test()'"
+    agents:
+      queue: Oceananigans
+      architecture: GPU
+    retry:
+      automatic:
+        - exit_status: 1      
+          limit: 1
+    depends_on: "init_gpu"
+
+  - label: "🥒 cpu vertical coordinate"
+    env:
+      JULIA_DEPOT_PATH: "$TARTARUS_HOME/.julia-$BUILDKITE_BUILD_NUMBER"
+      TEST_GROUP: "vertical_coordinate"
+      CUDA_VISIBLE_DEVICES: "-1"
+    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"
+
 #####
 ##### Enzyme extension tests
 #####

.github/workflows/ci.yml

@@ -15,12 +15,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.95.5"
+version = "0.95.8"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -12,6 +12,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -44,7 +45,7 @@ OceananigansMakieExt = ["MakieCore", "Makie"]
 
 [compat]
 Adapt = "4.1.1"
-CUDA = "4.1.1, 5"
+CUDA = "5"
 Crayons = "4"
 CubedSphere = "0.2, 0.3"
 Dates = "1.9"
@@ -53,6 +54,7 @@ DocStringExtensions = "0.8, 0.9"
 Enzyme = "0.13.14"
 FFTW = "1"
 Glob = "1.3"
+GPUArrays = "10.3"
 IncompleteLU = "0.2"
 InteractiveUtils = "1.9"
 IterativeSolvers = "0.9"
@@ -61,8 +63,8 @@ KernelAbstractions = "0.9.21"
 LinearAlgebra = "1.9"
 Logging = "1.9"
 MPI = "0.16, 0.17, 0.18, 0.19, 0.20"
-Makie = "0.21"
-MakieCore = "0.7, 0.8"
+Makie = "0.21, 0.22"
+MakieCore = "0.7, 0.8, 0.9"
 NCDatasets = "0.12.10, 0.13.1, 0.14"
 OffsetArrays = "1.4"
 OrderedCollections = "1.1"

docs/src/index.md

@@ -24,8 +24,11 @@ julia> using Pkg
 julia> Pkg.add("Oceananigans")
 ```
 
-!!! compat "Julia 1.9 is required"
+!!! compat "Julia 1.9 or later is required"
     Oceananigans requires Julia 1.9 or later.
+
+!!! info "Tested Julia versions"
+    Oceananigans is currently tested on Julia 1.10.
 
 If you're [new to Julia](https://docs.julialang.org/en/v1/manual/getting-started/) and its [wonderful `Pkg` manager](https://docs.julialang.org/en/v1/stdlib/Pkg/), the [Oceananigans wiki](https://github.com/CliMA/Oceananigans.jl/wiki) provides [more detailed installation instructions](https://github.com/CliMA/Oceananigans.jl/wiki/Installation-and-getting-started-with-Oceananigans).
 
@@ -79,7 +82,9 @@ cite our work and mention Oceananigans by name.
 
 If you have work using Oceananigans that you would like to have listed here, please open a pull request to add it or let us know!
 
-1. Bisits, J. I., Zika, J. D., and Evans, D. G. (2024) [Does cabbeling shape the thermohaline structure of high-latitude oceans?](https://doi.org/10.1175/JPO-D-24-0061.1) _Journal of Physical Oceanography_, in press. DOI: [10.1175/JPO-D-24-0061.1](https://doi.org/10.1175/JPO-D-24-0061.1)
+1. Abbott, K. and Mahadevan, A. (2024). [Why is the monsoon coastal upwelling signal subdued in the Bay of Bengal?](https://doi.org/10.1029/2024JC022023), _Journal of Geophysical Research: Oceans_, **129**, e2024JC022023. DOI: [10.1029/2024JC022023](https://doi.org/10.1029/2024JC022023)
+
+1. Bisits, J. I., Zika, J. D., and Evans, D. G. (2024) [Does cabbeling shape the thermohaline structure of high-latitude oceans?](https://doi.org/10.1175/JPO-D-24-0061.1), _Journal of Physical Oceanography_, in press. DOI: [10.1175/JPO-D-24-0061.1](https://doi.org/10.1175/JPO-D-24-0061.1)
 
 1. Strong-Wright J. and Taylor, J. R. (2024) [A model of tidal flow and tracer release in a giant kelp forest](https://doi.org/10.1017/flo.2024.13), _Flow_, **4**, E21. DOI: [10.1017/flo.2024.13](https://doi.org/10.1017/flo.2024.13)
 
@@ -91,10 +96,10 @@ If you have work using Oceananigans that you would like to have listed here, ple
 
 1. Silvestri, S., Wagner, G. L., Constantinou, N. C., Hill, C., Campin, J.-M., Souza, A., Bishnu, S., Churavy, V., Marshall, J., and Ferrari, R. (2024) [A GPU-based ocean dynamical core for routine mesoscale-resolving climate simulations](https://doi.org/10.22541/essoar.171708158.82342448/v1), _ESS Open Archive_. DOI: [10.22541/essoar.171708158.82342448/v1](https://doi.org/10.22541/essoar.171708158.82342448/v1)
 
-1. Silvestri, S., Wagner, G. L., Campin, J.-M., Constantinou, N. C., Hill, C., Souza, A., and Ferrari, R. (2024). [A new WENO-based momentum advection scheme for simulations of ocean mesoscale turbulence](https://doi.org/10.22541/essoar.170110657.76489860/v2), _ESS Open Archive_. DOI: [10.22541/essoar.170110657.76489860/v2](https://doi.org/10.22541/essoar.170110657.76489860/v2)
-
 1. Whitley V. and Wenegrat, J. O. (2024) [Breaking internal waves on sloping topography: connecting parcel displacements to overturn size, interior-boundary exchanges, and mixing](https://doi.org/10.31223/X5PM5Q), _Earth Arxiv_. DOI: [10.31223/X5PM5Q](https://doi.org/10.31223/X5PM5Q)
 
+1. Silvestri, S., Wagner, G. L., Campin, J.-M., Constantinou, N. C., Hill, C., Souza, A., and Ferrari, R. (2024). [A new WENO-based momentum advection scheme for simulations of ocean mesoscale turbulence](https://doi.org/10.1029/2023MS004130), _Journal of Advances in Modeling Earth Systems_, **16(7)**, e2023MS004130. DOI: [10.1029/2023MS004130](https://doi.org/10.1029/2023MS004130)
+
 1. Chen S., Strong-Wright J., and Taylor, J. R. (2024) [Modeling carbon dioxide removal via sinking of particulate organic carbon from macroalgae cultivation](https://doi.org/10.3389/fmars.2024.1359614), _Frontiers in Marine Science_, **11**, 1359614. DOI: [10.3389/fmars.2024.1359614](https://doi.org/10.3389/fmars.2024.1359614)
 
 1. Gupta, M., Gürcan, E., and Thompson, A. F. (2024). [Eddy-induced dispersion of sea ice floes at the marginal ice zone](https://doi.org/10.1029/2023GL105656), _Geophysical Research Letters_, **51**, e2023GL105656. DOI: [10.1029/2023GL105656](https://doi.org/10.1029/2023GL105656)

ext/OceananigansEnzymeExt.jl

@@ -17,6 +17,7 @@ EnzymeCore.EnzymeRules.inactive_noinl(::typeof(Base.:(==)), ::Oceananigans.Abstr
 EnzymeCore.EnzymeRules.inactive_noinl(::typeof(Oceananigans.AbstractOperations.validate_grid), x...) = nothing
 EnzymeCore.EnzymeRules.inactive_noinl(::typeof(Oceananigans.AbstractOperations.metric_function), x...) = nothing
 EnzymeCore.EnzymeRules.inactive_noinl(::typeof(Oceananigans.Utils.flatten_reduced_dimensions), x...) = nothing
+EnzymeCore.EnzymeRules.inactive_noinl(::typeof(Oceananigans.Utils.prettytime), x...) = nothing
 EnzymeCore.EnzymeRules.inactive(::typeof(Oceananigans.Grids.total_size), x...) = nothing
 EnzymeCore.EnzymeRules.inactive(::typeof(Oceananigans.BoundaryConditions.parent_size_and_offset), x...) = nothing
 @inline EnzymeCore.EnzymeRules.inactive_type(v::Type{Oceananigans.Utils.KernelParameters}) = true

ext/OceananigansMakieExt.jl

@@ -1,12 +1,13 @@
 module OceananigansMakieExt
 
 using Oceananigans
-using Oceananigans.Grids: OrthogonalSphericalShellGrid
+using Oceananigans.Grids: OrthogonalSphericalShellGrid, topology
 using Oceananigans.Fields: AbstractField
 using Oceananigans.AbstractOperations: AbstractOperation
 using Oceananigans.Architectures: on_architecture
 using Oceananigans.ImmersedBoundaries: mask_immersed_field!
 
+using Makie: Observable
 using MakieCore: AbstractPlot
 import MakieCore: convert_arguments, _create_plot
 import Makie: args_preferred_axis
@@ -46,11 +47,14 @@ end
 
 axis_str(::RectilinearGrid, dim) = ("x", "y", "z")[dim]
 axis_str(::LatitudeLongitudeGrid, dim) = ("Longitude (deg)", "Latitude (deg)", "z")[dim]
+axis_str(grid::ImmersedBoundaryGrid, dim) = axis_str(grid.underlying_grid, dim)
+
+const LLGOrIBLLG = Union{LatitudeLongitudeGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:LatitudeLongitudeGrid}}
 
 function _create_plot(F::Function, attributes::Dict, f::Field)
     converted_args = convert_field_argument(f)
 
-    if !(:axis ∈ keys(attributes)) # Let's try to add labels automatically
+    if !(:axis ∈ keys(attributes)) # Let's try to automatically add labels and ticks
         d1, d2, D = deduce_dimensionality(f) 
         grid = f.grid
 
@@ -69,6 +73,11 @@ function _create_plot(F::Function, attributes::Dict, f::Field)
             throw(ArgumentError("Cannot create axis labels for a 3D field!"))
         end
 
+        # if longitude wraps around the globe then adjust the longitude ticks
+        if grid isa LLGOrIBLLG && grid.Lx == 360 && topology(grid, 1) == Periodic
+            axis = merge(axis, (xticks = -360:60:360,))
+        end
+
         attributes[:axis] = axis
     end
 
@@ -81,6 +90,9 @@ function _create_plot(F::Function, attributes::Dict, op::AbstractOperation)
     return _create_plot(F::Function, attributes::Dict, f)
 end
 
+_create_plot(F::Function, attributes::Dict, f::Observable{<:Field}) =
+    _create_plot(F, attributes, f[])
+
 convert_arguments(pl::Type{<:AbstractPlot}, f::Field) =
     convert_arguments(pl, convert_field_argument(f)...)
 
@@ -183,4 +195,3 @@ function convert_arguments(pl::Type{<:AbstractPlot}, ξ1::AbstractArray, ξ2::Ab
 end
 
 end # module
-

src/AbstractOperations/at.jl

@@ -1,5 +1,5 @@
 using Oceananigans
-using Oceananigans.Fields: default_indices
+using Oceananigans.Fields: default_indices, compute_index_intersection
 
 """
     insert_location(ex::Expr, location)
@@ -55,8 +55,8 @@ macro at(location, abstract_operation)
 end
 
 # Numbers and functions do not have index restrictions
-indices(f::Function) = default_indices(3)
-indices(f::Number)   = default_indices(3)
+indices(::Function) = default_indices(3)
+indices(::Number)   = default_indices(3)
 
 """
     intersect_indices(loc, operands...)
@@ -65,66 +65,22 @@ Utility to compute the intersection of `operands' indices.
 """
 function intersect_indices(loc, operands...)
 
-    idx1 = compute_index_intersection(Colon(), loc[1], operands...; dim=1)
-    idx2 = compute_index_intersection(Colon(), loc[2], operands...; dim=2)
-    idx3 = compute_index_intersection(Colon(), loc[3], operands...; dim=3)
+    idx1 = compute_operand_intersection(Colon(), loc[1], operands...; dim=1)
+    idx2 = compute_operand_intersection(Colon(), loc[2], operands...; dim=2)
+    idx3 = compute_operand_intersection(Colon(), loc[3], operands...; dim=3)
 
     return (idx1, idx2, idx3)
 end
 
 # Fallback for `KernelFunctionOperation`s with no argument 
-compute_index_intersection(::Colon, to_loc; kw...) = Colon()
+compute_operand_intersection(::Colon, to_loc; kw...) = Colon()
 
-compute_index_intersection(to_idx, to_loc, op; dim) =
-    _compute_index_intersection(to_idx, indices(op)[dim],
-                                to_loc, location(op, dim))
+compute_operand_intersection(to_idx, to_loc, op; dim) =
+    compute_index_intersection(to_idx, indices(op)[dim],
+                               to_loc, location(op, dim))
 
 """Compute index intersection recursively for `dim`ension ∈ (1, 2, 3)."""
-function compute_index_intersection(to_idx, to_loc, op1, op2, more_ops...; dim)
-    new_to_idx = _compute_index_intersection(to_idx, indices(op1)[dim], to_loc, location(op1, dim))
-    return compute_index_intersection(new_to_idx, to_loc, op2, more_ops...; dim)
-end
-
-# Life is pretty simple in this case.
-_compute_index_intersection(to_idx::Colon, from_idx::Colon, args...) = Colon()
-
-# Because `from_idx` imposes no restrictions, we just return `to_idx`.
-_compute_index_intersection(to_idx::UnitRange, from_idx::Colon, args...) = to_idx
-
-# This time we account for the possible range-reducing effect of interpolation on `from_idx`.
-function _compute_index_intersection(to_idx::Colon, from_idx::UnitRange, to_loc, from_loc)
-    shifted_idx = restrict_index_for_interpolation(from_idx, from_loc, to_loc)
-    validate_shifted_index(shifted_idx)
-    return shifted_idx
-end
-
-# Compute the intersection of two index ranges
-function _compute_index_intersection(to_idx::UnitRange, from_idx::UnitRange, to_loc, from_loc)
-    shifted_idx = restrict_index_for_interpolation(from_idx, from_loc, to_loc)
-    validate_shifted_index(shifted_idx)
-
-    range_intersection = UnitRange(max(first(shifted_idx), first(to_idx)), min(last(shifted_idx), last(to_idx)))
-
-    # Check validity of the intersection index range
-    first(range_intersection) > last(range_intersection) &&
-        throw(ArgumentError("Indices $(from_idx) and $(to_idx) interpolated from $(from_loc) to $(to_loc) do not intersect!"))
-
-    return range_intersection
-end
-
-validate_shifted_index(shifted_idx) = first(shifted_idx) > last(shifted_idx) &&
-    throw(ArgumentError("Cannot compute index intersection for indices $(from_idx) interpolating from $(from_loc) to $(to_loc)!"))
-
-"""
-    restrict_index_for_interpolation(from_idx, from_loc, to_loc)
-
-Return a "restricted" index range for the result of interpolating from
-`from_loc` to `to_loc`, over the index range `from_idx`:
-
-* Windowed fields interpolated from `Center`s to `Face`s lose the first index.
-* Conversely, windowed fields interpolated from `Face`s to `Center`s lose the last index
-"""
-restrict_index_for_interpolation(from_idx, ::Type{Face},   ::Type{Face})   = UnitRange(first(from_idx),   last(from_idx))
-restrict_index_for_interpolation(from_idx, ::Type{Center}, ::Type{Center}) = UnitRange(first(from_idx),   last(from_idx))
-restrict_index_for_interpolation(from_idx, ::Type{Face},   ::Type{Center}) = UnitRange(first(from_idx),   last(from_idx)-1)
-restrict_index_for_interpolation(from_idx, ::Type{Center}, ::Type{Face})   = UnitRange(first(from_idx)+1, last(from_idx))
+function compute_operand_intersection(to_idx, to_loc, op1, op2, more_ops...; dim)
+    new_to_idx = compute_index_intersection(to_idx, indices(op1)[dim], to_loc, location(op1, dim))
+    return compute_operand_intersection(new_to_idx, to_loc, op2, more_ops...; dim)
+end

src/Advection/Advection.jl

@@ -28,12 +28,13 @@ using DocStringExtensions
 using Base: @propagate_inbounds
 using Adapt 
 using OffsetArrays
-
 using Oceananigans.Grids
 using Oceananigans.Grids: with_halo
 using Oceananigans.Architectures: architecture, CPU
+using Oceananigans: supported_float_types
 
 using Oceananigans.Operators
+using Oceananigans.Operators: flux_div_xyᶜᶜᶜ, Γᶠᶠᶜ, ∂t_σ
 
 import Base: show, summary
 import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z

src/Advection/centered_reconstruction.jl

@@ -98,16 +98,16 @@ const ACAS = AbstractCenteredAdvectionScheme
 @inline biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::ACAS, bias, c, args...) = symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, c, args...)
 
 # uniform centered reconstruction
-for buffer in advection_buffers
+for buffer in advection_buffers, FT in supported_float_types
     @eval begin
-        @inline inner_symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::Centered{$buffer, FT, <:Nothing}, ψ, idx, loc, args...)           where FT = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :x, false))
-        @inline inner_symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::Centered{$buffer, FT, <:Nothing}, ψ::Function, idx, loc, args...) where FT = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :x,  true))
+        @inline inner_symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::Centered{$buffer, $FT, <:Nothing}, ψ, idx, loc, args...)           = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :x, false))
+        @inline inner_symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::Centered{$buffer, $FT, <:Nothing}, ψ::Function, idx, loc, args...) = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :x,  true))
 
-        @inline inner_symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::Centered{$buffer, FT, XT, <:Nothing}, ψ, idx, loc, args...)           where {FT, XT} = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :y, false))
-        @inline inner_symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::Centered{$buffer, FT, XT, <:Nothing}, ψ::Function, idx, loc, args...) where {FT, XT} = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :y,  true))
+        @inline inner_symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::Centered{$buffer, $FT, XT, <:Nothing}, ψ, idx, loc, args...)           where XT = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :y, false))
+        @inline inner_symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::Centered{$buffer, $FT, XT, <:Nothing}, ψ::Function, idx, loc, args...) where XT = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :y,  true))
 
-        @inline inner_symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::Centered{$buffer, FT, XT, YT, <:Nothing}, ψ, idx, loc, args...)           where {FT, XT, YT} = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :z, false))
-        @inline inner_symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::Centered{$buffer, FT, XT, YT, <:Nothing}, ψ::Function, idx, loc, args...) where {FT, XT, YT} = @inbounds $(calc_reconstruction_stencil(buffer, :symmetric, :z,  true))
+        @inline inner_symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::Centered{$buffer, $FT, XT, YT, <:Nothing}, ψ, idx, loc, args...)           where {XT, YT} = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :z, false))
+        @inline inner_symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::Centered{$buffer, $FT, XT, YT, <:Nothing}, ψ::Function, idx, loc, args...) where {XT, YT} = @inbounds $(calc_reconstruction_stencil(FT, buffer, :symmetric, :z,  true))
     end
 end