Fix RotatedLatitudeLongitudeGrid and add more tests (#4879)

* Fix RotatedLatitudeLongitudeGrid and add more tests

* fix error comment

* fix docstring

* infer topology as for ordinary lat lon grid

Author: glwagner

src/OrthogonalSphericalShellGrids/rotated_latitude_longitude_grid.jl

@@ -1,5 +1,5 @@
 using CubedSphere.SphericalGeometry
-using Oceananigans.Grids: LatitudeLongitudeGrid, Bounded
+using Oceananigans.Grids: LatitudeLongitudeGrid, Bounded, validate_lat_lon_grid_args
 using Oceananigans.Utils: KernelParameters
 using StaticArrays
 using LinearAlgebra
@@ -59,24 +59,26 @@ z = (0, 1)
 grid = RotatedLatitudeLongitudeGrid(; size, longitude, latitude, z, north_pole=(70, 55))
 
 # output
-90×40×1 OrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
-├── centered at (λ, φ) = (146.656, 11.3134)
-├── longitude: Bounded  extent 360.0 degrees variably spaced with min(Δλ)=0.694593, max(Δλ)=4.0
-├── latitude:  Bounded  extent 160.0 degrees variably spaced with min(Δφ)=4.0, max(Δφ)=4.0
-└── z:         Bounded  z ∈ [0.0, 1.0]       regularly spaced with Δz=1.0
+90×40×1 OrthogonalSphericalShellGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
+├── centered at (λ, φ) = (-110.0, 35.0)
+├── longitude: Periodic  extent 360.0 degrees variably spaced with min(Δλ)=0.694593, max(Δλ)=4.0
+├── latitude:  Bounded  extent 160.0 degrees  variably spaced with min(Δφ)=4.0, max(Δφ)=4.0
+└── z:         Bounded  z ∈ [0.0, 1.0]        regularly spaced with Δz=1.0
 ```
 
-We can also make an ordinary LatitudeLongitudeGrid using `north_polar = (0, 90)`:
+Note that the center latitude ``λ = -110`` follows from ``180 + 70 - 360 = -110``:
+a clockwise rotation of 70 degrees modulo 360 degrees.
+We can also make an ordinary LatitudeLongitudeGrid using `north_pole = (0, 90)`:
 
 ```jldoctest rllg
 grid = RotatedLatitudeLongitudeGrid(; size, longitude, latitude, z, north_pole=(0, 90))
 
 # output
-90×40×1 OrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
+90×40×1 OrthogonalSphericalShellGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
 ├── centered at (λ, φ) = (180.0, 0.0)
-├── longitude: Bounded  extent 360.0 degrees variably spaced with min(Δλ)=0.694593, max(Δλ)=4.0
-├── latitude:  Bounded  extent 160.0 degrees variably spaced with min(Δφ)=4.0, max(Δφ)=4.0
-└── z:         Bounded  z ∈ [0.0, 1.0]       regularly spaced with Δz=1.0
+├── longitude: Periodic  extent 360.0 degrees variably spaced with min(Δλ)=0.694593, max(Δλ)=4.0
+├── latitude:  Bounded  extent 160.0 degrees  variably spaced with min(Δφ)=4.0, max(Δφ)=4.0
+└── z:         Bounded  z ∈ [0.0, 1.0]        regularly spaced with Δz=1.0
 ```
 """
 function RotatedLatitudeLongitudeGrid(arch::AbstractArchitecture = CPU(),
@@ -88,15 +90,23 @@ function RotatedLatitudeLongitudeGrid(arch::AbstractArchitecture = CPU(),
                                       z,
                                       halo = (3, 3, 3),
                                       radius = Oceananigans.defaults.planet_radius,
-                                      topology = (Bounded, Bounded, Bounded))
+                                      topology = nothing)
 
-    shifted_halo = halo .+ 1
-    source_grid = LatitudeLongitudeGrid(arch, FT; size, z, topology, radius,
-                                        latitude, longitude, halo = shifted_halo)
+    precompute_metrics = true
+    topology, size, halo, latitude, longitude, z, precompute_metrics =
+        validate_lat_lon_grid_args(topology, size, halo, FT, latitude, longitude, z, precompute_metrics)
+
+    _, φ₀ = north_pole
+
+    if φ₀ < 0
+        throw(ArgumentError("North pole latitude must be >= 0."))
+    elseif φ₀ > 90
+        throw(ArgumentError("North pole latitude must be <= 90 degrees."))
+    end
 
-    Nx, Ny, Nz = size
-    Hx, Hy, Hz = halo_size(source_grid)
-    Lz = source_grid.Lz
+    shifted_halo = halo .+ 1
+    source_grid = LatitudeLongitudeGrid(arch, FT; halo=shifted_halo, precompute_metrics,
+                                        size, z, topology, radius, latitude, longitude)
 
     conformal_mapping = LatitudeLongitudeRotation(north_pole)
     grid = OrthogonalSphericalShellGrid(arch, FT; size, z, radius, halo, topology, conformal_mapping)
@@ -136,44 +146,50 @@ function cartesian_to_spherical(X)
 end
 
 # Rotation about x-axis by dλ (Change in Longitude)
-x_rotation(dλ) = @SMatrix [1        0       0
-                           0  cos(dλ) -sin(dλ)
-                           0  sin(dλ)  cos(dλ)]
+x_rotation_matrix(dλ) = @SMatrix [1        0       0
+                                  0  cos(dλ) -sin(dλ)
+                                  0  sin(dλ)  cos(dλ)]
 
 # Rotation about y-axis by dφ (Change in Latitude)
-y_rotation(dφ) = @SMatrix [ cos(dφ)  0  sin(dφ)
-                            0        1  0
-                           -sin(dφ)  0  cos(dφ)]
+y_rotation_matrix(dφ) = @SMatrix [ cos(dφ)  0  sin(dφ)
+                                   0        1  0
+                                  -sin(dφ)  0  cos(dφ)]
 
 # Rotation about z-axis by dλ (Change in Longitude)
-z_rotation(dλ) = @SMatrix [cos(dλ) -sin(dλ) 0
-                           sin(dλ)  cos(dλ) 0
-                           0        0       1]
+z_rotation_matrix(dλ) = @SMatrix [cos(dλ) -sin(dλ) 0
+                                  sin(dλ)  cos(dλ) 0
+                                  0        0       1]
+
+"""
+    rotate_coordinates(λ′, φ′, λ₀, φ₀)
 
-# Perform the rotation
+Return the geographic longitude and latitude `(λ, φ)` corresponding to the rotated
+coordinates `(λ′, φ′)` on a grid whose north pole is located at `(λ₀, φ₀)`. All
+arguments are interpreted in degrees. The rotation aligns the grid pole with the
+geographic pole, then maps the rotated point back to geographic coordinates.
+"""
 function rotate_coordinates(λ′, φ′, λ₀, φ₀)
-    λ′ *= π/180
-    φ′ *= π/180
-    λ₀ *= π/180
-    φ₀ *= π/180
+    λ′ *= π / 180
+    φ′ *= π / 180
+    λ₀ *= π / 180
+    φ₀ *= π / 180
 
-    dλ = - λ₀
+    dλ = λ₀
     dφ = π/2 - φ₀
 
     # Convert to Cartesian
     X′ = SVector(spherical_to_cartesian(φ′, λ′; check_latitude_bounds = false)...)
 
     # Rotate Cartesian coordinates
-    Rx = x_rotation(dλ)
-    Ry = y_rotation(dφ)
-    Rz = z_rotation(dλ)
-    X = Rx * Ry * X′
+    Ry = y_rotation_matrix(dφ)
+    Rz = z_rotation_matrix(dλ)
+    X = Rz * Ry * X′
 
     # Convert back to Spherical
     φ, λ = cartesian_to_spherical(X)
 
-    λ *= 180/π
-    φ *= 180/π
+    λ *= 180 / π
+    φ *= 180 / π
 
     return λ, φ
 end

test/test_grids.jl

@@ -6,7 +6,7 @@ using Oceananigans.Grids: total_extent, ColumnEnsembleSize,
                           xnode, ynode, znode, λnode, φnode,
                           λspacings, φspacings
 
-using Oceananigans.OrthogonalSphericalShellGrids: RotatedLatitudeLongitudeGrid, ConformalCubedSpherePanelGrid
+using Oceananigans.OrthogonalSphericalShellGrids: RotatedLatitudeLongitudeGrid, ConformalCubedSpherePanelGrid, rotate_coordinates
 
 using Oceananigans.Operators: Δx, Δy, Δz, Δλ, Δφ, Ax, Ay, Az, volume
 using Oceananigans.Operators: Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶠᶠᵃ, Δxᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ
@@ -1243,13 +1243,38 @@ end
                                             north_pole = (0, 0),
                                             topology = (Bounded, Bounded, Bounded))
 
-        @show grid.conformal_mapping
-
         @test grid isa OrthogonalSphericalShellGrid
         @test grid isa RotatedLatitudeLongitudeGrid
         @test grid.Lz == 1000
         @test size(grid) == (10, 10, 1)
 
+        @testset "RotatedLatitudeLongitudeGrid: rotate_coordinates utility" begin
+            for φ₀ = (90, 30, 10, 0)
+                λ, φ = rotate_coordinates(0, 0, 0, φ₀)
+                @test λ ≈ 0
+                @test φ ≈ φ₀ - 90
+            end
+
+            λᴺ, φᴺ = rotate_coordinates(0, 90, 70, 55)
+            @test λᴺ ≈ 70
+            @test φᴺ ≈ 55
+        end
+
+        @testset "RotatedLatitudeLongitudeGrid respects north_pole argument" begin
+            test_kw = (size = (6, 6, 1),
+                       latitude = (-80, 80),
+                       longitude = (-120, 120),
+                       z = (-100, 0),
+                       topology = (Bounded, Bounded, Bounded))
+
+            tilted_pole = (70, 55)
+            tilted_grid = RotatedLatitudeLongitudeGrid(north_pole = tilted_pole; test_kw...)
+            @test tilted_grid.conformal_mapping.north_pole == tilted_pole
+
+            @test_throws ArgumentError RotatedLatitudeLongitudeGrid(north_pole=(0, 91); test_kw...)
+            @test_throws ArgumentError RotatedLatitudeLongitudeGrid(north_pole=(0, -1); test_kw...)
+        end
+
         for arch in archs
             for FT in float_types
                 z = (0, 1)