Skip to content

(0.96.30) Fix vector rotation operation #4560

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 21 commits into from
May 30, 2025
Merged

(0.96.30) Fix vector rotation operation #4560

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
d62bad6
fix vector rotation
simone-silvestri May 29, 2025
ed02bf8
add a more comprehensive test
simone-silvestri May 29, 2025
50a478d
Update vector_rotation_operators.jl
simone-silvestri May 29, 2025
f3b26be
add more tests
simone-silvestri May 29, 2025
8f4ae83
Merge branch 'ss/fix-vector-rotation' of github.com:CliMA/Oceananigan…
simone-silvestri May 29, 2025
fc442a9
fix vector rotation in case of zero spacing
simone-silvestri May 29, 2025
6248cfe
test also the return
simone-silvestri May 29, 2025
57bb82a
this is the only test we need
simone-silvestri May 29, 2025
32a841b
clean up testing
simone-silvestri May 29, 2025
67a5031
Merge branch 'main' into ss/fix-vector-rotation
simone-silvestri May 29, 2025
33203b6
is this screwing up with the tests?
simone-silvestri May 30, 2025
9cf7842
Merge branch 'ss/fix-vector-rotation' of github.com:CliMA/Oceananigan…
simone-silvestri May 30, 2025
98d7d2f
maybe it's a random problem
simone-silvestri May 30, 2025
2f618c1
bump the project
simone-silvestri May 30, 2025
5afa53a
add the rotation
simone-silvestri May 30, 2025
36d366c
improve comment
simone-silvestri May 30, 2025
26950af
correction
simone-silvestri May 30, 2025
6a13ba5
new test
simone-silvestri May 30, 2025
531e594
rotation_matrix -> rotation_angle
navidcy May 30, 2025
f9f99d6
add rotation_angle docs
navidcy May 30, 2025
300c7af
cleanup
navidcy May 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 13 additions & 7 deletions src/Operators/vector_rotation_operators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,11 @@ _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)

cosθ₁ = ifelse(Δyᶠᶜᵃ⁺ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺)
cosθ₂ = ifelse(Δyᶠᶜᵃ⁻ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻)

# θᵢ is the rotation angle between intrinsic and extrinsic reference frame
Rcosθ = (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2
Rcosθ = (cosθ₁ + cosθ₂) / 2
Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2

# Normalization for the rotation angles
Expand All @@ -80,8 +83,8 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f
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θ
Comment on lines +103 to +104
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

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

Copy link
Member

@navidcy navidcy May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.


return uᵢ, vᵢ
end
Expand All @@ -108,8 +111,11 @@ end
Δxᶜᶠᵃ⁺ = Δxᶜᶠᶜ(i, j+1, 1, grid)
Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i, j, 1, grid)

cosθ₁ = ifelse(Δyᶠᶜᵃ⁺ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺)
cosθ₂ = ifelse(Δyᶠᶜᵃ⁻ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻)

# θᵢ is the rotation angle between intrinsic and extrinsic reference frame
Rcosθ = (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺ + deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) / 2
Rcosθ = (cosθ₁ + cosθ₂) / 2
Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2

# Normalization for the rotation angles
Expand All @@ -121,8 +127,8 @@ end
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θ
Comment on lines +128 to +129
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

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

Copy link
Member

@navidcy navidcy May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Collaborator Author

@simone-silvestri simone-silvestri May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

@siddharthabishnu siddharthabishnu May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Collaborator Author

@simone-silvestri simone-silvestri May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Contributor

@siddharthabishnu siddharthabishnu May 30, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cool

I'll merge if/when tests pass


return uₑ, vₑ
end
Expand All @@ -134,4 +140,4 @@ end
wₑ = getvalue(wᵢ, i, j, k, grid)

return uₑ, vₑ, wₑ
end
end
206 changes: 58 additions & 148 deletions test/test_vector_rotation_operators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ include("dependencies_for_runtests.jl")

using Statistics: mean
using Oceananigans.Operators
using Random

# To be used in the test below as `KernelFunctionOperation`s
@inline intrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) =
Expand All @@ -16,163 +17,72 @@ using Oceananigans.Operators
@inline extrinsic_vector_y_component(i, j, k, grid, uᵢ, vᵢ) =
@inbounds extrinsic_vector(i, j, k, grid, uᵢ, vᵢ)[2]

@inline function kinetic_energyᶜᶜᶜ(i, j, k, grid, uᶜᶜᶜ, vᶜᶜᶜ)
@inbounds u² = uᶜᶜᶜ[i, j, k]^2
@inbounds v² = vᶜᶜᶜ[i, j, k]^2
return (u² + v²) / 2
end

function kinetic_energy(u, v)
grid = u.grid
ke_op = KernelFunctionOperation{Center, Center, Center}(kinetic_energyᶜᶜᶜ, grid, u, v)
ke = Field(ke_op)
return compute!(ke)
end

function pointwise_approximate_equal(field, val)
CPU_field = on_architecture(CPU(), field)
@test all(interior(CPU_field) .≈ val)
end

# A purely zonal flow with an west-east velocity > 0
# on a cubed sphere in an intrinsic coordinate system
# has the following properties:
function test_purely_zonal_flow(uᵢ, vᵢ, grid)
c1 = maximum(uᵢ) ≈ - minimum(vᵢ)
c2 = minimum(uᵢ) ≈ - maximum(vᵢ)
c3 = mean(uᵢ) ≈ - mean(vᵢ)
c4 = mean(uᵢ) > 0 # The mean value should be positive)

return c1 & c2 & c3 & c4
end

# A purely meridional flow with a south-north velocity > 0
# on a cubed sphere in an intrinsic coordinate system
# has the following properties:
function test_purely_meridional_flow(uᵢ, vᵢ, grid)
c1 = maximum(uᵢ) ≈ maximum(vᵢ)
c2 = minimum(uᵢ) ≈ minimum(vᵢ)
c3 = mean(uᵢ) ≈ mean(vᵢ)
c4 = mean(vᵢ) > 0 # The mean value should be positive

return c1 & c2 & c3 & c4
function pointwise_approximate_equal(field1, field2)
CPU_field1 = on_architecture(CPU(), field1)
CPU_field2 = on_architecture(CPU(), field2)
@test all(interior(CPU_field1) .≈ interior(CPU_field2))
end

# Test vector invariants i.e.
# -> dot product of two vectors
# -> cross product of two vectors
function test_vector_rotation(grid)
u = CenterField(grid)
Copy link
Collaborator Author

@simone-silvestri simone-silvestri May 29, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

v = CenterField(grid)

# Purely longitudinal flow in the extrinsic coordinate system
fill!(u, 1)
fill!(v, 0)

# Convert it to an "Instrinsic" reference frame
uᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_x_component, grid, u, v)
vᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_y_component, grid, u, v)

uᵢ = compute!(Field(uᵢ))
vᵢ = compute!(Field(vᵢ))

# The extrema of u and v, as well as their mean value should
# be equivalent on an "intrinsic" frame
@test test_purely_zonal_flow(uᵢ, vᵢ, grid)

# Kinetic energy should remain the same
KE = kinetic_energy(uᵢ, vᵢ)
@apply_regionally pointwise_approximate_equal(KE, 0.5)

# Convert it back to a purely zonal velocity (vₑ == 0)
uₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ)
vₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_y_component, grid, uᵢ, vᵢ)

uₑ = compute!(Field(uₑ))
vₑ = compute!(Field(vₑ))

# Make sure that the flow was converted back to a
# purely zonal flow in the extrensic frame (v ≈ 0)
if architecture(grid) isa CPU
# Note that on the GPU, there are (apparently?) larger numerical errors
# which lead to -1e-17 < vₑ < 1e-17 for which this test fails.
@apply_regionally pointwise_approximate_equal(vₑ, 0)
end

@apply_regionally pointwise_approximate_equal(uₑ, 1)

# Purely meridional flow in the extrinsic coordinate system
fill!(u, 0)
fill!(v, 1)

# Convert it to an "Instrinsic" reference frame
uᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_x_component, grid, u, v)
vᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_y_component, grid, u, v)

uᵢ = compute!(Field(uᵢ))
vᵢ = compute!(Field(vᵢ))

# The extrema of u and v, as well as their mean value should
# be equivalent on an "intrinsic" frame
@test test_purely_meridional_flow(uᵢ, vᵢ, grid)

# Kinetic energy should remain the same
KE = kinetic_energy(uᵢ, vᵢ)
@apply_regionally pointwise_approximate_equal(KE, 0.5)

# Convert it back to a purely zonal velocity (vₑ == 0)
uₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ)
vₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_y_component, grid, uᵢ, vᵢ)

uₑ = compute!(Field(uₑ))
vₑ = compute!(Field(vₑ))

# Make sure that the flow was converted back to a
# purely zonal flow in the extrensic frame (v ≈ 0)
@apply_regionally pointwise_approximate_equal(vₑ, 1)

if architecture(grid) isa CPU
# Note that on the GPU, there are (apparently?) larger numerical errors
# which lead to - 4e-17 < uₑ < 4e-17 for which this test fails.
@apply_regionally pointwise_approximate_equal(uₑ, 0)
end

# Mixed zonal and meridional flow.
fill!(u, 0.5)
fill!(v, 0.5)

# Convert it to an "Instrinsic" reference frame
uᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_x_component, grid, u, v)
vᵢ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_y_component, grid, u, v)

uᵢ = compute!(Field(uᵢ))
vᵢ = compute!(Field(vᵢ))

# The extrema of u and v, should be equivalent on an "intrinsic" frame
# when u == v on an extrinsic frame
@test maximum(uᵢ) ≈ maximum(vᵢ)
@test minimum(uᵢ) ≈ minimum(vᵢ)

# Kinetic energy should remain the same
KE = kinetic_energy(uᵢ, vᵢ)
@apply_regionally pointwise_approximate_equal(KE, 0.25)

# Convert it back to a purely zonal velocity (vₑ == 0)
uₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ)
vₑ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_y_component, grid, uᵢ, vᵢ)

uₑ = compute!(Field(uₑ))
vₑ = compute!(Field(vₑ))

# Make sure that the flow was converted back to a
# purely zonal flow in the extrensic frame (v ≈ 0)
@apply_regionally pointwise_approximate_equal(vₑ, 0.5)
@apply_regionally pointwise_approximate_equal(uₑ, 0.5)
x₁ = CenterField(grid)
y₁ = CenterField(grid)
x₂ = CenterField(grid)
y₂ = CenterField(grid)

Random.seed!(1234)
set!(x₁, (x, y, z) -> rand())
set!(y₁, (x, y, z) -> rand())
set!(x₂, (x, y, z) -> rand())
set!(y₂, (x, y, z) -> rand())

fill_halo_regions!((x₁, y₁, x₂, y₂))

d = compute!(Field(x₁ * x₂ + y₁ * y₂))
c = compute!(Field(x₁ * y₂ - y₁ * x₂))

xᵢ₁ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_x_component, grid, x₁, y₁)
yᵢ₁ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_y_component, grid, x₁, y₁)
xᵢ₂ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_x_component, grid, x₂, y₂)
yᵢ₂ = KernelFunctionOperation{Center, Center, Center}(intrinsic_vector_y_component, grid, x₂, y₂)

xᵢ₁ = compute!(Field(xᵢ₁))
yᵢ₁ = compute!(Field(yᵢ₁))
xᵢ₂ = compute!(Field(xᵢ₂))
yᵢ₂ = compute!(Field(yᵢ₂))

dᵢ = compute!(Field(xᵢ₁ * xᵢ₂ + yᵢ₁ * yᵢ₂))
cᵢ = compute!(Field(xᵢ₁ * yᵢ₂ - yᵢ₁ * xᵢ₂))

@apply_regionally pointwise_approximate_equal(dᵢ, d)
@apply_regionally pointwise_approximate_equal(cᵢ, c)

xₑ₁ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_x_component, grid, xᵢ₁, yᵢ₁)
yₑ₁ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_y_component, grid, xᵢ₁, yᵢ₁)
xₑ₂ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_x_component, grid, xᵢ₂, yᵢ₂)
yₑ₂ = KernelFunctionOperation{Center, Center, Center}(extrinsic_vector_y_component, grid, xᵢ₂, yᵢ₂)

xₑ₁ = compute!(Field(xₑ₁))
yₑ₁ = compute!(Field(yₑ₁))
xₑ₂ = compute!(Field(xₑ₂))
yₑ₂ = compute!(Field(yₑ₂))

@apply_regionally pointwise_approximate_equal(xₑ₁, xᵢ₁)
@apply_regionally pointwise_approximate_equal(xₑ₁, xᵢ₁)
@apply_regionally pointwise_approximate_equal(yₑ₂, yᵢ₂)
@apply_regionally pointwise_approximate_equal(yₑ₂, yᵢ₂)
end

@testset "Vector rotation" begin
for arch in archs
@testset "Conversion from Intrinsic to Extrinsic reference frame [$(typeof(arch))]" begin
@info " Testing the conversion of a vector between the Intrinsic and Extrinsic reference frame"
grid = ConformalCubedSphereGrid(arch; panel_size=(10, 10, 1), z=(-1, 0))
test_vector_rotation(grid)
cubed_sphere_grid = ConformalCubedSphereGrid(arch; panel_size=(10, 10, 1), z=(-1, 0))
tripolar_grid = TripolarGrid(arch; size = (40, 40, 1), z=(-1, 0))
test_vector_rotation(cubed_sphere_grid)
test_vector_rotation(tripolar_grid)
end
end
end
Loading