diff --git a/docs/src/fields.md b/docs/src/fields.md index 6d39014d78..e37215ee48 100644 --- a/docs/src/fields.md +++ b/docs/src/fields.md @@ -3,7 +3,7 @@ `Field`s and its relatives are core Oceananigans data structures. `Field`s are more or less arrays of `data` located on a `grid`, whose entries correspond to the average value of some quantity over some finite-sized volume. -`Field`s also may contain `boundary_conditions`, may be computed from an `operand` +`Field`s also may contain `boundary_conditions`, may be computed from an `operand` or expression involving other fields, and may cover only a portion of the total `indices` spanned by the grid. @@ -43,10 +43,10 @@ primary mesh. For example, the primary or `Center` cell spacings in ``z`` are ```jldoctest fields -zspacings(grid, Center()) +zspacings(grid, Center())[:, :, 1:4] # output -4-element view(OffsetArray(::Vector{Float64}, 0:5), 1:4) with eltype Float64: +4-element Vector{Float64}: 0.1 0.19999999999999998 0.3 @@ -57,10 +57,10 @@ corresponding to cell interfaces located at `z = [0, 0.1, 0.3, 0.6, 1]`. But then for the grid which is staggered in `z` relative to the primary mesh, ```jldoctest fields -zspacings(grid, Face()) +zspacings(grid, Face())[:, :, 1:5] # output -5-element view(OffsetArray(::Vector{Float64}, -1:5), 1:5) with eltype Float64: +5-element Vector{Float64}: 0.1 0.15000000000000002 0.24999999999999994 @@ -277,7 +277,7 @@ Note that indexing into `c` is the same as indexing into `c.data`. ```jldoctest fields c[:, :, :] == c.data - + # output true ``` @@ -333,7 +333,7 @@ non-`Flat` directions must be included. For example, to `set!` on a one-dimensio one_d_grid = RectilinearGrid(size=7, x=(0, 7), topology=(Periodic, Flat, Flat)) one_d_c = CenterField(one_d_grid) -# The one-dimensional grid varies only in `x` +# The one-dimensional grid varies only in `x` still_pretty_fun(x) = 3x set!(one_d_c, still_pretty_fun) diff --git a/docs/src/grids.md b/docs/src/grids.md index 71ec02adc0..b02a045d82 100644 --- a/docs/src/grids.md +++ b/docs/src/grids.md @@ -219,7 +219,7 @@ current_figure() ## Once more with feeling -In summary, making a grid requires +In summary, making a grid requires * The machine architecture, or whether data is stored on the CPU, GPU, or distributed across multiple devices or nodes. * Information about the domain geometry. Domains can take a variety of shapes, including @@ -413,7 +413,7 @@ hidespines!(axy) axΔy = Axis(fig[2, 1]; xlabel = "y (m)", ylabel = "y-spacing (m)") scatter!(axΔy, yc, Δy) -hidespines!(axΔy, :t, :r) +hidespines!(axΔy, :t, :r) axz = Axis(fig[3, 1], title="z-grid") lines!(axz, [-Lz, 0], [0, 0], color=:gray) @@ -440,18 +440,17 @@ using Oceananigans ```@example latlon_nodes grid = LatitudeLongitudeGrid(size = (1, 44), - longitude = (0, 1), + longitude = (0, 1), latitude = (0, 88), topology = (Bounded, Bounded, Flat)) -φ = φnodes(grid, Center()) Δx = xspacings(grid, Center(), Center()) using CairoMakie fig = Figure(size=(600, 400)) ax = Axis(fig[1, 1], xlabel="Zonal spacing on 2 degree grid (km)", ylabel="Latitude (degrees)") -scatter!(ax, Δx ./ 1e3, φ) +scatter!(ax, Δx / 1e3) current_figure() ``` @@ -473,7 +472,7 @@ m = 2 # spacing at the equator in degrees function latitude_faces(j) if j == 1 # equator return 0 - else # crudely estimate the location of the jth face + else # crudely estimate the location of the jth face φ₋ = latitude_faces(j-1) φ′ = φ₋ + m * scale_factor(φ₋) / 2 return φ₋ + m * scale_factor(φ′) @@ -493,7 +492,7 @@ grid = LatitudeLongitudeGrid(size = (Nx, Ny), 180×28×1 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo and with precomputed metrics ├── longitude: Bounded λ ∈ [0.0, 360.0] regularly spaced with Δλ=2.0 ├── latitude: Bounded φ ∈ [0.0, 77.2679] variably spaced with min(Δφ)=2.0003, max(Δφ)=6.95319 -└── z: Flat z +└── z: Flat z ``` We've also illustrated the construction of a grid that is `Flat` in the vertical direction. @@ -508,7 +507,7 @@ m = 2 # spacing at the equator in degrees function latitude_faces(j) if j == 1 # equator return 0 - else # crudely estimate the location of the jth face + else # crudely estimate the location of the jth face φ₋ = latitude_faces(j-1) φ′ = φ₋ + m * scale_factor(φ₋) / 2 return φ₋ + m * scale_factor(φ′) @@ -530,17 +529,17 @@ grid = LatitudeLongitudeGrid(size = (Nx, Ny), ```@example plot φ = φnodes(grid, Center()) -Δx = xspacings(grid, Center(), Center(), with_halos=true)[1:Ny] -Δy = yspacings(grid, Center())[1:Ny] +Δx = xspacings(grid, Center(), Center())[1, 1:Ny] +Δy = yspacings(grid, Center(), Center())[1, 1:Ny] using CairoMakie fig = Figure(size=(800, 400), title="Spacings on a Mercator grid") axx = Axis(fig[1, 1], xlabel="Zonal spacing (km)", ylabel="Latitude (degrees)") -scatter!(axx, Δx ./ 1e3, φ) +scatter!(axx, Δx / 1e3, φ) axy = Axis(fig[1, 2], xlabel="Meridional spacing (km)") -scatter!(axy, Δy ./ 1e3, φ) +scatter!(axy, Δy / 1e3, φ) hidespines!(axx, :t, :r) hidespines!(axy, :t, :l, :r) diff --git a/examples/ocean_wind_mixing_and_convection.jl b/examples/ocean_wind_mixing_and_convection.jl index 01629ef4ff..6daa4c0dc0 100644 --- a/examples/ocean_wind_mixing_and_convection.jl +++ b/examples/ocean_wind_mixing_and_convection.jl @@ -50,13 +50,13 @@ h(k) = (k - 1) / Nz ## Linear near-surface generator ζ₀(k) = 1 + (h(k) - 1) / refinement -## Bottom-intensified stretching function +## Bottom-intensified stretching function Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching)) ## Generating function z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1) -grid = RectilinearGrid(size = (Nx, Nx, Nz), +grid = RectilinearGrid(size = (Nx, Nx, Nz), x = (0, Lx), y = (0, Ly), z = z_faces) @@ -66,8 +66,8 @@ grid = RectilinearGrid(size = (Nx, Nx, Nz), fig = Figure(size=(1200, 800)) ax = Axis(fig[1, 1], ylabel = "Depth (m)", xlabel = "Vertical spacing (m)") -lines!(ax, zspacings(grid, Center()), znodes(grid, Center())) -scatter!(ax, zspacings(grid, Center()), znodes(grid, Center())) +lines!(ax, zspacings(grid, Center())) +scatter!(ax, zspacings(grid, Center())) current_figure() #hide fig diff --git a/examples/tilted_bottom_boundary_layer.jl b/examples/tilted_bottom_boundary_layer.jl index d351e3f2bf..c6b0b15eb3 100644 --- a/examples/tilted_bottom_boundary_layer.jl +++ b/examples/tilted_bottom_boundary_layer.jl @@ -3,7 +3,7 @@ # away from a constant along-slope (y-direction) velocity constant density stratification. # This perturbation develops into a turbulent bottom boundary layer due to momentum # loss at the bottom boundary modeled with a quadratic drag law. -# +# # This example illustrates # # * changing the direction of gravitational acceleration in the buoyancy model; @@ -33,7 +33,7 @@ Nz = 64 ## Creates a grid with near-constant spacing `refinement * Lz / Nz` ## near the bottom: refinement = 1.8 # controls spacing near surface (higher means finer spaced) -stretching = 10 # controls rate of stretching at bottom +stretching = 10 # controls rate of stretching at bottom ## "Warped" height coordinate h(k) = (Nz + 1 - k) / Nz @@ -41,7 +41,7 @@ h(k) = (Nz + 1 - k) / Nz ## Linear near-surface generator ζ(k) = 1 + (h(k) - 1) / refinement -## Bottom-intensified stretching function +## Bottom-intensified stretching function Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching)) ## Generating function @@ -56,11 +56,9 @@ grid = RectilinearGrid(topology = (Periodic, Flat, Bounded), using CairoMakie -lines(zspacings(grid, Center()), znodes(grid, Center()), - axis = (ylabel = "Depth (m)", - xlabel = "Vertical spacing (m)")) - -scatter!(zspacings(grid, Center()), znodes(grid, Center())) +scatterlines(zspacings(grid, Center()), + axis = (ylabel = "Depth (m)", + xlabel = "Vertical spacing (m)")) current_figure() #hide diff --git a/ext/OceananigansMakieExt.jl b/ext/OceananigansMakieExt.jl index e88cb8f0a7..6627ce8491 100644 --- a/ext/OceananigansMakieExt.jl +++ b/ext/OceananigansMakieExt.jl @@ -2,6 +2,7 @@ module OceananigansMakieExt using Oceananigans using Oceananigans.Grids: OrthogonalSphericalShellGrid +using Oceananigans.Fields: AbstractField using Oceananigans.AbstractOperations: AbstractOperation using Oceananigans.Architectures: on_architecture using Oceananigans.ImmersedBoundaries: mask_immersed_field! @@ -14,7 +15,7 @@ import Makie: args_preferred_axis # do not overstate a preference for being plotted in a 3D LScene. # Because often we are trying to plot 1D and 2D Field, even though # (perhaps incorrectly) all Field are AbstractArray{3}. -args_preferred_axis(::Field) = nothing +args_preferred_axis(::AbstractField) = nothing function drop_singleton_indices(N) if N == 1 @@ -146,3 +147,4 @@ function convert_arguments(pl::Type{<:AbstractPlot}, ξ1::AbstractArray, ξ2::Ab end end # module + diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index f873242aa1..4205135ef0 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -1,26 +1,29 @@ using Adapt using Oceananigans.Operators -using Oceananigans.Fields: default_indices +using Oceananigans.Grids: AbstractGrid +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 +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 +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 +struct VolumeMetric <: AbstractGridMetric end metric_function_prefix(::VolumeMetric) = :V @@ -33,7 +36,7 @@ const Δy = YSpacingMetric() Instance of `ZSpacingMetric` that generates `BinaryOperation`s between `AbstractField`s and the vertical grid spacing evaluated -at the same location as the `AbstractField`. +at the same location as the `AbstractField`. `Δx` and `Δy` play a similar role for horizontal grid spacings. @@ -133,7 +136,7 @@ Adapt.adapt_structure(to, gm::GridMetricOperation{LX, LY, LZ}) where {LX, LY, LZ 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) @@ -141,3 +144,156 @@ indices(::GridMetricOperation) = default_indices(3) # Special constructor for BinaryOperation GridMetricOperation(L, metric, grid) = GridMetricOperation{L[1], L[2], L[3]}(metric_function(L, metric), grid) + +##### +##### Spacings +##### + +""" + xspacings(grid, ℓx, ℓy, ℓz) + +Return a `KernelFunctionOperation` that computes the grid spacings for `grid` +in the ``x`` direction at location `ℓx, ℓy, ℓz`. + +Examples +======== +```jldoctest +julia> using Oceananigans + +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) +└── 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) + return Δx_op +end + +""" + yspacings(grid, ℓx, ℓy, ℓz) + +Return a `KernelFunctionOperation` that computes the grid spacings for `grid` +in the ``y`` direction at location `ℓx, ℓy, ℓz`. + +Examples +======== +```jldoctest +julia> using Oceananigans + +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) +└── 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) + return Δy_op +end + +""" + zspacings(grid, ℓx, ℓy, ℓz) + +Return a `KernelFunctionOperation` that computes the grid spacings for `grid` +in the ``z`` direction at location `ℓx, ℓy, ℓz`. + +Examples +======== +```jldoctest +julia> using Oceananigans + +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) +└── 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) + return Δz_op +end + +""" + λspacings(grid, ℓx, ℓy, ℓz) + +Return a `KernelFunctionOperation` that computes the grid spacings for `grid` +in the ``z`` direction at location `ℓx, ℓy, ℓz`. + +Examples +======== +```jldoctest +julia> using Oceananigans + +julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25), + longitude = (-180, 180), + latitude = (-85, 85), + z = (-1000, 0)); + +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) +└── 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) + return Δλ_op +end + +""" + φspacings(grid, ℓx, ℓy, ℓz) + +Return a `KernelFunctionOperation` that computes the grid spacings for `grid` +in the ``z`` direction at location `ℓx, ℓy, ℓz`. + +Examples +======== +```jldoctest +julia> using Oceananigans + +julia> grid = LatitudeLongitudeGrid(size=(36, 34, 25), + longitude = (-180, 180), + latitude = (-85, 85), + z = (-1000, 0)); + +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) +└── 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) + return Δφ_op +end + +@inline xspacings(field::AbstractField) = xspacings(field.grid, location(field)...) +@inline yspacings(field::AbstractField) = yspacings(field.grid, location(field)...) +@inline zspacings(field::AbstractField) = zspacings(field.grid, location(field)...) +@inline λspacings(field::AbstractField) = λspacings(field.grid, location(field)...) +@inline φspacings(field::AbstractField) = φspacings(field.grid, location(field)...) + +# Some defaults for e.g. easy CFL computations. +@inline xspacings(grid::AbstractGrid) = xspacings(grid, Center(), Center(), Center()) +@inline yspacings(grid::AbstractGrid) = yspacings(grid, Center(), Center(), Center()) +@inline zspacings(grid::AbstractGrid) = zspacings(grid, Center(), Center(), Center()) +@inline λspacings(grid::AbstractGrid) = λspacings(grid, Center(), Center(), Center()) +@inline φspacings(grid::AbstractGrid) = φspacings(grid, Center(), Center(), Center()) diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index 6acd5b7ad9..47eb75be5b 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -66,14 +66,14 @@ end @inline function fractional_x_index(x, locs, grid::XRegularRG) x₀ = xnode(1, 1, 1, grid, locs...) - Δx = xspacings(grid, locs...) + Δx = @inbounds first(xspacings(grid, locs...)) FT = eltype(grid) return convert(FT, (x - x₀) / Δx) end @inline function fractional_x_index(λ, locs, grid::XRegularLLG) λ₀ = λnode(1, 1, 1, grid, locs...) - Δλ = λspacings(grid, locs...) + Δλ = @inbounds first(λspacings(grid, locs...)) FT = eltype(grid) return convert(FT, (λ - λ₀) / Δλ) end @@ -98,14 +98,14 @@ end @inline function fractional_y_index(y, locs, grid::YRegularRG) y₀ = ynode(1, 1, 1, grid, locs...) - Δy = yspacings(grid, locs...) + Δy = @inbounds first(yspacings(grid, locs...)) FT = eltype(grid) return convert(FT, (y - y₀) / Δy) end @inline function fractional_y_index(φ, locs, grid::YRegularLLG) φ₀ = φnode(1, 1, 1, grid, locs...) - Δφ = φspacings(grid, locs...) + Δφ = @inbounds first(φspacings(grid, locs...)) FT = eltype(grid) return convert(FT, (φ - φ₀) / Δφ) end @@ -132,7 +132,7 @@ ZRegGrid = Union{ZRegularRG, ZRegularLLG, ZRegOrthogonalSphericalShellGrid} @inline function fractional_z_index(z::FT, locs, grid::ZRegGrid) where FT z₀ = znode(1, 1, 1, grid, locs...) - Δz = zspacings(grid, locs...) + Δz = @inbounds first(zspacings(grid, locs...)) return convert(FT, (z - z₀) / Δz) end diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl index 27bd8f3a4a..519233561b 100644 --- a/src/Grids/Grids.jl +++ b/src/Grids/Grids.jl @@ -16,7 +16,8 @@ export ξnode, ηnode, rnode export xnode, ynode, znode, λnode, φnode export xnodes, ynodes, znodes, λnodes, φnodes export spacings -export xspacings, yspacings, zspacings, xspacing, yspacing, zspacing +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ᶠᶠᵃ export offset_data, new_data diff --git a/src/Grids/latitude_longitude_grid.jl b/src/Grids/latitude_longitude_grid.jl index 205b02aa6f..bd411ea5d1 100644 --- a/src/Grids/latitude_longitude_grid.jl +++ b/src/Grids/latitude_longitude_grid.jl @@ -655,6 +655,11 @@ end @inline xnodes(grid::LLG, ℓx, ℓy, ℓz; with_halos=false) = xnodes(grid, ℓx, ℓy; with_halos) @inline ynodes(grid::LLG, ℓx, ℓy, ℓz; with_halos=false) = ynodes(grid, ℓy; with_halos) +@inline λnodes(grid::LLG, ℓx::F; with_halos=false) = _property(grid.λᶠᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) +@inline λnodes(grid::LLG, ℓx::C; with_halos=false) = _property(grid.λᶜᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) +@inline φnodes(grid::LLG, ℓy::F; with_halos=false) = _property(grid.φᵃᶠᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) +@inline φnodes(grid::LLG, ℓy::C; with_halos=false) = _property(grid.φᵃᶜᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) + # Generalized coordinates @inline ξnodes(grid::LLG, ℓx; kwargs...) = λnodes(grid, ℓx; kwargs...) @inline ηnodes(grid::LLG, ℓy; kwargs...) = φnodes(grid, ℓy; kwargs...) @@ -664,76 +669,10 @@ end @inline ηnodes(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = φnodes(grid, ℓy; kwargs...) @inline rnodes(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...) -##### -##### Grid spacings in x, y, z (in meters) -##### - -@inline xspacings(grid::LLG, ℓx::C, ℓy::C; with_halos=false) = _property(grid.Δxᶜᶜᵃ, ℓx, ℓy, - topology(grid, 1), topology(grid, 2), - size(grid, 1), size(grid, 2), with_halos) - -@inline xspacings(grid::LLG, ℓx::C, ℓy::F; with_halos=false) = _property(grid.Δxᶜᶠᵃ, ℓx, ℓy, - topology(grid, 1), topology(grid, 2), - size(grid, 1), size(grid, 2), with_halos) - -@inline xspacings(grid::LLG, ℓx::F, ℓy::C; with_halos=false) = _property(grid.Δxᶠᶜᵃ, ℓx, ℓy, - topology(grid, 1), topology(grid, 2), - size(grid, 1), size(grid, 2), with_halos) - -@inline xspacings(grid::LLG, ℓx::F, ℓy::F; with_halos=false) = _property(grid.Δxᶠᶠᵃ, ℓx, ℓy, - topology(grid, 1), topology(grid, 2), - size(grid, 1), size(grid, 2), with_halos) - -@inline xspacings(grid::HRegularLLG, ℓx::C, ℓy::C; with_halos=false) = _property(grid.Δxᶜᶜᵃ, ℓy, topology(grid, 2), - size(grid, 2), with_halos) - -@inline xspacings(grid::HRegularLLG, ℓx::C, ℓy::F; with_halos=false) = _property(grid.Δxᶜᶠᵃ, ℓy, topology(grid, 2), - size(grid, 2), with_halos) - -@inline xspacings(grid::HRegularLLG, ℓx::F, ℓy::C; with_halos=false) = _property(grid.Δxᶠᶜᵃ, ℓy, topology(grid, 2), - size(grid, 2), with_halos) - -@inline xspacings(grid::HRegularLLG, ℓx::F, ℓy::F; with_halos=false) = _property(grid.Δxᶠᶠᵃ, ℓy, topology(grid, 2), - size(grid, 2), with_halos) - - -@inline yspacings(grid::YNonRegularLLG, ℓx::C, ℓy::F; with_halos=false) = _property(grid.Δyᶜᶠᵃ, ℓy, topoloy(grid, 2), - size(grid, 2), with_halos) - -@inline yspacings(grid::YNonRegularLLG, ℓx::F, ℓy::C; with_halos=false) = _property(grid.Δyᶠᶜᵃ, ℓy, topoloy(grid, 2), - size(grid, 2), with_halos) - -@inline yspacings(grid::YRegularLLG, ℓx, ℓy; with_halos=false) = yspacings(grid, ℓy; with_halos) -@inline yspacings(grid, ℓy::C; kwargs...) = grid.Δyᶠᶜᵃ -@inline yspacings(grid, ℓy::F; kwargs...) = grid.Δyᶜᶠᵃ - -@inline xspacings(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = xspacings(grid, ℓx, ℓy; kwargs...) -@inline yspacings(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = yspacings(grid, ℓx, ℓy; kwargs...) -@inline zspacings(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = zspacings(grid, ℓz; kwargs...) - -##### -##### Grid spacings in λ, φ (in degrees) -##### - -@inline λnodes(grid::LLG, ℓx::F; with_halos=false) = _property(grid.λᶠᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline λnodes(grid::LLG, ℓx::C; with_halos=false) = _property(grid.λᶜᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline φnodes(grid::LLG, ℓy::F; with_halos=false) = _property(grid.φᵃᶠᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) -@inline φnodes(grid::LLG, ℓy::C; with_halos=false) = _property(grid.φᵃᶜᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) - ##### ##### Grid spacings ##### -@inline λspacings(grid::LLG, ℓx::C; with_halos=false) = _property(grid.Δλᶜᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline λspacings(grid::LLG, ℓx::F; with_halos=false) = _property(grid.Δλᶠᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline φspacings(grid::LLG, ℓy::C; with_halos=false) = _property(grid.Δφᵃᶜᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) -@inline φspacings(grid::LLG, ℓy::F; with_halos=false) = _property(grid.Δφᵃᶠᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) -@inline zspacings(grid::LLG, ℓz::C; with_halos=false) = _property(grid.Δzᵃᵃᶜ, ℓz, topology(grid, 3), size(grid, 3), with_halos) -@inline zspacings(grid::LLG, ℓz::F; with_halos=false) = _property(grid.Δzᵃᵃᶠ, ℓz, topology(grid, 3), size(grid, 3), with_halos) - -@inline λspacings(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = λspacings(grid, ℓx; kwargs...) -@inline φspacings(grid::LLG, ℓx, ℓy, ℓz; kwargs...) = φspacings(grid, ℓy; kwargs...) - @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.Δλᶜᵃᵃ @@ -746,3 +685,10 @@ end @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) + +@inline λspacings(grid::LLG, ℓx) = λspacings(grid, ℓx, nothing, nothing) +@inline φspacings(grid::LLG, ℓy) = φspacings(grid, nothing, ℓy, nothing) diff --git a/src/Grids/nodes_and_spacings.jl b/src/Grids/nodes_and_spacings.jl index e348c0cb47..fd814d94e5 100644 --- a/src/Grids/nodes_and_spacings.jl +++ b/src/Grids/nodes_and_spacings.jl @@ -161,96 +161,36 @@ nodes(grid::AbstractGrid, (ℓx, ℓy, ℓz); reshape=false, with_halos=false) = ##### << Spacings >> ##### -# placeholders; see Oceananigans.Operators for x/y/zspacing definitions +# placeholders +# see Oceananigans/AbstractOperations/grid_metrics.jl for definitions function xspacing end function yspacing end function zspacing end +function λspacing end +function φspacing end -""" - xspacings(grid, ℓx, ℓy, ℓz; with_halos=true) - -Return the spacings over the interior nodes on `grid` in the ``x``-direction for the location `ℓx`, -`ℓy`, `ℓz`. For `Bounded` directions, `Face` nodes include the boundary points. - -```jldoctest xspacings -julia> using Oceananigans - -julia> grid = LatitudeLongitudeGrid(size=(8, 15, 10), longitude=(-20, 60), latitude=(-10, 50), z=(-100, 0)); - -julia> xspacings(grid, Center(), Face(), Center()) -16-element view(OffsetArray(::Vector{Float64}, -2:18), 1:16) with eltype Float64: - 1.0950562585518518e6 - 1.1058578920188267e6 - 1.1112718969963323e6 - 1.1112718969963323e6 - 1.1058578920188267e6 - 1.0950562585518518e6 - 1.0789196210678827e6 - 1.0575265956426917e6 - 1.0309814069457315e6 - 999413.38046802 - 962976.3124613502 - 921847.720658409 - 876227.979424229 - 826339.3435524226 - 772424.8654621692 - 714747.2110712599 -``` -""" -@inline xspacings(grid, ℓx, ℓy, ℓz; with_halos=true) = xspacings(grid, ℓx; with_halos) - -""" - yspacings(grid, ℓx, ℓy, ℓz; with_halos=true) - -Return the spacings over the interior nodes on `grid` in the ``y``-direction for the location `ℓx`, -`ℓy`, `ℓz`. For `Bounded` directions, `Face` nodes include the boundary points. - -```jldoctest yspacings -julia> using Oceananigans - -julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0)); - -julia> yspacings(grid, Center(), Center(), Center()) -222389.85328911748 -``` -""" -@inline yspacings(grid, ℓx, ℓy, ℓz; with_halos=true) = yspacings(grid, ℓy; with_halos) - -""" - zspacings(grid, ℓx, ℓy, ℓz; with_halos=true) - -Return the spacings over the interior nodes on `grid` in the ``z``-direction for the location `ℓx`, -`ℓy`, `ℓz`. For `Bounded` directions, `Face` nodes include the boundary points. - -```jldoctest zspacings -julia> using Oceananigans - -julia> grid = LatitudeLongitudeGrid(size=(20, 15, 10), longitude=(0, 20), latitude=(-15, 15), z=(-100, 0)); - -julia> zspacings(grid, Center(), Center(), Center()) -10.0 -``` -""" -@inline zspacings(grid, ℓx, ℓy, ℓz; with_halos=true) = zspacings(grid, ℓz; with_halos) +function xspacings end +function yspacings end +function zspacings end +function λspacings end +function φspacings end destantiate(::Face) = Face destantiate(::Center) = Center -spacing_function(::Val{:x}) = xspacing -spacing_function(::Val{:y}) = yspacing -spacing_function(::Val{:z}) = zspacing +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) - spacing = spacing_function(s) - LX, LY, LZ = map(destantiate, (ℓx, ℓy, ℓz)) - Δ = KernelFunctionOperation{LX, LY, LZ}(spacing, grid, ℓx, ℓy, ℓz) - - return minimum(Δ) + spacings = spacings_function(s) + return minimum(spacings(grid, ℓx, ℓy, ℓz)) end """ minimum_xspacing(grid, ℓx, ℓy, ℓz) - minimum_xspacing(grid) = minimum_xspacing(grid, Center(), Center(), Center()) Return the minimum spacing for `grid` in ``x`` direction at location `ℓx, ℓy, ℓz`. @@ -265,11 +205,11 @@ julia> minimum_xspacing(grid, Center(), Center(), Center()) 0.5 ``` """ -minimum_xspacing(grid, ℓx, ℓy, ℓz) = minimum_spacing(Val(:x), grid, ℓx, ℓy, ℓz) -minimum_xspacing(grid) = minimum_spacing(Val(:x), grid, Center(), Center(), Center()) +minimum_xspacing(grid, loc...) = minimum(xspacings(grid, loc...)) +minimum_xspacing(grid) = minimum(xspacings(grid)) + """ minimum_yspacing(grid, ℓx, ℓy, ℓz) - minimum_yspacing(grid) = minimum_yspacing(grid, Center(), Center(), Center()) Return the minimum spacing for `grid` in ``y`` direction at location `ℓx, ℓy, ℓz`. @@ -284,8 +224,8 @@ julia> minimum_yspacing(grid, Center(), Center(), Center()) 0.25 ``` """ -minimum_yspacing(grid, ℓx, ℓy, ℓz) = minimum_spacing(Val(:y), grid, ℓx, ℓy, ℓz) -minimum_yspacing(grid) = minimum_spacing(Val(:y), grid, Center(), Center(), Center()) +minimum_yspacing(grid, loc...) = minimum(yspacings(grid, loc...)) +minimum_yspacing(grid) = minimum(yspacings(grid)) """ minimum_zspacing(grid, ℓx, ℓy, ℓz) @@ -304,6 +244,5 @@ julia> minimum_zspacing(grid, Center(), Center(), Center()) 0.125 ``` """ -minimum_zspacing(grid, ℓx, ℓy, ℓz) = minimum_spacing(Val(:z), grid, ℓx, ℓy, ℓz) -minimum_zspacing(grid) = minimum_spacing(Val(:z), grid, Center(), Center(), Center()) - +minimum_zspacing(grid, loc...) = minimum(zspacings(grid, loc...)) +minimum_zspacing(grid) = minimum(zspacings(grid)) diff --git a/src/Grids/orthogonal_spherical_shell_grid.jl b/src/Grids/orthogonal_spherical_shell_grid.jl index 7d6edd4f5b..7a30c2b7a4 100644 --- a/src/Grids/orthogonal_spherical_shell_grid.jl +++ b/src/Grids/orthogonal_spherical_shell_grid.jl @@ -100,7 +100,7 @@ OrthogonalSphericalShellGrid(architecture, Nx, Ny, Nz, Hx, Hy, Hz, Lz, halo = (1, 1, 1), rotation = nothing) -Create a `OrthogonalSphericalShellGrid` that represents a section of a sphere after it has been +Create a `OrthogonalSphericalShellGrid` that represents a section of a sphere after it has been conformally mapped from the face of a cube. The cube's coordinates are `ξ` and `η` (which, by default, both take values in the range ``[-1, 1]``. @@ -176,7 +176,7 @@ function conformal_cubed_sphere_panel(architecture::AbstractArchitecture = CPU() halo = (1, 1, 1), rotation = nothing) - if architecture == GPU() && !has_cuda() + if architecture == GPU() && !has_cuda() throw(ArgumentError("Cannot create a GPU grid. No CUDA-enabled GPU was detected!")) end @@ -364,9 +364,9 @@ function conformal_cubed_sphere_panel(architecture::AbstractArchitecture = CPU() j = 1 Δyᶠᶠᵃ[i, j] = 2haversine((λᶠᶜᵃ[i, j], φᶠᶜᵃ[i, j]), (λᶠᶠᵃ[i, j ], φᶠᶠᵃ[i, j ]), radius) end - + for i in 1:Nξ+1 - j = Nη+1 + j = Nη+1 Δyᶠᶠᵃ[i, j] = 2haversine((λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), (λᶠᶜᵃ[i, j-1], φᶠᶜᵃ[i, j-1]), radius) end end @@ -1201,31 +1201,9 @@ rname(::OSSG) = :z ##### Grid spacings in x, y, z (in meters) ##### -@inline xspacings(grid::OSSG, ℓx::Center, ℓy::Center; with_halos=false) = - with_halos ? grid.Δxᶜᶜᵃ : view(grid.Δxᶜᶜᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline xspacings(grid::OSSG, ℓx::Face , ℓy::Center; with_halos=false) = - with_halos ? grid.Δxᶠᶜᵃ : view(grid.Δxᶠᶜᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline xspacings(grid::OSSG, ℓx::Center, ℓy::Face ; with_halos=false) = - with_halos ? grid.Δxᶜᶠᵃ : view(grid.Δxᶜᶠᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline xspacings(grid::OSSG, ℓx::Face , ℓy::Face ; with_halos=false) = - with_halos ? grid.Δxᶠᶠᵃ : view(grid.Δxᶠᶠᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) - -@inline yspacings(grid::OSSG, ℓx::Center, ℓy::Center; with_halos=false) = - with_halos ? grid.Δyᶜᶜᵃ : view(grid.Δyᶜᶜᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline yspacings(grid::OSSG, ℓx::Face , ℓy::Center; with_halos=false) = - with_halos ? grid.Δyᶠᶜᵃ : view(grid.Δyᶠᶜᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline yspacings(grid::OSSG, ℓx::Center, ℓy::Face ; with_halos=false) = - with_halos ? grid.Δyᶜᶠᵃ : view(grid.Δyᶜᶠᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) -@inline yspacings(grid::OSSG, ℓx::Face , ℓy::Face ; with_halos=false) = - with_halos ? grid.Δyᶠᶠᵃ : view(grid.Δyᶠᶠᵃ, interior_indices(ℓx, topology(grid, 1)(), grid.Nx), interior_indices(ℓy, topology(grid, 2)(), grid.Ny)) - -@inline zspacings(grid::OSSG, ℓz::Center; with_halos=false) = with_halos ? grid.Δzᵃᵃᶜ : - view(grid.Δzᵃᵃᶜ, interior_indices(ℓz, topology(grid, 3)(), grid.Nz)) -@inline zspacings(grid::ZRegOSSG, ℓz::Center; with_halos=false) = grid.Δzᵃᵃᶜ -@inline zspacings(grid::OSSG, ℓz::Face; with_halos=false) = with_halos ? grid.Δzᵃᵃᶠ : - view(grid.Δzᵃᵃᶠ, interior_indices(ℓz, topology(grid, 3)(), grid.Nz)) -@inline zspacings(grid::ZRegOSSG, ℓz::Face; with_halos=false) = grid.Δzᵃᵃᶠ - -@inline xspacings(grid::OSSG, ℓx, ℓy, ℓz; with_halos=false) = xspacings(grid, ℓx, ℓy; with_halos) -@inline yspacings(grid::OSSG, ℓx, ℓy, ℓz; with_halos=false) = yspacings(grid, ℓx, ℓy; with_halos) -@inline zspacings(grid::OSSG, ℓx, ℓy, ℓz; with_halos=false) = zspacings(grid, ℓz; with_halos) +@inline xspacings(grid::OSSG, ℓx, ℓy) = xspacings(grid, ℓx, ℓy, nothing) +@inline yspacings(grid::OSSG, ℓx, ℓy) = yspacings(grid, ℓx, ℓy, nothing) +@inline zspacings(grid::OSSG, ℓz) = zspacings(grid, nothing, nothing, ℓz) + +@inline λspacings(grid::OSSG, ℓx, ℓy) = λspacings(grid, ℓx, ℓy, nothing) +@inline φspacings(grid::OSSG, ℓx, ℓy) = φspacings(grid, ℓx, ℓy, nothing) diff --git a/src/Grids/rectilinear_grid.jl b/src/Grids/rectilinear_grid.jl index 25a30e7cbc..f148625f15 100644 --- a/src/Grids/rectilinear_grid.jl +++ b/src/Grids/rectilinear_grid.jl @@ -518,19 +518,12 @@ const C = Center @inline ηnodes(grid::RG, ℓx, ℓy, ℓz; kwargs...) = ynodes(grid, ℓy; kwargs...) @inline rnodes(grid::RG, ℓx, ℓy, ℓz; kwargs...) = znodes(grid, ℓz; kwargs...) +@inline isrectilinear(::RG) = true + ##### -##### Grid spacings +##### Grid-specific grid spacings ##### -@inline xspacings(grid::RG, ℓx::C; with_halos=false) = _property(grid.Δxᶜᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline xspacings(grid::RG, ℓx::F; with_halos=false) = _property(grid.Δxᶠᵃᵃ, ℓx, topology(grid, 1), size(grid, 1), with_halos) -@inline yspacings(grid::RG, ℓy::C; with_halos=false) = _property(grid.Δyᵃᶜᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) -@inline yspacings(grid::RG, ℓy::F; with_halos=false) = _property(grid.Δyᵃᶠᵃ, ℓy, topology(grid, 2), size(grid, 2), with_halos) -@inline zspacings(grid::RG, ℓz::C; with_halos=false) = _property(grid.Δzᵃᵃᶜ, ℓz, topology(grid, 3), size(grid, 3), with_halos) -@inline zspacings(grid::RG, ℓz::F; with_halos=false) = _property(grid.Δzᵃᵃᶠ, ℓz, topology(grid, 3), size(grid, 3), with_halos) - -@inline xspacings(grid::RG, ℓx, ℓy, ℓz; kwargs...) = xspacings(grid, ℓx; kwargs...) -@inline yspacings(grid::RG, ℓx, ℓy, ℓz; kwargs...) = yspacings(grid, ℓy; kwargs...) -@inline zspacings(grid::RG, ℓx, ℓy, ℓz; kwargs...) = zspacings(grid, ℓz; kwargs...) - -@inline isrectilinear(::RG) = true +@inline xspacings(grid::RG, ℓx) = xspacings(grid, ℓx, nothing, nothing) +@inline yspacings(grid::RG, ℓy) = yspacings(grid, nothing, ℓy, nothing) +@inline zspacings(grid::RG, ℓz) = zspacings(grid, nothing, nothing, ℓz) diff --git a/src/ImmersedBoundaries/ImmersedBoundaries.jl b/src/ImmersedBoundaries/ImmersedBoundaries.jl index 7047948e31..2b45f4af9a 100644 --- a/src/ImmersedBoundaries/ImmersedBoundaries.jl +++ b/src/ImmersedBoundaries/ImmersedBoundaries.jl @@ -31,216 +31,9 @@ import Oceananigans.Architectures: on_architecture import Oceananigans.Fields: fractional_x_index, fractional_y_index, fractional_z_index -""" - abstract type AbstractImmersedBoundary - -Abstract supertype for immersed boundary grids. -""" -abstract type AbstractImmersedBoundary end - -##### -##### ImmersedBoundaryGrid -##### - -struct ImmersedBoundaryGrid{FT, TX, TY, TZ, G, I, M, S, Arch} <: AbstractGrid{FT, TX, TY, TZ, Arch} - architecture :: Arch - underlying_grid :: G - immersed_boundary :: I - interior_active_cells :: M - active_z_columns :: S - - # Internal interface - function ImmersedBoundaryGrid{TX, TY, TZ}(grid::G, ib::I, mi::M, ms::S) where {TX, TY, TZ, G <: AbstractUnderlyingGrid, I, M, S} - FT = eltype(grid) - arch = architecture(grid) - Arch = typeof(arch) - return new{FT, TX, TY, TZ, G, I, M, S, Arch}(arch, grid, ib, mi, ms) - end - - # Constructor with no active map - function ImmersedBoundaryGrid{TX, TY, TZ}(grid::G, ib::I) where {TX, TY, TZ, G <: AbstractUnderlyingGrid, I} - FT = eltype(grid) - arch = architecture(grid) - Arch = typeof(arch) - return new{FT, TX, TY, TZ, G, I, Nothing, Nothing, Arch}(arch, grid, ib, nothing, nothing) - end -end - -const IBG = ImmersedBoundaryGrid - -@inline Base.getproperty(ibg::IBG, property::Symbol) = get_ibg_property(ibg, Val(property)) -@inline get_ibg_property(ibg::IBG, ::Val{property}) where property = getfield(getfield(ibg, :underlying_grid), property) -@inline get_ibg_property(ibg::IBG, ::Val{:immersed_boundary}) = getfield(ibg, :immersed_boundary) -@inline get_ibg_property(ibg::IBG, ::Val{:underlying_grid}) = getfield(ibg, :underlying_grid) -@inline get_ibg_property(ibg::IBG, ::Val{:interior_active_cells}) = getfield(ibg, :interior_active_cells) -@inline get_ibg_property(ibg::IBG, ::Val{:active_z_columns}) = getfield(ibg, :active_z_columns) - -@inline architecture(ibg::IBG) = architecture(ibg.underlying_grid) - -@inline x_domain(ibg::IBG) = x_domain(ibg.underlying_grid) -@inline y_domain(ibg::IBG) = y_domain(ibg.underlying_grid) -@inline z_domain(ibg::IBG) = z_domain(ibg.underlying_grid) - -Adapt.adapt_structure(to, ibg::IBG{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = - ImmersedBoundaryGrid{TX, TY, TZ}(adapt(to, ibg.underlying_grid), - adapt(to, ibg.immersed_boundary), - nothing, - nothing) - -with_halo(halo, ibg::ImmersedBoundaryGrid) = - ImmersedBoundaryGrid(with_halo(halo, ibg.underlying_grid), ibg.immersed_boundary) - -# ImmersedBoundaryGrids require an extra halo point to check the "inactivity" of a `Face` node at N + H -# (which requires checking `Center` nodes at N + H and N + H + 1) -inflate_halo_size_one_dimension(req_H, old_H, _, ::IBG) = max(req_H + 1, old_H) -inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, ::IBG) = 0 - -function Base.summary(grid::ImmersedBoundaryGrid) - FT = eltype(grid) - TX, TY, TZ = topology(grid) - - return string(size_summary(size(grid)), - " ImmersedBoundaryGrid{$FT, $TX, $TY, $TZ} on ", summary(architecture(grid)), - " with ", size_summary(halo_size(grid)), " halo") -end - -function show(io::IO, ibg::ImmersedBoundaryGrid) - print(io, summary(ibg), ":", "\n", - "├── immersed_boundary: ", summary(ibg.immersed_boundary), "\n", - "├── underlying_grid: ", summary(ibg.underlying_grid), "\n") - - return show(io, ibg.underlying_grid, false) -end - -##### -##### Interface for immersed_boundary -##### - -""" - immersed_cell(i, j, k, grid) - -Return true if a `cell` is "completely" immersed, and thus -is not part of the prognostic state. -""" -@inline immersed_cell(i, j, k, grid) = false - -# Unpack to make defining new immersed boundaries more convenient -@inline immersed_cell(i, j, k, grid::ImmersedBoundaryGrid) = - immersed_cell(i, j, k, grid.underlying_grid, grid.immersed_boundary) - -""" - inactive_cell(i, j, k, grid::ImmersedBoundaryGrid) - -Return `true` if the tracer cell at `i, j, k` either (i) lies outside the `Bounded` domain -or (ii) lies within the immersed region of `ImmersedBoundaryGrid`. - -Example -======= - -Consider the configuration - -``` - Immersed Fluid - =========== ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ - - c c - i-1 i - - | ========= | | - × === ∘ === × ∘ × - | ========= | | - -i-1 i - f f f -``` - -We then have - -* `inactive_node(i, 1, 1, grid, f, c, c) = false` - -As well as - -* `inactive_node(i, 1, 1, grid, c, c, c) = false` -* `inactive_node(i-1, 1, 1, grid, c, c, c) = true` -* `inactive_node(i-1, 1, 1, grid, f, c, c) = true` -""" -@inline inactive_cell(i, j, k, ibg::IBG) = immersed_cell(i, j, k, ibg) | inactive_cell(i, j, k, ibg.underlying_grid) - -# Isolate periphery of the immersed boundary -@inline immersed_peripheral_node(i, j, k, ibg::IBG, LX, LY, LZ) = peripheral_node(i, j, k, ibg, LX, LY, LZ) & - !peripheral_node(i, j, k, ibg.underlying_grid, LX, LY, LZ) - -@inline immersed_inactive_node(i, j, k, ibg::IBG, LX, LY, LZ) = inactive_node(i, j, k, ibg, LX, LY, LZ) & - !inactive_node(i, j, k, ibg.underlying_grid, LX, LY, LZ) - -##### -##### Utilities -##### - -const c = Center() -const f = Face() - -@inline Base.zero(ibg::IBG) = zero(ibg.underlying_grid) - -@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) - -function on_architecture(arch, ibg::IBG) - underlying_grid = on_architecture(arch, ibg.underlying_grid) - immersed_boundary = on_architecture(arch, ibg.immersed_boundary) - return ImmersedBoundaryGrid(underlying_grid, immersed_boundary) -end - -isrectilinear(ibg::IBG) = isrectilinear(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) - +include("immersed_boundary_grid.jl") +include("immersed_boundary_interface.jl") +include("immersed_boundary_nodes.jl") include("active_cells_map.jl") include("immersed_grid_metrics.jl") include("abstract_grid_fitted_boundary.jl") diff --git a/src/ImmersedBoundaries/immersed_boundary_grid.jl b/src/ImmersedBoundaries/immersed_boundary_grid.jl new file mode 100644 index 0000000000..664dfb8ecf --- /dev/null +++ b/src/ImmersedBoundaries/immersed_boundary_grid.jl @@ -0,0 +1,90 @@ +""" + abstract type AbstractImmersedBoundary + +Abstract supertype for immersed boundary grids. +""" +abstract type AbstractImmersedBoundary end + +struct ImmersedBoundaryGrid{FT, TX, TY, TZ, G, I, M, S, Arch} <: AbstractGrid{FT, TX, TY, TZ, Arch} + architecture :: Arch + underlying_grid :: G + immersed_boundary :: I + interior_active_cells :: M + active_z_columns :: S + + # Internal interface + function ImmersedBoundaryGrid{TX, TY, TZ}(grid::G, ib::I, mi::M, ms::S) where {TX, TY, TZ, G <: AbstractUnderlyingGrid, I, M, S} + FT = eltype(grid) + arch = architecture(grid) + Arch = typeof(arch) + return new{FT, TX, TY, TZ, G, I, M, S, Arch}(arch, grid, ib, mi, ms) + end + + # Constructor with no active map + function ImmersedBoundaryGrid{TX, TY, TZ}(grid::G, ib::I) where {TX, TY, TZ, G <: AbstractUnderlyingGrid, I} + FT = eltype(grid) + arch = architecture(grid) + Arch = typeof(arch) + return new{FT, TX, TY, TZ, G, I, Nothing, Nothing, Arch}(arch, grid, ib, nothing, nothing) + end +end + +const IBG = ImmersedBoundaryGrid + +@inline Base.getproperty(ibg::IBG, property::Symbol) = get_ibg_property(ibg, Val(property)) +@inline get_ibg_property(ibg::IBG, ::Val{property}) where property = getfield(getfield(ibg, :underlying_grid), property) +@inline get_ibg_property(ibg::IBG, ::Val{:immersed_boundary}) = getfield(ibg, :immersed_boundary) +@inline get_ibg_property(ibg::IBG, ::Val{:underlying_grid}) = getfield(ibg, :underlying_grid) +@inline get_ibg_property(ibg::IBG, ::Val{:interior_active_cells}) = getfield(ibg, :interior_active_cells) +@inline get_ibg_property(ibg::IBG, ::Val{:active_z_columns}) = getfield(ibg, :active_z_columns) + +@inline architecture(ibg::IBG) = architecture(ibg.underlying_grid) + +@inline x_domain(ibg::IBG) = x_domain(ibg.underlying_grid) +@inline y_domain(ibg::IBG) = y_domain(ibg.underlying_grid) +@inline z_domain(ibg::IBG) = z_domain(ibg.underlying_grid) + +Adapt.adapt_structure(to, ibg::IBG{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = + ImmersedBoundaryGrid{TX, TY, TZ}(adapt(to, ibg.underlying_grid), + adapt(to, ibg.immersed_boundary), + nothing, + nothing) + +with_halo(halo, ibg::ImmersedBoundaryGrid) = + ImmersedBoundaryGrid(with_halo(halo, ibg.underlying_grid), ibg.immersed_boundary) + +# ImmersedBoundaryGrids require an extra halo point to check the "inactivity" of a `Face` node at N + H +# (which requires checking `Center` nodes at N + H and N + H + 1) +inflate_halo_size_one_dimension(req_H, old_H, _, ::IBG) = max(req_H + 1, old_H) +inflate_halo_size_one_dimension(req_H, old_H, ::Type{Flat}, ::IBG) = 0 + +# Defining the bottom +@inline z_bottom(i, j, grid) = znode(i, j, 1, grid, c, c, f) +@inline z_bottom(i, j, ibg::IBG) = error("The function `bottom` has not been defined for $(summary(ibg))!") + +function Base.summary(grid::ImmersedBoundaryGrid) + FT = eltype(grid) + TX, TY, TZ = topology(grid) + + return string(size_summary(size(grid)), + " ImmersedBoundaryGrid{$FT, $TX, $TY, $TZ} on ", summary(architecture(grid)), + " with ", size_summary(halo_size(grid)), " halo") +end + +function show(io::IO, ibg::ImmersedBoundaryGrid) + print(io, summary(ibg), ":", "\n", + "├── immersed_boundary: ", summary(ibg.immersed_boundary), "\n", + "├── underlying_grid: ", summary(ibg.underlying_grid), "\n") + + return show(io, ibg.underlying_grid, false) +end + +@inline Base.zero(ibg::IBG) = zero(ibg.underlying_grid) + +function on_architecture(arch, ibg::IBG) + underlying_grid = on_architecture(arch, ibg.underlying_grid) + immersed_boundary = on_architecture(arch, ibg.immersed_boundary) + return ImmersedBoundaryGrid(underlying_grid, immersed_boundary) +end + +isrectilinear(ibg::IBG) = isrectilinear(ibg.underlying_grid) diff --git a/src/ImmersedBoundaries/immersed_boundary_interface.jl b/src/ImmersedBoundaries/immersed_boundary_interface.jl new file mode 100644 index 0000000000..8f6e62153b --- /dev/null +++ b/src/ImmersedBoundaries/immersed_boundary_interface.jl @@ -0,0 +1,56 @@ +""" + immersed_cell(i, j, k, grid) + +Return true if a `cell` is "completely" immersed, and thus +is not part of the prognostic state. +""" +@inline immersed_cell(i, j, k, grid) = false + +# Unpack to make defining new immersed boundaries more convenient +@inline immersed_cell(i, j, k, grid::ImmersedBoundaryGrid) = + immersed_cell(i, j, k, grid.underlying_grid, grid.immersed_boundary) + +""" + inactive_cell(i, j, k, grid::ImmersedBoundaryGrid) + +Return `true` if the tracer cell at `i, j, k` either (i) lies outside the `Bounded` domain +or (ii) lies within the immersed region of `ImmersedBoundaryGrid`. + +Example +======= + +Consider the configuration + +``` + Immersed Fluid + =========== ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ + + c c + i-1 i + + | ========= | | + × === ∘ === × ∘ × + | ========= | | + +i-1 i + f f f +``` + +We then have + +* `inactive_node(i, 1, 1, grid, f, c, c) = false` + +As well as + +* `inactive_node(i, 1, 1, grid, c, c, c) = false` +* `inactive_node(i-1, 1, 1, grid, c, c, c) = true` +* `inactive_node(i-1, 1, 1, grid, f, c, c) = true` +""" +@inline inactive_cell(i, j, k, ibg::IBG) = immersed_cell(i, j, k, ibg) | inactive_cell(i, j, k, ibg.underlying_grid) + +# Isolate periphery of the immersed boundary +@inline immersed_peripheral_node(i, j, k, ibg::IBG, LX, LY, LZ) = peripheral_node(i, j, k, ibg, LX, LY, LZ) & + !peripheral_node(i, j, k, ibg.underlying_grid, LX, LY, LZ) + +@inline immersed_inactive_node(i, j, k, ibg::IBG, LX, LY, LZ) = inactive_node(i, j, k, ibg, LX, LY, LZ) & + !inactive_node(i, j, k, ibg.underlying_grid, LX, LY, LZ) diff --git a/src/ImmersedBoundaries/immersed_boundary_nodes.jl b/src/ImmersedBoundaries/immersed_boundary_nodes.jl new file mode 100644 index 0000000000..62f39894f8 --- /dev/null +++ b/src/ImmersedBoundaries/immersed_boundary_nodes.jl @@ -0,0 +1,69 @@ +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) + +xspacings(ibg::ImmersedBoundaryGrid, ℓ...) = xspacings(ibg.underlying_grid, ℓ...) +yspacings(ibg::ImmersedBoundaryGrid, ℓ...) = yspacings(ibg.underlying_grid, ℓ...) +zspacings(ibg::ImmersedBoundaryGrid, ℓ...) = zspacings(ibg.underlying_grid, ℓ...) diff --git a/src/ImmersedBoundaries/immersed_grid_metrics.jl b/src/ImmersedBoundaries/immersed_grid_metrics.jl index a7dfbcbe10..466ad4fa2f 100644 --- a/src/ImmersedBoundaries/immersed_grid_metrics.jl +++ b/src/ImmersedBoundaries/immersed_grid_metrics.jl @@ -1,10 +1,7 @@ using Oceananigans.AbstractOperations: GridMetricOperation import Oceananigans.Grids: coordinates - -const c = Center() -const f = Face() -const IBG = ImmersedBoundaryGrid +import Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ # Grid metrics for ImmersedBoundaryGrid # @@ -12,10 +9,10 @@ const IBG = ImmersedBoundaryGrid # # 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 @@ -30,10 +27,7 @@ for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) end end -@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) - coordinates(grid::IBG) = coordinates(grid.underlying_grid) -xspacings(X, grid::IBG) = xspacings(X, grid.underlying_grid) -yspacings(Y, grid::IBG) = yspacings(Y, grid.underlying_grid) -zspacings(Z, grid::IBG) = zspacings(Z, 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) diff --git a/src/ImmersedBoundaries/partial_cell_bottom.jl b/src/ImmersedBoundaries/partial_cell_bottom.jl index 9e7b552df5..36cc32487b 100644 --- a/src/ImmersedBoundaries/partial_cell_bottom.jl +++ b/src/ImmersedBoundaries/partial_cell_bottom.jl @@ -91,7 +91,7 @@ end bottom_cell = (z⁻ ≤ zb) & (z⁺ ≥ zb) capped_zb = min(z⁺ - ϵ * Δz, zb) - # If the size of the bottom cell is less than ϵ Δz, + # 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) end @@ -108,20 +108,20 @@ function on_architecture(arch, ib::PartialCellBottom{<:Field}) end Adapt.adapt_structure(to, ib::PartialCellBottom) = PartialCellBottom(adapt(to, ib.bottom_height), - ib.minimum_fractional_cell_height) + ib.minimum_fractional_cell_height) on_architecture(to, ib::PartialCellBottom) = PartialCellBottom(on_architecture(to, ib.bottom_height), - on_architecture(to, ib.minimum_fractional_cell_height)) + on_architecture(to, ib.minimum_fractional_cell_height)) """ immersed underlying --x-- --x-- - - + + ∘ ↑ ∘ k+1 | - | + | k+1 --x-- | k+1 --x-- ↑ <- node z ∘ ↓ | zb ⋅⋅x⋅⋅ | @@ -130,7 +130,7 @@ on_architecture(to, ib::PartialCellBottom) = PartialCellBottom(on_architecture(t | | k --x-- ↓ - + Criterion is zb ≥ z - ϵ Δz """ @@ -185,9 +185,9 @@ end @inline Δzᶠᶜᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i-1, j, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶜᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶜ(i, j-1, k, ibg), Δzᶜᶜᶜ(i, j, k, ibg)) @inline Δzᶠᶠᶜ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶜ(i, j-1, k, ibg), Δzᶠᶜᶜ(i, j, k, ibg)) - + @inline Δzᶠᶜᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶠ(i-1, j, k, ibg), Δzᶜᶜᶠ(i, j, k, ibg)) -@inline Δzᶜᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶠ(i, j-1, k, ibg), Δzᶜᶜᶠ(i, j, k, ibg)) +@inline Δzᶜᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶜᶜᶠ(i, j-1, k, ibg), Δzᶜᶜᶠ(i, j, k, ibg)) @inline Δzᶠᶠᶠ(i, j, k, ibg::PCBIBG) = min(Δzᶠᶜᶠ(i, j-1, k, ibg), Δzᶠᶜᶠ(i, j, k, ibg)) # Make sure Δz works for horizontally-Flat topologies. diff --git a/src/MultiRegion/multi_region_grid.jl b/src/MultiRegion/multi_region_grid.jl index 3e6ddb1cd7..264b6a29be 100644 --- a/src/MultiRegion/multi_region_grid.jl +++ b/src/MultiRegion/multi_region_grid.jl @@ -3,9 +3,10 @@ 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: minimum_xspacing, minimum_yspacing, minimum_zspacing import Oceananigans.Models.HydrostaticFreeSurfaceModels: default_free_surface import Oceananigans.DistributedComputations: reconstruct_global_grid -import Oceananigans.Grids: minimum_spacing, destantiate struct MultiRegionGrid{FT, TX, TY, TZ, P, C, G, D, Arch} <: AbstractMultiRegionGrid{FT, TX, TY, TZ, Arch} architecture :: Arch @@ -19,7 +20,7 @@ struct MultiRegionGrid{FT, TX, TY, TZ, P, C, G, D, Arch} <: AbstractMultiRegionG new{FT, TX, TY, TZ, P, C, G, D, A}(arch, partition, connectivity, region_grids, devices) end -const ImmersedMultiRegionGrid = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:MultiRegionGrid} +const ImmersedMultiRegionGrid = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:MultiRegionGrid} const MultiRegionGrids = Union{MultiRegionGrid, ImmersedMultiRegionGrid} @@ -41,6 +42,15 @@ 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)) + +minimum_yspacing(grid::MultiRegionGrid, ℓx, ℓy, ℓz) = + minimum(minimum_yspacing(grid[r], ℓx, ℓy, ℓz) for r in 1:number_of_regions(grid)) + +minimum_zspacing(grid::MultiRegionGrid, ℓx, ℓy, ℓz) = + minimum(minimum_zspacing(grid[r], ℓx, ℓy, ℓz) for r in 1:number_of_regions(grid)) + @inline getdevice(mrg::ImmersedMultiRegionGrid, i) = getdevice(mrg.underlying_grid.region_grids, i) @inline switch_device!(mrg::ImmersedMultiRegionGrid, i) = switch_device!(getdevice(mrg.underlying_grid, i)) @inline devices(mrg::ImmersedMultiRegionGrid) = devices(mrg.underlying_grid.region_grids) @@ -68,7 +78,7 @@ Positional Arguments Keyword Arguments ================= -- `partition`: the partitioning required. The implemented partitioning are `XPartition` +- `partition`: the partitioning required. The implemented partitioning are `XPartition` (division along the ``x`` direction) and `YPartition` (division along the ``y`` direction). @@ -97,7 +107,7 @@ function MultiRegionGrid(global_grid; partition = XPartition(2), validate = true) @warn "MultiRegion functionalities are experimental: help the development by reporting bugs or non-implemented features!" - + if length(partition) == 1 return global_grid end @@ -172,7 +182,7 @@ end ##### `ImmersedMultiRegionGrid` functionalities ##### -function reconstruct_global_grid(mrg::ImmersedMultiRegionGrid) +function reconstruct_global_grid(mrg::ImmersedMultiRegionGrid) global_grid = reconstruct_global_grid(mrg.underlying_grid) global_immersed_boundary = reconstruct_global_immersed_boundary(mrg.immersed_boundary) global_immersed_boundary = on_architecture(architecture(mrg), global_immersed_boundary) @@ -203,14 +213,14 @@ end # Fallback! multi_region_object_from_array(a::AbstractArray, grid) = on_architecture(architecture(grid), a) -#### +#### #### Utilities for MultiRegionGrid #### new_data(FT::DataType, mrg::MultiRegionGrids, args...) = construct_regionally(new_data, FT, mrg, args...) # This is kind of annoying but it is necessary to have compatible MultiRegion and Distributed -function with_halo(new_halo, mrg::MultiRegionGrid) +function with_halo(new_halo, mrg::MultiRegionGrid) devices = mrg.devices partition = mrg.partition cpu_mrg = on_architecture(CPU(), mrg) @@ -224,11 +234,11 @@ end function on_architecture(::CPU, mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} new_grids = construct_regionally(on_architecture, CPU(), mrg) - devices = Tuple(CPU() for i in 1:length(mrg)) + devices = Tuple(CPU() for i in 1:length(mrg)) return MultiRegionGrid{FT, TX, TY, TZ}(CPU(), mrg.partition, mrg.connectivity, new_grids, devices) end -Base.summary(mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = +Base.summary(mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = "MultiRegionGrid{$FT, $TX, $TY, $TZ} with $(summary(mrg.partition)) on $(string(typeof(mrg.region_grids[1]).name.wrapper))" Base.show(io::IO, mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = @@ -243,13 +253,13 @@ function Base.:(==)(mrg₁::MultiRegionGrid, mrg₂::MultiRegionGrid) vals = construct_regionally(Base.:(==), mrg₁, mrg₂) return all(vals.regional_objects) end - + #### #### This works only for homogenous partitioning #### -size(mrg::MultiRegionGrids) = size(getregion(mrg, 1)) -halo_size(mrg::MultiRegionGrids) = halo_size(getregion(mrg, 1)) +size(mrg::MultiRegionGrids) = size(getregion(mrg, 1)) +halo_size(mrg::MultiRegionGrids) = halo_size(getregion(mrg, 1)) #### #### Get property for `MultiRegionGrid` (gets the properties of region 1) diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl index d7d6234fc6..39701328c5 100644 --- a/src/Oceananigans.jl +++ b/src/Oceananigans.jl @@ -6,14 +6,14 @@ module Oceananigans export # Architectures - CPU, GPU, + CPU, GPU, # Grids Center, Face, - Periodic, Bounded, Flat, + Periodic, Bounded, Flat, RectilinearGrid, LatitudeLongitudeGrid, OrthogonalSphericalShellGrid, nodes, xnodes, ynodes, znodes, λnodes, φnodes, - xspacings, yspacings, zspacings, + xspacings, yspacings, zspacings, λspacings, φspacings, minimum_xspacing, minimum_yspacing, minimum_zspacing, # Immersed boundaries @@ -247,4 +247,3 @@ using .AbstractOperations using .MultiRegion end # module - diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index 0a80bc183e..70c295681f 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -30,6 +30,8 @@ export δyᶠᶠᶠ, δyᶠᶠᶜ, δyᶠᶜᶠ, δyᶠᶜᶜ, δyᶜᶠᶠ, δy export δzᶠᶠᶠ, δzᶠᶠᶜ, δzᶠᶜᶠ, δzᶠᶜᶜ, δzᶜᶠᶠ, δzᶜᶠᶜ, δzᶜᶜᶠ, δzᶜᶜᶜ # Derivatives +export ∂xᶜᵃᵃ, ∂xᶠᵃᵃ, ∂yᵃᶜᵃ, ∂yᵃᶠᵃ, ∂zᵃᵃᶜ, ∂zᵃᵃᶠ + export ∂xᶠᶠᶠ, ∂xᶠᶠᶜ, ∂xᶠᶜᶠ, ∂xᶠᶜᶜ, ∂xᶜᶠᶠ, ∂xᶜᶠᶜ, ∂xᶜᶜᶠ, ∂xᶜᶜᶜ export ∂yᶠᶠᶠ, ∂yᶠᶠᶜ, ∂yᶠᶜᶠ, ∂yᶠᶜᶜ, ∂yᶜᶠᶠ, ∂yᶜᶠᶜ, ∂yᶜᶜᶠ, ∂yᶜᶜᶜ export ∂zᶠᶠᶠ, ∂zᶠᶠᶜ, ∂zᶠᶜᶠ, ∂zᶠᶜᶜ, ∂zᶜᶠᶠ, ∂zᶜᶠᶜ, ∂zᶜᶜᶠ, ∂zᶜᶜᶜ @@ -75,12 +77,30 @@ import Oceananigans.Grids: xspacing, yspacing, zspacing ##### Convenient aliases ##### -const AG = AbstractGrid +const AG = AbstractGrid const Δx = xspacing const Δy = yspacing const Δz = zspacing +const RG = RectilinearGrid +const RGX = XRegularRG +const RGY = YRegularRG +const RGZ = ZRegularRG + +const OSSG = OrthogonalSphericalShellGrid +const OSSGZ = ZRegOrthogonalSphericalShellGrid + +const LLG = LatitudeLongitudeGrid +const LLGX = XRegularLLG +const LLGY = YRegularLLG +const LLGZ = ZRegularLLG + +# On the fly calculations of metrics +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} + include("difference_operators.jl") include("interpolation_operators.jl") include("interpolation_utils.jl") diff --git a/src/Operators/derivative_operators.jl b/src/Operators/derivative_operators.jl index 67f3287df3..a78dfef074 100644 --- a/src/Operators/derivative_operators.jl +++ b/src/Operators/derivative_operators.jl @@ -2,8 +2,8 @@ ##### First derivative operators ##### -for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) - +for LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) + x_derivative = Symbol(:∂x, LX, LY, LZ) x_spacing = Symbol(:Δx, LX, LY, LZ) x_difference = Symbol(:δx, LX, LY, LZ) @@ -35,7 +35,7 @@ end ##### Second, Third, and Fourth derivatives ##### -@inline insert_symbol(dir, L, L1, L2) = +@inline insert_symbol(dir, L, L1, L2) = dir == :x ? (L, L1, L2) : dir == :y ? @@ -43,7 +43,7 @@ end (L1, L2, L) -for dir in (:x, :y, :z), L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ) +for dir in (:x, :y, :z), L1 in (:ᶜ, :ᶠ, :ᵃ), L2 in (:ᶜ, :ᶠ, :ᵃ) first_order_face = Symbol(:∂, dir, insert_symbol(dir, :ᶠ, L1, L2)...) second_order_face = Symbol(:∂², dir, insert_symbol(dir, :ᶠ, L1, L2)...) @@ -78,8 +78,8 @@ end ##### Operators of the form A*∂(q) where A is an area and q is some quantity. ##### -for dir in (:x, :y, :z), LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ) - +for dir in (:x, :y, :z), LX in (:ᶜ, :ᶠ, :ᵃ), LY in (:ᶜ, :ᶠ, :ᵃ), LZ in (:ᶜ, :ᶠ, :ᵃ) + operator = Symbol(:A, dir, :_∂, dir, LX, LY, LZ) area = Symbol(:A, dir, LX, LY, LZ) derivative = Symbol(:∂, dir, LX, LY, LZ) diff --git a/src/Operators/interpolation_utils.jl b/src/Operators/interpolation_utils.jl index 748eacf242..a2ff28e4d5 100644 --- a/src/Operators/interpolation_utils.jl +++ b/src/Operators/interpolation_utils.jl @@ -8,10 +8,10 @@ interpolation_code(from, to) = interpolation_code(to) interpolation_code(::Type{Face}) = :ᶠ interpolation_code(::Type{Center}) = :ᶜ -interpolation_code(::Type{Nothing}) = :ᶜ +interpolation_code(::Type{Nothing}) = :ᵃ interpolation_code(::Face) = :ᶠ interpolation_code(::Center) = :ᶜ -interpolation_code(::Nothing) = :ᶜ +interpolation_code(::Nothing) = :ᵃ # Intercept non-interpolations interpolation_code(from::L, to::L) where L = :ᵃ @@ -41,7 +41,7 @@ for i = 1:number_of_identities @inline $identity(i, j, k, grid, F::TF, args...) where TF<:Function = F(i, j, k, grid, args...) end end - + torus(x, lower, upper) = lower + rem(x - lower, upper - lower, RoundDown) identify_an_identity(number) = Symbol(:identity, torus(number, 1, number_of_identities)) identity_counter = 0 @@ -118,10 +118,10 @@ for LX in (:Center, :Face), LY in (:Center, :Face), LZ in (:Center, :Face) to = (eval(IX), eval(IY), eval(IZ)) interp_func = Symbol(interpolation_operator(from, to)) @eval begin - @inline ℑxyz(i, j, k, grid, from::F, to::T, c) where {F<:Tuple{<:$LX, <:$LY, <:$LZ}, T<:Tuple{<:$IX, <:$IY, <:$IZ}} = + @inline ℑxyz(i, j, k, grid, from::F, to::T, c) where {F<:Tuple{<:$LX, <:$LY, <:$LZ}, T<:Tuple{<:$IX, <:$IY, <:$IZ}} = $interp_func(i, j, k, grid, c) - - @inline ℑxyz(i, j, k, grid, from::F, to::T, f, args...) where {F<:Tuple{<:$LX, <:$LY, <:$LZ}, T<:Tuple{<:$IX, <:$IY, <:$IZ}} = + + @inline ℑxyz(i, j, k, grid, from::F, to::T, f, args...) where {F<:Tuple{<:$LX, <:$LY, <:$LZ}, T<:Tuple{<:$IX, <:$IY, <:$IZ}} = $interp_func(i, j, k, grid, f, args...) end end diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 5735444c63..f82f49d83a 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -1,23 +1,5 @@ using Oceananigans.Grids: Center, Face -const RG = RectilinearGrid -const RGX = XRegularRG -const RGY = YRegularRG -const RGZ = ZRegularRG - -const OSSG = OrthogonalSphericalShellGrid -const OSSGZ = ZRegOrthogonalSphericalShellGrid - -const LLG = LatitudeLongitudeGrid -const LLGX = XRegularLLG -const LLGY = YRegularLLG -const LLGZ = ZRegularLLG - -# On the fly calculations of metrics -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} - @inline hack_cosd(φ) = cos(π * φ / 180) @inline hack_sind(φ) = sin(π * φ / 180) @@ -53,25 +35,6 @@ The operators in this file fall into three categories: ##### ##### -@inline Δxᶠᵃᵃ(i, j, k, grid) = nothing -@inline Δxᶜᵃᵃ(i, j, k, grid) = nothing -@inline Δyᵃᶠᵃ(i, j, k, grid) = nothing -@inline Δyᵃᶜᵃ(i, j, k, grid) = nothing - -const ZRG = Union{LLGZ, RGZ} - -@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::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ᵃᵃᶠ - # Convenience Functions for all grids for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ) @@ -105,51 +68,106 @@ end ##### Rectilinear Grids (Flat grids already have Δ = 1) ##### -@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 Δ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] +@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ᵃᶜᵃ + +## 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 ##### ## 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 Δ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) + +@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] @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 Δyᶜᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶜᶠᵃ[j] -@inline Δyᶠᶜᵃ(i, j, k, grid::LLG) = @inbounds grid.Δyᶠᶜᵃ[j] +### YRegularLLG with pre-computed 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) = Δyᶠᶜᵃ(i, j, k, grid) -@inline Δyᶠᶠᵃ(i, j, k, grid::LLG) = Δyᶜᶠᵃ(i, j, k, grid) + +### 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]) -@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::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]) +@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]) + +### 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.Δφᵃᶜᵃ) @@ -167,6 +185,12 @@ 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!! @@ -229,14 +253,16 @@ end ##### ##### -@inline Vᶜᶜᶜ(i, j, k, grid) = Azᶜᶜᶜ(i, j, k, grid) * Δzᶜᶜᶜ(i, j, k, grid) -@inline Vᶠᶜᶜ(i, j, k, grid) = Azᶠᶜᶜ(i, j, k, grid) * Δzᶠᶜᶜ(i, j, k, grid) -@inline Vᶜᶠᶜ(i, j, k, grid) = Azᶜᶠᶜ(i, j, k, grid) * Δzᶜᶠᶜ(i, j, k, grid) -@inline Vᶜᶜᶠ(i, j, k, grid) = Azᶜᶜᶠ(i, j, k, grid) * Δzᶜᶜᶠ(i, j, k, grid) -@inline Vᶠᶠᶜ(i, j, k, grid) = Azᶠᶠᶜ(i, j, k, grid) * Δzᶠᶠᶜ(i, j, k, grid) -@inline Vᶠᶜᶠ(i, j, k, grid) = Azᶠᶜᶠ(i, j, k, grid) * Δzᶠᶜᶠ(i, j, k, grid) -@inline Vᶜᶠᶠ(i, j, k, grid) = Azᶜᶠᶠ(i, j, k, grid) * Δzᶜᶠᶠ(i, j, k, grid) -@inline Vᶠᶠᶠ(i, j, k, grid) = Azᶠᶠᶠ(i, j, k, grid) * Δzᶠᶠᶠ(i, j, k, grid) +for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ) + + volume = Symbol(:V, LX, LY, LZ) + z_area = Symbol(:Az, LX, LY, LZ) + z_spacing = Symbol(:Δz, LX, LY, LZ) + + @eval begin + @inline $volume(i, j, k, grid) = $z_area(i, j, k, grid) * $z_spacing(i, j, k, grid) + end +end ##### ##### Generic functions for specified locations @@ -248,9 +274,9 @@ end location_code(LX, LY, LZ) = Symbol(interpolation_code(LX), interpolation_code(LY), interpolation_code(LZ)) -for LX in (:Center, :Face) - for LY in (:Center, :Face) - for LZ in (:Center, :Face) +for LX in (:Center, :Face, :Nothing) + for LY in (:Center, :Face, :Nothing) + for LZ in (:Center, :Face, :Nothing) LXe = @eval $LX LYe = @eval $LY LZe = @eval $LZ @@ -271,4 +297,3 @@ for LX in (:Center, :Face) end end end - diff --git a/test/test_grids.jl b/test/test_grids.jl index 0d94ee3ee1..b30b48c1bb 100644 --- a/test/test_grids.jl +++ b/test/test_grids.jl @@ -2,9 +2,9 @@ include("dependencies_for_runtests.jl") include("data_dependencies.jl") using Oceananigans.Grids: total_extent, - xspacings, yspacings, zspacings, + xspacings, yspacings, zspacings, xnode, ynode, znode, λnode, φnode, - λspacings, φspacings, λspacing, φspacing + λspacing, φspacing, λspacings, φspacings using Oceananigans.Operators: Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶠᶠᵃ, Δxᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ @@ -182,8 +182,12 @@ function test_regular_rectilinear_xnode_ynode_znode_and_spacings(arch, FT) variably_spaced_grid = RectilinearGrid(arch, FT; size, topology, x=domain, y=domain, z=domain) - grids_types = ["regularly spaced", "variably spaced"] - grids = [regular_spaced_grid, variably_spaced_grid] + ibg_regular_spaced_grid = ImmersedBoundaryGrid(regular_spaced_grid, GridFittedBottom((x, y) -> 0)) + + ibg_variably_spaced_grid = ImmersedBoundaryGrid(variably_spaced_grid, GridFittedBottom((x, y) -> 0)) + + grids_types = ["regularly spaced", "variably spaced", "IBG regularly spaced", "IBG variably spaced"] + grids = [regular_spaced_grid, variably_spaced_grid, ibg_regular_spaced_grid, ibg_variably_spaced_grid] for (grid_type, grid) in zip(grids_types, grids) @info " Testing grid utils on $grid_type grid...." @@ -208,9 +212,9 @@ function test_regular_rectilinear_xnode_ynode_znode_and_spacings(arch, FT) @test all(y ≈ FT(π/N) for y in yspacings(grid, Face())) @test all(z ≈ FT(π/N) for z in zspacings(grid, Face())) - @test xspacings(grid, Face()) == xspacings(grid, Face(), Center(), Center()) - @test yspacings(grid, Face()) == yspacings(grid, Center(), Face(), Center()) - @test zspacings(grid, Face()) == zspacings(grid, Center(), Center(), Face()) + @test all(xspacings(grid, Face()) .== xspacings(grid, Face(), Center(), Center())) + @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) @@ -407,8 +411,9 @@ function test_rectilinear_grid_correct_spacings(FT, N) @test all(isapprox.( grid.zᵃᵃᶜ[1:N], zᵃᵃᶜ.(1:N) )) @test all(isapprox.( grid.Δzᵃᵃᶜ[1:N], Δzᵃᵃᶜ.(1:N) )) - @test all(isapprox.(zspacings(grid, Face(), with_halos=true), grid.Δzᵃᵃᶠ)) - @test all(isapprox.(zspacings(grid, Center(), with_halos=true), grid.Δzᵃᵃᶜ)) + @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 minimum_zspacing(grid, Center(), Center(), Center()) ≈ minimum(grid.Δzᵃᵃᶜ[1:grid.Nz]) @@ -539,21 +544,21 @@ function test_basic_lat_lon_general_grid(FT) @test typeof(grid_reg.Δzᵃᵃᶜ) == typeof(grid_reg.Δzᵃᵃᶠ) == FT - @test xspacings(grid_reg, Center(), Center(), with_halos=true) == grid_reg.Δxᶜᶜᵃ - @test xspacings(grid_reg, Center(), Face(), with_halos=true) == grid_reg.Δxᶜᶠᵃ - @test xspacings(grid_reg, Face(), Center(), with_halos=true) == grid_reg.Δxᶠᶜᵃ - @test xspacings(grid_reg, Face(), Face(), with_halos=true) == grid_reg.Δxᶠᶠᵃ - @test yspacings(grid_reg, Center(), Face(), with_halos=true) == grid_reg.Δyᶜᶠᵃ - @test yspacings(grid_reg, Face(), Center(), with_halos=true) == grid_reg.Δyᶠᶜᵃ - @test zspacings(grid_reg, Center(), with_halos=true) == grid_reg.Δzᵃᵃᶜ - @test zspacings(grid_reg, Face(), with_halos=true) == grid_reg.Δzᵃᵃᶠ - - @test xspacings(grid_reg, Center(), Center(), Center()) == xspacings(grid_reg, Center(), Center()) - @test xspacings(grid_reg, Face(), Face(), Center()) == xspacings(grid_reg, Face(), Face()) - @test yspacings(grid_reg, Center(), Face(), Center()) == yspacings(grid_reg, Center(), Face()) - @test yspacings(grid_reg, Face(), Center(), Center()) == yspacings(grid_reg, Face(), Center()) - @test zspacings(grid_reg, Face(), Face(), Center()) == zspacings(grid_reg, Center()) - @test zspacings(grid_reg, Face(), Center(), Face() ) == zspacings(grid_reg, Face()) + @test all(xspacings(grid_reg, Center(), Center()) .== reshape(grid_reg.Δxᶜᶜᵃ[1:Nφ], 1, Nφ, 1)) + @test all(xspacings(grid_reg, Center(), Face() ) .== reshape(grid_reg.Δxᶜᶠᵃ[1:Nφ+1], 1, Nφ+1, 1)) + @test all(xspacings(grid_reg, Face(), Center()) .== reshape(grid_reg.Δxᶠᶜᵃ[1:Nφ], 1, Nφ, 1)) + @test all(xspacings(grid_reg, Face(), Face()) .== reshape(grid_reg.Δxᶠᶠᵃ[1:Nφ+1], 1, Nφ+1, 1)) + @test all(yspacings(grid_reg, Center(), Face()) .== grid_reg.Δyᶜᶠᵃ) + @test all(yspacings(grid_reg, Face(), Center()) .== grid_reg.Δyᶠᶜᵃ) + @test all(zspacings(grid_reg, Center()) .== grid_reg.Δzᵃᵃᶜ) + @test all(zspacings(grid_reg, Face()) .== grid_reg.Δzᵃᵃᶠ) + + @test all(xspacings(grid_reg, Center(), Center(), Center()) .== xspacings(grid_reg, Center(), Center())) + @test all(xspacings(grid_reg, Face(), Face(), Center()) .== xspacings(grid_reg, Face(), Face())) + @test all(yspacings(grid_reg, Center(), Face(), Center()) .== yspacings(grid_reg, Center(), Face())) + @test all(yspacings(grid_reg, Face(), Center(), Center()) .== yspacings(grid_reg, Face(), Center())) + @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] @@ -562,10 +567,10 @@ function test_basic_lat_lon_general_grid(FT) @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 λspacings(grid_reg, Center(), with_halos=true) == grid_reg.Δλᶜᵃᵃ - @test λspacings(grid_reg, Face(), with_halos=true) == grid_reg.Δλᶠᵃᵃ - @test φspacings(grid_reg, Center(), with_halos=true) == grid_reg.Δφᵃᶜᵃ - @test φspacings(grid_reg, Face(), with_halos=true) == grid_reg.Δφᵃᶠᵃ + @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.Δφᵃᶠᵃ @@ -580,13 +585,13 @@ function test_basic_lat_lon_general_grid(FT) @test length(grid_str.λᶠᵃᵃ) == length(grid_reg.λᶠᵃᵃ) == Nλ + 2Hλ @test length(grid_str.λᶜᵃᵃ) == length(grid_reg.λᶜᵃᵃ) == Nλ + 2Hλ - + @test length(grid_str.φᵃᶠᵃ) == length(grid_reg.φᵃᶠᵃ) == Nφ + 2Hφ + 1 @test length(grid_str.φᵃᶜᵃ) == length(grid_reg.φᵃᶜᵃ) == Nφ + 2Hφ - + @test length(grid_str.zᵃᵃᶠ) == length(grid_reg.zᵃᵃᶠ) == Nz + 2Hz + 1 @test length(grid_str.zᵃᵃᶜ) == length(grid_reg.zᵃᵃᶜ) == Nz + 2Hz - + @test length(grid_str.Δzᵃᵃᶠ) == Nz + 2Hz + 1 @test length(grid_str.Δzᵃᵃᶜ) == Nz + 2Hz @@ -600,22 +605,19 @@ function test_basic_lat_lon_general_grid(FT) @test sum(grid_str.Δzᵃᵃᶜ) == grid_reg.Δzᵃᵃᶜ * length(grid_str.Δzᵃᵃᶜ) @test sum(grid_str.Δzᵃᵃᶠ) == grid_reg.Δzᵃᵃᶠ * length(grid_str.Δzᵃᵃᶠ) - @test xspacings(grid_str, Center(), Center(), with_halos=true) == grid_str.Δxᶜᶜᵃ - @test xspacings(grid_str, Center(), Face(), with_halos=true) == grid_str.Δxᶜᶠᵃ - @test xspacings(grid_str, Face(), Center(), with_halos=true) == grid_str.Δxᶠᶜᵃ - @test xspacings(grid_str, Face(), Face(), with_halos=true) == grid_str.Δxᶠᶠᵃ - @test yspacings(grid_str, Center(), Face(), with_halos=true) == grid_str.Δyᶜᶠᵃ - @test yspacings(grid_str, Face(), Center(), with_halos=true) == grid_str.Δyᶠᶜᵃ - @test zspacings(grid_str, Center(), with_halos=true) == grid_str.Δzᵃᵃᶜ - @test zspacings(grid_str, Face(), with_halos=true) == grid_str.Δzᵃᵃᶠ + @test all(xspacings(grid_str, Center(), Center()) .== reshape(grid_str.Δxᶜᶜᵃ[1:Nλ, 1:Nφ], Nλ, Nφ, 1)) + @test all(xspacings(grid_str, Center(), Face()) .== reshape(grid_str.Δxᶜᶠᵃ[1:Nλ, 1:Nφ+1], Nλ, Nφ+1, 1)) + @test all(xspacings(grid_str, Face(), Center()) .== reshape(grid_str.Δxᶠᶜᵃ[1:Nλ, 1:Nφ], Nλ, Nφ, 1)) + @test all(xspacings(grid_str, Face(), Face()) .== reshape(grid_str.Δxᶠᶠᵃ[1:Nλ, 1:Nφ+1], Nλ, Nφ+1, 1)) + + @test all(yspacings(grid_str, Center(), Face()) .== grid_str.Δyᶜᶠᵃ) + @test all(yspacings(grid_str, Face(), Center()) .== grid_str.Δyᶠᶜᵃ) - @test xspacings(grid_str, Center(), Center()) == grid_str.Δxᶜᶜᵃ[1:grid_str.Nx, 1:grid_str.Ny] - @test xspacings(grid_str, Center(), Face()) == grid_str.Δxᶜᶠᵃ[1:grid_str.Nx, 1:grid_str.Ny+1] - @test zspacings(grid_str, Center()) == grid_str.Δzᵃᵃᶜ[1:grid_str.Nz] - @test zspacings(grid_str, Face()) == grid_str.Δzᵃᵃᶠ[1:grid_str.Nz+1] + @test all(zspacings(grid_str, Center()) .== reshape(grid_str.Δzᵃᵃᶜ[1:Nz], 1, 1, Nz)) + @test all(zspacings(grid_str, Face()) .== reshape(grid_str.Δzᵃᵃᶠ[1:Nz+1], 1, 1, Nz+1)) - @test zspacings(grid_str, Face(), Face(), Center()) == zspacings(grid_str, Center()) - @test zspacings(grid_str, Face(), Center(), Face() ) == zspacings(grid_str, Face()) + @test all(zspacings(grid_str, Center()) .== zspacings(grid_str, Center(), Center(), Center())) + @test all(zspacings(grid_str, Face()) .== zspacings(grid_str, Face(), Center(), Face())) return nothing end @@ -689,7 +691,7 @@ function test_lat_lon_precomputed_metrics(FT, arch) println("$lat, $lon, $z") grid_pre = LatitudeLongitudeGrid(arch, FT, size=N, halo=H, latitude=lat, longitude=lon, z=z, precompute_metrics=true) grid_fly = LatitudeLongitudeGrid(arch, FT, size=N, halo=H, latitude=lat, longitude=lon, z=z) - + @test all(Array([all(Array([Δxᶠᶜᵃ(i, j, 1, grid_pre) ≈ Δxᶠᶜᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) @test all(Array([all(Array([Δxᶜᶠᵃ(i, j, 1, grid_pre) ≈ Δxᶜᶠᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) @test all(Array([all(Array([Δxᶠᶠᵃ(i, j, 1, grid_pre) ≈ Δxᶠᶠᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) @@ -699,7 +701,7 @@ function test_lat_lon_precomputed_metrics(FT, arch) @test all(Array([all(Array([Azᶜᶠᵃ(i, j, 1, grid_pre) ≈ Azᶜᶠᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) @test all(Array([all(Array([Azᶠᶠᵃ(i, j, 1, grid_pre) ≈ Azᶠᶠᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) @test all(Array([all(Array([Azᶜᶜᵃ(i, j, 1, grid_pre) ≈ Azᶜᶜᵃ(i, j, 1, grid_fly) for i in 1-Hλ+1:Nλ+Hλ-1])) for j in 1-Hφ+1:Nφ+Hφ-1])) - end + end end end @@ -737,18 +739,18 @@ function test_orthogonal_shell_grid_array_sizes_and_spacings(FT) @test size(grid.φᶜᶠᵃ) == (Nx + 2Hx, Ny + 2Hy + 1) @test size(grid.φᶠᶠᵃ) == (Nx + 2Hx + 1, Ny + 2Hy + 1) - @test xspacings(grid, Center(), Center(), Face(), with_halos=true) == xspacings(grid, Center(), Center(), with_halos=true) == grid.Δxᶜᶜᵃ - @test xspacings(grid, Center(), Face(), Face(), with_halos=true) == xspacings(grid, Center(), Face(), with_halos=true) == grid.Δxᶜᶠᵃ - @test xspacings(grid, Face(), Center(), Face()) == xspacings(grid, Face(), Center()) == grid.Δxᶠᶜᵃ[1:grid.Nx+1, 1:grid.Ny] - @test xspacings(grid, Face(), Face(), Face()) == xspacings(grid, Face(), Face()) == grid.Δxᶠᶠᵃ[1:grid.Nx+1, 1:grid.Ny+1] + @test all(xspacings(grid, Center(), Center(), Face()) .== xspacings(grid, Center(), Center()) .== grid.Δxᶜᶜᵃ[1:Nx, 1:Ny]) + @test all(xspacings(grid, Center(), Face(), Face()) .== xspacings(grid, Center(), Face() ) .== grid.Δxᶜᶠᵃ[1:Nx, 1:Ny+1]) + @test all(xspacings(grid, Face(), Center(), Face()) .== xspacings(grid, Face(), Center()) .== grid.Δxᶠᶜᵃ[1:Nx+1, 1:Ny]) + @test all(xspacings(grid, Face(), Face(), Face()) .== xspacings(grid, Face(), Face() ) .== grid.Δxᶠᶠᵃ[1:Nx+1, 1:Ny+1]) - @test yspacings(grid, Center(), Center(), Face(), with_halos=true) == yspacings(grid, Center(), Center(), with_halos=true) == grid.Δyᶜᶜᵃ - @test yspacings(grid, Center(), Face(), Face(), with_halos=true) == yspacings(grid, Center(), Face(), with_halos=true) == grid.Δyᶜᶠᵃ - @test yspacings(grid, Face(), Center(), Face()) == yspacings(grid, Face(), Center()) == grid.Δyᶠᶜᵃ[1:grid.Nx+1, 1:grid.Ny] - @test yspacings(grid, Face(), Face(), Face()) == yspacings(grid, Face(), Face()) == grid.Δyᶠᶠᵃ[1:grid.Nx+1, 1:grid.Ny+1] + @test all(yspacings(grid, Center(), Center(), Face()) .== yspacings(grid, Center(), Center()) .== grid.Δyᶜᶜᵃ[1:Nx, 1:Ny]) + @test all(yspacings(grid, Center(), Face(), Face()) .== yspacings(grid, Center(), Face() ) .== grid.Δyᶜᶠᵃ[1:Nx, 1:Ny+1]) + @test all(yspacings(grid, Face(), Center(), Face()) .== yspacings(grid, Face(), Center()) .== grid.Δyᶠᶜᵃ[1:Nx+1, 1:Ny]) + @test all(yspacings(grid, Face(), Face(), Face()) .== yspacings(grid, Face(), Face() ) .== grid.Δyᶠᶠᵃ[1:Nx+1, 1:Ny+1]) - @test zspacings(grid, Center(), Face(), Face(), with_halos=true) == zspacings(grid, Face(), with_halos=true) == grid.Δzᵃᵃᶠ - @test zspacings(grid, Center(), Face(), Center()) == zspacings(grid, Center()) == grid.Δzᵃᵃᶜ + @test all(zspacings(grid, Center(), Center(), Face() ) .== zspacings(grid, Face() ) .== grid.Δzᵃᵃᶠ) + @test all(zspacings(grid, Center(), Center(), Center()) .== zspacings(grid, Center()) .== grid.Δzᵃᵃᶜ) return nothing end @@ -806,7 +808,7 @@ end @testset "Grid equality" begin @info " Testing grid equality operator (==)..." - + for arch in archs test_grid_equality(arch) end @@ -818,7 +820,7 @@ end # Testing show function topo = (Periodic, Periodic, Periodic) - + grid = RectilinearGrid(CPU(), topology=topo, size=(3, 7, 9), x=(0, 1), y=(-π, π), z=(0, 2π)) @test try @@ -829,7 +831,7 @@ end println(sprint(showerror, err)) false end - + @test grid isa RectilinearGrid end @@ -859,7 +861,7 @@ end # Testing show function Nz = 20 grid = RectilinearGrid(arch, size=(1, 1, Nz), x=(0, 1), y=(0, 1), z=collect(0:Nz).^2) - + @test try show(grid); println() true @@ -868,7 +870,7 @@ end println(sprint(showerror, err)) false end - + @test grid isa RectilinearGrid end @@ -882,7 +884,7 @@ end @test cpu_grid_again == cpu_grid end end - + @testset "Latitude-longitude grid" begin @info " Testing general latitude-longitude grid..." @@ -901,7 +903,7 @@ end # Testing show function for regular grid grid = LatitudeLongitudeGrid(CPU(), size=(36, 32, 1), longitude=(-180, 180), latitude=(-80, 80), z=(0, 1)) - + @test try show(grid); println() true @@ -990,7 +992,7 @@ end end end end - + @testset "Conformal cubed sphere face grid" begin @info " Testing OrthogonalSphericalShellGrid grid..." @@ -1000,7 +1002,7 @@ end # Testing show function grid = conformal_cubed_sphere_panel(CPU(), size=(10, 10, 1), z=(0, 1)) - + @test try show(grid); println() true @@ -1043,4 +1045,3 @@ end end end end - diff --git a/test/utils_for_runtests.jl b/test/utils_for_runtests.jl index a027daf4de..11dc6ecd5e 100644 --- a/test/utils_for_runtests.jl +++ b/test/utils_for_runtests.jl @@ -3,23 +3,23 @@ using Oceananigans.DistributedComputations: Distributed, Partition, child_archit import Oceananigans.Fields: interior -# Are the test running on the GPUs? +# Are the test running on the GPUs? # Are the test running in parallel? child_arch = get(ENV, "GPU_TEST", nothing) == "true" ? GPU() : CPU() mpi_test = get(ENV, "MPI_TEST", nothing) == "true" # Sometimes when running tests in parallel, the CUDA.jl package is not loaded correctly. -# This function is a failsafe to re-load CUDA.jl using the suggested cach compilation from +# This function is a failsafe to re-load CUDA.jl using the suggested cach compilation from # https://github.com/JuliaGPU/CUDA.jl/blob/a085bbb3d7856dfa929e6cdae04a146a259a2044/src/initialization.jl#L105 # To make sure Julia restarts, an error is thrown. function reset_cuda_if_necessary() - + # Do nothing if we are on the CPU if child_arch isa CPU return end - - try + + try c = CUDA.zeros(10) # This will fail if CUDA is not available catch err @@ -37,9 +37,9 @@ function reset_cuda_if_necessary() end end -function test_architectures() +function test_architectures() # If MPI is initialized with MPI.Comm_size > 0, we are running in parallel. - # We test several different configurations: `Partition(x = 4)`, `Partition(y = 4)`, + # We test several different configurations: `Partition(x = 4)`, `Partition(y = 4)`, # `Partition(x = 2, y = 2)`, and different fractional subdivisions in x, y and xy if mpi_test if MPI.Initialized() && MPI.Comm_size(MPI.COMM_WORLD) == 4 @@ -48,7 +48,7 @@ function test_architectures() Distributed(child_arch; partition = Partition(2, 2)), Distributed(child_arch; partition = Partition(x = Fractional(1, 2, 3, 4))), Distributed(child_arch; partition = Partition(y = Fractional(1, 2, 3, 4))), - Distributed(child_arch; partition = Partition(x = Fractional(1, 2), y = Equal()))) + Distributed(child_arch; partition = Partition(x = Fractional(1, 2), y = Equal()))) else return throw("The MPI partitioning is not correctly configured.") end @@ -59,9 +59,9 @@ end # For nonhydrostatic simulations we cannot use `Fractional` at the moment (requirements # for the tranpose are more stringent than for hydrostatic simulations). -function nonhydrostatic_regression_test_architectures() +function nonhydrostatic_regression_test_architectures() # If MPI is initialized with MPI.Comm_size > 0, we are running in parallel. - # We test 3 different configurations: `Partition(x = 4)`, `Partition(y = 4)` + # We test 3 different configurations: `Partition(x = 4)`, `Partition(y = 4)` # and `Partition(x = 2, y = 2)` if mpi_test if MPI.Initialized() && MPI.Comm_size(MPI.COMM_WORLD) == 4 @@ -70,7 +70,7 @@ function nonhydrostatic_regression_test_architectures() Distributed(child_arch; partition = Partition(2, 2))) else return throw("The MPI partitioning is not correctly configured.") - end + end else return tuple(child_arch) end @@ -98,8 +98,8 @@ end # TODO: docstring? function center_clustered_coord(N, L, x₀) - Δz(k) = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1 - z_faces = zeros(N+1) + Δz(k) = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1 + z_faces = zeros(N+1) for k = 2:N+1 z_faces[k] = z_faces[k-1] + 3 - Δz(k-1) end @@ -109,12 +109,12 @@ end # TODO: docstring? function boundary_clustered_coord(N, L, x₀) - Δz(k) = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1 - z_faces = zeros(N+1) + Δz(k) = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1 + z_faces = zeros(N+1) for k = 2:N+1 z_faces[k] = z_faces[k-1] + Δz(k-1) end - z_faces = z_faces ./ z_faces[end] .* L .+ x₀ + z_faces = z_faces ./ z_faces[end] .* L .+ x₀ return z_faces end