(0.96.30) Fix vector rotation operation

[simone-silvestri]

The rotation and inverse rotation matrix seem to be reversed, leading to the issue encountered in CliMA/ClimaOcean.jl#549

I think I will add a couple of more tests here, since this operation seems to be quite elusive and a big source of error.

[simone-silvestri]

I have dropped this test because it was weirdly concocted, and most likely grid dependent. The new test is grid independent and tests only that the invariants of a vector remain invariant when moving to the intrinsic frame and return equal to the extrinsic components once transformed back. I think this is the only objective test we need for vector rotation.

[simone-silvestri]

Should be ready to merge when tests pass

[siddharthabishnu]

    uᵢ =   u * cosθ + v * sinθ
    vᵢ = - u * sinθ + v * cosθ

[navidcy]

@siddharthabishnu are you suggesting that there is a mistake in what @simone-silvestri proposed? I mean, is what you point out what you reckon is the right expression? (I'm only asking because you didn't include it as a suggestion so I just want to ensure that I understand properly your comment.)

I'll double check the calculation here.

[siddharthabishnu]

@navidcy, apologies for not providing a more detailed explanation earlier. Instead of using a rotation matrix, I sketched the setup by hand and resolved the components of a vector along the new directions of the rotated frame. I realized that my suggestion holds only if the frame containing the intrinsic vector is rotated counterclockwise relative to the frame with the extrinsic vector by an angle θ. If the rotation is in the opposite direction, then @simone-silvestri’s implementation is correct.

[siddharthabishnu]

    uₑ = + u * cosθ - v * sinθ
    vₑ =   u * sinθ + v * cosθ

[navidcy]

@siddharthabishnu are you suggesting that there is a mistake in what @simone-silvestri proposed? I mean, is what you point out what you reckon is the right expression? (I'm only asking because you didn't include it as a suggestion so I just want to ensure that I understand properly your comment.)

I'll double check the calculation here.

[simone-silvestri]

I think the expression was like that before. The change in this PR is indeed changing

uₑ = + u * cosθ - v * sinθ
vₑ =   u * sinθ + v * cosθ

to

uₑ = + u * cosθ + v * sinθ
vₑ = - u * sinθ + v * cosθ

This is because in the intrisic operation we are rotating the vector to comply with the grid we pass as an argument to the function. Therefore, once the angle is computed we need to use a simple rotation matrix.

In the extrinsic operation we are rotating against the grid passed in as an argument, so we are counter-rotating and we need to use the inverse of the rotation matrix, which is:
$$R^{-1} = [cosθ, sinθ; -sinθ, cosθ]$$
complying with the changes made in this PR

[siddharthabishnu]

@simone-silvestri, can you clarify if the frame of reference containing the intrinsic vector is rotated counterclockwise relative to the frame containing the extrinsic vector by an angle θ? If so, my suggestion is correct else your implementation is.

[simone-silvestri]

The frame of reference containing the intrinsic vector is rotated clockwise relative to frame containing the extrinsic vector. Or better put, the angle that we compute is positive if the intrinsic frame is clockwise rotated with respect to the extrinsic frame, and negative in the case that the intrinsic frame is rotated counterclockwise with respect to the extrinsic frame

[navidcy]

@siddharthabishnu, have a look at the docstring of the rotation_angle method at

Oceananigans.jl/src/Operators/vector_rotation_operators.jl

Lines 52 to 85 in 300c7af

	"""
	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())
	φᶠᶠᵃ⁻⁺ = φnode(i, j+1, 1, grid, Face(), Face(), Center())
	φᶠᶠᵃ⁻⁻ = φnode(i, j, 1, grid, Face(), Face(), Center())
	Δ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)
	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)
	cosθ, sinθ = Rcosθ / R, Rsinθ / R
	θ_degrees = atand(sinθ / cosθ)
	return θ_degrees
	end

Does this clarify your concerns and do you agree now that the PR changes in the rotation are correct (or at least consistent with the definition of the angle as per the docstring)?

We added also a test at

Oceananigans.jl/test/test_vector_rotation_operators.jl

Lines 85 to 125 in 300c7af

	@info " Testing the calculation of the rotation angle between the Intrinsic and Extrinsic reference frame"
	# Build a custom grid that is rotated by θᵢ degrees _clockwise_ from the vertical
	# and test that the rotation_angle is computed correctly.
	# Since we want to test the rotation angle, we build a grid with coordinates
	# that are uniformly spaced Δ << 1 apart. Small coordinate angle spacing ensures
	# that geometric factors related with spherical geometry don't come into play.
	Δ = 1e-3
	angles = [-22.5, 30, 45, 60]
	for θᵢ in angles
	radius = 1
	Nx, Ny = 4, 4
	# allocate an empty grid
	grid = OrthogonalSphericalShellGrid(; size=(Nx, Ny, 1), z=(0, 1), radius, conformal_mapping=nothing)
	λᶠᶠᵃ = grid.λᶠᶠᵃ
	φᶠᶠᵃ = grid.φᶠᶠᵃ
	# fill in coordinates
	sinθ = sind(θᵢ)
	cosθ = cosd(θᵢ)
	for i in 1:Nx+1, j in 1:Ny+1
	λᶠᶠᵃ[i, j, 1] = (i-1) * Δ * cosθ + (j-1) * Δ * sinθ
	φᶠᶠᵃ[i, j, 1] = - (i-1) * Δ * sinθ + (j-1) * Δ * cosθ
	end
	# fill in metrics
	for i in 1:Nx+1, j in 1:Ny+1
	grid.Δxᶜᶠᵃ[i, j] = haversine((λᶠᶠᵃ[i+1, j], φᶠᶠᵃ[i+1, j]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), radius)
	grid.Δyᶠᶜᵃ[i, j] = haversine((λᶠᶠᵃ[i, j+1], φᶠᶠᵃ[i, j+1]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), radius)
	end
	# ensure that rotation_angle returns θᵢ
	for i in 1:Nx, j in 1:Ny
	θ = rotation_angle(i, j, grid)
	@test θ ≈ θᵢ
	end
	end

[navidcy]

@simone-silvestri, can you clarify if the frame of reference containing the intrinsic vector is rotated counterclockwise relative to the frame containing the extrinsic vector by an angle θ? If so, my suggestion is correct else your implementation is.

the angle θ is the angle we need to rotate the intrinsic system to match the extrinsic

that is, the intrinsic system is rotated by an angle -θ compared to the extrinsic one

that's why I think it the PR is correct. But @siddharthabishnu, I admit it wasn't as clear before and I was also confused by that; I added the docstring on rotation_angle to clarify this.

[siddharthabishnu]

@simone-silvestri and @navidcy, thank you both for the detailed explanations, which along with the docstring and the new test, helped clarify my confusion. I was unsure about the relative orientation of the intrinsic and extrinsic reference frames, and I always associate a positive angle with counterclockwise rotation by default, which is what led to my misunderstanding. It's all clear now, and I agree that the PR implementation is correct and consistent with the documentation.

[navidcy]

cool

I'll merge if/when tests pass