Add unit tests for background field flux divergence

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.94.1"
+version = "0.94.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -72,13 +72,13 @@ Rotations = "1.0"
 SeawaterPolynomials = "0.3.5"
 SparseArrays = "1.9"
 Statistics = "1.9"
-StructArrays = "0.4, 0.5, 0.6"
+StructArrays = "0.4, 0.5, 0.6, 0.7"
 julia = "1.9"
 
 [extras]
+CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
-CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/oceananigans.bib

@@ -979,3 +979,13 @@ @article{BouZeid05
     year = {2005},
     doi = {10.1063/1.1839152},
 }
+
+@article{Lan2022,
+	author = {Lan, Rihui and Ju, Lili and Wanh, Zhu and Gunzburger, Max and Jones, Philip},
+	title = {High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations},
+	journal = {Journal of Computational Physics},
+	volume = {457},
+	pages = {111050},
+	year = {2022},
+	doi = {https://doi.org/10.1016/j.jcp.2022.111050},
+}

docs/src/grids.md

@@ -174,7 +174,7 @@ mountain_grid = ImmersedBoundaryGrid(grid, GridFittedBottom(mountain))
 
 # output
 20×20×20 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo:
-├── immersed_boundary: GridFittedBottom(mean(z)=6.28318, min(z)=1.58939e-8, max(z)=93.9413)
+├── immersed_boundary: GridFittedBottom(mean(z)=4.5, min(z)=0.0, max(z)=100.0)
 ├── underlying_grid: 20×20×20 RectilinearGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo
 ├── Bounded  x ∈ [-5000.0, 5000.0] regularly spaced with Δx=500.0
 ├── Bounded  y ∈ [-5000.0, 5000.0] regularly spaced with Δy=500.0

src/AbstractOperations/binary_operations.jl

@@ -1,3 +1,5 @@
+using Oceananigans.Operators
+
 const binary_operators = Set()
 
 struct BinaryOperation{LX, LY, LZ, O, A, B, IA, IB, G, T} <: AbstractOperation{LX, LY, LZ, G, T}

src/AbstractOperations/grid_metrics.jl

@@ -2,121 +2,32 @@ using Adapt
 using Oceananigans.Operators
 using Oceananigans.Grids: AbstractGrid
 using Oceananigans.Fields: AbstractField, default_indices, location
+using Oceananigans.Operators: Δx, Δy, Δz, Ax, Δλ, Δφ, Ay, Az, volume
 
 import Oceananigans.Grids: xspacings, yspacings, zspacings, λspacings, φspacings
 
-abstract type AbstractGridMetric end
-
-struct XSpacingMetric <: AbstractGridMetric end
-struct YSpacingMetric <: AbstractGridMetric end
-struct ZSpacingMetric <: AbstractGridMetric end
-
-metric_function_prefix(::XSpacingMetric) = :Δx
-metric_function_prefix(::YSpacingMetric) = :Δy
-metric_function_prefix(::ZSpacingMetric) = :Δz
-
-struct XAreaMetric <: AbstractGridMetric end
-struct YAreaMetric <: AbstractGridMetric end
-struct ZAreaMetric <: AbstractGridMetric end
-
-metric_function_prefix(::XAreaMetric) = :Ax
-metric_function_prefix(::YAreaMetric) = :Ay
-metric_function_prefix(::ZAreaMetric) = :Az
-
-struct VolumeMetric <: AbstractGridMetric end
-
-metric_function_prefix(::VolumeMetric) = :V
-
-# Convenient instances for users
-const Δx = XSpacingMetric()
-const Δy = YSpacingMetric()
-
-"""
-    Δz = ZSpacingMetric()
-
-Instance of `ZSpacingMetric` that generates `BinaryOperation`s
-between `AbstractField`s and the vertical grid spacing evaluated
-at the same location as the `AbstractField`.
-
-`Δx` and `Δy` play a similar role for horizontal grid spacings.
-
-Example
-=======
-
-```jldoctest
-julia> using Oceananigans
-
-julia> using Oceananigans.AbstractOperations: Δz
-
-julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3)));
-
-julia> c_dz = c * Δz # returns BinaryOperation between Field and GridMetricOperation
-BinaryOperation at (Center, Center, Center)
-├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo
-└── tree:
-    * at (Center, Center, Center)
-    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
-    └── Δzᶜᶜᶜ at (Center, Center, Center)
-
-julia> c .= 1;
-
-julia> c_dz[1, 1, 1]
-3.0
-```
-"""
-const Δz = ZSpacingMetric()
-
-const Ax = XAreaMetric()
-const Ay = YAreaMetric()
-const Az = ZAreaMetric()
-
-"""
-    volume = VolumeMetric()
-
-Instance of `VolumeMetric` that generates `BinaryOperation`s
-between `AbstractField`s and their cell volumes. Summing
-this `BinaryOperation` yields an integral of `AbstractField`
-over the domain.
-
-Example
-=======
-
-```jldoctest
-julia> using Oceananigans
-
-julia> using Oceananigans.AbstractOperations: volume
-
-julia> c = CenterField(RectilinearGrid(size=(2, 2, 2), extent=(1, 2, 3)));
-
-julia> c .= 1;
-
-julia> c_dV = c * volume
-BinaryOperation at (Center, Center, Center)
-├── grid: 2×2×2 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×2×2 halo
-└── tree:
-    * at (Center, Center, Center)
-    ├── 2×2×2 Field{Center, Center, Center} on RectilinearGrid on CPU
-    └── Vᶜᶜᶜ at (Center, Center, Center)
-
-julia> c_dV[1, 1, 1]
-0.75
-
-julia> sum(c_dV)
-6.0
-```
-"""
-const volume = VolumeMetric()
+const AbstractGridMetric = Union{typeof(Δx), 
+                                 typeof(Δy), 
+                                 typeof(Δz), 
+                                 typeof(Δλ),
+                                 typeof(Δφ),
+                                 typeof(Ax), 
+                                 typeof(Ay), 
+                                 typeof(Az), 
+                                 typeof(volume)} # Do we want it to be `volume` or just `V` like in the Operators module?
 
 """
     metric_function(loc, metric::AbstractGridMetric)
 
-Return the function associated with `metric::AbstractGridMetric`
-at `loc`ation.
+Return the function associated with `metric::AbstractGridMetric` at `loc`ation.
 """
-function metric_function(loc, metric::AbstractGridMetric)
+function metric_function(loc, metric)
     code = Tuple(interpolation_code(ℓ) for ℓ in loc)
-    prefix = metric_function_prefix(metric)
-    metric_function_symbol = Symbol(prefix, code...)
+    if metric isa typeof(volume)
+        metric_function_symbol = Symbol(:V, code...)
+    else
+        metric_function_symbol = Symbol(metric, code...)
+    end
     return getglobal(@__MODULE__, metric_function_symbol)
 end
 
@@ -137,12 +48,33 @@ on_architecture(to, gm::GridMetricOperation{LX, LY, LZ}) where {LX, LY, LZ} =
     GridMetricOperation{LX, LY, LZ}(on_architecture(to, gm.metric),
                                     on_architecture(to, gm.grid))
 
-
 @inline Base.getindex(gm::GridMetricOperation, i, j, k) = gm.metric(i, j, k, gm.grid)
 
 indices(::GridMetricOperation) = default_indices(3)
 
-# Special constructor for BinaryOperation
+"""
+    GridMetricOperation(L, metric, grid)
+
+Instance of `GridMetricOperation` that generates `BinaryOperation`s between `AbstractField`s and the metric `metric`
+at the same location as the `AbstractField`.
+
+Example
+=======
+```jldoctest
+julia> using Oceananigans
+
+julia> using Oceananigans.Operators: Δz
+
+julia> c = CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3)));
+
+julia> c_dz = c * Δz; # returns BinaryOperation between Field and GridMetricOperation
+
+julia> c .= 1;
+
+julia> c_dz[1, 1, 1]
+3.0
+```
+"""
 GridMetricOperation(L, metric, grid) = GridMetricOperation{L[1], L[2], L[3]}(metric_function(L, metric), grid)
 
 #####
@@ -165,13 +97,13 @@ julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
 julia> xspacings(grid, Center(), Center(), Center())
 KernelFunctionOperation at (Center, Center, Center)
 ├── grid: 2×4×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×3×3 halo
-├── kernel_function: xspacing (generic function with 27 methods)
+├── kernel_function: Δx (generic function with 29 methods)
 └── arguments: ("Center", "Center", "Center")
 ```
 """
 function xspacings(grid, ℓx, ℓy, ℓz)
     LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
-    Δx_op = KernelFunctionOperation{LX, LY, LZ}(xspacing, grid, ℓx, ℓy, ℓz)
+    Δx_op = KernelFunctionOperation{LX, LY, LZ}(Δx, grid, ℓx, ℓy, ℓz)
     return Δx_op
 end
 
@@ -191,13 +123,13 @@ julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
 julia> yspacings(grid, Center(), Face(), Center())
 KernelFunctionOperation at (Center, Face, Center)
 ├── grid: 2×4×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×3×3 halo
-├── kernel_function: yspacing (generic function with 27 methods)
+├── kernel_function: Δy (generic function with 29 methods)
 └── arguments: ("Center", "Face", "Center")
 ```
 """
 function yspacings(grid, ℓx, ℓy, ℓz)
     LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
-    Δy_op = KernelFunctionOperation{LX, LY, LZ}(yspacing, grid, ℓx, ℓy, ℓz)
+    Δy_op = KernelFunctionOperation{LX, LY, LZ}(Δy, grid, ℓx, ℓy, ℓz)
     return Δy_op
 end
 
@@ -217,21 +149,21 @@ julia> grid = RectilinearGrid(size=(2, 4, 8), extent=(1, 1, 1));
 julia> zspacings(grid, Center(), Center(), Face())
 KernelFunctionOperation at (Center, Center, Face)
 ├── grid: 2×4×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 2×3×3 halo
-├── kernel_function: zspacing (generic function with 27 methods)
+├── kernel_function: Δz (generic function with 28 methods)
 └── arguments: ("Center", "Center", "Face")
 ```
 """
 function zspacings(grid, ℓx, ℓy, ℓz)
     LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
-    Δz_op = KernelFunctionOperation{LX, LY, LZ}(zspacing, grid, ℓx, ℓy, ℓz)
+    Δz_op = KernelFunctionOperation{LX, LY, LZ}(Δz, grid, ℓx, ℓy, ℓz)
     return Δz_op
 end
 
 """
     λspacings(grid, ℓx, ℓy, ℓz)
 
 Return a `KernelFunctionOperation` that computes the grid spacings for `grid`
-in the ``z`` direction at location `ℓx, ℓy, ℓz`.
+in the ``λ`` direction at location `ℓx, ℓy, ℓz`.
 
 Examples
 ========
@@ -246,21 +178,21 @@ julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25),
 julia> λspacings(grid, Center(), Face(), Center())
 KernelFunctionOperation at (Center, Face, Center)
 ├── grid: 36×34×25 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
-├── kernel_function: λspacing (generic function with 5 methods)
+├── kernel_function: Δλ (generic function with 29 methods)
 └── arguments: ("Center", "Face", "Center")
 ```
 """
 function λspacings(grid, ℓx, ℓy, ℓz)
     LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
-    Δλ_op = KernelFunctionOperation{LX, LY, LZ}(λspacing, grid, ℓx, ℓy, ℓz)
+    Δλ_op = KernelFunctionOperation{LX, LY, LZ}(Δλ, grid, ℓx, ℓy, ℓz)
     return Δλ_op
 end
 
 """
     φspacings(grid, ℓx, ℓy, ℓz)
 
 Return a `KernelFunctionOperation` that computes the grid spacings for `grid`
-in the ``z`` direction at location `ℓx, ℓy, ℓz`.
+in the ``φ`` direction at location `ℓx, ℓy, ℓz`.
 
 Examples
 ========
@@ -275,13 +207,13 @@ julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25),
 julia> φspacings(grid, Center(), Face(), Center())
 KernelFunctionOperation at (Center, Face, Center)
 ├── grid: 36×34×25 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
-├── kernel_function: φspacing (generic function with 5 methods)
+├── kernel_function: Δφ (generic function with 29 methods)
 └── arguments: ("Center", "Face", "Center")
 ```
 """
 function φspacings(grid, ℓx, ℓy, ℓz)
     LX, LY, LZ = map(typeof, (ℓx, ℓy, ℓz))
-    Δφ_op = KernelFunctionOperation{LX, LY, LZ}(φspacing, grid, ℓx, ℓy, ℓz)
+    Δφ_op = KernelFunctionOperation{LX, LY, LZ}(Δφ, grid, ℓx, ℓy, ℓz)
     return Δφ_op
 end
 

src/Advection/positivity_preserving_tracer_advection_operators.jl

@@ -7,6 +7,8 @@ const ε₂ = 1e-20
 # Here in the future we can easily add UpwindBiased
 const BoundPreservingScheme = PositiveWENO
 
+@inline div_Uc(i, j, k, grid, advection::BoundPreservingScheme, U, ::ZeroField) = zero(grid)
+
 # Is this immersed-boundary safe without having to extend it in ImmersedBoundaries.jl? I think so... (velocity on immmersed boundaries is masked to 0)
 @inline function div_Uc(i, j, k, grid, advection::BoundPreservingScheme, U, c)
 

src/Advection/weno_reconstruction.jl

@@ -104,6 +104,10 @@ function WENO(FT::DataType=Float64;
 
     mod(order, 2) == 0 && throw(ArgumentError("WENO reconstruction scheme is defined only for odd orders"))
 
+    if !isnothing(bounds)
+        @warn "Bounds preserving WENO is experimental."
+    end
+
     if order < 3
         # WENO(order=1) is equivalent to UpwindBiased(order=1)
         return UpwindBiased(FT; order=1)

src/Biogeochemistry.jl

@@ -27,7 +27,7 @@ Update biogeochemical state variables. Called at the end of update_state!.
 """
 update_biogeochemical_state!(bgc, model) = nothing
 
-@inline biogeochemical_drift_velocity(bgc, val_tracer_name) = (u = ZeroField(), v = ZeroField(), w = ZeroField())
+@inline biogeochemical_drift_velocity(bgc, val_tracer_name) = nothing
 @inline biogeochemical_auxiliary_fields(bgc) = NamedTuple()
 
 """

src/BoundaryConditions/BoundaryConditions.jl

@@ -37,6 +37,7 @@ include("fill_halo_regions_nothing.jl")
 include("apply_flux_bcs.jl")
 
 include("update_boundary_conditions.jl")
+include("polar_boundary_condition.jl")
 
 include("flat_extrapolation_open_boundary_matching_scheme.jl")
 end # module