Just updating this branch with recent changes to main. (#4707)

* (0.97.6) Move the computation of flux BC before the time-stepping (#4696)

* go ahead

* ad project

* fix tests

* start with compute

* more changes

* more changes

* more fixes

* Update cubed_sphere_exchange_bcs.jl

* Update compute_hydrostatic_flux_bcs.jl

* Update hydrostatic_free_surface_tendency_kernel_functions.jl

* Update explicit_free_surface.jl

---------

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

* Create methods for equality of ImmersedBoundaryGrids (#4700)

* Move storage space to from `MutableVerticalDiscretization` to `ZStarCoordinate` (#4686)

* first change

* comment

* bugfix

* Update hydrostatic_free_surface_model.jl

* bugfix

* bugfix

* bugfix

* fix tests

* export ZCoordinate, ZStarCoordinate

* add summary/show method for ZStarCoordinate

* bit of rephrasing

* add vertical_coordinate in show method

* Update src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

* Update validation/z_star_coordinate/lock_release.jl

* Update validation/implicit_free_surface/geostrophic_adjustment_test.jl

* fix doctests

* bump patch release

* slight rephrase

* Update lock_release.jl

---------

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

* Regularize field boundary conditions with a location rather than a name (#4704)

* Update field_boundary_conditions.jl

* add a field_name

---------

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>
Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>
Co-authored-by: Tomás Chor <tomaschor@gmail.com>

Author: 4 people

Project.toml

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

docs/src/model_setup/buoyancy_and_equation_of_state.md

@@ -65,6 +65,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration =
 │   └── solver: FFTImplicitFreeSurfaceSolver
 ├── advection scheme:
 │   └── momentum: VectorInvariant
+├── vertical_coordinate: ZCoordinate
 └── coriolis: Nothing
 ```
 
@@ -100,6 +101,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration =
 ├── advection scheme:
 │   ├── momentum: VectorInvariant
 │   └── b: Centered(order=2)
+├── vertical_coordinate: ZCoordinate
 └── coriolis: Nothing
 ```
 
@@ -142,6 +144,7 @@ HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration =
 │   ├── momentum: VectorInvariant
 │   ├── T: Centered(order=2)
 │   └── S: Centered(order=2)
+├── vertical_coordinate: ZCoordinate
 └── coriolis: Nothing
 ```
 
@@ -246,4 +249,3 @@ NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
 ├── buoyancy: BuoyancyTracer with ĝ = (0.0, 0.707107, 0.707107)
 └── coriolis: Nothing
 ```
-

docs/src/physics/hydrostatic_free_surface_model.md

@@ -93,11 +93,11 @@ equation are
 
 ## Vertical coordinates
 
-We can use either `ZCoordinate`, that is height-coordinate, or the
-`ZStar` [generalized vertical coordinate](@ref generalized_vertical_coordinates).
+We can use either `ZCoordinate`, that is height coordinate, or the
+`ZStarCoordinate` [generalized vertical coordinate](@ref generalized_vertical_coordinates).
 
-The `ZStar` vertical coordinate conserves tracers and volume with the grid following the evolution of the
-free surface in the domain [adcroft2004rescaled](@citep).
+The `ZStarCoordinate` vertical coordinate conserves tracers and volume with the grid following
+ the evolution of the free surface in the domain [adcroft2004rescaled](@citep).
 
 In terms of the notation in the [Generalized vertical coordinates](@ref generalized_vertical_coordinates)
 section, for a `ZCoordinate` we have that
@@ -106,39 +106,40 @@ r(x, y, z, t) = z
 ```
 and the specific thickness is ``\sigma = \partial z / \partial r = 1``.
 
-For the `ZStar` generalized vertical coordinate is often denoted as ``z^*`` (zee-star), i.e.,
+The `ZStarCoordinate` generalized vertical coordinate is often denoted as ``z^*`` (zee-star), i.e.,
 ```math
 \begin{equation}
     r(x, y, z, t) = z^*(x, y, z, t) = \frac{H(x, y)}{H(x, y) + \eta(x, y, t)}[z - \eta(x, y, t)] \label{zstardef}
 \end{equation}
 ```
 where ``\eta`` is the free surface and ``z = -H(x, y)`` is the bottom of the domain.
 
-![Schematic of the quantities involved in the ZStar generalized vertical coordinate](../assets/zstar_schematic.png)
+![Schematic of the quantities involved in the ZStarCoordinate generalized vertical coordinate](../assets/zstar_schematic.png)
 
-Note, that in both depth and ``z^*`` coordinates the bottom boundary is the same ``z = z^* = - H(x, y)``.
+Note that in both depth and ``z^*`` coordinates, the bottom boundary is the same ``z = z^* = - H(x, y)``.
 On the other hand, while in depth coordinates the upper boundary ``z = \eta(x, y, t)`` changes with time,
-but in ``z^*`` coordinates it's fixed to ``z^* = 0``.
+but in ``z^*`` coordinates is fixed to ``z^* = 0``.
 
-The `ZStar` coordinate definition \eqref{zstardef} implies a specific thickness
+The `ZStarCoordinate` definition \eqref{zstardef} implies a specific thickness
 
 ```math
 \sigma = 1 + \frac{\eta}{H}
 ```
 
-All the equations transformed in ``r``-coordinates are described in the [Generalized vertical coordinates](@ref generalized_vertical_coordinates)
-section.
+All the equations transformed in ``r``-coordinates are described in the
+[Generalized vertical coordinates](@ref generalized_vertical_coordinates) section.
 
-For the specific choice of `ZStar` coordinate \eqref{zstardef}, the ``\partial \eta/\partial r`` identically vanishes and
-thus the horizontal gradient of the free surface remain unchanged under vertical coordinate transformation, i.e.,
+For the specific choice of `ZStarCoordinate` coordinate \eqref{zstardef}, the ``\partial \eta/\partial r``
+identically vanishes and thus the horizontal gradient of the free surface remain unchanged under vertical
+coordinate transformation, i.e.,
 ```math
 \begin{align}
     \frac{\partial \eta}{\partial x} \bigg\rvert_z & = \frac{\partial \eta}{\partial x} \bigg\rvert_r \\
     \frac{\partial \eta}{\partial y} \bigg\rvert_z & = \frac{\partial \eta}{\partial y} \bigg\rvert_r
 \end{align}
 ```
 
-An example of how the vertical coordinate surfaces differ for `ZCoordinate` and the time-varying `ZStar` coordinate is shown below.
+An example of how the vertical coordinate surfaces differ for `ZCoordinate` and the time-varying `ZStarCoordinate` coordinate is shown below.
 
 ```@example
 using CairoMakie
@@ -165,7 +166,7 @@ t = Observable(0.0)
 fig = Figure(size=(1000, 400))
 axis_kwargs = (titlesize = 20, xlabel = "x", ygridvisible = false)
 ax1 = Axis(fig[1, 1]; title="ZCoordinate", ylabel="z", axis_kwargs...)
-ax2 = Axis(fig[1, 2]; title="ZStar", axis_kwargs...)
+ax2 = Axis(fig[1, 2]; title="ZStarCoordinate", axis_kwargs...)
 
 for ax in (ax1, ax2)
     band!(ax, x, -H, η, color = (:dodgerblue, 0.5))
@@ -179,7 +180,7 @@ for r in range(-Lz, stop=0, length=6)
     z = r * ones(size(x))
     lines!(ax1, x, z, color=:crimson, linestyle=:dash)
 
-    # ZStar
+    # ZStarCoordinate
     z = lift(η) do η_val
         @. r * (H + η_val) / H + η_val
     end
@@ -197,6 +198,6 @@ nothing #hide
 
 ![](z-zstar.gif)
 
-Near the top, the surfaces of `ZStar` mimic that of the free surface.
-As we move away from the fluid's surface, the surfaces of `ZStar` resemble more surfaces
+Near the top, the surfaces of `ZStarCoordinate` mimic the free surface.
+Further away from the fluid's surface, the surfaces of `ZStarCoordinate` resemble more surfaces
 of constant depth `ZCoordinate`.

src/BoundaryConditions/BoundaryConditions.jl

@@ -7,7 +7,7 @@ export
     FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, DistributedCommunicationBoundaryCondition,
     validate_boundary_condition_topology, validate_boundary_condition_architecture,
     FieldBoundaryConditions,
-    apply_x_bcs!, apply_y_bcs!, apply_z_bcs!,
+    compute_x_bcs!, compute_y_bcs!, compute_z_bcs!,
     fill_halo_regions!
 
 using Adapt
@@ -35,7 +35,7 @@ include("fill_halo_regions_flux.jl")
 include("fill_halo_regions_nothing.jl")
 include("fill_halo_regions_zipper.jl")
 
-include("apply_flux_bcs.jl")
+include("compute_flux_bcs.jl")
 
 include("update_boundary_conditions.jl")
 include("polar_boundary_condition.jl")

src/BoundaryConditions/apply_flux_bcs.jl renamed to src/BoundaryConditions/compute_flux_bcs.jl

@@ -9,95 +9,95 @@ using Oceananigans.Grids: AbstractGrid
 #####
 
 # Unpack
-apply_x_bcs!(Gc, c, args...) = apply_x_bcs!(Gc, Gc.grid, c, c.boundary_conditions.west,   c.boundary_conditions.east, args...)
-apply_y_bcs!(Gc, c, args...) = apply_y_bcs!(Gc, Gc.grid, c, c.boundary_conditions.south,  c.boundary_conditions.north, args...)
-apply_z_bcs!(Gc, c, args...) = apply_z_bcs!(Gc, Gc.grid, c, c.boundary_conditions.bottom, c.boundary_conditions.top, args...)
+compute_x_bcs!(Gc, c, args...) = compute_x_bcs!(Gc, Gc.grid, c, c.boundary_conditions.west,   c.boundary_conditions.east, args...)
+compute_y_bcs!(Gc, c, args...) = compute_y_bcs!(Gc, Gc.grid, c, c.boundary_conditions.south,  c.boundary_conditions.north, args...)
+compute_z_bcs!(Gc, c, args...) = compute_z_bcs!(Gc, Gc.grid, c, c.boundary_conditions.bottom, c.boundary_conditions.top, args...)
 
 # Shortcuts for...
 #
 # Nothing tendencies.
-apply_x_bcs!(::Nothing, args...) = nothing
-apply_y_bcs!(::Nothing, args...) = nothing
-apply_z_bcs!(::Nothing, args...) = nothing
+compute_x_bcs!(::Nothing, args...) = nothing
+compute_y_bcs!(::Nothing, args...) = nothing
+compute_z_bcs!(::Nothing, args...) = nothing
 
 # Not-flux boundary conditions
 const NotFluxBC = Union{PBC, MCBC, DCBC, VBC, GBC, OBC, ZFBC, Nothing}
 
-apply_x_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
-apply_y_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
-apply_z_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
+compute_x_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
+compute_y_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
+compute_z_bcs!(Gc, ::AbstractGrid, c, ::NotFluxBC, ::NotFluxBC, ::AbstractArchitecture, args...) = nothing
 
 # The real deal
 """
 Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
 the source term `Gc` at the left and right.
 """
-apply_x_bcs!(Gc, grid::AbstractGrid, c, west_bc, east_bc, arch::AbstractArchitecture, args...) =
-    launch!(arch, grid, :yz, _apply_x_bcs!, Gc, instantiated_location(Gc), grid, west_bc, east_bc, Tuple(args))
+compute_x_bcs!(Gc, grid::AbstractGrid, c, west_bc, east_bc, arch::AbstractArchitecture, args...) =
+    launch!(arch, grid, :yz, _compute_x_bcs!, Gc, instantiated_location(Gc), grid, west_bc, east_bc, Tuple(args))
 
 """
 Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
 the source term `Gc` at the left and right.
 """
-apply_y_bcs!(Gc, grid::AbstractGrid, c, south_bc, north_bc, arch::AbstractArchitecture, args...) =
-    launch!(arch, grid, :xz, _apply_y_bcs!, Gc, instantiated_location(Gc), grid, south_bc, north_bc, Tuple(args))
+compute_y_bcs!(Gc, grid::AbstractGrid, c, south_bc, north_bc, arch::AbstractArchitecture, args...) =
+    launch!(arch, grid, :xz, _compute_y_bcs!, Gc, instantiated_location(Gc), grid, south_bc, north_bc, Tuple(args))
 
 """
 Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
 the source term `Gc` at the top and bottom.
 """
-apply_z_bcs!(Gc, grid::AbstractGrid, c, bottom_bc, top_bc, arch::AbstractArchitecture, args...) =
-    launch!(arch, grid, :xy, _apply_z_bcs!, Gc, instantiated_location(Gc), grid, bottom_bc, top_bc, Tuple(args))
+compute_z_bcs!(Gc, grid::AbstractGrid, c, bottom_bc, top_bc, arch::AbstractArchitecture, args...) =
+    launch!(arch, grid, :xy, _compute_z_bcs!, Gc, instantiated_location(Gc), grid, bottom_bc, top_bc, Tuple(args))
 
 """
-    _apply_x_bcs!(Gc, grid, west_bc, east_bc, args...)
+    _compute_x_bcs!(Gc, grid, west_bc, east_bc, args...)
 
 Apply a west and/or east boundary condition to variable `c`.
 """
-@kernel function _apply_x_bcs!(Gc, loc, grid, west_bc, east_bc, args)
+@kernel function _compute_x_bcs!(Gc, loc, grid, west_bc, east_bc, args)
     j, k = @index(Global, NTuple)
-    apply_x_west_bc!(Gc, loc, west_bc, j, k, grid, args...)
-    apply_x_east_bc!(Gc, loc, east_bc, j, k, grid, args...)
+    compute_x_west_bc!(Gc, loc, west_bc, j, k, grid, args...)
+    compute_x_east_bc!(Gc, loc, east_bc, j, k, grid, args...)
 end
 
 """
-    _apply_y_bcs!(Gc, grid, south_bc, north_bc, args...)
+    _compute_y_bcs!(Gc, grid, south_bc, north_bc, args...)
 
 Apply a south and/or north boundary condition to variable `c`.
 """
-@kernel function _apply_y_bcs!(Gc, loc, grid, south_bc, north_bc, args)
+@kernel function _compute_y_bcs!(Gc, loc, grid, south_bc, north_bc, args)
     i, k = @index(Global, NTuple)
-    apply_y_south_bc!(Gc, loc, south_bc, i, k, grid, args...)
-    apply_y_north_bc!(Gc, loc, north_bc, i, k, grid, args...)
+    compute_y_south_bc!(Gc, loc, south_bc, i, k, grid, args...)
+    compute_y_north_bc!(Gc, loc, north_bc, i, k, grid, args...)
 end
 
 """
-    _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
+    _compute_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
 
 Apply a top and/or bottom boundary condition to variable `c`.
 """
-@kernel function _apply_z_bcs!(Gc, loc, grid, bottom_bc, top_bc, args)
+@kernel function _compute_z_bcs!(Gc, loc, grid, bottom_bc, top_bc, args)
     i, j = @index(Global, NTuple)
-    apply_z_bottom_bc!(Gc, loc, bottom_bc, i, j, grid, args...)
-       apply_z_top_bc!(Gc, loc, top_bc,    i, j, grid, args...)
+    compute_z_bottom_bc!(Gc, loc, bottom_bc, i, j, grid, args...)
+       compute_z_top_bc!(Gc, loc, top_bc,    i, j, grid, args...)
 end
 
 # Shortcuts for zero flux or non-flux boundary conditions
-@inline   apply_x_east_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
-@inline   apply_x_west_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
-@inline  apply_y_north_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
-@inline  apply_y_south_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
-@inline    apply_z_top_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
-@inline apply_z_bottom_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline   compute_x_east_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline   compute_x_west_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline  compute_y_north_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline  compute_y_south_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline    compute_z_top_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
+@inline compute_z_bottom_bc!(Gc, loc, ::NotFluxBC, args...) = nothing
 
 # shortcut for the zipper BC
-@inline apply_y_north_bc!(Gc, loc, ::ZBC, args...) = nothing
+@inline compute_y_north_bc!(Gc, loc, ::ZBC, args...) = nothing
 
 @inline flip(::Center) = Face()
 @inline flip(::Face) = Center()
 
 """
-    apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
+    compute_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
 
 Add the flux divergence associated with a west flux boundary condition on `c`.
 Note that because
@@ -111,26 +111,26 @@ If `west_bc.condition` is a function, the function must have the signature
 
 The same logic holds for south and bottom boundary conditions in `y`, and `z`, respectively.
 """
-@inline function apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
+@inline function compute_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[1, j, k] += getbc(west_flux, j, k, grid, args...) * Ax(1, j, k, grid, flip(LX), LY, LZ) / volume(1, j, k, grid, LX, LY, LZ)
     return nothing
 end
 
-@inline function apply_y_south_bc!(Gc, loc, south_flux::BC{<:Flux}, i, k, grid, args...)
+@inline function compute_y_south_bc!(Gc, loc, south_flux::BC{<:Flux}, i, k, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[i, 1, k] += getbc(south_flux, i, k, grid, args...) * Ay(i, 1, k, grid, LX, flip(LY), LZ) / volume(i, 1, k, grid, LX, LY, LZ)
     return nothing
 end
 
-@inline function apply_z_bottom_bc!(Gc, loc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
+@inline function compute_z_bottom_bc!(Gc, loc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) * Az(i, j, 1, grid, LX, LY, flip(LZ)) / volume(i, j, 1, grid, LX, LY, LZ)
     return nothing
 end
 
 """
-    apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
+    compute_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
 
 Add the part of flux divergence associated with a east boundary condition on `c`.
 Note that because
@@ -144,19 +144,19 @@ If `east_bc.condition` is a function, the function must have the signature
 
 The same logic holds for north and top boundary conditions in `y`, and `z`, respectively.
 """
-@inline function apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
+@inline function compute_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[grid.Nx, j, k] -= getbc(east_flux, j, k, grid, args...) * Ax(grid.Nx+1, j, k, grid, flip(LX), LY, LZ) / volume(grid.Nx, j, k, grid, LX, LY, LZ)
     return nothing
 end
 
-@inline function apply_y_north_bc!(Gc, loc, north_flux::BC{<:Flux}, i, k, grid, args...)
+@inline function compute_y_north_bc!(Gc, loc, north_flux::BC{<:Flux}, i, k, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[i, grid.Ny, k] -= getbc(north_flux, i, k, grid, args...) * Ay(i, grid.Ny+1, k, grid, LX, flip(LY), LZ) / volume(i, grid.Ny, k, grid, LX, LY, LZ)
     return nothing
 end
 
-@inline function apply_z_top_bc!(Gc, loc, top_flux::BC{<:Flux}, i, j, grid, args...)
+@inline function compute_z_top_bc!(Gc, loc, top_flux::BC{<:Flux}, i, j, grid, args...)
     LX, LY, LZ = loc
     @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) * Az(i, j, grid.Nz+1, grid, LX, LY, flip(LZ)) / volume(i, j, grid.Nz, grid, LX, LY, LZ)
     return nothing

src/BoundaryConditions/field_boundary_conditions.jl

@@ -229,7 +229,15 @@ function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions,
                                               prognostic_names=nothing)
 
     loc = assumed_field_location(field_name)
+    return regularize_field_boundary_conditions(bcs, grid, loc, prognostic_names, field_name)
+end
 
+function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions,
+                                              grid::AbstractGrid,
+                                              loc::Tuple,
+                                              prognostic_names=nothing,
+                                              field_name=nothing)
+    
     west   = regularize_west_boundary_condition(bcs.west,     grid, loc, 1, LeftBoundary,  prognostic_names)
     east   = regularize_east_boundary_condition(bcs.east,     grid, loc, 1, RightBoundary, prognostic_names)
     south  = regularize_south_boundary_condition(bcs.south,   grid, loc, 2, LeftBoundary,  prognostic_names)