From b19461beb02845cd182daf9acc29e3c0ddcf4071 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 16:23:21 -0800 Subject: [PATCH 01/47] small correction --- src/Grids/nodes_and_spacings.jl | 11 --- .../immersed_boundary_nodes.jl | 93 ++++++++++++++++++- 2 files changed, 90 insertions(+), 14 deletions(-) diff --git a/src/Grids/nodes_and_spacings.jl b/src/Grids/nodes_and_spacings.jl index fd814d94e5..21561e80ab 100644 --- a/src/Grids/nodes_and_spacings.jl +++ b/src/Grids/nodes_and_spacings.jl @@ -178,17 +178,6 @@ function φspacings end destantiate(::Face) = Face destantiate(::Center) = Center -spacings_function(::Val{:x}) = xspacings -spacings_function(::Val{:y}) = yspacings -spacings_function(::Val{:z}) = zspacings -spacings_function(::Val{:λ}) = λspacings -spacings_function(::Val{:φ}) = φspacings - -function minimum_spacing(s, grid, ℓx, ℓy, ℓz) - spacings = spacings_function(s) - return minimum(spacings(grid, ℓx, ℓy, ℓz)) -end - """ minimum_xspacing(grid, ℓx, ℓy, ℓz) diff --git a/src/ImmersedBoundaries/immersed_boundary_nodes.jl b/src/ImmersedBoundaries/immersed_boundary_nodes.jl index 62f39894f8..71b56f0aee 100644 --- a/src/ImmersedBoundaries/immersed_boundary_nodes.jl +++ b/src/ImmersedBoundaries/immersed_boundary_nodes.jl @@ -63,7 +63,94 @@ rname(ibg::IBG) = rname(ibg.underlying_grid) @inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid) @inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid) @inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid) +import Oceananigans.Grids: xspacings, yspacings, zspacings + +const c = Center() +const f = Face() + +@inline xnode(i, ibg::IBG, ℓx) = xnode(i, ibg.underlying_grid, ℓx) +@inline ynode(j, ibg::IBG, ℓy) = ynode(j, ibg.underlying_grid, ℓy) +@inline znode(k, ibg::IBG, ℓz) = znode(k, ibg.underlying_grid, ℓz) + +@inline λnode(i, ibg::IBG, ℓx) = λnode(i, ibg.underlying_grid, ℓx) +@inline φnode(j, ibg::IBG, ℓy) = φnode(i, ibg.underlying_grid, ℓy) + +@inline xnode(i, j, ibg::IBG, ℓx, ℓy) = xnode(i, j, ibg.underlying_grid, ℓx, ℓy) + +@inline xnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = xnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) +@inline ynode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ynode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) +@inline znode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = znode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) + +@inline λnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = λnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) +@inline φnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = φnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) + +@inline ξnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ξnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) +@inline ηnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ηnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) +@inline rnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = rnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) + +@inline node(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = node(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) + +nodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = nodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +nodes(ibg::IBG, (ℓx, ℓy, ℓz); kwargs...) = nodes(ibg, ℓx, ℓy, ℓz; kwargs...) + +xnodes(ibg::IBG, loc; kwargs...) = xnodes(ibg.underlying_grid, loc; kwargs...) +ynodes(ibg::IBG, loc; kwargs...) = ynodes(ibg.underlying_grid, loc; kwargs...) +znodes(ibg::IBG, loc; kwargs...) = znodes(ibg.underlying_grid, loc; kwargs...) + +λnodes(ibg::IBG, loc; kwargs...) = λnodes(ibg.underlying_grid, loc; kwargs...) +φnodes(ibg::IBG, loc; kwargs...) = φnodes(ibg.underlying_grid, loc; kwargs...) + +ξnodes(ibg::IBG, loc; kwargs...) = ξnodes(ibg.underlying_grid, loc; kwargs...) +ηnodes(ibg::IBG, loc; kwargs...) = ηnodes(ibg.underlying_grid, loc; kwargs...) +rnodes(ibg::IBG, loc; kwargs...) = rnodes(ibg.underlying_grid, loc; kwargs...) + +xnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = xnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +ynodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ynodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +znodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = znodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) + +λnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = λnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +φnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = φnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) + +ξnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ξnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +ηnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ηnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) +rnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = rnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) + +@inline cpu_face_constructor_x(ibg::IBG) = cpu_face_constructor_x(ibg.underlying_grid) +@inline cpu_face_constructor_y(ibg::IBG) = cpu_face_constructor_y(ibg.underlying_grid) +@inline cpu_face_constructor_z(ibg::IBG) = cpu_face_constructor_z(ibg.underlying_grid) + +node_names(ibg::IBG, ℓx, ℓy, ℓz) = node_names(ibg.underlying_grid, ℓx, ℓy, ℓz) + +ξname(ibg::IBG) = ξname(ibg.underlying_grid) +ηname(ibg::IBG) = ηname(ibg.underlying_grid) +rname(ibg::IBG) = rname(ibg.underlying_grid) + +@inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid) +@inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid) +@inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid) + +##### +##### Grid-specific grid spacings +##### + +const RGIBG{F, X, Y, Z} = ImmersedBoundaryGrid{F, X, Y, Z, <:RectilinearGrid} +const LLIBG{F, X, Y, Z} = ImmersedBoundaryGrid{F, X, Y, Z, <:LatitudeLongitudeGrid} +const OSIBG{F, X, Y, Z} = ImmersedBoundaryGrid{F, X, Y, Z, <:OrthogonalSphericalShellGrid} + +@inline xspacings(grid::RGIBG, ℓx) = xspacings(grid, ℓx, nothing, nothing) +@inline yspacings(grid::RGIBG, ℓy) = yspacings(grid, nothing, ℓy, nothing) +@inline zspacings(grid::RGIBG, ℓz) = zspacings(grid, nothing, nothing, ℓz) + +@inline xspacings(grid::LLIBG, ℓx, ℓy) = xspacings(grid, ℓx, ℓy, nothing) +@inline yspacings(grid::LLIBG, ℓx, ℓy) = yspacings(grid, ℓx, ℓy, nothing) +@inline zspacings(grid::LLIBG, ℓz) = zspacings(grid, nothing, nothing, ℓz) + +@inline xspacings(grid::OSIBG, ℓx, ℓy) = xspacings(grid, ℓx, ℓy, nothing) +@inline yspacings(grid::OSIBG, ℓx, ℓy) = yspacings(grid, ℓx, ℓy, nothing) +@inline zspacings(grid::OSIBG, ℓz) = zspacings(grid, nothing, nothing, ℓz) + +@inline λspacings(grid::LLIBG, ℓx) = λspacings(grid, ℓx, nothing, nothing) +@inline φspacings(grid::LLIBG, ℓy) = φspacings(grid, nothing, ℓy, nothing) -xspacings(ibg::ImmersedBoundaryGrid, ℓ...) = xspacings(ibg.underlying_grid, ℓ...) -yspacings(ibg::ImmersedBoundaryGrid, ℓ...) = yspacings(ibg.underlying_grid, ℓ...) -zspacings(ibg::ImmersedBoundaryGrid, ℓ...) = zspacings(ibg.underlying_grid, ℓ...) +@inline λspacings(grid::OSIBG, ℓx, ℓy) = λspacings(grid, ℓx, ℓy, nothing) +@inline φspacings(grid::OSIBG, ℓx, ℓy) = φspacings(grid, ℓx, ℓy, nothing) From 73b9b4a7ac626d176d6a5827ee8dadbd1ff06373 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 16:23:59 -0800 Subject: [PATCH 02/47] correction --- src/MultiRegion/multi_region_grid.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/MultiRegion/multi_region_grid.jl b/src/MultiRegion/multi_region_grid.jl index 264b6a29be..69d0b9410d 100644 --- a/src/MultiRegion/multi_region_grid.jl +++ b/src/MultiRegion/multi_region_grid.jl @@ -3,7 +3,7 @@ using Oceananigans.ImmersedBoundaries: GridFittedBottom, PartialCellBottom, Grid import Oceananigans.Grids: architecture, size, new_data, halo_size import Oceananigans.Grids: with_halo, on_architecture -import Oceananigans.Grids: minimum_spacing, destantiate +import Oceananigans.Grids: destantiate import Oceananigans.Grids: minimum_xspacing, minimum_yspacing, minimum_zspacing import Oceananigans.Models.HydrostaticFreeSurfaceModels: default_free_surface import Oceananigans.DistributedComputations: reconstruct_global_grid @@ -39,9 +39,6 @@ const MultiRegionGrids = Union{MultiRegionGrid, ImmersedMultiRegionGrid} @inline Base.lastindex(mrg::MultiRegionGrids) = length(mrg) number_of_regions(mrg::MultiRegionGrids) = lastindex(mrg) -minimum_spacing(dir, grid::MultiRegionGrid, ℓx, ℓy, ℓz) = - minimum(minimum_spacing(dir, grid[r], ℓx, ℓy, ℓz) for r in 1:number_of_regions(grid)) - minimum_xspacing(grid::MultiRegionGrid, ℓx, ℓy, ℓz) = minimum(minimum_xspacing(grid[r], ℓx, ℓy, ℓz) for r in 1:number_of_regions(grid)) From 1c1dea377e9b53bae666364269563b51efc6615f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 16:31:23 -0800 Subject: [PATCH 03/47] add a test --- test/test_grids.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_grids.jl b/test/test_grids.jl index b30b48c1bb..2cf965895b 100644 --- a/test/test_grids.jl +++ b/test/test_grids.jl @@ -660,6 +660,11 @@ function test_lat_lon_xyzλφ_node_nodes(FT, arch) @test minimum_yspacing(grid) / grid.radius ≈ FT(π/6) @test minimum_zspacing(grid) ≈ 5 + grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> y < 20 && y > -20)) + + @test minimum_xspacing(grid, Face(), Face(), Face()) / grid.radius ≈ π/6 * cosd(15) + @test minimum_xspacing(grid) / grid.radius ≈ π/6 * cosd(15) + return nothing end From 7e5b287e13750cdc0b81cdf36bdb81c69b30ed55 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 16:32:57 -0800 Subject: [PATCH 04/47] weird --- .../immersed_boundary_nodes.jl | 67 +------------------ 1 file changed, 1 insertion(+), 66 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_boundary_nodes.jl b/src/ImmersedBoundaries/immersed_boundary_nodes.jl index 71b56f0aee..f5f11877f3 100644 --- a/src/ImmersedBoundaries/immersed_boundary_nodes.jl +++ b/src/ImmersedBoundaries/immersed_boundary_nodes.jl @@ -60,71 +60,6 @@ node_names(ibg::IBG, ℓx, ℓy, ℓz) = node_names(ibg.underlying_grid, ℓx, ηname(ibg::IBG) = ηname(ibg.underlying_grid) rname(ibg::IBG) = rname(ibg.underlying_grid) -@inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid) -@inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid) -@inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid) -import Oceananigans.Grids: xspacings, yspacings, zspacings - -const c = Center() -const f = Face() - -@inline xnode(i, ibg::IBG, ℓx) = xnode(i, ibg.underlying_grid, ℓx) -@inline ynode(j, ibg::IBG, ℓy) = ynode(j, ibg.underlying_grid, ℓy) -@inline znode(k, ibg::IBG, ℓz) = znode(k, ibg.underlying_grid, ℓz) - -@inline λnode(i, ibg::IBG, ℓx) = λnode(i, ibg.underlying_grid, ℓx) -@inline φnode(j, ibg::IBG, ℓy) = φnode(i, ibg.underlying_grid, ℓy) - -@inline xnode(i, j, ibg::IBG, ℓx, ℓy) = xnode(i, j, ibg.underlying_grid, ℓx, ℓy) - -@inline xnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = xnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) -@inline ynode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ynode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) -@inline znode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = znode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) - -@inline λnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = λnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) -@inline φnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = φnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) - -@inline ξnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ξnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) -@inline ηnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = ηnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) -@inline rnode(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = rnode(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) - -@inline node(i, j, k, ibg::IBG, ℓx, ℓy, ℓz) = node(i, j, k, ibg.underlying_grid, ℓx, ℓy, ℓz) - -nodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = nodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -nodes(ibg::IBG, (ℓx, ℓy, ℓz); kwargs...) = nodes(ibg, ℓx, ℓy, ℓz; kwargs...) - -xnodes(ibg::IBG, loc; kwargs...) = xnodes(ibg.underlying_grid, loc; kwargs...) -ynodes(ibg::IBG, loc; kwargs...) = ynodes(ibg.underlying_grid, loc; kwargs...) -znodes(ibg::IBG, loc; kwargs...) = znodes(ibg.underlying_grid, loc; kwargs...) - -λnodes(ibg::IBG, loc; kwargs...) = λnodes(ibg.underlying_grid, loc; kwargs...) -φnodes(ibg::IBG, loc; kwargs...) = φnodes(ibg.underlying_grid, loc; kwargs...) - -ξnodes(ibg::IBG, loc; kwargs...) = ξnodes(ibg.underlying_grid, loc; kwargs...) -ηnodes(ibg::IBG, loc; kwargs...) = ηnodes(ibg.underlying_grid, loc; kwargs...) -rnodes(ibg::IBG, loc; kwargs...) = rnodes(ibg.underlying_grid, loc; kwargs...) - -xnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = xnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -ynodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ynodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -znodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = znodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) - -λnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = λnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -φnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = φnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) - -ξnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ξnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -ηnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = ηnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) -rnodes(ibg::IBG, ℓx, ℓy, ℓz; kwargs...) = rnodes(ibg.underlying_grid, ℓx, ℓy, ℓz; kwargs...) - -@inline cpu_face_constructor_x(ibg::IBG) = cpu_face_constructor_x(ibg.underlying_grid) -@inline cpu_face_constructor_y(ibg::IBG) = cpu_face_constructor_y(ibg.underlying_grid) -@inline cpu_face_constructor_z(ibg::IBG) = cpu_face_constructor_z(ibg.underlying_grid) - -node_names(ibg::IBG, ℓx, ℓy, ℓz) = node_names(ibg.underlying_grid, ℓx, ℓy, ℓz) - -ξname(ibg::IBG) = ξname(ibg.underlying_grid) -ηname(ibg::IBG) = ηname(ibg.underlying_grid) -rname(ibg::IBG) = rname(ibg.underlying_grid) - @inline fractional_x_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_x_index(x, locs, grid.underlying_grid) @inline fractional_y_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_y_index(x, locs, grid.underlying_grid) @inline fractional_z_index(x, locs, grid::ImmersedBoundaryGrid) = fractional_z_index(x, locs, grid.underlying_grid) @@ -153,4 +88,4 @@ const OSIBG{F, X, Y, Z} = ImmersedBoundaryGrid{F, X, Y, Z, <:OrthogonalSpherical @inline φspacings(grid::LLIBG, ℓy) = φspacings(grid, nothing, ℓy, nothing) @inline λspacings(grid::OSIBG, ℓx, ℓy) = λspacings(grid, ℓx, ℓy, nothing) -@inline φspacings(grid::OSIBG, ℓx, ℓy) = φspacings(grid, ℓx, ℓy, nothing) +@inline φspacings(grid::OSIBG, ℓx, ℓy) = φspacings(grid, ℓx, ℓy, nothing) \ No newline at end of file From b578440079163e06d45b0975bccd8c861b6a39e6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 17:21:51 -0800 Subject: [PATCH 05/47] introduce an error for dz --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 2 +- src/ImmersedBoundaries/partial_cell_bottom.jl | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 48cb239768..fd968a0c3a 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,7 +10,7 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) +for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) for dir in (:x, :y, :z), operator in (:Δ, :A) metric = Symbol(operator, dir, LX, LY, LZ) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 36cc32487b..3cab9e6183 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -202,3 +202,10 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial @inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) @inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg) @inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) + +# Make sure a PartialCellBottom grid cannot use :ᵃ for the z-spacings +# in the z-direction +for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) + metric = Symbol(operator, dir, LX, LY, :ᵃ) + @eval @inline $metric(i, j, k, ibg::PCBIBG) = throw(ArgumentError("PartialCellBottom grids cannot use `:ᵃ` for spacings in the z-direction")) +end From 63c374f48f801aa03ea85b631a90cd3af74b4793 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 20:52:46 -0800 Subject: [PATCH 06/47] bugfix --- .../immersed_grid_metrics.jl | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index fd968a0c3a..2157e89ba9 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,20 +10,22 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) - for dir in (:x, :y, :z), operator in (:Δ, :A) - - metric = Symbol(operator, dir, LX, LY, LZ) +for dir in (:x, :y, :z) + for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) + spacing = Symbol(:Δ, dir, LX, LY, LZ) @eval begin - import Oceananigans.Operators: $metric - @inline $metric(i, j, k, ibg::IBG) = $metric(i, j, k, ibg.underlying_grid) + import Oceananigans.Operators: $spacing + @inline $spacing(i, j, k, ibg::IBG) = $spacing(i, j, k, ibg.underlying_grid) end end - - volume = Symbol(:V, LX, LY, LZ) - @eval begin - import Oceananigans.Operators: $volume - @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid) + for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) + area = Symbol(:A, dir, LX, LY, LZ) + volume = Symbol(:V, LX, LY, LZ) + @eval begin + import Oceananigans.Operators: $area, $volume + @inline $area(i, j, k, ibg::IBG) = $area(i, j, k, ibg.underlying_grid) + @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid) + end end end From 93e6305e132ade777f3107274de695faf0cff80e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 10:55:47 +0100 Subject: [PATCH 07/47] improve grid metrics --- .../immersed_grid_metrics.jl | 6 ++- src/ImmersedBoundaries/partial_cell_bottom.jl | 9 +--- .../spacings_and_areas_and_volumes.jl | 44 +++++++++---------- 3 files changed, 26 insertions(+), 33 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 2157e89ba9..da3adbc9c6 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,8 +10,12 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. +spacing_x_loc(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +spacing_y_loc(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +spacing_z_loc(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) + for dir in (:x, :y, :z) - for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) + for LX in spacing_x_loc(dir), LY in spacing_y_loc(dir), LZ in spacing_z_loc(dir) spacing = Symbol(:Δ, dir, LX, LY, LZ) @eval begin import Oceananigans.Operators: $spacing diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 3cab9e6183..214ab40619 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -201,11 +201,4 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial @inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) @inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg) -@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) - -# Make sure a PartialCellBottom grid cannot use :ᵃ for the z-spacings -# in the z-direction -for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) - metric = Symbol(operator, dir, LX, LY, :ᵃ) - @eval @inline $metric(i, j, k, ibg::PCBIBG) = throw(ArgumentError("PartialCellBottom grids cannot use `:ᵃ` for spacings in the z-direction")) -end +@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) \ No newline at end of file diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 649be89975..613070bbe7 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -35,32 +35,28 @@ The operators in this file fall into three categories: ##### ##### -# Convenience Functions for all grids -for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) +spacing_x_loc(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +spacing_y_loc(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +spacing_z_loc(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) - x_spacing_1D = Symbol(:Δx, LX, :ᵃ, :ᵃ) - x_spacing_2D = Symbol(:Δx, LX, LY, :ᵃ) - - y_spacing_1D = Symbol(:Δy, :ᵃ, LY, :ᵃ) - y_spacing_2D = Symbol(:Δy, LX, LY, :ᵃ) +spacing_1D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, :ᵃ, :ᵃ) +spacing_2D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, LY, :ᵃ) +spacing_1D(::Val{:y}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, LY, :ᵃ) +spacing_2D(::Val{:y}, LX, LY, LZ) = Symbol(:Δx, LX, LY, :ᵃ) +spacing_1D(::Val{:z}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, :ᵃ, LZ) +spacing_2D(::Val{:z}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, LY, LZ) - @eval begin - @inline $x_spacing_2D(i, j, k, grid) = $x_spacing_1D(i, j, k, grid) - @inline $y_spacing_2D(i, j, k, grid) = $y_spacing_1D(i, j, k, grid) - end - - for LZ in (:ᶜ, :ᶠ) - x_spacing_3D = Symbol(:Δx, LX, LY, LZ) - y_spacing_3D = Symbol(:Δy, LX, LY, LZ) - - z_spacing_1D = Symbol(:Δz, :ᵃ, :ᵃ, LZ) - z_spacing_3D = Symbol(:Δz, LX, LY, LZ) - - @eval begin - @inline $x_spacing_3D(i, j, k, grid) = $x_spacing_2D(i, j, k, grid) - @inline $y_spacing_3D(i, j, k, grid) = $y_spacing_2D(i, j, k, grid) - @inline $z_spacing_3D(i, j, k, grid) = $z_spacing_1D(i, j, k, grid) - end +# Convenience Functions for all grids +# This metaprogramming loop defines all the allowed combinations of Δx, Δy, and Δz +# Note `:ᵃ` is not allowed for the location associated with the spacing +for dir in (:x, :y, :z) + for LX in spacing_x_loc(dir), LY in spacing_y_loc(dir), LZ in spacing_z_loc(dir) + spacing1D = spacing_1D(Val(dir), LX, LY, LZ) + spacing2D = spacing_2D(Val(dir), LX, LY, LZ) + spacing_3D = Symbol(:Δ, dir, LX, LY, LZ) + + @eval @inline $spacing2D(i, j, k, grid) = $spacing1D(i, j, k, grid) + @eval @inline $spacing3D(i, j, k, grid) = $spacing2D(i, j, k, grid) end end From da88ed29e86e78c8602814881694f8c1ef15501d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 14:49:42 +0100 Subject: [PATCH 08/47] =?UTF-8?q?remove=20=CF=88spacing...?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/AbstractOperations/grid_metrics.jl | 132 ++---------------- ...apolation_open_boundary_matching_scheme.jl | 38 ++--- src/Grids/Grids.jl | 1 - src/Grids/nodes_and_spacings.jl | 6 - src/Operators/Operators.jl | 9 +- .../spacings_and_areas_and_volumes.jl | 20 ++- 6 files changed, 50 insertions(+), 156 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 4205135ef0..469aac68d1 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -5,118 +5,9 @@ using Oceananigans.Fields: AbstractField, default_indices, location 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() - -""" - metric_function(loc, metric::AbstractGridMetric) - -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...) + metric_function_symbol = Symbol(metric, code...) return getglobal(@__MODULE__, metric_function_symbol) end @@ -137,7 +28,6 @@ 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) @@ -165,13 +55,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 27 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 +81,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 27 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,13 +107,13 @@ 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 27 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 @@ -246,13 +136,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 5 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 @@ -275,7 +165,7 @@ 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 5 methods) └── arguments: ("Center", "Face", "Center") ``` """ diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index b91ac36249..fbbcb748f3 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -1,4 +1,4 @@ -using Oceananigans.Grids: xspacing, yspacing, zspacing +using Oceananigans.Operators: Δx, Δy, Δz """ FlatExtrapolation @@ -62,9 +62,9 @@ end const c = Center() @inline function _fill_west_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δx₁ = xspacing(1, j, k, grid, c, c, c) - Δx₂ = xspacing(2, j, k, grid, c, c, c) - Δx₃ = xspacing(3, j, k, grid, c, c, c) + Δx₁ = Δx(1, j, k, grid, c, c, c) + Δx₂ = Δx(2, j, k, grid, c, c, c) + Δx₃ = Δx(3, j, k, grid, c, c, c) spacing_factor = Δx₁ / (Δx₂ + Δx₃) @@ -78,9 +78,9 @@ end @inline function _fill_east_open_halo!(j, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) i = grid.Nx + 1 - Δx₁ = xspacing(i-1, j, k, grid, c, c, c) - Δx₂ = xspacing(i-2, j, k, grid, c, c, c) - Δx₃ = xspacing(i-3, j, k, grid, c, c, c) + Δx₁ = Δx(i-1, j, k, grid, c, c, c) + Δx₂ = Δx(i-2, j, k, grid, c, c, c) + Δx₃ = Δx(i-3, j, k, grid, c, c, c) spacing_factor = Δx₁ / (Δx₂ + Δx₃) @@ -92,9 +92,9 @@ end end @inline function _fill_south_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δy₁ = yspacing(i, 1, k, grid, c, c, c) - Δy₂ = yspacing(i, 2, k, grid, c, c, c) - Δy₃ = yspacing(i, 3, k, grid, c, c, c) + Δy₁ = Δy(i, 1, k, grid, c, c, c) + Δy₂ = Δy(i, 2, k, grid, c, c, c) + Δy₃ = Δy(i, 3, k, grid, c, c, c) spacing_factor = Δy₁ / (Δy₂ + Δy₃) @@ -108,9 +108,9 @@ end @inline function _fill_north_open_halo!(i, k, grid, ϕ, bc::FEOBC, loc, clock, model_fields) j = grid.Ny + 1 - Δy₁ = yspacing(i, j-1, k, grid, c, c, c) - Δy₂ = yspacing(i, j-2, k, grid, c, c, c) - Δy₃ = yspacing(i, j-3, k, grid, c, c, c) + Δy₁ = Δy(i, j-1, k, grid, c, c, c) + Δy₂ = Δy(i, j-2, k, grid, c, c, c) + Δy₃ = Δy(i, j-3, k, grid, c, c, c) spacing_factor = Δy₁ / (Δy₂ + Δy₃) @@ -122,9 +122,9 @@ end end @inline function _fill_bottom_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) - Δz₁ = zspacing(i, j, 1, grid, c, c, c) - Δz₂ = zspacing(i, j, 2, grid, c, c, c) - Δz₃ = zspacing(i, j, 3, grid, c, c, c) + Δz₁ = Δz(i, j, 1, grid, c, c, c) + Δz₂ = Δz(i, j, 2, grid, c, c, c) + Δz₃ = Δz(i, j, 3, grid, c, c, c) spacing_factor = Δz₁ / (Δz₂ + Δz₃) @@ -138,9 +138,9 @@ end @inline function _fill_top_open_halo!(i, j, grid, ϕ, bc::FEOBC, loc, clock, model_fields) k = grid.Nz + 1 - Δz₁ = zspacing(i, j, k-1, grid, c, c, c) - Δz₂ = zspacing(i, j, k-2, grid, c, c, c) - Δz₃ = zspacing(i, j, k-3, grid, c, c, c) + Δz₁ = Δz(i, j, k-1, grid, c, c, c) + Δz₂ = Δz(i, j, k-2, grid, c, c, c) + Δz₃ = Δz(i, j, k-3, grid, c, c, c) spacing_factor = Δz₁ / (Δz₂ + Δz₃) diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl index 519233561b..8f43052461 100644 --- a/src/Grids/Grids.jl +++ b/src/Grids/Grids.jl @@ -16,7 +16,6 @@ export ξnode, ηnode, rnode export xnode, ynode, znode, λnode, φnode export xnodes, ynodes, znodes, λnodes, φnodes export spacings -export xspacing, yspacing, zspacing, λspacing, φspacing export xspacings, yspacings, zspacings, λspacings, φspacings export minimum_xspacing, minimum_yspacing, minimum_zspacing export static_column_depthᶜᶜᵃ, static_column_depthᶠᶜᵃ, static_column_depthᶜᶠᵃ, static_column_depthᶠᶠᵃ diff --git a/src/Grids/nodes_and_spacings.jl b/src/Grids/nodes_and_spacings.jl index 21561e80ab..21f93cd729 100644 --- a/src/Grids/nodes_and_spacings.jl +++ b/src/Grids/nodes_and_spacings.jl @@ -163,12 +163,6 @@ nodes(grid::AbstractGrid, (ℓx, ℓy, ℓz); reshape=false, with_halos=false) = # placeholders # see Oceananigans/AbstractOperations/grid_metrics.jl for definitions -function xspacing end -function yspacing end -function zspacing end -function λspacing end -function φspacing end - function xspacings end function yspacings end function zspacings end diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index a033f62d9e..4d18f10f51 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -1,5 +1,8 @@ module Operators +# General spacings +export Δx, Δy, Δz, Δλ, Δφ + # Spacings export Δxᶠᶠᶠ, Δxᶠᶠᶜ, Δxᶠᶜᶠ, Δxᶠᶜᶜ, Δxᶜᶠᶠ, Δxᶜᶠᶜ, Δxᶜᶜᶠ, Δxᶜᶜᶜ export Δyᶠᶠᶠ, Δyᶠᶠᶜ, Δyᶠᶜᶠ, Δyᶠᶜᶜ, Δyᶜᶠᶠ, Δyᶜᶠᶜ, Δyᶜᶜᶠ, Δyᶜᶜᶜ @@ -75,18 +78,12 @@ export intrinsic_vector, extrinsic_vector using Oceananigans.Grids -import Oceananigans.Grids: xspacing, yspacing, zspacing - ##### ##### Convenient aliases ##### const AG = AbstractGrid -const Δx = xspacing -const Δy = yspacing -const Δz = zspacing - const RG = RectilinearGrid const RGX = XRegularRG const RGY = YRegularRG diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 613070bbe7..6ca30a9248 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -51,9 +51,9 @@ spacing_2D(::Val{:z}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, LY, LZ) # Note `:ᵃ` is not allowed for the location associated with the spacing for dir in (:x, :y, :z) for LX in spacing_x_loc(dir), LY in spacing_y_loc(dir), LZ in spacing_z_loc(dir) - spacing1D = spacing_1D(Val(dir), LX, LY, LZ) - spacing2D = spacing_2D(Val(dir), LX, LY, LZ) - spacing_3D = Symbol(:Δ, dir, LX, LY, LZ) + spacing1D = spacing_1D(Val(dir), LX, LY, LZ) + spacing2D = spacing_2D(Val(dir), LX, LY, LZ) + spacing3D = Symbol(:Δ, dir, LX, LY, LZ) @eval @inline $spacing2D(i, j, k, grid) = $spacing1D(i, j, k, grid) @eval @inline $spacing3D(i, j, k, grid) = $spacing2D(i, j, k, grid) @@ -293,3 +293,17 @@ for LX in (:Center, :Face, :Nothing) end end end + +# Special curvilinear spacings for curvilinear grids +@inline Δλ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δλᶜᵃᵃ[i] +@inline Δλ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δλᶠᵃᵃ[i] +@inline Δλ(i, j, k, grid::LLGX, ::Center, ℓy, ℓz) = @inbounds grid.Δλᶜᵃᵃ +@inline Δλ(i, j, k, grid::LLGX, ::Face, ℓy, ℓz) = @inbounds grid.Δλᶠᵃᵃ + +@inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ[j] +@inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ[j] +@inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ +@inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ + +@inline Δλ(i, j, k, grid::OSSG, ℓx, ℓy, ℓz) = δxᶠᵃᵃ(i, j, k, grid, λnode, flip(ℓx), ℓy, ℓz) +@inline Δφ(i, j, k, grid::OSSG, ℓx, ℓy, ℓz) = δxᶠᵃᵃ(i, j, k, grid, λnode, ℓx, flip(ℓy), ℓz) From e129e212b2cb0e7e77f7b4d2b5f438481aaf4b8e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 15:05:46 +0100 Subject: [PATCH 09/47] remove `AbstractMetric` --- src/AbstractOperations/binary_operations.jl | 2 + src/AbstractOperations/grid_metrics.jl | 42 ++++++++++++++++++- src/Operators/Operators.jl | 4 +- .../spacings_and_areas_and_volumes.jl | 17 ++++---- 4 files changed, 55 insertions(+), 10 deletions(-) diff --git a/src/AbstractOperations/binary_operations.jl b/src/AbstractOperations/binary_operations.jl index 45956d239f..780a2e6c5e 100644 --- a/src/AbstractOperations/binary_operations.jl +++ b/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} diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 469aac68d1..8913de1e36 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -5,6 +5,17 @@ using Oceananigans.Fields: AbstractField, default_indices, location import Oceananigans.Grids: xspacings, yspacings, zspacings, λspacings, φspacings +const AbstractGridMetric = Union{typeof(Δx), typeof(Δy), typeof(Δz), + typeof(Ax), typeof(Ay), typeof(Az), + typeof(volume), + typeof(Δλ), typeof(Δφ)} + + +""" + metric_function(loc, metric::AbstractGridMetric) + +Return the function associated with `metric::AbstractGridMetric` at `loc`ation. +""" function metric_function(loc, metric) code = Tuple(interpolation_code(ℓ) for ℓ in loc) metric_function_symbol = Symbol(metric, code...) @@ -32,7 +43,36 @@ on_architecture(to, gm::GridMetricOperation{LX, LY, LZ}) where {LX, LY, LZ} = indices(::GridMetricOperation) = default_indices(3) -# Special constructor for BinaryOperation +""" + gm = 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 +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 +``` +""" GridMetricOperation(L, metric, grid) = GridMetricOperation{L[1], L[2], L[3]}(metric_function(L, metric), grid) ##### diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index 4d18f10f51..aed8286c6f 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -1,7 +1,7 @@ module Operators -# General spacings -export Δx, Δy, Δz, Δλ, Δφ +# General metric operators +export Δx, Δy, Δz, Δλ, Δφ, Ax, Ay, Az, volume # Spacings export Δxᶠᶠᶠ, Δxᶠᶠᶜ, Δxᶠᶜᶠ, Δxᶠᶜᶜ, Δxᶜᶠᶠ, Δxᶜᶠᶜ, Δxᶜᶜᶠ, Δxᶜᶜᶜ diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 6ca30a9248..9a5135b3de 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -41,10 +41,10 @@ spacing_z_loc(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) spacing_1D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, :ᵃ, :ᵃ) spacing_2D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, LY, :ᵃ) -spacing_1D(::Val{:y}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, LY, :ᵃ) -spacing_2D(::Val{:y}, LX, LY, LZ) = Symbol(:Δx, LX, LY, :ᵃ) -spacing_1D(::Val{:z}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, :ᵃ, LZ) -spacing_2D(::Val{:z}, LX, LY, LZ) = Symbol(:Δx, :ᵃ, LY, LZ) +spacing_1D(::Val{:y}, LX, LY, LZ) = Symbol(:Δy, :ᵃ, LY, :ᵃ) +spacing_2D(::Val{:y}, LX, LY, LZ) = Symbol(:Δy, LX, LY, :ᵃ) +spacing_1D(::Val{:z}, LX, LY, LZ) = Symbol(:Δz, :ᵃ, :ᵃ, LZ) +spacing_2D(::Val{:z}, LX, LY, LZ) = Symbol(:Δz, :ᵃ, LY, LZ) # Convenience Functions for all grids # This metaprogramming loop defines all the allowed combinations of Δx, Δy, and Δz @@ -54,9 +54,12 @@ for dir in (:x, :y, :z) spacing1D = spacing_1D(Val(dir), LX, LY, LZ) spacing2D = spacing_2D(Val(dir), LX, LY, LZ) spacing3D = Symbol(:Δ, dir, LX, LY, LZ) - - @eval @inline $spacing2D(i, j, k, grid) = $spacing1D(i, j, k, grid) - @eval @inline $spacing3D(i, j, k, grid) = $spacing2D(i, j, k, grid) + if spacing2D != spacing1D + @eval @inline $spacing2D(i, j, k, grid) = $spacing1D(i, j, k, grid) + end + if spacing3D != spacing2D + @eval @inline $spacing3D(i, j, k, grid) = $spacing2D(i, j, k, grid) + end end end From 11ac809e086fc59978d20dca60d3b1bf32532adf Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 15:11:09 +0100 Subject: [PATCH 10/47] better formatting --- src/AbstractOperations/grid_metrics.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 8913de1e36..db5d6625c9 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -5,11 +5,15 @@ using Oceananigans.Fields: AbstractField, default_indices, location import Oceananigans.Grids: xspacings, yspacings, zspacings, λspacings, φspacings -const AbstractGridMetric = Union{typeof(Δx), typeof(Δy), typeof(Δz), - typeof(Ax), typeof(Ay), typeof(Az), - typeof(volume), - typeof(Δλ), typeof(Δφ)} - +const AbstractGridMetric = Union{typeof(Δx), + typeof(Δy), + typeof(Δz), + typeof(Δλ), + typeof(Δφ), + typeof(Ax), + typeof(Ay), + typeof(Az), + typeof(volume)} """ metric_function(loc, metric::AbstractGridMetric) From ed636bc52d0e7bcde4a260d1beea660c3802e364 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 15:23:59 +0100 Subject: [PATCH 11/47] continue with the spacings --- src/Operators/Operators.jl | 3 ++ .../spacings_and_areas_and_volumes.jl | 53 +++++++++---------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index aed8286c6f..15a3a1ad8d 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -97,6 +97,9 @@ const LLGX = XRegularLLG const LLGY = YRegularLLG const LLGZ = ZRegularLLG +# Vertically regular grids +const ZRG = Union{RGZ, OSSGZ, LLGZ} + # On the fly calculations of metrics const LLGF = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing} const LLGFX = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Number} diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 9a5135b3de..225787ea71 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -63,6 +63,26 @@ for dir in (:x, :y, :z) end end +##### +##### Vertical spacings (same for all grids) +##### + +@inline Δzᵃᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] +@inline Δzᵃᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] + +@inline Δzᶜᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] +@inline Δzᶠᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] +@inline Δzᶜᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] +@inline Δzᶠᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] + +@inline Δzᵃᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ +@inline Δzᵃᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ + +@inline Δzᶜᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ +@inline Δzᶠᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ +@inline Δzᶜᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ +@inline Δzᶠᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ + ##### ##### Rectilinear Grids (Flat grids already have Δ = 1) ##### @@ -82,8 +102,6 @@ end @inline Δzᵃᵃᶠ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶠ[k] @inline Δzᵃᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶜ[k] @inline Δzᶜᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᶠᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᶜᵃᶠ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶠ[k] ## XRegularRG @@ -103,15 +121,6 @@ end @inline Δyᶠᵃᶜ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ @inline Δyᶜᵃᶠ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ -## ZRegularRG - -@inline Δzᵃᵃᶠ(i, j, k, grid::RGZ) = grid.Δzᵃᵃᶠ -@inline Δzᵃᵃᶜ(i, j, k, grid::RGZ) = grid.Δzᵃᵃᶜ - -@inline Δzᶜᵃᶜ(i, j, k, grid::RGZ) = grid.Δzᵃᵃᶜ -@inline Δzᶠᵃᶜ(i, j, k, grid::RGZ) = grid.Δzᵃᵃᶜ -@inline Δzᶜᵃᶠ(i, j, k, grid::RGZ) = grid.Δzᵃᵃᶠ - ##### ##### LatitudeLongitudeGrid ##### @@ -128,9 +137,6 @@ end @inline Δyᶜᶜᵃ(i, j, k, grid::LLG) = Δyᶠᶜᵃ(i, j, k, grid) @inline Δyᶠᶠᵃ(i, j, k, grid::LLG) = Δyᶜᶠᵃ(i, j, k, grid) -@inline Δzᵃᵃᶠ(i, j, k, grid::LLG) = @inbounds grid.Δzᵃᵃᶠ[k] -@inline Δzᵃᵃᶜ(i, j, k, grid::LLG) = @inbounds grid.Δzᵃᵃᶜ[k] - ### XRegularLLG with pre-computed metrics @inline Δxᶠᶜᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶠᶜᵃ[j] @@ -143,11 +149,6 @@ end @inline Δyᶜᶠᵃ(i, j, k, grid::LLGY) = grid.Δyᶜᶠᵃ @inline Δyᶠᶜᵃ(i, j, k, grid::LLGY) = grid.Δyᶠᶜᵃ -### ZRegularLLG with pre-computed metrics - -@inline Δzᵃᵃᶠ(i, j, k, grid::LLGZ) = grid.Δzᵃᵃᶠ -@inline Δzᵃᵃᶜ(i, j, k, grid::LLGZ) = grid.Δzᵃᵃᶜ - ## On the fly metrics @inline Δxᶠᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j]) @@ -184,12 +185,6 @@ end @inline Δyᶜᶠᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δyᶜᶠᵃ[i, j] @inline Δyᶠᶠᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δyᶠᶠᵃ[i, j] -@inline Δzᵃᵃᶜ(i, j, k, grid::OSSG) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᵃᵃᶠ(i, j, k, grid::OSSG) = @inbounds grid.Δzᵃᵃᶠ[k] - -@inline Δzᵃᵃᶜ(i, j, k, grid::OSSGZ) = grid.Δzᵃᵃᶜ -@inline Δzᵃᵃᶠ(i, j, k, grid::OSSGZ) = grid.Δzᵃᵃᶠ - ##### ##### ##### Areas!! @@ -221,7 +216,6 @@ for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) end end - #### #### Special 2D z Areas for LatitudeLongitudeGrid and OrthogonalSphericalShellGrid #### @@ -308,5 +302,8 @@ end @inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ @inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ -@inline Δλ(i, j, k, grid::OSSG, ℓx, ℓy, ℓz) = δxᶠᵃᵃ(i, j, k, grid, λnode, flip(ℓx), ℓy, ℓz) -@inline Δφ(i, j, k, grid::OSSG, ℓx, ℓy, ℓz) = δxᶠᵃᵃ(i, j, k, grid, λnode, ℓx, flip(ℓy), ℓz) +@inline Δλ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), ℓy, ℓz) +@inline Δλ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Center(), ℓy, ℓz) + +@inline Δφ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Face(), ℓz) +@inline Δφ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Center(), ℓz) From 20c94d6d451579d0f458681cdb2b19827c19591a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 15:54:40 +0100 Subject: [PATCH 12/47] a bit of cleanup in the operators --- .../spacings_and_areas_and_volumes.jl | 166 +++++++++--------- 1 file changed, 83 insertions(+), 83 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 225787ea71..50e2338e14 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -35,53 +35,60 @@ The operators in this file fall into three categories: ##### ##### -spacing_x_loc(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -spacing_y_loc(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -spacing_z_loc(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) - -spacing_1D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, :ᵃ, :ᵃ) -spacing_2D(::Val{:x}, LX, LY, LZ) = Symbol(:Δx, LX, LY, :ᵃ) -spacing_1D(::Val{:y}, LX, LY, LZ) = Symbol(:Δy, :ᵃ, LY, :ᵃ) -spacing_2D(::Val{:y}, LX, LY, LZ) = Symbol(:Δy, LX, LY, :ᵃ) -spacing_1D(::Val{:z}, LX, LY, LZ) = Symbol(:Δz, :ᵃ, :ᵃ, LZ) -spacing_2D(::Val{:z}, LX, LY, LZ) = Symbol(:Δz, :ᵃ, LY, LZ) +# This metaprogramming loop defines all possible combinations of locations for spacings +# All the spacings are reconducted to their one - dimensional counterparts. +# Grids that do not have a specific one - dimensional spacing for a given location need to +# extend these functions (for example, LatitudeLongitudeGrid). + +# Calling a non existing function (for example Δxᶜᵃᶜ on an OrthogonalSphericalShellGrid) will throw an error because +# the associated one - dimensional function is not defined. This is a feature, not a bug. + +spacing_1D(::Val{:x}, L1) = Symbol(:Δx, L1, :ᵃ, :ᵃ) +spacing_1D(::Val{:y}, L1) = Symbol(:Δy, :ᵃ, L1, :ᵃ) +spacing_1D(::Val{:z}, L1) = Symbol(:Δz, :ᵃ, :ᵃ, L1) +spacing_2D1(::Val{:x}, L1, L2) = Symbol(:Δx, L1, L2, :ᵃ) +spacing_2D1(::Val{:y}, L1, L2) = Symbol(:Δy, L2, L1, :ᵃ) +spacing_2D1(::Val{:z}, L1, L2) = Symbol(:Δz, :ᵃ, L2, L1) +spacing_2D2(::Val{:x}, L1, L2) = Symbol(:Δx, L1, :ᵃ, L2) +spacing_2D2(::Val{:y}, L1, L2) = Symbol(:Δy, :ᵃ, L1, L2) +spacing_2D2(::Val{:z}, L1, L2) = Symbol(:Δz, L2, :ᵃ, L1) +spacing_3D(::Val{:x}, L1, L2, L3) = Symbol(:Δx, L1, L2, L3) +spacing_3D(::Val{:y}, L1, L2, L3) = Symbol(:Δy, L2, L1, L3) +spacing_3D(::Val{:z}, L1, L2, L3) = Symbol(:Δz, L3, L2, L1) # Convenience Functions for all grids # This metaprogramming loop defines all the allowed combinations of Δx, Δy, and Δz # Note `:ᵃ` is not allowed for the location associated with the spacing for dir in (:x, :y, :z) - for LX in spacing_x_loc(dir), LY in spacing_y_loc(dir), LZ in spacing_z_loc(dir) - spacing1D = spacing_1D(Val(dir), LX, LY, LZ) - spacing2D = spacing_2D(Val(dir), LX, LY, LZ) - spacing3D = Symbol(:Δ, dir, LX, LY, LZ) - if spacing2D != spacing1D - @eval @inline $spacing2D(i, j, k, grid) = $spacing1D(i, j, k, grid) - end - if spacing3D != spacing2D - @eval @inline $spacing3D(i, j, k, grid) = $spacing2D(i, j, k, grid) + for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) + spacing1D = spacing_1D(Val(dir), L1) + spacing2D1 = spacing_2D(Val(dir), L1, L2) + spacing2D2 = spacing_2D(Val(dir), L1, L2) + @eval @inline $spacing2D1(i, j, k, grid) = $spacing1D(i, j, k, grid) + @eval @inline $spacing2D2(i, j, k, grid) = $spacing1D(i, j, k, grid) + + for L3 in (:ᶜ, :ᶠ) + spacing3D = spacing_3D(Val(dir), L1, L2, L3) + @eval @inline $spacing3D(i, j, k, grid) = $spacing2D1(i, j, k, grid) end end end ##### -##### Vertical spacings (same for all grids) +##### One - dimensional Vertical spacing (same for all grids) ##### @inline Δzᵃᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] @inline Δzᵃᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] -@inline Δzᶜᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᶠᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᶜᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] -@inline Δzᶠᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k] - @inline Δzᵃᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ @inline Δzᵃᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ -@inline Δzᶜᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ -@inline Δzᶠᵃᶜ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶜ -@inline Δzᶜᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ -@inline Δzᶠᵃᶠ(i, j, k, grid::ZRG) = grid.Δzᵃᵃᶠ +##### +##### +##### One - Dimensional Horizontal Spacings +##### +##### ##### ##### Rectilinear Grids (Flat grids already have Δ = 1) @@ -89,90 +96,79 @@ end @inline Δxᶠᵃᵃ(i, j, k, grid::RG) = @inbounds grid.Δxᶠᵃᵃ[i] @inline Δxᶜᵃᵃ(i, j, k, grid::RG) = @inbounds grid.Δxᶜᵃᵃ[i] -@inline Δxᶜᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δxᶜᵃᵃ[i] -@inline Δxᶠᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δxᶠᵃᵃ[i] -@inline Δxᶜᵃᶠ(i, j, k, grid::RG) = @inbounds grid.Δxᶜᵃᵃ[i] @inline Δyᵃᶠᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶠᵃ[j] @inline Δyᵃᶜᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j] -@inline Δyᶜᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j] -@inline Δyᶠᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j] -@inline Δyᶜᵃᶠ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j] - -@inline Δzᵃᵃᶠ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶠ[k] -@inline Δzᵃᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶜ[k] -@inline Δzᶜᵃᶜ(i, j, k, grid::RG) = @inbounds grid.Δzᵃᵃᶜ[k] ## XRegularRG @inline Δxᶠᵃᵃ(i, j, k, grid::RGX) = grid.Δxᶠᵃᵃ @inline Δxᶜᵃᵃ(i, j, k, grid::RGX) = grid.Δxᶜᵃᵃ -@inline Δxᶜᵃᶜ(i, j, k, grid::RGX) = grid.Δxᶜᵃᵃ -@inline Δxᶠᵃᶜ(i, j, k, grid::RGX) = grid.Δxᶠᵃᵃ -@inline Δxᶜᵃᶠ(i, j, k, grid::RGX) = grid.Δxᶜᵃᵃ ## YRegularRG @inline Δyᵃᶠᵃ(i, j, k, grid::RGY) = grid.Δyᵃᶠᵃ @inline Δyᵃᶜᵃ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ -@inline Δyᶜᵃᶜ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ -@inline Δyᶠᵃᶜ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ -@inline Δyᶜᵃᶠ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ +##### +##### LatitudeLongitude Grids (define both precomputed and non-precomputed metrics) +##### + +# Precomputed metrics + +@inline Δyᵃᶜᵃ(i, j, k, grid::LLGY) = grid.Δyᶠᶜᵃ +@inline Δyᵃᶠᵃ(i, j, k, grid::LLGY) = grid.Δyᶜᶠᵃ +@inline Δyᵃᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶠᶜᵃ[j] +@inline Δyᵃᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶜᶠᵃ[j] + +# On-the-fly metrics + +@inline Δyᵃᶠᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶠᵃ) +@inline Δyᵃᶜᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶜᵃ) +@inline Δyᵃᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶠᵃ[j]) +@inline Δyᵃᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶜᵃ[j]) ##### -##### LatitudeLongitudeGrid +##### +##### Two - Dimensional Horizontal Spacings +##### ##### -## Pre computed metrics +##### +##### LatitudeLongitudeGrid (only the Δx are required, Δy are 1D) +##### + +### Pre computed metrics @inline Δxᶜᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δxᶜᶠᵃ[i, j] @inline Δxᶠᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δxᶠᶜᵃ[i, j] @inline Δxᶠᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δxᶠᶠᵃ[i, j] @inline Δxᶜᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δxᶜᶜᵃ[i, j] -@inline Δyᶜᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶜᶠᵃ[j] -@inline Δyᶠᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶠᶜᵃ[j] -@inline Δyᶜᶜᵃ(i, j, k, grid::LLG) = Δyᶠᶜᵃ(i, j, k, grid) -@inline Δyᶠᶠᵃ(i, j, k, grid::LLG) = Δyᶜᶠᵃ(i, j, k, grid) - -### XRegularLLG with pre-computed metrics +### XRegularLLG with pre computed metrics @inline Δxᶠᶜᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶠᶜᵃ[j] @inline Δxᶜᶠᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶜᶠᵃ[j] @inline Δxᶠᶠᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶠᶠᵃ[j] @inline Δxᶜᶜᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶜᶜᵃ[j] -### YRegularLLG with pre-computed metrics - -@inline Δyᶜᶠᵃ(i, j, k, grid::LLGY) = grid.Δyᶜᶠᵃ -@inline Δyᶠᶜᵃ(i, j, k, grid::LLGY) = grid.Δyᶠᶜᵃ +### On-the-fly metrics -## On the fly metrics - -@inline Δxᶠᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j]) -@inline Δxᶜᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j]) -@inline Δxᶠᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j]) -@inline Δxᶜᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j]) - -@inline Δyᶜᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶠᵃ[j]) -@inline Δyᶠᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶜᵃ[j]) +@inline Δxᶠᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j]) +@inline Δxᶜᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j]) +@inline Δxᶠᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j]) +@inline Δxᶜᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j]) ### XRegularLLG with on-the-fly metrics -@inline Δxᶠᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ) * hack_cosd(grid.φᵃᶜᵃ[j]) -@inline Δxᶜᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ) * hack_cosd(grid.φᵃᶠᵃ[j]) -@inline Δxᶠᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ) * hack_cosd(grid.φᵃᶠᵃ[j]) -@inline Δxᶜᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ) * hack_cosd(grid.φᵃᶜᵃ[j]) - -### YRegularLLG with on-the-fly metrics - -@inline Δyᶜᶠᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶠᵃ) -@inline Δyᶠᶜᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶜᵃ) +@inline Δxᶠᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ) * hack_cosd(grid.φᵃᶜᵃ[j]) +@inline Δxᶜᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ) * hack_cosd(grid.φᵃᶠᵃ[j]) +@inline Δxᶠᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶠᵃᵃ) * hack_cosd(grid.φᵃᶠᵃ[j]) +@inline Δxᶜᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius * deg2rad(grid.Δλᶜᵃᵃ) * hack_cosd(grid.φᵃᶜᵃ[j]) ##### -##### OrthogonalSphericalShellGrid +##### OrthogonalSphericalShellGrid (does not have one-dimensional spacings) ##### @inline Δxᶜᶜᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δxᶜᶜᵃ[i, j] @@ -191,6 +187,9 @@ end ##### ##### +# We do the same thing as for the spacings: define general areas and then specialize for each grid. +# Areas need to be at least 2D + for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) x_spacing_2D = Symbol(:Δx, LX, LY, :ᵃ) @@ -224,10 +223,11 @@ end @inline Azᶜᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) @inline Azᶠᶠᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) @inline Azᶜᶜᵃ(i, j, k, grid::LLGF) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j])) -@inline Azᶠᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j])) -@inline Azᶜᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) -@inline Azᶠᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) -@inline Azᶜᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j])) + +@inline Azᶠᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j])) +@inline Azᶜᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) +@inline Azᶠᶠᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j]) - hack_sind(grid.φᵃᶜᵃ[j-1])) +@inline Azᶜᶜᵃ(i, j, k, grid::LLGFX) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j])) for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ) @@ -242,14 +242,14 @@ end ##### ##### -##### Volumes!! +##### Volumes!! (quite unambiguous) ##### ##### for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ) - volume = Symbol(:V, LX, LY, LZ) - z_area = Symbol(:Az, LX, LY, LZ) + volume = Symbol(:V, LX, LY, LZ) + z_area = Symbol(:Az, LX, LY, LZ) z_spacing = Symbol(:Δz, LX, LY, LZ) @eval begin From c495dd08eb2373825fa1f8c6fb9cb52d388c6ab2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 15:55:19 +0100 Subject: [PATCH 13/47] comment --- src/Operators/spacings_and_areas_and_volumes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 50e2338e14..7965581cde 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -35,8 +35,8 @@ The operators in this file fall into three categories: ##### ##### -# This metaprogramming loop defines all possible combinations of locations for spacings -# All the spacings are reconducted to their one - dimensional counterparts. +# This metaprogramming loop defines all possible combinations of locations for spacings. +# The general 2D and 3D spacings are reconducted to their one - dimensional counterparts. # Grids that do not have a specific one - dimensional spacing for a given location need to # extend these functions (for example, LatitudeLongitudeGrid). From 5f04304d3aaad85ef427ca16e172aac26f46781c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:01:19 +0100 Subject: [PATCH 14/47] fixes --- src/Operators/spacings_and_areas_and_volumes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 7965581cde..bddac04dd4 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -62,8 +62,8 @@ spacing_3D(::Val{:z}, L1, L2, L3) = Symbol(:Δz, L3, L2, L1) for dir in (:x, :y, :z) for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) spacing1D = spacing_1D(Val(dir), L1) - spacing2D1 = spacing_2D(Val(dir), L1, L2) - spacing2D2 = spacing_2D(Val(dir), L1, L2) + spacing2D1 = spacing_2D1(Val(dir), L1, L2) + spacing2D2 = spacing_2D2(Val(dir), L1, L2) @eval @inline $spacing2D1(i, j, k, grid) = $spacing1D(i, j, k, grid) @eval @inline $spacing2D2(i, j, k, grid) = $spacing1D(i, j, k, grid) From 81db44b2db991eef7dbed06ad73edbd1db7d5e3b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:03:41 +0100 Subject: [PATCH 15/47] last one --- src/Operators/spacings_and_areas_and_volumes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index bddac04dd4..2170be63f5 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -299,8 +299,8 @@ end @inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ[j] @inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ[j] -@inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ -@inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ +@inline Δφ(i, j, k, grid::LLGY, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ +@inline Δφ(i, j, k, grid::LLGY, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ @inline Δλ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), ℓy, ℓz) @inline Δλ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Center(), ℓy, ℓz) From 3c37ee4d60cc2bdc38c52bd84c1f06ebaea43bd6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:09:24 +0100 Subject: [PATCH 16/47] simplification --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index da3adbc9c6..2fb885e42b 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -35,9 +35,6 @@ end coordinates(grid::IBG) = coordinates(grid.underlying_grid) -@inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid) -@inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid) - # Extend both 2D and 3D methods @inline intrinsic_vector(i, j, k, ibg::IBG, u, v) = intrinsic_vector(i, j, k, ibg.underlying_grid, u, v) @inline extrinsic_vector(i, j, k, ibg::IBG, u, v) = extrinsic_vector(i, j, k, ibg.underlying_grid, u, v) From ae6ddebde5ad8efa4acb89cc49d7f2dd26f7580a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:31:51 +0100 Subject: [PATCH 17/47] coorect import --- .../immersed_grid_metrics.jl | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 2fb885e42b..640fe6fd6c 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,25 +10,26 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -spacing_x_loc(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -spacing_y_loc(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -spacing_z_loc(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +x_superscript(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +y_superscript(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +z_superscript(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) for dir in (:x, :y, :z) - for LX in spacing_x_loc(dir), LY in spacing_y_loc(dir), LZ in spacing_z_loc(dir) + for LX in x_superscript(dir), LY in y_superscript(dir), LZ in z_superscript(dir) spacing = Symbol(:Δ, dir, LX, LY, LZ) @eval begin import Oceananigans.Operators: $spacing - @inline $spacing(i, j, k, ibg::IBG) = $spacing(i, j, k, ibg.underlying_grid) + $spacing(i, j, k, ibg::IBG) = $spacing(i, j, k, ibg.underlying_grid) end end - for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) - area = Symbol(:A, dir, LX, LY, LZ) - volume = Symbol(:V, LX, LY, LZ) +end + +for L1 in (:ᶜ, :ᶠ) + for L2 in (:ᶜ, :ᶠ) + zarea2D = Symbol(:Az, L1, L2, :ᵃ) @eval begin - import Oceananigans.Operators: $area, $volume - @inline $area(i, j, k, ibg::IBG) = $area(i, j, k, ibg.underlying_grid) - @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid) + import Oceananigans.Operators: $zarea2D + $zarea2D(i, j, k, ibg::IBG) = $zarea2D(i, j, k, ibg.underlying_grid) end end end From ed5ef2720286cb716d6e420a289ed2c193692a18 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:56:53 +0100 Subject: [PATCH 18/47] correct areas --- .../spacings_and_areas_and_volumes.jl | 52 +++++++++++-------- 1 file changed, 31 insertions(+), 21 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 2170be63f5..1d15c1679e 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -188,33 +188,43 @@ end ##### # We do the same thing as for the spacings: define general areas and then specialize for each grid. -# Areas need to be at least 2D +# Areas need to be at least 2D so we use the respective 2D spacings to define them. +for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) + + Δxˡᵃ = Symbol(:Δx, L1, L2, :ᵃ) + Δxᵃˡ = Symbol(:Δx, L1, :ᵃ, L2) + Δyˡᵃ = Symbol(:Δy, L1, L2, :ᵃ) + Δyᵃˡ = Symbol(:Δy, :ᵃ, L1, L2) + Δzˡᵃ = Symbol(:Δz, L1, :ᵃ, L2) + Δzᵃˡ = Symbol(:Δz, :ᵃ, L1, L2) + + # 2D areas + Axᵃˡˡ = Symbol(:Ax, :ᵃ, L1, L2) + Ayˡᵃˡ = Symbol(:Ay, L1, :ᵃ, L2) + Azˡˡᵃ = Symbol(:Az, L1, L2, :ᵃ) -for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ) - - x_spacing_2D = Symbol(:Δx, LX, LY, :ᵃ) - y_spacing_2D = Symbol(:Δy, LX, LY, :ᵃ) - z_area_2D = Symbol(:Az, LX, LY, :ᵃ) - - @eval $z_area_2D(i, j, k, grid) = $x_spacing_2D(i, j, k, grid) * $y_spacing_2D(i, j, k, grid) - - for LZ in (:ᶜ, :ᶠ) - x_spacing_3D = Symbol(:Δx, LX, LY, LZ) - y_spacing_3D = Symbol(:Δy, LX, LY, LZ) - z_spacing_3D = Symbol(:Δz, LX, LY, LZ) - - x_area_3D = Symbol(:Ax, LX, LY, LZ) - y_area_3D = Symbol(:Ay, LX, LY, LZ) - z_area_3D = Symbol(:Az, LX, LY, LZ) + @eval begin + @inline $Axᵃˡˡ(i, j, k, grid) = $Δyᵃˡ(i, j, k, grid) * $Δzᵃˡ(i, j, k, grid) + @inline $Ayˡᵃˡ(i, j, k, grid) = $Δxᵃˡ(i, j, k, grid) * $Δzˡᵃ(i, j, k, grid) + @inline $Azˡˡᵃ(i, j, k, grid) = $Δxˡᵃ(i, j, k, grid) * $Δyˡᵃ(i, j, k, grid) + end + for L3 in (:ᶜ, :ᶠ) + # 3D areas + Axˡˡˡ = Symbol(:Ax, L3, L1, L2) + Ayˡˡˡ = Symbol(:Ay, L1, L3, L2) + Azˡˡˡ = Symbol(:Az, L1, L2, L3) + @eval begin - @inline $x_area_3D(i, j, k, grid) = $y_spacing_3D(i, j, k, grid) * $z_spacing_3D(i, j, k, grid) - @inline $y_area_3D(i, j, k, grid) = $x_spacing_3D(i, j, k, grid) * $z_spacing_3D(i, j, k, grid) - @inline $z_area_3D(i, j, k, grid) = $z_area_2D(i, j, k, grid) + # For the moment 3D areas equal their 2D counterpart. This might change if + # we want to implement deep atmospheres where Az is a function of z + @inline $Axˡˡˡ(i, j, k, grid) = $Axᵃˡˡ(i, j, k, grid) + @inline $Ayˡˡˡ(i, j, k, grid) = $Axᵃˡˡ(i, j, k, grid) + @inline $Azˡˡˡ(i, j, k, grid) = $Azˡˡᵃ(i, j, k, grid) end end end - + #### #### Special 2D z Areas for LatitudeLongitudeGrid and OrthogonalSphericalShellGrid #### From 57890ec3320a2b0e89adf49b091125608ee2735b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 16:58:16 +0100 Subject: [PATCH 19/47] better naming --- .../spacings_and_areas_and_volumes.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 1d15c1679e..fbb59b4bfc 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -191,12 +191,12 @@ end # Areas need to be at least 2D so we use the respective 2D spacings to define them. for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) - Δxˡᵃ = Symbol(:Δx, L1, L2, :ᵃ) - Δxᵃˡ = Symbol(:Δx, L1, :ᵃ, L2) - Δyˡᵃ = Symbol(:Δy, L1, L2, :ᵃ) - Δyᵃˡ = Symbol(:Δy, :ᵃ, L1, L2) - Δzˡᵃ = Symbol(:Δz, L1, :ᵃ, L2) - Δzᵃˡ = Symbol(:Δz, :ᵃ, L1, L2) + Δxˡˡᵃ = Symbol(:Δx, L1, L2, :ᵃ) + Δxˡᵃˡ = Symbol(:Δx, L1, :ᵃ, L2) + Δyˡˡᵃ = Symbol(:Δy, L1, L2, :ᵃ) + Δyᵃˡˡ = Symbol(:Δy, :ᵃ, L1, L2) + Δzˡᵃˡ = Symbol(:Δz, L1, :ᵃ, L2) + Δzᵃˡˡ = Symbol(:Δz, :ᵃ, L1, L2) # 2D areas Axᵃˡˡ = Symbol(:Ax, :ᵃ, L1, L2) @@ -204,9 +204,9 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) Azˡˡᵃ = Symbol(:Az, L1, L2, :ᵃ) @eval begin - @inline $Axᵃˡˡ(i, j, k, grid) = $Δyᵃˡ(i, j, k, grid) * $Δzᵃˡ(i, j, k, grid) - @inline $Ayˡᵃˡ(i, j, k, grid) = $Δxᵃˡ(i, j, k, grid) * $Δzˡᵃ(i, j, k, grid) - @inline $Azˡˡᵃ(i, j, k, grid) = $Δxˡᵃ(i, j, k, grid) * $Δyˡᵃ(i, j, k, grid) + @inline $Axᵃˡˡ(i, j, k, grid) = $Δyᵃˡˡ(i, j, k, grid) * $Δzᵃˡˡ(i, j, k, grid) + @inline $Ayˡᵃˡ(i, j, k, grid) = $Δxˡᵃˡ(i, j, k, grid) * $Δzˡᵃˡ(i, j, k, grid) + @inline $Azˡˡᵃ(i, j, k, grid) = $Δxˡˡᵃ(i, j, k, grid) * $Δyˡˡᵃ(i, j, k, grid) end for L3 in (:ᶜ, :ᶠ) From 64235b9d09d2256be4a3c11143b546544ca4ef19 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 17:00:51 +0100 Subject: [PATCH 20/47] correct grid metric --- src/AbstractOperations/grid_metrics.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index db5d6625c9..75b9fba620 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -13,7 +13,7 @@ const AbstractGridMetric = Union{typeof(Δx), typeof(Ax), typeof(Ay), typeof(Az), - typeof(volume)} + typeof(volume)} # Do we want it to be `volume` or just `V` like in the Operators module? """ metric_function(loc, metric::AbstractGridMetric) @@ -22,7 +22,11 @@ Return the function associated with `metric::AbstractGridMetric` at `loc`ation. """ function metric_function(loc, metric) code = Tuple(interpolation_code(ℓ) for ℓ in loc) - metric_function_symbol = Symbol(metric, 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 From 5d030f1e3595691f8f29aa3057fb6578ca7b1d3e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 17:01:19 +0100 Subject: [PATCH 21/47] space --- src/Operators/spacings_and_areas_and_volumes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index fbb59b4bfc..2547a729e4 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -258,7 +258,7 @@ end for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ) - volume = Symbol(:V, LX, LY, LZ) + volume = Symbol(:V, LX, LY, LZ) z_area = Symbol(:Az, LX, LY, LZ) z_spacing = Symbol(:Δz, LX, LY, LZ) From 33ad6005ecc477d7f1c671e127f6c3cbce943c01 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 17:32:10 +0100 Subject: [PATCH 22/47] bugfix --- src/Operators/spacings_and_areas_and_volumes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 2547a729e4..50b271d64b 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -219,7 +219,7 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) # For the moment 3D areas equal their 2D counterpart. This might change if # we want to implement deep atmospheres where Az is a function of z @inline $Axˡˡˡ(i, j, k, grid) = $Axᵃˡˡ(i, j, k, grid) - @inline $Ayˡˡˡ(i, j, k, grid) = $Axᵃˡˡ(i, j, k, grid) + @inline $Ayˡˡˡ(i, j, k, grid) = $Ayˡᵃˡ(i, j, k, grid) @inline $Azˡˡˡ(i, j, k, grid) = $Azˡˡᵃ(i, j, k, grid) end end From 5941eb6128fde610f637bee6721c73bfb7ab5990 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 17:41:28 +0100 Subject: [PATCH 23/47] correct 3D vertical areas --- src/Operators/spacings_and_areas_and_volumes.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 50b271d64b..ae6f2f7d12 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -210,16 +210,22 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) end for L3 in (:ᶜ, :ᶠ) + # 3D spacings + Δxˡˡˡ = Symbol(:Δx, L1, L2, L3) + Δyˡˡˡ = Symbol(:Δy, L1, L2, L3) + Δzˡˡˡ = Symbol(:Δz, L1, L2, L3) + # 3D areas - Axˡˡˡ = Symbol(:Ax, L3, L1, L2) - Ayˡˡˡ = Symbol(:Ay, L1, L3, L2) + Axˡˡˡ = Symbol(:Ax, L1, L2, L3) + Ayˡˡˡ = Symbol(:Ay, L1, L2, L3) Azˡˡˡ = Symbol(:Az, L1, L2, L3) @eval begin - # For the moment 3D areas equal their 2D counterpart. This might change if + @inline $Axˡˡˡ(i, j, k, grid) = $Δyˡˡˡ(i, j, k, grid) * $Δzˡˡˡ(i, j, k, grid) + @inline $Ayˡˡˡ(i, j, k, grid) = $Axˡˡˡ(i, j, k, grid) * $Δzˡˡˡ(i, j, k, grid) + + # For the moment the horizontal area is independent of `z`. This might change if # we want to implement deep atmospheres where Az is a function of z - @inline $Axˡˡˡ(i, j, k, grid) = $Axᵃˡˡ(i, j, k, grid) - @inline $Ayˡˡˡ(i, j, k, grid) = $Ayˡᵃˡ(i, j, k, grid) @inline $Azˡˡˡ(i, j, k, grid) = $Azˡˡᵃ(i, j, k, grid) end end From 35a85045e3832011bab9fe499f6e8a2b1c1096fe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 18:06:43 +0100 Subject: [PATCH 24/47] another bugfix --- src/ImmersedBoundaries/partial_cell_bottom.jl | 2 +- src/Operators/spacings_and_areas_and_volumes.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 214ab40619..36cc32487b 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -201,4 +201,4 @@ YFlatPCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Flat, <:Any, <:Any, <:Partial @inline Δzᶜᶠᶠ(i, j, k, ibg::YFlatPCBIBG) = Δzᶜᶜᶠ(i, j, k, ibg) @inline Δzᶠᶠᶜ(i, j, k, ibg::XFlatPCBIBG) = Δzᶜᶠᶜ(i, j, k, ibg) -@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) \ No newline at end of file +@inline Δzᶠᶠᶜ(i, j, k, ibg::YFlatPCBIBG) = Δzᶠᶜᶜ(i, j, k, ibg) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index ae6f2f7d12..9c8fc0ca01 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -222,7 +222,7 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) @eval begin @inline $Axˡˡˡ(i, j, k, grid) = $Δyˡˡˡ(i, j, k, grid) * $Δzˡˡˡ(i, j, k, grid) - @inline $Ayˡˡˡ(i, j, k, grid) = $Axˡˡˡ(i, j, k, grid) * $Δzˡˡˡ(i, j, k, grid) + @inline $Ayˡˡˡ(i, j, k, grid) = $Δxˡˡˡ(i, j, k, grid) * $Δzˡˡˡ(i, j, k, grid) # For the moment the horizontal area is independent of `z`. This might change if # we want to implement deep atmospheres where Az is a function of z From 75b5a138297c2716ab5791ee5d61b5bae46804e2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 18:20:48 +0100 Subject: [PATCH 25/47] remove all spacing functions --- src/Grids/latitude_longitude_grid.jl | 13 ------------- test/test_grids.jl | 27 ++++++++++++++------------- 2 files changed, 14 insertions(+), 26 deletions(-) diff --git a/src/Grids/latitude_longitude_grid.jl b/src/Grids/latitude_longitude_grid.jl index bd411ea5d1..a7d06621fc 100644 --- a/src/Grids/latitude_longitude_grid.jl +++ b/src/Grids/latitude_longitude_grid.jl @@ -673,19 +673,6 @@ end ##### Grid spacings ##### -@inline λspacing(i, grid::LLG, ::C) = @inbounds grid.Δλᶜᵃᵃ[i] -@inline λspacing(i, grid::LLG, ::F) = @inbounds grid.Δλᶠᵃᵃ[i] -@inline λspacing(i, grid::XRegularLLG, ::C) = grid.Δλᶜᵃᵃ -@inline λspacing(i, grid::XRegularLLG, ::F) = grid.Δλᶠᵃᵃ - -@inline φspacing(j, grid::LLG, ::C) = @inbounds grid.Δφᵃᶜᵃ[j] -@inline φspacing(j, grid::LLG, ::F) = @inbounds grid.Δφᵃᶠᵃ[j] -@inline φspacing(j, grid::YRegularLLG, ::C) = grid.Δφᵃᶜᵃ -@inline φspacing(j, grid::YRegularLLG, ::F) = grid.Δφᵃᶠᵃ - -@inline λspacing(i, j, k, grid::LLG, ℓx, ℓy, ℓz) = λspacing(i, grid, ℓx) -@inline φspacing(i, j, k, grid::LLG, ℓx, ℓy, ℓz) = φspacing(j, grid, ℓy) - @inline xspacings(grid::LLG, ℓx, ℓy) = xspacings(grid, ℓx, ℓy, nothing) @inline yspacings(grid::LLG, ℓx, ℓy) = yspacings(grid, ℓx, ℓy, nothing) @inline zspacings(grid::LLG, ℓz) = zspacings(grid, nothing, nothing, ℓz) diff --git a/test/test_grids.jl b/test/test_grids.jl index 2cf965895b..50a0a3f7f1 100644 --- a/test/test_grids.jl +++ b/test/test_grids.jl @@ -4,8 +4,9 @@ include("data_dependencies.jl") using Oceananigans.Grids: total_extent, xspacings, yspacings, zspacings, xnode, ynode, znode, λnode, φnode, - λspacing, φspacing, λspacings, φspacings + λspacings, φspacings +using Oceananigans.Operators using Oceananigans.Operators: Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶠᶠᵃ, Δxᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ ##### @@ -216,9 +217,9 @@ function test_regular_rectilinear_xnode_ynode_znode_and_spacings(arch, FT) @test all(yspacings(grid, Face()) .== yspacings(grid, Center(), Face(), Center())) @test all(zspacings(grid, Face()) .== zspacings(grid, Center(), Center(), Face())) - @test xspacing(1, 1, 1, grid, Face(), Center(), Center()) ≈ FT(π/N) - @test yspacing(1, 1, 1, grid, Center(), Face(), Center()) ≈ FT(π/N) - @test zspacing(1, 1, 1, grid, Center(), Center(), Face()) ≈ FT(π/N) + @test Δx(1, 1, 1, grid, Face(), Center(), Center()) ≈ FT(π/N) + @test Δy(1, 1, 1, grid, Center(), Face(), Center()) ≈ FT(π/N) + @test Δz(1, 1, 1, grid, Center(), Center(), Face()) ≈ FT(π/N) end return nothing @@ -414,7 +415,7 @@ function test_rectilinear_grid_correct_spacings(FT, N) @test all(isapprox.(zspacings(grid, Face()), reshape(grid.Δzᵃᵃᶠ[1:N+1], 1, 1, N+1))) @test all(isapprox.(zspacings(grid, Center()), reshape(grid.Δzᵃᵃᶜ[1:N], 1, 1, N))) - @test zspacing(1, 1, 2, grid, Center(), Center(), Face()) == grid.Δzᵃᵃᶠ[2] + @test Δz(1, 1, 2, grid, Center(), Center(), Face()) == grid.Δzᵃᵃᶠ[2] @test minimum_zspacing(grid, Center(), Center(), Center()) ≈ minimum(grid.Δzᵃᵃᶜ[1:grid.Nz]) @@ -560,20 +561,20 @@ function test_basic_lat_lon_general_grid(FT) @test all(zspacings(grid_reg, Face(), Center(), Center()) .== zspacings(grid_reg, Center())) @test all(zspacings(grid_reg, Face(), Center(), Face() ) .== zspacings(grid_reg, Face())) - @test xspacing(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δxᶜᶜᵃ[2] - @test xspacing(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δxᶜᶠᵃ[2] - @test yspacing(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δyᶜᶠᵃ - @test yspacing(1, 2, 3, grid_reg, Face(), Center(), Center()) == grid_reg.Δyᶠᶜᵃ - @test zspacing(1, 2, 3, grid_reg, Center(), Center(), Face() ) == grid_reg.Δzᵃᵃᶠ - @test zspacing(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δzᵃᵃᶜ + @test Δx(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δxᶜᶜᵃ[2] + @test Δx(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δxᶜᶠᵃ[2] + @test Δy(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δyᶜᶠᵃ + @test Δy(1, 2, 3, grid_reg, Face(), Center(), Center()) == grid_reg.Δyᶠᶜᵃ + @test Δz(1, 2, 3, grid_reg, Center(), Center(), Face() ) == grid_reg.Δzᵃᵃᶠ + @test Δz(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δzᵃᵃᶜ @test all(λspacings(grid_reg, Center()) .== grid_reg.Δλᶜᵃᵃ) @test all(λspacings(grid_reg, Face()) .== grid_reg.Δλᶠᵃᵃ) @test all(φspacings(grid_reg, Center()) .== grid_reg.Δφᵃᶜᵃ) @test all(φspacings(grid_reg, Face()) .== grid_reg.Δφᵃᶠᵃ) - @test λspacing(1, 2, 3, grid_reg, Face(), Center(), Face()) == grid_reg.Δλᶠᵃᵃ - @test φspacing(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δφᵃᶠᵃ + @test Δλ(1, 2, 3, grid_reg, Face(), Center(), Face()) == grid_reg.Δλᶠᵃᵃ + @test Δφ(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δφᵃᶠᵃ Δλ = grid_reg.Δλᶠᵃᵃ λₛ = (-grid_reg.Lx/2):Δλ:(grid_reg.Lx/2) From cdaf207fac4ccbde1a76d39ade16a36bee5f61d5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 18:56:53 +0100 Subject: [PATCH 26/47] tests should pass now --- src/AbstractOperations/grid_metrics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 75b9fba620..32386aaaa6 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -219,7 +219,7 @@ KernelFunctionOperation at (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 From 2f924f8b44c49c02de7f4d7ae260318e40352c7f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 21:21:34 +0100 Subject: [PATCH 27/47] correct bottom height --- src/ImmersedBoundaries/partial_cell_bottom.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 36cc32487b..93be418142 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -80,8 +80,6 @@ end domain_bottom = znode(i, j, 1, grid, c, c, f) domain_top = znode(i, j, grid.Nz+1, grid, c, c, f) @inbounds bottom_field[i, j, 1] = clamp(zb, domain_bottom, domain_top) - adjusted_zb = bottom_field[i, j, 1] - ϵ = ib.minimum_fractional_cell_height for k in 1:grid.Nz @@ -93,10 +91,8 @@ end # If the size of the bottom cell is less than ϵ Δz, # we enforce a minimum size of ϵ Δz. - adjusted_zb = ifelse(bottom_cell, capped_zb, zb) + @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, capped_zb, bottom_field[i, j, 1]) end - - @inbounds bottom_field[i, j, 1] = adjusted_zb end function on_architecture(arch, ib::PartialCellBottom{<:Field}) @@ -159,7 +155,7 @@ end # Get bottom z-coordinate and fractional Δz parameter zb = @inbounds ib.bottom_height[i, j, 1] - + # Are we in a bottom cell? at_the_bottom = bottom_cell(i, j, k, ibg) From b7a797196ca391a62fb707fe0b6588abf8c4558d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 25 Nov 2024 23:15:13 +0100 Subject: [PATCH 28/47] fix tests --- .../spacings_and_areas_and_volumes.jl | 12 +++++----- test/test_grids.jl | 24 +++++++++---------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 9c8fc0ca01..f643227f89 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -313,13 +313,13 @@ end @inline Δλ(i, j, k, grid::LLGX, ::Center, ℓy, ℓz) = @inbounds grid.Δλᶜᵃᵃ @inline Δλ(i, j, k, grid::LLGX, ::Face, ℓy, ℓz) = @inbounds grid.Δλᶠᵃᵃ -@inline Δφ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ[j] -@inline Δφ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ[j] -@inline Δφ(i, j, k, grid::LLGY, ::Center, ℓy, ℓz) = @inbounds grid.Δφᵃᶜᵃ -@inline Δφ(i, j, k, grid::LLGY, ::Face, ℓy, ℓz) = @inbounds grid.Δφᵃᶠᵃ +@inline Δφ(i, j, k, grid::LLG, ℓx, ::Center, ℓz) = @inbounds grid.Δφᵃᶜᵃ[j] +@inline Δφ(i, j, k, grid::LLG, ℓx, ::Face, ℓz) = @inbounds grid.Δφᵃᶠᵃ[j] +@inline Δφ(i, j, k, grid::LLGY, ℓx, ::Center, ℓz) = @inbounds grid.Δφᵃᶜᵃ +@inline Δφ(i, j, k, grid::LLGY, ℓx, ::Face, ℓz) = @inbounds grid.Δφᵃᶠᵃ @inline Δλ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), ℓy, ℓz) @inline Δλ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Center(), ℓy, ℓz) -@inline Δφ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Face(), ℓz) -@inline Δφ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Center(), ℓz) +@inline Δφ(i, j, k, grid::OSSG, ℓy, ::Center, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Face(), ℓz) +@inline Δφ(i, j, k, grid::OSSG, ℓy, ::Face, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Center(), ℓz) diff --git a/test/test_grids.jl b/test/test_grids.jl index 50a0a3f7f1..7a39d3b64b 100644 --- a/test/test_grids.jl +++ b/test/test_grids.jl @@ -6,7 +6,7 @@ using Oceananigans.Grids: total_extent, xnode, ynode, znode, λnode, φnode, λspacings, φspacings -using Oceananigans.Operators +using Oceananigans.Operators: Δx, Δy, Δz, Δλ, Δφ, Ax, Ay, Az, volume using Oceananigans.Operators: Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶠᶠᵃ, Δxᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ ##### @@ -561,20 +561,20 @@ function test_basic_lat_lon_general_grid(FT) @test all(zspacings(grid_reg, Face(), Center(), Center()) .== zspacings(grid_reg, Center())) @test all(zspacings(grid_reg, Face(), Center(), Face() ) .== zspacings(grid_reg, Face())) - @test Δx(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δxᶜᶜᵃ[2] - @test Δx(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δxᶜᶠᵃ[2] - @test Δy(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δyᶜᶠᵃ - @test Δy(1, 2, 3, grid_reg, Face(), Center(), Center()) == grid_reg.Δyᶠᶜᵃ - @test Δz(1, 2, 3, grid_reg, Center(), Center(), Face() ) == grid_reg.Δzᵃᵃᶠ - @test Δz(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δzᵃᵃᶜ + @test Operators.Δx(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δxᶜᶜᵃ[2] + @test Operators.Δx(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δxᶜᶠᵃ[2] + @test Operators.Δy(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δyᶜᶠᵃ + @test Operators.Δy(1, 2, 3, grid_reg, Face(), Center(), Center()) == grid_reg.Δyᶠᶜᵃ + @test Operators.Δz(1, 2, 3, grid_reg, Center(), Center(), Face() ) == grid_reg.Δzᵃᵃᶠ + @test Operators.Δz(1, 2, 3, grid_reg, Center(), Center(), Center()) == grid_reg.Δzᵃᵃᶜ @test all(λspacings(grid_reg, Center()) .== grid_reg.Δλᶜᵃᵃ) @test all(λspacings(grid_reg, Face()) .== grid_reg.Δλᶠᵃᵃ) @test all(φspacings(grid_reg, Center()) .== grid_reg.Δφᵃᶜᵃ) @test all(φspacings(grid_reg, Face()) .== grid_reg.Δφᵃᶠᵃ) - @test Δλ(1, 2, 3, grid_reg, Face(), Center(), Face()) == grid_reg.Δλᶠᵃᵃ - @test Δφ(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δφᵃᶠᵃ + @test Operators.Δλ(1, 2, 3, grid_reg, Face(), Center(), Face()) == grid_reg.Δλᶠᵃᵃ + @test Operators.Δφ(1, 2, 3, grid_reg, Center(), Face(), Center()) == grid_reg.Δφᵃᶠᵃ Δλ = grid_reg.Δλᶠᵃᵃ λₛ = (-grid_reg.Lx/2):Δλ:(grid_reg.Lx/2) @@ -661,10 +661,10 @@ function test_lat_lon_xyzλφ_node_nodes(FT, arch) @test minimum_yspacing(grid) / grid.radius ≈ FT(π/6) @test minimum_zspacing(grid) ≈ 5 - grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> y < 20 && y > -20)) + grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> y < 20 && y > -20 ? -50 : -0)) - @test minimum_xspacing(grid, Face(), Face(), Face()) / grid.radius ≈ π/6 * cosd(15) - @test minimum_xspacing(grid) / grid.radius ≈ π/6 * cosd(15) + @test minimum_xspacing(grid, Face(), Face(), Face()) / grid.radius ≈ FT(π/6) * cosd(30) + @test minimum_xspacing(grid) / grid.radius ≈ FT(π/6) * cosd(15) return nothing end From 9d2df2ff80ffc77a7e15b92844b88e883eccc857 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 09:45:29 +0100 Subject: [PATCH 29/47] fix example --- examples/internal_tide.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 08d2ef4db3..0e39bb4f49 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -46,7 +46,7 @@ width = 20kilometers hill(x) = h₀ * exp(-x^2 / 2width^2) bottom(x) = - H + hill(x) -grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom)) +grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom, minimum_fractional_cell_height=1.0)) # Let's see how the domain with the bathymetry is. @@ -148,9 +148,10 @@ wall_clock = Ref(time_ns()) function progress(sim) elapsed = 1e-9 * (time_ns() - wall_clock[]) - msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹\n", + cfl = Δt / Oceananigans.Advection.cell_advection_timescale(sim.model.grid, sim.model.velocities) + msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹, CFL %6.3e\n", iteration(sim), prettytime(sim), prettytime(elapsed), - maximum(abs, sim.model.velocities.w)) + maximum(abs, sim.model.velocities.w), cfl) wall_clock[] = time_ns() From 70f0bbf4b278d4ffb651a1bd92abfd7f030a3bc6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 08:50:57 +0100 Subject: [PATCH 30/47] export the verbose versions --- src/AbstractOperations/grid_metrics.jl | 1 + src/Oceananigans.jl | 1 + src/Operators/Operators.jl | 15 ++++++++++++--- 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 32386aaaa6..d58957137c 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -2,6 +2,7 @@ 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 diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index 39701328c5..52d9df483d 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -13,6 +13,7 @@ export Periodic, Bounded, Flat, RectilinearGrid, LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, nodes, xnodes, ynodes, znodes, λnodes, φnodes, + xspacing, yspacing, zspacing, xarea, yarea, zarea, volume, xspacings, yspacings, zspacings, λspacings, φspacings, minimum_xspacing, minimum_yspacing, minimum_zspacing, diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index 15a3a1ad8d..d9a7941622 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -1,8 +1,5 @@ module Operators -# General metric operators -export Δx, Δy, Δz, Δλ, Δφ, Ax, Ay, Az, volume - # Spacings export Δxᶠᶠᶠ, Δxᶠᶠᶜ, Δxᶠᶜᶠ, Δxᶠᶜᶜ, Δxᶜᶠᶠ, Δxᶜᶠᶜ, Δxᶜᶜᶠ, Δxᶜᶜᶜ export Δyᶠᶠᶠ, Δyᶠᶠᶜ, Δyᶠᶜᶠ, Δyᶠᶜᶜ, Δyᶜᶠᶠ, Δyᶜᶠᶜ, Δyᶜᶜᶠ, Δyᶜᶜᶜ @@ -105,6 +102,9 @@ const LLGF = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing} const LLGFX = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Number} const LLGFY = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Any, <:Number} +# General metric operators +export xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume + include("difference_operators.jl") include("interpolation_operators.jl") include("interpolation_utils.jl") @@ -120,4 +120,13 @@ include("laplacian_operators.jl") include("vector_rotation_operators.jl") +@inline xspacing(args...) = Δx(args...) +@inline yspacing(args...) = Δy(args...) +@inline zspacing(args...) = Δz(args...) +@inline λspacing(abs...) = Δλ(abs...) +@inline φspacing(abs...) = Δφ(abs...) +@inline xarea(args...) = Ax(args...) +@inline yarea(args...) = Ay(args...) +@inline zarea(args...) = Az(args...) + end # module From da22227fb3e899f16e0b56d05e6ccf168ecf28e7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 09:00:23 +0100 Subject: [PATCH 31/47] revert partial cell bottom --- examples/internal_tide.jl | 7 ++--- .../immersed_grid_metrics.jl | 28 ++++++++----------- src/ImmersedBoundaries/partial_cell_bottom.jl | 8 ++++-- src/Oceananigans.jl | 4 ++- 4 files changed, 24 insertions(+), 23 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 0e39bb4f49..08d2ef4db3 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -46,7 +46,7 @@ width = 20kilometers hill(x) = h₀ * exp(-x^2 / 2width^2) bottom(x) = - H + hill(x) -grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom, minimum_fractional_cell_height=1.0)) +grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom)) # Let's see how the domain with the bathymetry is. @@ -148,10 +148,9 @@ wall_clock = Ref(time_ns()) function progress(sim) elapsed = 1e-9 * (time_ns() - wall_clock[]) - cfl = Δt / Oceananigans.Advection.cell_advection_timescale(sim.model.grid, sim.model.velocities) - msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹, CFL %6.3e\n", + msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹\n", iteration(sim), prettytime(sim), prettytime(elapsed), - maximum(abs, sim.model.velocities.w), cfl) + maximum(abs, sim.model.velocities.w)) wall_clock[] = time_ns() diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 640fe6fd6c..48cb239768 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,32 +10,28 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -x_superscript(dir) = dir == :x ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -y_superscript(dir) = dir == :y ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) -z_superscript(dir) = dir == :z ? (:ᶜ, :ᶠ) : (:ᶜ, :ᶠ, :ᵃ) +for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) + for dir in (:x, :y, :z), operator in (:Δ, :A) -for dir in (:x, :y, :z) - for LX in x_superscript(dir), LY in y_superscript(dir), LZ in z_superscript(dir) - spacing = Symbol(:Δ, dir, LX, LY, LZ) + metric = Symbol(operator, dir, LX, LY, LZ) @eval begin - import Oceananigans.Operators: $spacing - $spacing(i, j, k, ibg::IBG) = $spacing(i, j, k, ibg.underlying_grid) + import Oceananigans.Operators: $metric + @inline $metric(i, j, k, ibg::IBG) = $metric(i, j, k, ibg.underlying_grid) end end -end -for L1 in (:ᶜ, :ᶠ) - for L2 in (:ᶜ, :ᶠ) - zarea2D = Symbol(:Az, L1, L2, :ᵃ) - @eval begin - import Oceananigans.Operators: $zarea2D - $zarea2D(i, j, k, ibg::IBG) = $zarea2D(i, j, k, ibg.underlying_grid) - end + volume = Symbol(:V, LX, LY, LZ) + @eval begin + import Oceananigans.Operators: $volume + @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid) end end coordinates(grid::IBG) = coordinates(grid.underlying_grid) +@inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid) +@inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid) + # Extend both 2D and 3D methods @inline intrinsic_vector(i, j, k, ibg::IBG, u, v) = intrinsic_vector(i, j, k, ibg.underlying_grid, u, v) @inline extrinsic_vector(i, j, k, ibg::IBG, u, v) = extrinsic_vector(i, j, k, ibg.underlying_grid, u, v) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 93be418142..36cc32487b 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -80,6 +80,8 @@ end domain_bottom = znode(i, j, 1, grid, c, c, f) domain_top = znode(i, j, grid.Nz+1, grid, c, c, f) @inbounds bottom_field[i, j, 1] = clamp(zb, domain_bottom, domain_top) + adjusted_zb = bottom_field[i, j, 1] + ϵ = ib.minimum_fractional_cell_height for k in 1:grid.Nz @@ -91,8 +93,10 @@ end # If the size of the bottom cell is less than ϵ Δz, # we enforce a minimum size of ϵ Δz. - @inbounds bottom_field[i, j, 1] = ifelse(bottom_cell, capped_zb, bottom_field[i, j, 1]) + adjusted_zb = ifelse(bottom_cell, capped_zb, zb) end + + @inbounds bottom_field[i, j, 1] = adjusted_zb end function on_architecture(arch, ib::PartialCellBottom{<:Field}) @@ -155,7 +159,7 @@ end # Get bottom z-coordinate and fractional Δz parameter zb = @inbounds ib.bottom_height[i, j, 1] - + # Are we in a bottom cell? at_the_bottom = bottom_cell(i, j, k, ibg) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index 52d9df483d..d49accb100 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -13,10 +13,12 @@ export Periodic, Bounded, Flat, RectilinearGrid, LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, nodes, xnodes, ynodes, znodes, λnodes, φnodes, - xspacing, yspacing, zspacing, xarea, yarea, zarea, volume, xspacings, yspacings, zspacings, λspacings, φspacings, minimum_xspacing, minimum_yspacing, minimum_zspacing, + # Pointwise spacing, area, and volume operators + xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume, + # Immersed boundaries ImmersedBoundaryGrid, GridFittedBoundary, GridFittedBottom, PartialCellBottom, From 8cb1fc1ddf8313964ddaea656418507514efadbe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 17:37:01 +0100 Subject: [PATCH 32/47] bugfix --- src/AbstractOperations/grid_metrics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index d58957137c..a1efb0054d 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -2,7 +2,7 @@ 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 +using Oceananigans.Operators: Δx, Δy, Δz, Ax, Δλ, Δφ, Ay, Az, volume import Oceananigans.Grids: xspacings, yspacings, zspacings, λspacings, φspacings From 6d7f4352f401dbdc02f49c7cc8b3c42eb297e673 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 08:27:33 +0100 Subject: [PATCH 33/47] fix tests --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 48cb239768..fd968a0c3a 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,7 +10,7 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) +for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) for dir in (:x, :y, :z), operator in (:Δ, :A) metric = Symbol(operator, dir, LX, LY, LZ) From 653b672fed62bdc8ac14ecd81d6784925cfaf50e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 29 Nov 2024 10:37:15 +0100 Subject: [PATCH 34/47] Update spacings_and_areas_and_volumes.jl --- src/Operators/spacings_and_areas_and_volumes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index f643227f89..d1063a7b49 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -41,7 +41,7 @@ The operators in this file fall into three categories: # extend these functions (for example, LatitudeLongitudeGrid). # Calling a non existing function (for example Δxᶜᵃᶜ on an OrthogonalSphericalShellGrid) will throw an error because -# the associated one - dimensional function is not defined. This is a feature, not a bug. +# the associated one - dimensional function is not defined. spacing_1D(::Val{:x}, L1) = Symbol(:Δx, L1, :ᵃ, :ᵃ) spacing_1D(::Val{:y}, L1) = Symbol(:Δy, :ᵃ, L1, :ᵃ) From 01c0e4a28a3ef6000fdb26615a10ed3cebd526e2 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 1 Dec 2024 08:13:27 +1100 Subject: [PATCH 35/47] fix docstrings --- src/AbstractOperations/grid_metrics.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index a1efb0054d..b825cba9f8 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -53,7 +53,7 @@ on_architecture(to, gm::GridMetricOperation{LX, LY, LZ}) where {LX, LY, LZ} = indices(::GridMetricOperation) = default_indices(3) """ - gm = GridMetricOperation(L, metric, grid) + 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`. @@ -170,7 +170,7 @@ 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 ======== @@ -199,7 +199,7 @@ 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 ======== From 5ee71bf60ead5e09b6915d39c59d32439ac8a6b5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 12:25:28 +0100 Subject: [PATCH 36/47] better organization --- .../spacings_and_areas_and_volumes.jl | 159 ++++++++++++------ 1 file changed, 107 insertions(+), 52 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index f643227f89..7bb95f21de 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -42,35 +42,52 @@ The operators in this file fall into three categories: # Calling a non existing function (for example Δxᶜᵃᶜ on an OrthogonalSphericalShellGrid) will throw an error because # the associated one - dimensional function is not defined. This is a feature, not a bug. +for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) + Δxˡᵃᵃ = Symbol(:Δx, L1, :ᵃ, :ᵃ) + Δyᵃˡᵃ = Symbol(:Δy, :ᵃ, L1, :ᵃ) + Δzᵃᵃˡ = Symbol(:Δz, :ᵃ, :ᵃ, L1) + Δλˡᵃᵃ = Symbol(:Δλ, L1, :ᵃ, :ᵃ) + Δφᵃˡᵃ = Symbol(:Δφ, :ᵃ, L1, :ᵃ) + + Δxˡˡᵃ = Symbol(:Δx, L1, L2, :ᵃ) + Δyˡˡᵃ = Symbol(:Δy, L2, L1, :ᵃ) + Δzˡᵃˡ = Symbol(:Δz, L2, :ᵃ, L1) + Δλˡˡᵃ = Symbol(:Δλ, L1, L2, :ᵃ) + Δφˡˡᵃ = Symbol(:Δφ, L2, L1, :ᵃ) -spacing_1D(::Val{:x}, L1) = Symbol(:Δx, L1, :ᵃ, :ᵃ) -spacing_1D(::Val{:y}, L1) = Symbol(:Δy, :ᵃ, L1, :ᵃ) -spacing_1D(::Val{:z}, L1) = Symbol(:Δz, :ᵃ, :ᵃ, L1) -spacing_2D1(::Val{:x}, L1, L2) = Symbol(:Δx, L1, L2, :ᵃ) -spacing_2D1(::Val{:y}, L1, L2) = Symbol(:Δy, L2, L1, :ᵃ) -spacing_2D1(::Val{:z}, L1, L2) = Symbol(:Δz, :ᵃ, L2, L1) -spacing_2D2(::Val{:x}, L1, L2) = Symbol(:Δx, L1, :ᵃ, L2) -spacing_2D2(::Val{:y}, L1, L2) = Symbol(:Δy, :ᵃ, L1, L2) -spacing_2D2(::Val{:z}, L1, L2) = Symbol(:Δz, L2, :ᵃ, L1) -spacing_3D(::Val{:x}, L1, L2, L3) = Symbol(:Δx, L1, L2, L3) -spacing_3D(::Val{:y}, L1, L2, L3) = Symbol(:Δy, L2, L1, L3) -spacing_3D(::Val{:z}, L1, L2, L3) = Symbol(:Δz, L3, L2, L1) - -# Convenience Functions for all grids -# This metaprogramming loop defines all the allowed combinations of Δx, Δy, and Δz -# Note `:ᵃ` is not allowed for the location associated with the spacing -for dir in (:x, :y, :z) - for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) - spacing1D = spacing_1D(Val(dir), L1) - spacing2D1 = spacing_2D1(Val(dir), L1, L2) - spacing2D2 = spacing_2D2(Val(dir), L1, L2) - @eval @inline $spacing2D1(i, j, k, grid) = $spacing1D(i, j, k, grid) - @eval @inline $spacing2D2(i, j, k, grid) = $spacing1D(i, j, k, grid) + Δxˡᵃˡ = Symbol(:Δx, L1, :ᵃ, L2) + Δyᵃˡˡ = Symbol(:Δy, :ᵃ, L1, L2) + Δzᵃˡˡ = Symbol(:Δz, :ᵃ, L2, L1) + Δλˡᵃˡ = Symbol(:Δλ, L1, :ᵃ, L2) + Δφᵃˡˡ = Symbol(:Δφ, :ᵃ, L1, L2) + + @eval @inline $Δxˡˡᵃ(i, j, k, grid) = $Δxˡᵃᵃ(i, j, k, grid) + @eval @inline $Δxˡᵃˡ(i, j, k, grid) = $Δxˡᵃᵃ(i, j, k, grid) + + @eval @inline $Δyˡˡᵃ(i, j, k, grid) = $Δyᵃˡᵃ(i, j, k, grid) + @eval @inline $Δyᵃˡˡ(i, j, k, grid) = $Δyᵃˡᵃ(i, j, k, grid) + + @eval @inline $Δzˡᵃˡ(i, j, k, grid) = $Δzᵃᵃˡ(i, j, k, grid) + @eval @inline $Δzᵃˡˡ(i, j, k, grid) = $Δzᵃᵃˡ(i, j, k, grid) - for L3 in (:ᶜ, :ᶠ) - spacing3D = spacing_3D(Val(dir), L1, L2, L3) - @eval @inline $spacing3D(i, j, k, grid) = $spacing2D1(i, j, k, grid) - end + @eval @inline $Δλˡˡᵃ(i, j, k, grid) = $Δλˡᵃᵃ(i, j, k, grid) + @eval @inline $Δλˡᵃˡ(i, j, k, grid) = $Δλˡᵃᵃ(i, j, k, grid) + + @eval @inline $Δφˡˡᵃ(i, j, k, grid) = $Δφᵃˡᵃ(i, j, k, grid) + @eval @inline $Δφᵃˡˡ(i, j, k, grid) = $Δφᵃˡᵃ(i, j, k, grid) + + for L3 in (:ᶜ, :ᶠ) + Δxˡˡˡ = Symbol(:Δx, L1, L2, L3) + Δyˡˡˡ = Symbol(:Δy, L2, L1, L3) + Δzˡˡˡ = Symbol(:Δz, L2, L3, L1) + Δλˡˡˡ = Symbol(:Δλ, L1, L2, L3) + Δφˡˡˡ = Symbol(:Δφ, L2, L1, L3) + + @eval @inline $Δxˡˡˡ(i, j, k, grid) = $Δxˡˡᵃ(i, j, k, grid) + @eval @inline $Δyˡˡˡ(i, j, k, grid) = $Δyˡˡᵃ(i, j, k, grid) + @eval @inline $Δzˡˡˡ(i, j, k, grid) = $Δzˡᵃˡ(i, j, k, grid) + @eval @inline $Δλˡˡˡ(i, j, k, grid) = $Δλˡˡᵃ(i, j, k, grid) + @eval @inline $Δφˡˡˡ(i, j, k, grid) = $Δφˡˡᵃ(i, j, k, grid) end end @@ -100,13 +117,13 @@ end @inline Δyᵃᶠᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶠᵃ[j] @inline Δyᵃᶜᵃ(i, j, k, grid::RG) = @inbounds grid.Δyᵃᶜᵃ[j] -## XRegularRG +### XRegularRG @inline Δxᶠᵃᵃ(i, j, k, grid::RGX) = grid.Δxᶠᵃᵃ @inline Δxᶜᵃᵃ(i, j, k, grid::RGX) = grid.Δxᶜᵃᵃ -## YRegularRG +### YRegularRG @inline Δyᵃᶠᵃ(i, j, k, grid::RGY) = grid.Δyᵃᶠᵃ @inline Δyᵃᶜᵃ(i, j, k, grid::RGY) = grid.Δyᵃᶜᵃ @@ -115,14 +132,28 @@ end ##### LatitudeLongitude Grids (define both precomputed and non-precomputed metrics) ##### -# Precomputed metrics +### Curvilinear spacings + +@inline Δλᶜᵃᵃ(i, j, k, grid::LLG) = @inbounds grid.Δλᶜᵃᵃ[i] +@inline Δλᶠᵃᵃ(i, j, k, grid::LLG) = @inbounds grid.Δλᶠᵃᵃ[i] +@inline Δλᶜᵃᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δλᶜᵃᵃ +@inline Δλᶠᵃᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δλᶠᵃᵃ + +@inline Δφᵃᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δφᵃᶜᵃ[j] +@inline Δφᵃᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δφᵃᶠᵃ[j] +@inline Δφᵃᶜᵃ(i, j, k, grid::LLGY) = @inbounds grid.Δφᵃᶜᵃ +@inline Δφᵃᶠᵃ(i, j, k, grid::LLGY) = @inbounds grid.Δφᵃᶠᵃ + +### Linear spacings + +### Precomputed metrics @inline Δyᵃᶜᵃ(i, j, k, grid::LLGY) = grid.Δyᶠᶜᵃ @inline Δyᵃᶠᵃ(i, j, k, grid::LLGY) = grid.Δyᶜᶠᵃ @inline Δyᵃᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶠᶜᵃ[j] @inline Δyᵃᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶜᶠᵃ[j] -# On-the-fly metrics +### On-the-fly metrics @inline Δyᵃᶠᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶠᵃ) @inline Δyᵃᶜᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶜᵃ) @@ -136,7 +167,7 @@ end ##### ##### -##### LatitudeLongitudeGrid (only the Δx are required, Δy are 1D) +##### LatitudeLongitudeGrid (only the Δx are required, Δy, Δλ, and Δφ are 1D) ##### ### Pre computed metrics @@ -171,6 +202,20 @@ end ##### OrthogonalSphericalShellGrid (does not have one-dimensional spacings) ##### +### Curvilinear spacings + +@inline Δλᶜᶜᵃ(i, j, k, grid::OSSG) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), Center(), nothing) +@inline Δλᶠᶜᵃ(i, j, k, grid::OSSG) = δxᶠᵃᵃ(i, j, k, grid, λnode, Center(), Center(), nothing) +@inline Δλᶜᶠᵃ(i, j, k, grid::OSSG) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), Face(), nothing) +@inline Δλᶠᶠᵃ(i, j, k, grid::OSSG) = δxᶠᵃᵃ(i, j, k, grid, λnode, Center(), Face(), nothing) + +@inline Δφᶜᶜᵃ(i, j, k, grid::OSSG) = δyᵃᶠᵃ(i, j, k, grid, λnode, Center(), Face(), nothing) +@inline Δφᶠᶜᵃ(i, j, k, grid::OSSG) = δyᵃᶜᵃ(i, j, k, grid, λnode, Face(), Face(), nothing) +@inline Δφᶜᶠᵃ(i, j, k, grid::OSSG) = δyᵃᶠᵃ(i, j, k, grid, λnode, Center(), Center(), nothing) +@inline Δφᶠᶠᵃ(i, j, k, grid::OSSG) = δyᵃᶜᵃ(i, j, k, grid, λnode, Face(), Center(), nothing) + +### Linear spacings + @inline Δxᶜᶜᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δxᶜᶜᵃ[i, j] @inline Δxᶠᶜᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δxᶠᶜᵃ[i, j] @inline Δxᶜᶠᵃ(i, j, k, grid::OSSG) = @inbounds grid.Δxᶜᶠᵃ[i, j] @@ -258,7 +303,7 @@ end ##### ##### -##### Volumes!! (quite unambiguous) +##### Volumes!! (always 3D) ##### ##### @@ -290,36 +335,46 @@ for LX in (:Center, :Face, :Nothing) LYe = @eval $LY LZe = @eval $LZ - volume_function = Symbol(:V, location_code(LXe, LYe, LZe)) - @eval begin - @inline volume(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $volume_function(i, j, k, grid) + # General spacing functions + for dir in (:x, :y, :λ, :φ, :z) + func = Symbol(:Δ, dir) + metric = Symbol(:Δ, dir, location_code(LXe, LYe, LZe)) + + @eval begin + @inline $func(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $metric(i, j, k, grid) + end end - for op in (:Δ, :A), dir in (:x, :y, :z) - func = Symbol(op, dir) - metric = Symbol(op, dir, location_code(LXe, LYe, LZe)) + # General area functions + for dir in (:x, :y, :z) + func = Symbol(:A, dir) + metric = Symbol(:A, dir, location_code(LXe, LYe, LZe)) @eval begin @inline $func(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $metric(i, j, k, grid) end end + + # General volume function + volume_function = Symbol(:V, location_code(LXe, LYe, LZe)) + @eval begin + @inline volume(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $volume_function(i, j, k, grid) + end end end end -# Special curvilinear spacings for curvilinear grids -@inline Δλ(i, j, k, grid::LLG, ::Center, ℓy, ℓz) = @inbounds grid.Δλᶜᵃᵃ[i] -@inline Δλ(i, j, k, grid::LLG, ::Face, ℓy, ℓz) = @inbounds grid.Δλᶠᵃᵃ[i] -@inline Δλ(i, j, k, grid::LLGX, ::Center, ℓy, ℓz) = @inbounds grid.Δλᶜᵃᵃ -@inline Δλ(i, j, k, grid::LLGX, ::Face, ℓy, ℓz) = @inbounds grid.Δλᶠᵃᵃ +# One-dimensional convenience spacings (for grids that support them) -@inline Δφ(i, j, k, grid::LLG, ℓx, ::Center, ℓz) = @inbounds grid.Δφᵃᶜᵃ[j] -@inline Δφ(i, j, k, grid::LLG, ℓx, ::Face, ℓz) = @inbounds grid.Δφᵃᶠᵃ[j] -@inline Δφ(i, j, k, grid::LLGY, ℓx, ::Center, ℓz) = @inbounds grid.Δφᵃᶜᵃ -@inline Δφ(i, j, k, grid::LLGY, ℓx, ::Face, ℓz) = @inbounds grid.Δφᵃᶠᵃ +Δx(i, grid, ℓx) = Δx(i, 1, 1, grid, ℓx, nothing, nothing) +Δy(j, grid, ℓy) = Δy(1, j, 1, grid, nothing, ℓy, nothing) +Δz(k, grid, ℓz) = Δz(1, 1, k, grid, nothing, nothing, ℓz) +Δλ(i, grid, ℓx) = Δλ(i, 1, 1, grid, ℓx, nothing, nothing) +Δφ(j, grid, ℓy) = Δφ(1, j, 1, grid, nothing, ℓy, nothing) -@inline Δλ(i, j, k, grid::OSSG, ::Center, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), ℓy, ℓz) -@inline Δλ(i, j, k, grid::OSSG, ::Face, ℓy, ℓz) = δxᶜᵃᵃ(i, j, k, grid, λnode, Center(), ℓy, ℓz) +# Two-dimensional horizontal convenience spacings (for grids that support them) -@inline Δφ(i, j, k, grid::OSSG, ℓy, ::Center, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Face(), ℓz) -@inline Δφ(i, j, k, grid::OSSG, ℓy, ::Face, ℓz) = δyᵃᶜᵃ(i, j, k, grid, λnode, ℓx, Center(), ℓz) +Δx(i, j, grid, ℓx, ℓy) = Δx(i, j, 1, grid, ℓx, ℓy, nothing) +Δy(i, j, grid, ℓx, ℓy) = Δy(i, j, 1, grid, ℓx, ℓy, nothing) +Δλ(i, j, grid, ℓx, ℓy) = Δλ(i, j, 1, grid, ℓx, ℓy, nothing) +Δφ(i, j, grid, ℓx, ℓy) = Δφ(i, j, 1, grid, ℓx, ℓy, nothing) From 91083db225f33a9fc534a4a02a687d84a6695624 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 13:45:46 +0100 Subject: [PATCH 37/47] remove comment --- src/Operators/spacings_and_areas_and_volumes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 7bb95f21de..7aca0d6a04 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -41,7 +41,7 @@ The operators in this file fall into three categories: # extend these functions (for example, LatitudeLongitudeGrid). # Calling a non existing function (for example Δxᶜᵃᶜ on an OrthogonalSphericalShellGrid) will throw an error because -# the associated one - dimensional function is not defined. This is a feature, not a bug. +# the associated one - dimensional function is not defined. for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) Δxˡᵃᵃ = Symbol(:Δx, L1, :ᵃ, :ᵃ) Δyᵃˡᵃ = Symbol(:Δy, :ᵃ, L1, :ᵃ) From 6c11e7363041cb9ceb98ba0d5aed96143e09cd31 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 13:54:11 +0100 Subject: [PATCH 38/47] correct docstring --- src/AbstractOperations/grid_metrics.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index b825cba9f8..0b516e1fa0 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -73,8 +73,8 @@ 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) + ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + └── Δzᶜᶜᶜ at (Center, Center, Center) julia> c .= 1; From e7b38ec625815592ecd2135682c7ddd798c4a655 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:09:15 +0100 Subject: [PATCH 39/47] Update Operators.jl --- src/Operators/Operators.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index d9a7941622..d6384dd5ca 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -17,6 +17,9 @@ export Azᵃᶜᶜ, Azᵃᶠᶠ, Azᶜᵃᶜ, Azᶠᵃᶠ, Azᶜᶜᵃ, Azᶠᶠ # Volumes export Vᶠᶠᶠ, Vᶠᶠᶜ, Vᶠᶜᶠ, Vᶠᶜᶜ, Vᶜᶠᶠ, Vᶜᶠᶜ, Vᶜᶜᶠ, Vᶜᶜᶜ +# General metric operators +export xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume + # Product between spacings and fields export Δx_qᶠᶠᶠ, Δx_qᶠᶠᶜ, Δx_qᶠᶜᶠ, Δx_qᶠᶜᶜ, Δx_qᶜᶠᶠ, Δx_qᶜᶠᶜ, Δx_qᶜᶜᶠ, Δx_qᶜᶜᶜ export Δy_qᶠᶠᶠ, Δy_qᶠᶠᶜ, Δy_qᶠᶜᶠ, Δy_qᶠᶜᶜ, Δy_qᶜᶠᶠ, Δy_qᶜᶠᶜ, Δy_qᶜᶜᶠ, Δy_qᶜᶜᶜ @@ -102,9 +105,6 @@ const LLGF = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing} const LLGFX = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Number} const LLGFY = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Any, <:Number} -# General metric operators -export xspacing, yspacing, zspacing, λspacing, φspacing, xarea, yarea, zarea, volume - include("difference_operators.jl") include("interpolation_operators.jl") include("interpolation_utils.jl") From 6236cf349c616d8b7af7b57d255ddc526a3431d3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:11:26 +0100 Subject: [PATCH 40/47] correction --- src/Operators/spacings_and_areas_and_volumes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 7aca0d6a04..8e14d55f99 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -209,10 +209,10 @@ end @inline Δλᶜᶠᵃ(i, j, k, grid::OSSG) = δxᶜᵃᵃ(i, j, k, grid, λnode, Face(), Face(), nothing) @inline Δλᶠᶠᵃ(i, j, k, grid::OSSG) = δxᶠᵃᵃ(i, j, k, grid, λnode, Center(), Face(), nothing) -@inline Δφᶜᶜᵃ(i, j, k, grid::OSSG) = δyᵃᶠᵃ(i, j, k, grid, λnode, Center(), Face(), nothing) +@inline Δφᶜᶜᵃ(i, j, k, grid::OSSG) = δyᵃᶜᵃ(i, j, k, grid, λnode, Center(), Face(), nothing) @inline Δφᶠᶜᵃ(i, j, k, grid::OSSG) = δyᵃᶜᵃ(i, j, k, grid, λnode, Face(), Face(), nothing) @inline Δφᶜᶠᵃ(i, j, k, grid::OSSG) = δyᵃᶠᵃ(i, j, k, grid, λnode, Center(), Center(), nothing) -@inline Δφᶠᶠᵃ(i, j, k, grid::OSSG) = δyᵃᶜᵃ(i, j, k, grid, λnode, Face(), Center(), nothing) +@inline Δφᶠᶠᵃ(i, j, k, grid::OSSG) = δyᵃᶠᵃ(i, j, k, grid, λnode, Face(), Center(), nothing) ### Linear spacings From 52fe6c8753f041becdf55201b510ad740424fe97 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 9 Dec 2024 21:40:13 +0100 Subject: [PATCH 41/47] correct this later --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index fd968a0c3a..48cb239768 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,7 +10,7 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) +for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) for dir in (:x, :y, :z), operator in (:Δ, :A) metric = Symbol(operator, dir, LX, LY, LZ) From a86183deba97bba77ab639d460d2519d45bf8815 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 9 Dec 2024 21:46:02 +0100 Subject: [PATCH 42/47] fix docs --- src/AbstractOperations/grid_metrics.jl | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 0b516e1fa0..9b9718485f 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -68,13 +68,7 @@ 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 -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_dz = c * Δz; # returns BinaryOperation between Field and GridMetricOperation julia> c .= 1; @@ -104,7 +98,7 @@ 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: Δx (generic function with 27 methods) +├── kernel_function: Δx (generic function with 29 methods) └── arguments: ("Center", "Center", "Center") ``` """ @@ -130,7 +124,7 @@ 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: Δy (generic function with 27 methods) +├── kernel_function: Δy (generic function with 29 methods) └── arguments: ("Center", "Face", "Center") ``` """ @@ -156,7 +150,7 @@ 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: Δz (generic function with 27 methods) +├── kernel_function: Δz (generic function with 28 methods) └── arguments: ("Center", "Center", "Face") ``` """ @@ -185,7 +179,7 @@ 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: Δλ (generic function with 5 methods) +├── kernel_function: Δλ (generic function with 29 methods) └── arguments: ("Center", "Face", "Center") ``` """ @@ -214,7 +208,7 @@ 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: Δφ (generic function with 5 methods) +├── kernel_function: Δφ (generic function with 29 methods) └── arguments: ("Center", "Face", "Center") ``` """ From 109756db85d0f8cb17f3ce6b0a2284a999079434 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 9 Dec 2024 22:41:06 +0100 Subject: [PATCH 43/47] correcting also immersed spacings --- .../immersed_grid_metrics.jl | 63 +++++++++++++------ 1 file changed, 45 insertions(+), 18 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 48cb239768..bb5b0f98b0 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -5,33 +5,60 @@ import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, ext # Grid metrics for ImmersedBoundaryGrid # -# All grid metrics are defined here. +# The "basic" grid metric functions are defined for the underlying grids. # # For non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. -for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) - for dir in (:x, :y, :z), operator in (:Δ, :A) - - metric = Symbol(operator, dir, LX, LY, LZ) - @eval begin - import Oceananigans.Operators: $metric - @inline $metric(i, j, k, ibg::IBG) = $metric(i, j, k, ibg.underlying_grid) - end - end - - volume = Symbol(:V, LX, LY, LZ) - @eval begin - import Oceananigans.Operators: $volume - @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid) - end -end - coordinates(grid::IBG) = coordinates(grid.underlying_grid) +# Vertical spacings + @inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid) @inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid) +# 1D Horizontal spacings + +@inline Δxᶠᵃᵃ(i, j, k, ibg::RGIBG) = Δxᶠᵃᵃ(i, j, k, ibg.underlying_grid) +@inline Δxᶜᵃᵃ(i, j, k, ibg::RGIBG) = Δxᶜᵃᵃ(i, j, k, ibg.underlying_grid) + +@inline Δyᵃᶠᵃ(i, j, k, ibg::Union{RGIBG, LLIBG}) = Δyᵃᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δyᵃᶜᵃ(i, j, k, ibg::Union{RGIBG, LLIBG}) = Δyᵃᶜᵃ(i, j, k, ibg.underlying_grid) + +@inline Δλᶜᵃᵃ(i, j, k, ibg::LLIBG) = Δλᶜᵃᵃ(i, j, k, ibg.underlying_grid) +@inline Δλᶠᵃᵃ(i, j, k, ibg::LLIBG) = Δλᶠᵃᵃ(i, j, k, ibg.underlying_grid) +@inline Δφᵃᶜᵃ(i, j, k, ibg::LLIBG) = Δφᵃᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δφᵃᶠᵃ(i, j, k, ibg::LLIBG) = Δφᵃᶠᵃ(i, j, k, ibg.underlying_grid) + +# 2D Horizontal spacings + +@inline Δxᶜᶠᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Δxᶜᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δxᶠᶜᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Δxᶠᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δxᶠᶠᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Δxᶠᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δxᶜᶜᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Δxᶜᶜᵃ(i, j, k, ibg.underlying_grid) + +@inline Δyᶜᶜᵃ(i, j, k, ibg::OSIBG) = Δyᶜᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δyᶠᶜᵃ(i, j, k, ibg::OSIBG) = Δyᶠᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δyᶜᶠᵃ(i, j, k, ibg::OSIBG) = Δyᶜᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δyᶠᶠᵃ(i, j, k, ibg::OSIBG) = Δyᶠᶠᵃ(i, j, k, ibg.underlying_grid) + +@inline Δλᶜᶜᵃ(i, j, k, ibg::OSIBG) = Δλᶜᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δλᶠᶜᵃ(i, j, k, ibg::OSIBG) = Δλᶠᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δλᶜᶠᵃ(i, j, k, ibg::OSIBG) = Δλᶜᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δλᶠᶠᵃ(i, j, k, ibg::OSIBG) = Δλᶠᶠᵃ(i, j, k, ibg.underlying_grid) + +@inline Δφᶜᶜᵃ(i, j, k, ibg::OSIBG) = Δφᶜᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δφᶠᶜᵃ(i, j, k, ibg::OSIBG) = Δφᶠᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Δφᶜᶠᵃ(i, j, k, ibg::OSIBG) = Δφᶜᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Δφᶠᶠᵃ(i, j, k, ibg::OSIBG) = Δφᶠᶠᵃ(i, j, k, ibg.underlying_grid) + +# Areas + +@inline Azᶠᶜᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Azᶠᶜᵃ(i, j, k, ibg.underlying_grid) +@inline Azᶜᶠᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Azᶜᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Azᶠᶠᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Azᶠᶠᵃ(i, j, k, ibg.underlying_grid) +@inline Azᶜᶜᵃ(i, j, k, ibg::Union{LLIBG, OSIBG}) = Azᶜᶜᵃ(i, j, k, ibg.underlying_grid) + # Extend both 2D and 3D methods @inline intrinsic_vector(i, j, k, ibg::IBG, u, v) = intrinsic_vector(i, j, k, ibg.underlying_grid, u, v) @inline extrinsic_vector(i, j, k, ibg::IBG, u, v) = extrinsic_vector(i, j, k, ibg.underlying_grid, u, v) From 1cd9f3213430fc54e27dc6c154937d0613990744 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 9 Dec 2024 23:04:22 +0100 Subject: [PATCH 44/47] import correct metrics --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index bb5b0f98b0..8389ac652d 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -1,7 +1,12 @@ using Oceananigans.AbstractOperations: GridMetricOperation import Oceananigans.Grids: coordinates -import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ, intrinsic_vector, extrinsic_vector +import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ +import Oceananigans.Operators: Δxᶠᵃᵃ, Δxᶜᵃᵃ, Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶠᶠᵃ, Δxᶜᶜᵃ +import Oceananigans.Operators: Δyᵃᶠᵃ, Δyᵃᶜᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Δyᶠᶠᵃ, Δyᶜᶜᵃ +import Oceananigans.Operators: Azᶠᶜᵃ, Azᶜᶠᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ + +import Oceananigans.Operators: intrinsic_vector, extrinsic_vector # Grid metrics for ImmersedBoundaryGrid # From b4cea75ca48718861ac1a6eae9679c309a4f9063 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 10 Dec 2024 09:57:09 +0100 Subject: [PATCH 45/47] import spacings --- src/ImmersedBoundaries/partial_cell_bottom.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 36cc32487b..87840bf28d 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -2,6 +2,8 @@ using Oceananigans.Utils: prettysummary using Oceananigans.Fields: fill_halo_regions! using Printf +import Oceananigans.Operators: Δzᶜᶜᶜ, Δzᶜᶜᶠ, Δzᶜᶠᶜ, Δzᶜᶠᶠ, Δzᶠᶜᶜ, Δzᶠᶜᶠ, Δzᶠᶠᶜ, Δzᶠᶠᶠ + ##### ##### PartialCellBottom ##### From f40627a1a966601d4d190b9b92a77b338773422b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 10 Dec 2024 10:11:46 +0100 Subject: [PATCH 46/47] add comment --- src/ImmersedBoundaries/immersed_grid_metrics.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index 8389ac652d..28018da192 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -10,9 +10,13 @@ import Oceananigans.Operators: intrinsic_vector, extrinsic_vector # Grid metrics for ImmersedBoundaryGrid # -# The "basic" grid metric functions are defined for the underlying grids. +# The grid metric functions "specialized" for the underlying grids in the +# Operators module are extended for immersed boundary grids here. # -# For non "full-cell" immersed boundaries, grid metric functions +# The other grid metric functions (for example 3D metrics) are general and do not need +# extension for "full-cell" immersed boundaries. +# +# However, for non "full-cell" immersed boundaries, grid metric functions # must be extended for the specific immersed boundary grid in question. coordinates(grid::IBG) = coordinates(grid.underlying_grid) From 960452ec1b67d652a7c6aa237a431355bbb60e97 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 10 Dec 2024 11:22:22 +0100 Subject: [PATCH 47/47] doctests --- src/AbstractOperations/grid_metrics.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 9b9718485f..0452e043e2 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -61,7 +61,6 @@ at the same location as the `AbstractField`. Example ======= ```jldoctest - julia> using Oceananigans julia> using Oceananigans.Operators: Δz