Insert updates from PRs #4724, #4743 and #4744

Author: siddharthabishnu

src/Grids/grid_utils.jl

@@ -329,81 +329,6 @@ coordinate_summary(topo, Δ::Union{AbstractVector, AbstractMatrix}, name) =
 @inline column_depthᶜᶠᵃ(i, j, k, grid, η) = static_column_depthᶜᶠᵃ(i, j, grid)
 @inline column_depthᶠᶠᵃ(i, j, k, grid, η) = static_column_depthᶠᶠᵃ(i, j, grid)
 
-#####
-##### Spherical geometry
-#####
-
-"""
-    spherical_area_triangle(a::Number, b::Number, c::Number)
-
-Return the area of a spherical triangle on the unit sphere with sides `a`, `b`, and `c`.
-
-The area of a spherical triangle on the unit sphere is ``E = A + B + C - π``, where ``A``, ``B``, and ``C``
-are the triangle's inner angles.
-
-It has been known since the time of Euler and Lagrange that
-``\\tan(E/2) = P / (1 + \\cos a + \\cos b + \\cos c)``, where
-``P = (1 - \\cos²a - \\cos²b - \\cos²c + 2 \\cos a \\cos b \\cos c)^{1/2}``.
-
-References
-==========
-
-* Euler, L. (1778) De mensura angulorum solidorum, Opera omnia, 26, 204-233 (Orig. in Acta adac. sc. Petrop. 1778)
-* Lagrange,  J.-L. (1798) Solutions de quilquies problèmes relatifs au triangles sphéruques, Oeuvres, 7, 331-359.
-"""
-function spherical_area_triangle(a::Number, b::Number, c::Number)
-    cosa = cos(a)
-    cosb = cos(b)
-    cosc = cos(c)
-
-    tan½E = sqrt(1 - cosa^2 - cosb^2 - cosc^2 + 2cosa * cosb * cosc)
-    tan½E /= 1 + cosa + cosb + cosc
-
-    return 2atan(tan½E)
-end
-
-"""
-    spherical_area_triangle(a::AbstractVector, b::AbstractVector, c::AbstractVector)
-
-Return the area of a spherical triangle on the unit sphere with vertices given by the 3-vectors
-`a`, `b`, and `c` whose origin is the the center of the sphere. The formula was first given by
-Eriksson (1990).
-
-If we denote with ``A``, ``B``, and ``C`` the inner angles of the spherical triangle and with
-``a``, ``b``, and ``c`` the side of the triangle then, it has been known since Euler and Lagrange
-that ``\\tan(E/2) = P / (1 + \\cos a + \\cos b + \\cos c)``, where ``E = A + B + C - π`` is the
-triangle's excess and ``P = (1 - \\cos²a - \\cos²b - \\cos²c + 2 \\cos a \\cos b \\cos c)^{1/2}``.
-On the unit sphere, ``E`` is precisely the area of the spherical triangle. Erikkson (1990) showed
-that ``P`` above is the same as the volume defined by the vectors `a`, `b`, and `c`, that is
-``P = |𝐚 \\cdot (𝐛 \\times 𝐜)|``.
-
-References
-==========
-
-* Eriksson, F. (1990) On the measure of solid angles, Mathematics Magazine, 63 (3), 184-187, doi:10.1080/0025570X.1990.11977515
-"""
-function spherical_area_triangle(a₁::AbstractVector, a₂::AbstractVector, a₃::AbstractVector)
-    (sum(a₁.^2) ≈ 1 && sum(a₂.^2) ≈ 1 && sum(a₃.^2) ≈ 1) || error("a₁, a₂, a₃ must be unit vectors")
-
-    tan½E = abs(dot(a₁, cross(a₂, a₃)))
-    tan½E /= 1 + dot(a₁, a₂) + dot(a₂, a₃) + dot(a₁, a₃)
-
-    return 2atan(tan½E)
-end
-
-"""
-    spherical_area_quadrilateral(a₁, a₂, a₃, a₄)
-
-Return the area of a spherical quadrilateral on the unit sphere whose points are given by 3-vectors,
-`a`, `b`, `c`, and `d`. The area of the quadrilateral is given as the sum of the ares of the two
-non-overlapping triangles. To avoid having to pick the triangles appropriately ensuring they are not
-overlapping, we compute the area of the quadrilateral as the half the sum of the areas of all four potential
-triangles formed by `a₁`, `a₂`, `a₃`, and `a₄`.
-"""
-spherical_area_quadrilateral(a::AbstractVector, b::AbstractVector, c::AbstractVector, d::AbstractVector) =
-    1/2 * (spherical_area_triangle(a, b, c) + spherical_area_triangle(a, b, d) +
-           spherical_area_triangle(a, c, d) + spherical_area_triangle(b, c, d))
-
 """
     add_halos(data, loc, topo, sz, halo_sz; warnings=true)
 

src/MultiRegion/multi_region_field.jl

@@ -9,7 +9,7 @@ using Base: @propagate_inbounds
 import Oceananigans.DistributedComputations: reconstruct_global_field, CommunicationBuffers
 import Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
 import Oceananigans.Grids: xnodes, ynodes
-import Oceananigans.Fields: set!, compute!, compute_at!, validate_field_data, validate_boundary_conditions
+import Oceananigans.Fields: set!, compute!, compute_at!, interior, validate_field_data, validate_boundary_conditions
 import Oceananigans.Fields: validate_indices, communication_buffers
 import Oceananigans.Diagnostics: hasnan
 
@@ -140,6 +140,16 @@ function compute!(comp::MultiRegionComputedField, time=nothing)
     return comp
 end
 
+function interior(mrf::MultiRegionField)
+    @apply_regionally interior_mrf = interior(mrf)
+    return interior_mrf
+end
+
+function interior(mrf::MultiRegionField, I...)
+    @apply_regionally interior_mrf = interior(mrf, I...)
+    return interior_mrf
+end
+
 @inline hasnan(field::MultiRegionField) = (&)(construct_regionally(hasnan, field).regional_objects...)
 
 validate_indices(indices, loc, mrg::MultiRegionGrids) =

src/MultiRegion/multi_region_grid.jl

@@ -39,6 +39,15 @@ const MultiRegionGrids{FT, TX, TY, TZ} = Union{MultiRegionGrid{FT, TX, TY, TZ},
 @inline Base.lastindex(mrg::MultiRegionGrids) = length(mrg)
 number_of_regions(mrg::MultiRegionGrids) = lastindex(mrg)
 
+minimum_xspacing(grid::MultiRegionGrid) =
+    minimum(minimum_xspacing(grid[r]) for r in 1:number_of_regions(grid))
+
+minimum_yspacing(grid::MultiRegionGrid) =
+    minimum(minimum_yspacing(grid[r]) for r in 1:number_of_regions(grid))
+
+minimum_zspacing(grid::MultiRegionGrid) =
+    minimum(minimum_zspacing(grid[r]) 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))
 

src/OrthogonalSphericalShellGrids/OrthogonalSphericalShellGrids.jl

@@ -6,6 +6,8 @@ export TripolarGrid, RotatedLatitudeLongitudeGrid, ConformalCubedSpherePanelGrid
 import Oceananigans
 import Oceananigans.Architectures: on_architecture
 
+using CubedSphere
+
 using Oceananigans
 using Oceananigans.Grids
 using Oceananigans.Architectures: device, on_architecture, AbstractArchitecture, CPU, GPU
@@ -15,9 +17,7 @@ using Oceananigans.Utils
 using Oceananigans.BoundaryConditions: Zipper
 using Oceananigans.Fields: index_binary_search, convert_to_0_360
 using Oceananigans.Grids: RightConnected
-using Oceananigans.Grids: R_Earth,
-                          halo_size, spherical_area_quadrilateral,
-                          lat_lon_to_cartesian, generate_coordinate, topology
+using Oceananigans.Grids: R_Earth, halo_size, lat_lon_to_cartesian, generate_coordinate, topology
 using Oceananigans.Grids: total_length, add_halos, fill_metric_halo_regions!
 
 using Oceananigans.Operators

src/OrthogonalSphericalShellGrids/rotated_latitude_longitude_grid.jl

@@ -1,3 +1,4 @@
+using CubedSphere
 using Oceananigans.Grids: LatitudeLongitudeGrid, Bounded
 using Oceananigans.Utils: KernelParameters
 using StaticArrays
@@ -58,23 +59,6 @@ function rotate_metrics!(grid, shifted_lat_lon_grid)
     return nothing
 end
 
-# Convert from Spherical to Cartesian
-function spherical_to_cartesian(λ, φ, r=1)
-    x = r * cos(φ) * cos(λ)
-    y = r * cos(φ) * sin(λ)
-    z = r * sin(φ)
-    return SVector(x, y, z)
-end
-
-# Convert from Cartesian to Spherical
-function cartesian_to_spherical(X)
-    x, y, z = X
-    r = norm(X)
-    φ = asin(z / r)
-    λ = atan(y, x)
-    return λ, φ
-end
-
 # Rotation about x-axis by dλ (Change in Longitude)
 x_rotation(dλ) = @SMatrix [1        0       0
                            0  cos(dλ) -sin(dλ)