(0.94.0) Overhaul spacings functions to return `KernelFunctionOperati…

…on`s (#3143) * fix spacings methods for immersed boundary grids * Update src/ImmersedBoundaries/immersed_grid_metrics.jl Co-authored-by: Navid C. Constantinou <[email protected]> * Update src/ImmersedBoundaries/immersed_grid_metrics.jl Co-authored-by: Navid C. Constantinou <[email protected]> * add tests for spacing(s) on IBG * don't redefine x/y/zspacing methods * Correct merge `main` * Disambiguate `yspacings` * Cleanup? * Add `*spacing` functions that rely on `KernelFunctionOperation`s * Declutter `ImmersedBoundaries` module * Changes needed to fix tests * Fix tests * Bump v0.94.0 * spacings functions now return a `KernelFunctionOperation` * Address some review comments and fix OSSG tests * Address more review comments, update and fix tests * Fix more tests * Fix bug in Makie ext * Update src/Grids/rectilinear_grid.jl Co-authored-by: Gregory L. Wagner <[email protected]> * Update src/Grids/latitude_longitude_grid.jl Co-authored-by: Gregory L. Wagner <[email protected]> * Update src/Grids/orthogonal_spherical_shell_grid.jl Co-authored-by: Gregory L. Wagner <[email protected]> * Update src/Grids/orthogonal_spherical_shell_grid.jl Co-authored-by: Gregory L. Wagner <[email protected]> * Update src/Grids/orthogonal_spherical_shell_grid.jl Co-authored-by: Gregory L. Wagner <[email protected]> * Use correct locations for OSSG spacings functions * Simpler `scatterlines` plotting for vertical grid spacing * Forgot to use `first` * Define spacing functions for `::Nothing` using superscript `ᵃ` * Define methods only once * Update some tests * Need to define more spacing functions for `RectilinearGrid` * Also define `:ᵃ` derivative operators * Need to export some derivatives * Also use `:ᵃ` for immersed grid metrics. * Update vertical grid plotting in ocean wind mixing example * Trying to define the complete set without duplicates * Revert "Trying to define the complete set without duplicates" This reverts commit 735a16f. * Revert "Also use `:ᵃ` for immersed grid metrics." This reverts commit bbe6967. * Revert "Need to export some derivatives" This reverts commit 000473f. * Revert "Also define `:ᵃ` derivative operators" This reverts commit c748494. * Let's try this again * Reorganize and define some specialized spacings * Some more spacings * Fix doctests * Fix docs * Fix plotting examples in docs --------- Co-authored-by: Navid C. Constantinou <[email protected]> Co-authored-by: Ali Ramadhan <[email protected]> Co-authored-by: Gregory Wagner <[email protected]>
CliMA · Nov 14, 2024 · cb9e6d2 · cb9e6d2 · glwagner · Nov 14, 2024
1 parent bb42ddd
commit cb9e6d2
Show file tree

Hide file tree

Showing 26 changed files with 724 additions and 655 deletions.
diff --git a/docs/src/fields.md b/docs/src/fields.md
@@ -3,7 +3,7 @@
 `Field`s and its relatives are core Oceananigans data structures.
 `Field`s are more or less arrays of `data` located on a `grid`, whose entries
 correspond to the average value of some quantity over some finite-sized volume.
-`Field`s also may contain `boundary_conditions`, may be computed from an `operand` 
+`Field`s also may contain `boundary_conditions`, may be computed from an `operand`
 or expression involving other fields, and may cover only a portion of the total
 `indices` spanned by the grid.
 
@@ -43,10 +43,10 @@ primary mesh.
 For example, the primary or `Center` cell spacings in ``z`` are
 
 ```jldoctest fields
-zspacings(grid, Center())
+zspacings(grid, Center())[:, :, 1:4]
 
 # output
-4-element view(OffsetArray(::Vector{Float64}, 0:5), 1:4) with eltype Float64:
+4-element Vector{Float64}:
  0.1
  0.19999999999999998
  0.3
@@ -57,10 +57,10 @@ corresponding to cell interfaces located at `z = [0, 0.1, 0.3, 0.6, 1]`.
 But then for the grid which is staggered in `z` relative to the primary mesh,
 
 ```jldoctest fields
-zspacings(grid, Face())
+zspacings(grid, Face())[:, :, 1:5]
 
 # output
-5-element view(OffsetArray(::Vector{Float64}, -1:5), 1:5) with eltype Float64:
+5-element Vector{Float64}:
  0.1
  0.15000000000000002
  0.24999999999999994
@@ -277,7 +277,7 @@ Note that indexing into `c` is the same as indexing into `c.data`.
 
 ```jldoctest fields
 c[:, :, :] == c.data
- 
+
 # output
 true
 ```
@@ -333,7 +333,7 @@ non-`Flat` directions must be included. For example, to `set!` on a one-dimensio
 one_d_grid = RectilinearGrid(size=7, x=(0, 7), topology=(Periodic, Flat, Flat))
 one_d_c = CenterField(one_d_grid)
 
-# The one-dimensional grid varies only in `x` 
+# The one-dimensional grid varies only in `x`
 still_pretty_fun(x) = 3x
 set!(one_d_c, still_pretty_fun)
 

diff --git a/docs/src/grids.md b/docs/src/grids.md
@@ -219,7 +219,7 @@ current_figure()
 
 ## Once more with feeling
 
-In summary, making a grid requires 
+In summary, making a grid requires
 
 * The machine architecture, or whether data is stored on the CPU, GPU, or distributed across multiple devices or nodes.
 * Information about the domain geometry. Domains can take a variety of shapes, including
@@ -413,7 +413,7 @@ hidespines!(axy)
 
 axΔy = Axis(fig[2, 1]; xlabel = "y (m)", ylabel = "y-spacing (m)")
 scatter!(axΔy, yc, Δy)
-hidespines!(axΔy, :t, :r) 
+hidespines!(axΔy, :t, :r)
 
 axz = Axis(fig[3, 1], title="z-grid")
 lines!(axz, [-Lz, 0], [0, 0], color=:gray)
@@ -440,18 +440,17 @@ using Oceananigans
 
 ```@example latlon_nodes
 grid = LatitudeLongitudeGrid(size = (1, 44),
-                             longitude = (0, 1),   
+                             longitude = (0, 1),
                              latitude = (0, 88),
                              topology = (Bounded, Bounded, Flat))
 
-φ = φnodes(grid, Center())
 Δx = xspacings(grid, Center(), Center())
 
 using CairoMakie
 
 fig = Figure(size=(600, 400))
 ax = Axis(fig[1, 1], xlabel="Zonal spacing on 2 degree grid (km)", ylabel="Latitude (degrees)")
-scatter!(ax, Δx ./ 1e3, φ)
+scatter!(ax, Δx / 1e3)
 
 current_figure()
 ```
@@ -473,7 +472,7 @@ m = 2 # spacing at the equator in degrees
 function latitude_faces(j)
     if j == 1 # equator
         return 0
-    else # crudely estimate the location of the jth face 
+    else # crudely estimate the location of the jth face
         φ₋ = latitude_faces(j-1)
         φ′ = φ₋ + m * scale_factor(φ₋) / 2
         return φ₋ + m * scale_factor(φ′)
@@ -493,7 +492,7 @@ grid = LatitudeLongitudeGrid(size = (Nx, Ny),
 180×28×1 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo and with precomputed metrics
 ├── longitude: Bounded  λ ∈ [0.0, 360.0]   regularly spaced with Δλ=2.0
 ├── latitude:  Bounded  φ ∈ [0.0, 77.2679] variably spaced with min(Δφ)=2.0003, max(Δφ)=6.95319
-└── z:         Flat z                      
+└── z:         Flat z
 ```
 
 We've also illustrated the construction of a grid that is `Flat` in the vertical direction.
@@ -508,7 +507,7 @@ m = 2 # spacing at the equator in degrees
 function latitude_faces(j)
     if j == 1 # equator
         return 0
-    else # crudely estimate the location of the jth face 
+    else # crudely estimate the location of the jth face
         φ₋ = latitude_faces(j-1)
         φ′ = φ₋ + m * scale_factor(φ₋) / 2
         return φ₋ + m * scale_factor(φ′)
@@ -530,17 +529,17 @@ grid = LatitudeLongitudeGrid(size = (Nx, Ny),
 
 ```@example plot
 φ = φnodes(grid, Center())
-Δx = xspacings(grid, Center(), Center(), with_halos=true)[1:Ny]
-Δy = yspacings(grid, Center())[1:Ny]
+Δx = xspacings(grid, Center(), Center())[1, 1:Ny]
+Δy = yspacings(grid, Center(), Center())[1, 1:Ny]
 
 using CairoMakie
 
 fig = Figure(size=(800, 400), title="Spacings on a Mercator grid")
 axx = Axis(fig[1, 1], xlabel="Zonal spacing (km)", ylabel="Latitude (degrees)")
-scatter!(axx, Δx ./ 1e3, φ)
+scatter!(axx, Δx / 1e3, φ)
 
 axy = Axis(fig[1, 2], xlabel="Meridional spacing (km)")
-scatter!(axy, Δy ./ 1e3, φ)
+scatter!(axy, Δy / 1e3, φ)
 
 hidespines!(axx, :t, :r)
 hidespines!(axy, :t, :l, :r)

diff --git a/examples/ocean_wind_mixing_and_convection.jl b/examples/ocean_wind_mixing_and_convection.jl
@@ -50,13 +50,13 @@ h(k) = (k - 1) / Nz
 ## Linear near-surface generator
 ζ₀(k) = 1 + (h(k) - 1) / refinement
 
-## Bottom-intensified stretching function 
+## Bottom-intensified stretching function
 Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
 
 ## Generating function
 z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)
 
-grid = RectilinearGrid(size = (Nx, Nx, Nz), 
+grid = RectilinearGrid(size = (Nx, Nx, Nz),
                           x = (0, Lx),
                           y = (0, Ly),
                           z = z_faces)
@@ -66,8 +66,8 @@ grid = RectilinearGrid(size = (Nx, Nx, Nz),
 fig = Figure(size=(1200, 800))
 ax = Axis(fig[1, 1], ylabel = "Depth (m)", xlabel = "Vertical spacing (m)")
 
-lines!(ax, zspacings(grid, Center()), znodes(grid, Center()))
-scatter!(ax, zspacings(grid, Center()), znodes(grid, Center()))
+lines!(ax, zspacings(grid, Center()))
+scatter!(ax, zspacings(grid, Center()))
 
 current_figure() #hide
 fig

diff --git a/examples/tilted_bottom_boundary_layer.jl b/examples/tilted_bottom_boundary_layer.jl
@@ -3,7 +3,7 @@
 # away from a constant along-slope (y-direction) velocity constant density stratification.
 # This perturbation develops into a turbulent bottom boundary layer due to momentum
 # loss at the bottom boundary modeled with a quadratic drag law.
-# 
+#
 # This example illustrates
 #
 #   * changing the direction of gravitational acceleration in the buoyancy model;
@@ -33,15 +33,15 @@ Nz = 64
 ## Creates a grid with near-constant spacing `refinement * Lz / Nz`
 ## near the bottom:
 refinement = 1.8 # controls spacing near surface (higher means finer spaced)
-stretching = 10  # controls rate of stretching at bottom 
+stretching = 10  # controls rate of stretching at bottom
 
 ## "Warped" height coordinate
 h(k) = (Nz + 1 - k) / Nz
 
 ## Linear near-surface generator
 ζ(k) = 1 + (h(k) - 1) / refinement
 
-## Bottom-intensified stretching function 
+## Bottom-intensified stretching function
 Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
 
 ## Generating function
@@ -56,11 +56,9 @@ grid = RectilinearGrid(topology = (Periodic, Flat, Bounded),
 
 using CairoMakie
 
-lines(zspacings(grid, Center()), znodes(grid, Center()),
-      axis = (ylabel = "Depth (m)",
-              xlabel = "Vertical spacing (m)"))
-
-scatter!(zspacings(grid, Center()), znodes(grid, Center()))
+scatterlines(zspacings(grid, Center()),
+             axis = (ylabel = "Depth (m)",
+                     xlabel = "Vertical spacing (m)"))
 
 current_figure() #hide
 

diff --git a/ext/OceananigansMakieExt.jl b/ext/OceananigansMakieExt.jl
@@ -2,6 +2,7 @@ module OceananigansMakieExt
 
 using Oceananigans
 using Oceananigans.Grids: OrthogonalSphericalShellGrid
+using Oceananigans.Fields: AbstractField
 using Oceananigans.AbstractOperations: AbstractOperation
 using Oceananigans.Architectures: on_architecture
 using Oceananigans.ImmersedBoundaries: mask_immersed_field!
@@ -14,7 +15,7 @@ import Makie: args_preferred_axis
 # do not overstate a preference for being plotted in a 3D LScene.
 # Because often we are trying to plot 1D and 2D Field, even though
 # (perhaps incorrectly) all Field are AbstractArray{3}.
-args_preferred_axis(::Field) = nothing
+args_preferred_axis(::AbstractField) = nothing
 
 function drop_singleton_indices(N)
     if N == 1
@@ -146,3 +147,4 @@ function convert_arguments(pl::Type{<:AbstractPlot}, ξ1::AbstractArray, ξ2::Ab
 end
 
 end # module
+