(0.96.30) Fix vector rotation operation (#4560)

* fix vector rotation

* add a more comprehensive test

* Update vector_rotation_operators.jl

* add more tests

* fix vector rotation in case of zero spacing

* test also the return

* this is the only test we need

* clean up testing

* is this screwing up with the tests?

* maybe it's a random problem

* bump the project

* add the rotation

* improve comment

* correction

* new test

* rotation_matrix -> rotation_angle

* add rotation_angle docs

* cleanup

---------

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>

Author: simone-silvestri

Project.toml

@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.96.29"
+version = "0.96.30"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/Operators/vector_rotation_operators.jl

@@ -13,8 +13,8 @@ coordinate system associated with the domain, to the coordinate system _intrinsi
 
 _extrinsic_ coordinate systems are:
 
-- Cartesian for any grid that discretizes a Cartesian domain (e.g. a `RectilinearGrid`)
-- Geographic coordinates for any grid that discretizes a Spherical domain (e.g. an `AbstractCurvilinearGrid`)
+- Cartesian coordinates for any grid that discretizes a cartesian domain (e.g. a `RectilinearGrid`)
+- Geographic coordinates for any grid that discretizes a spherical domain (e.g. an `AbstractCurvilinearGrid`)
 
 Therefore, for the [`RectilinearGrid`](@ref) and the [`LatitudeLongitudeGrid`](@ref), the _extrinsic_ and the
 _intrinsic_ coordinate system are equivalent. However, for other grids (e.g., for the
@@ -31,8 +31,8 @@ system of the grid, to the _extrinsic_ coordinate system associated with the dom
 
 _extrinsic_ coordinate systems are:
 
-- Cartesian for any grid that discretizes a Cartesian domain (e.g. a `RectilinearGrid`)
-- Geographic coordinates for any grid that discretizes a Spherical domain (e.g. an `AbstractCurvilinearGrid`)
+- Cartesian coordinates for any grid that discretizes a cartesian domain (e.g. a `RectilinearGrid`)
+- Geographic coordinates for any grid that discretizes a spherical domain (e.g. an `AbstractCurvilinearGrid`)
 
 Therefore, for the [`RectilinearGrid`](@ref) and the [`LatitudeLongitudeGrid`](@ref), the _extrinsic_ and the
 _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., for the
@@ -48,14 +48,15 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f
 @inline extrinsic_vector(i, j, k, grid::AbstractGrid, uᵢ, vᵢ) =
     getvalue(uᵢ, i, j, k, grid), getvalue(vᵢ, i, j, k, grid)
 
-# Intrinsic and extrinsic conversion for `OrthogonalSphericalShellGrid`s,
-# i.e. curvilinear grids defined on a sphere which are locally orthogonal.
-# If the coordinates match with the coordinates of a latitude-longitude grid
-# (i.e. globally orthogonal), these functions collapse to
-# uₑ, vₑ, wₑ = uᵢ, vᵢ, wᵢ
 
-# 2D vectors
-@inline function intrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uₑ, vₑ)
+"""
+    rotation_angle(i, j, grid::OrthogonalSphericalShellGrid)
+
+Return the rotation angle (in degrees) of the `i, j`-th point of the `grid`.
+The rotation angle is the angle (positive counter-clockwise) that we need to rotate
+the grid's intrinsic coordinates in order to match the grid's extrinsic coordinates.
+"""
+@inline function rotation_angle(i, j, grid::OrthogonalSphericalShellGrid)
 
     φᶠᶠᵃ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), Center())
     φᶠᶠᵃ⁺⁻ = φnode(i+1, j,   1, grid, Face(), Face(), Center())
@@ -67,21 +68,40 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f
     Δxᶜᶠᵃ⁺ = Δxᶜᶠᶜ(i,   j+1, 1, grid)
     Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i,   j,   1, grid)
 
-    # θᵢ is the rotation angle between intrinsic and extrinsic reference frame
-    Rcosθ =   (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2
+    Rcosθ₁ = ifelse(Δyᶠᶜᵃ⁺ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺)
+    Rcosθ₂ = ifelse(Δyᶠᶜᵃ⁻ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻)
+
+    # θ is the rotation angle between intrinsic and extrinsic reference frame
+    Rcosθ =   (Rcosθ₁ + Rcosθ₂) / 2
     Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2
 
     # Normalization for the rotation angles
     R = sqrt(Rcosθ^2 + Rsinθ^2)
 
-    u  = getvalue(uₑ, i, j, k, grid)
-    v  = getvalue(vₑ, i, j, k, grid)
+    cosθ, sinθ = Rcosθ / R, Rsinθ / R
+
+    θ_degrees = atand(sinθ / cosθ)
+    return θ_degrees
+end
+
+# Intrinsic and extrinsic conversion for `OrthogonalSphericalShellGrid`s,
+# i.e. curvilinear grids defined on a sphere which are locally orthogonal.
+# If the coordinates match with the coordinates of a latitude-longitude grid
+# (i.e. globally orthogonal), these functions collapse to
+# uₑ, vₑ, wₑ = uᵢ, vᵢ, wᵢ
+
+# 2D vectors
+@inline function intrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uₑ, vₑ)
 
-    cosθ = Rcosθ / R
-    sinθ = Rsinθ / R
+    u = getvalue(uₑ, i, j, k, grid)
+    v = getvalue(vₑ, i, j, k, grid)
 
-    uᵢ =   u * cosθ + v * sinθ
-    vᵢ = - u * sinθ + v * cosθ
+    θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid)
+    sinθ = sind(θ_degrees)
+    cosθ = cosd(θ_degrees)
+
+    uᵢ = u * cosθ - v * sinθ
+    vᵢ = u * sinθ + v * cosθ
 
     return uᵢ, vᵢ
 end
@@ -98,31 +118,15 @@ end
 # 2D vectors
 @inline function extrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uᵢ, vᵢ)
 
-    φᶠᶠᵃ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), Center())
-    φᶠᶠᵃ⁺⁻ = φnode(i+1, j,   1, grid, Face(), Face(), Center())
-    φᶠᶠᵃ⁻⁺ = φnode(i,   j+1, 1, grid, Face(), Face(), Center())
-    φᶠᶠᵃ⁻⁻ = φnode(i,   j,   1, grid, Face(), Face(), Center())
+    u = getvalue(uᵢ, i, j, k, grid)
+    v = getvalue(vᵢ, i, j, k, grid)
 
-    Δyᶠᶜᵃ⁺ = Δyᶠᶜᶜ(i+1, j,   1, grid)
-    Δyᶠᶜᵃ⁻ = Δyᶠᶜᶜ(i,   j,   1, grid)
-    Δxᶜᶠᵃ⁺ = Δxᶜᶠᶜ(i,   j+1, 1, grid)
-    Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i,   j,   1, grid)
+    θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid)
+    sinθ = sind(θ_degrees)
+    cosθ = cosd(θ_degrees)
 
-    # θᵢ is the rotation angle between intrinsic and extrinsic reference frame
-    Rcosθ =   (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2
-    Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2
-
-    # Normalization for the rotation angles
-    R  = sqrt(Rcosθ^2 + Rsinθ^2)
-
-    u  = getvalue(uᵢ, i, j, k, grid)
-    v  = getvalue(vᵢ, i, j, k, grid)
-
-    cosθ = Rcosθ / R
-    sinθ = Rsinθ / R
-
-    uₑ = u * cosθ - v * sinθ
-    vₑ = u * sinθ + v * cosθ
+    uₑ = + u * cosθ + v * sinθ
+    vₑ = - u * sinθ + v * cosθ
 
     return uₑ, vₑ
 end
@@ -134,4 +138,4 @@ end
     wₑ = getvalue(wᵢ, i, j, k, grid)
 
     return uₑ, vₑ, wₑ
-end
+end