From d62bad6854e86862fe72b90195f940de0ddb850a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:02:49 +0200 Subject: [PATCH 01/18] fix vector rotation --- src/Operators/vector_rotation_operators.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 8875faa829..66d123664c 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -80,8 +80,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θ return uᵢ, vᵢ end @@ -121,8 +121,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θ return uₑ, vₑ end From ed02bf8b2f86bf60f27fc1f0e288968089b045a6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:31:06 +0200 Subject: [PATCH 02/18] add a more comprehensive test --- test/test_vector_rotation_operators.jl | 48 ++++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index dbf8495466..1dda368e09 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -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ₑ) = @@ -29,9 +30,15 @@ function kinetic_energy(u, v) return compute!(ke) end -function pointwise_approximate_equal(field, val) - CPU_field = on_architecture(CPU(), field) - @test all(interior(CPU_field) .≈ val) +function pointwise_approximate_equal(field1, val::Number) + CPU_field1 = on_architecture(CPU(), field1) + @test all(interior(CPU_field1) .≈ val) +end + +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 # A purely zonal flow with an west-east velocity > 0 @@ -167,12 +174,47 @@ function test_vector_rotation(grid) @apply_regionally pointwise_approximate_equal(uₑ, 0.5) end +# Test vector invariants i.e. +# -> dot product of two vectors +# -> cross product of two vectors +function test_tripolar_vector_rotation(grid) + 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₂) + + dᵢ = compute!(Field(xᵢ₁ * xᵢ₂ + yᵢ₁ * yᵢ₂)) + cᵢ = compute!(Field(xᵢ₁ * yᵢ₂ - yᵢ₁ * xᵢ₂)) + + pointwise_approximate_equal(dᵢ, d) + pointwise_approximate_equal(cᵢ, c) +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) + grid = TripolarGrid(arch; size = (100, 100, 1), z=(-1, 0)) + TestModel_VerticallyStrectedRectGrid end end end From 50a478daaa80c2f9290bab7fba2a08a0ef9bd3f5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:32:36 +0200 Subject: [PATCH 03/18] Update vector_rotation_operators.jl --- src/Operators/vector_rotation_operators.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 66d123664c..223e630527 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -121,8 +121,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θ return uₑ, vₑ end @@ -134,4 +134,4 @@ end wₑ = getvalue(wᵢ, i, j, k, grid) return uₑ, vₑ, wₑ -end \ No newline at end of file +end From f3b26be3e560eb7ea8367a7c57a517f224384146 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:34:21 +0200 Subject: [PATCH 04/18] add more tests --- test/test_vector_rotation_operators.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 1dda368e09..3c70132e91 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -65,7 +65,7 @@ function test_purely_meridional_flow(uᵢ, vᵢ, grid) return c1 & c2 & c3 & c4 end -function test_vector_rotation(grid) +function test_cubed_sphere_vector_rotation(grid) u = CenterField(grid) v = CenterField(grid) @@ -177,7 +177,7 @@ end # Test vector invariants i.e. # -> dot product of two vectors # -> cross product of two vectors -function test_tripolar_vector_rotation(grid) +function test_vector_rotation(grid) x₁ = CenterField(grid) y₁ = CenterField(grid) @@ -203,18 +203,19 @@ function test_tripolar_vector_rotation(grid) dᵢ = compute!(Field(xᵢ₁ * xᵢ₂ + yᵢ₁ * yᵢ₂)) cᵢ = compute!(Field(xᵢ₁ * yᵢ₂ - yᵢ₁ * xᵢ₂)) - pointwise_approximate_equal(dᵢ, d) - pointwise_approximate_equal(cᵢ, c) + @apply_regionally pointwise_approximate_equal(dᵢ, d) + @apply_regionally pointwise_approximate_equal(cᵢ, c) 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) - grid = TripolarGrid(arch; size = (100, 100, 1), z=(-1, 0)) - TestModel_VerticallyStrectedRectGrid + cubed_sphere_grid = ConformalCubedSphereGrid(arch; panel_size=(10, 10, 1), z=(-1, 0)) + tripolar_grid = TripolarGrid(arch; size = (100, 100, 1), z=(-1, 0)) + test_vector_rotation(cubed_sphere_grid) + test_vector_rotation(tripolar_grid) + test_cubed_sphere_vector_rotation(cubed_sphere_grid) end end end From fc442a97aa149560cfd78f87e1b324a3c4b810c5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:51:53 +0200 Subject: [PATCH 05/18] fix vector rotation in case of zero spacing --- src/Operators/vector_rotation_operators.jl | 10 ++++++++-- test/test_vector_rotation_operators.jl | 8 ++++++-- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 223e630527..4fd06ae08f 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -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 @@ -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 diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 3c70132e91..2ba7281d35 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -180,7 +180,6 @@ end function test_vector_rotation(grid) x₁ = CenterField(grid) y₁ = CenterField(grid) - x₂ = CenterField(grid) y₂ = CenterField(grid) @@ -200,6 +199,11 @@ function test_vector_rotation(grid) 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ᵢ₂)) @@ -212,7 +216,7 @@ end @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" cubed_sphere_grid = ConformalCubedSphereGrid(arch; panel_size=(10, 10, 1), z=(-1, 0)) - tripolar_grid = TripolarGrid(arch; size = (100, 100, 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) test_cubed_sphere_vector_rotation(cubed_sphere_grid) From 6248cfe0c4e2400640c6d88baae101c24667fcc9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:53:44 +0200 Subject: [PATCH 06/18] test also the return --- test/test_vector_rotation_operators.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 2ba7281d35..7f2e6723f7 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -209,6 +209,21 @@ function test_vector_rotation(grid) @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 From 57bb82ab9e28b50fe69bb8906f80e0f93aad6878 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:54:23 +0200 Subject: [PATCH 07/18] this is the only test we need --- test/test_vector_rotation_operators.jl | 152 ------------------------- 1 file changed, 152 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 7f2e6723f7..96feb3c8a0 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -17,163 +17,12 @@ using Random @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(field1, val::Number) - CPU_field1 = on_architecture(CPU(), field1) - @test all(interior(CPU_field1) .≈ val) -end - 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 -# 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 -end - -function test_cubed_sphere_vector_rotation(grid) - u = CenterField(grid) - 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) -end - # Test vector invariants i.e. # -> dot product of two vectors # -> cross product of two vectors @@ -234,7 +83,6 @@ end tripolar_grid = TripolarGrid(arch; size = (40, 40, 1), z=(-1, 0)) test_vector_rotation(cubed_sphere_grid) test_vector_rotation(tripolar_grid) - test_cubed_sphere_vector_rotation(cubed_sphere_grid) end end end From 32a841bb5fb40856af90eb08b0c9a6422b5918e4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 29 May 2025 11:58:37 +0200 Subject: [PATCH 08/18] clean up testing --- test/test_vector_rotation_operators.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 96feb3c8a0..86b80f7bb3 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -43,6 +43,7 @@ function test_vector_rotation(grid) d = compute!(Field(x₁ * x₂ + y₁ * y₂)) c = compute!(Field(x₁ * y₂ - y₁ * x₂)) + @info " Testing the conversion of a vector between to the Intrinsice reference frame on $(summary(grid))" 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₂) @@ -69,10 +70,11 @@ function test_vector_rotation(grid) 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ᵢ₂) + @info " Testing the conversion of a vector between to the Extrinsic reference frame on $(summary(grid))" + @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 From 33203b6d88b1ca5d917598f0963b309387b8e61f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 09:11:24 +0200 Subject: [PATCH 09/18] is this screwing up with the tests? --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5c5ee49c99..487d511fa8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,7 +36,7 @@ CUDA.allowscalar() do include("test_grids.jl") include("test_immersed_boundary_grid.jl") include("test_operators.jl") - include("test_vector_rotation_operators.jl") + # include("test_vector_rotation_operators.jl") include("test_boundary_conditions.jl") include("test_field.jl") include("test_regrid.jl") From 98d7d2f64cdf037c272625a3198f5c59d2d26efa Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 09:12:36 +0200 Subject: [PATCH 10/18] maybe it's a random problem --- test/runtests.jl | 2 +- test/test_vector_rotation_operators.jl | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 487d511fa8..5c5ee49c99 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,7 +36,7 @@ CUDA.allowscalar() do include("test_grids.jl") include("test_immersed_boundary_grid.jl") include("test_operators.jl") - # include("test_vector_rotation_operators.jl") + include("test_vector_rotation_operators.jl") include("test_boundary_conditions.jl") include("test_field.jl") include("test_regrid.jl") diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 86b80f7bb3..138e55b32f 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -2,7 +2,6 @@ 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ₑ) = @@ -32,7 +31,6 @@ function test_vector_rotation(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()) From 2f618c142f1b7717b732e55dbc0c76fb07021feb Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 10:21:51 +0200 Subject: [PATCH 11/18] bump the project --- Project.toml | 2 +- src/Operators/vector_rotation_operators.jl | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 843953436d..5671e3ff1e 100644 --- a/Project.toml +++ b/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" diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 4fd06ae08f..be520eb987 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -67,11 +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ᶠᶜᵃ⁻) + 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θ = (cosθ₁ + cosθ₂) / 2 + Rcosθ = (Rcosθ₁ + Rcosθ₂) / 2 Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2 # Normalization for the rotation angles @@ -111,11 +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ᶠᶜᵃ⁻) + 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θ = (cosθ₁ + cosθ₂) / 2 + Rcosθ = (Rcosθ₁ + Rcosθ₂) / 2 Rsinθ = - (deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁻⁺) / Δxᶜᶠᵃ⁺ + deg2rad(φᶠᶠᵃ⁺⁻ - φᶠᶠᵃ⁻⁻) / Δxᶜᶠᵃ⁻) / 2 # Normalization for the rotation angles From 5afa53aacdbd380c68e82325163c578c20a89d8a Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 11:03:44 +0200 Subject: [PATCH 12/18] add the rotation --- src/Operators/vector_rotation_operators.jl | 47 +++++++++------------- test/test_vector_rotation_operators.jl | 37 +++++++++++++++++ 2 files changed, 55 insertions(+), 29 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index be520eb987..a43de724df 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -48,14 +48,7 @@ _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ₑ) +@inline function rotation_matrix(i, j, grid::OrthogonalSphericalShellGrid) φᶠᶠᵃ⁺⁺ = φnode(i+1, j+1, 1, grid, Face(), Face(), Center()) φᶠᶠᵃ⁺⁻ = φnode(i+1, j, 1, grid, Face(), Face(), Center()) @@ -65,7 +58,7 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f Δ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) + Δxᶜᶠᵃ⁻ = Δxᶜᶠᶜ(i, j, 1, grid) Rcosθ₁ = ifelse(Δyᶠᶜᵃ⁺ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁺⁺ - φᶠᶠᵃ⁺⁻) / Δyᶠᶜᵃ⁺) Rcosθ₂ = ifelse(Δyᶠᶜᵃ⁻ == 0, zero(grid), deg2rad(φᶠᶠᵃ⁻⁺ - φᶠᶠᵃ⁻⁻) / Δyᶠᶜᵃ⁻) @@ -77,6 +70,20 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f # Normalization for the rotation angles R = sqrt(Rcosθ^2 + Rsinθ^2) + return Rcosθ / R, Rsinθ / R +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θ, sinθ = rotation_matrix(i, j, grid) + u = getvalue(uₑ, i, j, k, grid) v = getvalue(vₑ, i, j, k, grid) @@ -101,26 +108,8 @@ 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()) - - Δ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θ = rotation_matrix(i, j, grid) + u = getvalue(uᵢ, i, j, k, grid) v = getvalue(vᵢ, i, j, k, grid) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 138e55b32f..fc14a93ee0 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -2,6 +2,8 @@ include("dependencies_for_runtests.jl") using Statistics: mean using Oceananigans.Operators +using Oceananigans.Operators: rotation_matrix +using Oceananigans.Grids: haversine # To be used in the test below as `KernelFunctionOperation`s @inline intrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) = @@ -76,6 +78,41 @@ function test_vector_rotation(grid) end @testset "Vector rotation" begin + # Build a custom grid that is rotated by 60 degrees + # and test that the rotation matrix is correct. + if arch isa CPU + grid = OrthogonalSphericalShellGrid(; size=(4, 4, 1), z=(0, 1), radius=1, conformal_mapping=nothing) + + # This is not really correct we are removing the lat-lon stretching factor, + # but computing the spacings like this leads to a grid with exactly 1m spacing every one degree. + # We use this to test that the grid is rotated by 60 degrees from the vertical -> i.e. a 60 degree + # clockwise rotation. + fill!(grid.Δxᶜᶠᵃ, 1) + fill!(grid.Δyᶜᶠᵃ, 1) + + angles = [-22.5, 30, 45, 60] + + # We build a grid that is rotated by θᵢ degrees clockwise from the vertical. + for θᵢ in angles + sinθ = sind(θᵢ) + cosθ = cosd(θᵢ) + + for i in 1:6, j in 1:6 + grid.λᶠᶠᵃ[1, j, 1] = (i-1) * cosθ + (j-1) * sinθ + grid.φᶠᶠᵃ[i, j, 1] = - (i-1) * sinθ + (j-1) * cosθ + end + + λᶠᶠᵃ = grid.λᶠᶠᵃ + φᶠᶠᵃ = grid.φᶠᶠᵃ + + for i in 1:3, j in 1:3 + cosθ, sinθ = rotation_matrix(i, j, grid) + θ = atand(sinθ / cosθ) + @test θ ≈ θᵢ + end + end + end + 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" From 36d366ce167bf94987197981f27cabcdb9b32b26 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 11:09:38 +0200 Subject: [PATCH 13/18] improve comment --- test/test_vector_rotation_operators.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index fc14a93ee0..f132dd6eb5 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -78,21 +78,21 @@ function test_vector_rotation(grid) end @testset "Vector rotation" begin - # Build a custom grid that is rotated by 60 degrees - # and test that the rotation matrix is correct. if arch isa CPU + + # Build a custom grid that is rotated and test that the rotation matrix is correct. + # We build a grid that is rotated by θᵢ degrees clockwise from the vertical. grid = OrthogonalSphericalShellGrid(; size=(4, 4, 1), z=(0, 1), radius=1, conformal_mapping=nothing) # This is not really correct we are removing the lat-lon stretching factor, # but computing the spacings like this leads to a grid with exactly 1m spacing every one degree. - # We use this to test that the grid is rotated by 60 degrees from the vertical -> i.e. a 60 degree - # clockwise rotation. + # We use this to test that the grid is rotated by θᵢ degrees from the vertical + # -> i.e. a θᵢ degree clockwise rotation. fill!(grid.Δxᶜᶠᵃ, 1) fill!(grid.Δyᶜᶠᵃ, 1) angles = [-22.5, 30, 45, 60] - # We build a grid that is rotated by θᵢ degrees clockwise from the vertical. for θᵢ in angles sinθ = sind(θᵢ) cosθ = cosd(θᵢ) From 26950afba61c4c21e4656d42ae2cb195836b8500 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 11:10:07 +0200 Subject: [PATCH 14/18] correction --- test/test_vector_rotation_operators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index f132dd6eb5..43efd5d80a 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -98,7 +98,7 @@ end cosθ = cosd(θᵢ) for i in 1:6, j in 1:6 - grid.λᶠᶠᵃ[1, j, 1] = (i-1) * cosθ + (j-1) * sinθ + grid.λᶠᶠᵃ[i, j, 1] = (i-1) * cosθ + (j-1) * sinθ grid.φᶠᶠᵃ[i, j, 1] = - (i-1) * sinθ + (j-1) * cosθ end From 6a13ba56bfdca1f7cbc92e4fb4d2326b8422114c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 30 May 2025 11:13:46 +0200 Subject: [PATCH 15/18] new test --- test/test_vector_rotation_operators.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 43efd5d80a..23bd83cb5e 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -84,13 +84,6 @@ end # We build a grid that is rotated by θᵢ degrees clockwise from the vertical. grid = OrthogonalSphericalShellGrid(; size=(4, 4, 1), z=(0, 1), radius=1, conformal_mapping=nothing) - # This is not really correct we are removing the lat-lon stretching factor, - # but computing the spacings like this leads to a grid with exactly 1m spacing every one degree. - # We use this to test that the grid is rotated by θᵢ degrees from the vertical - # -> i.e. a θᵢ degree clockwise rotation. - fill!(grid.Δxᶜᶠᵃ, 1) - fill!(grid.Δyᶜᶠᵃ, 1) - angles = [-22.5, 30, 45, 60] for θᵢ in angles @@ -98,13 +91,18 @@ end cosθ = cosd(θᵢ) for i in 1:6, j in 1:6 - grid.λᶠᶠᵃ[i, j, 1] = (i-1) * cosθ + (j-1) * sinθ - grid.φᶠᶠᵃ[i, j, 1] = - (i-1) * sinθ + (j-1) * cosθ + grid.λᶠᶠᵃ[i, j, 1] = ( (i-1) * cosθ + (j-1) * sinθ) * 1e-4 + grid.φᶠᶠᵃ[i, j, 1] = (- (i-1) * sinθ + (j-1) * cosθ) * 1e-4 end λᶠᶠᵃ = grid.λᶠᶠᵃ φᶠᶠᵃ = grid.φᶠᶠᵃ + for i in 1:6, j in 1:6 + grid.Δxᶜᶠᵃ[i, j] = haversine((λᶠᶠᵃ[i+1, j], φᶠᶠᵃ[i+1, j]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), 1) + grid.Δyᶠᶜᵃ[i, j] = haversine((λᶠᶠᵃ[i, j+1], φᶠᶠᵃ[i, j+1]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), 1) + end + for i in 1:3, j in 1:3 cosθ, sinθ = rotation_matrix(i, j, grid) θ = atand(sinθ / cosθ) From 531e594ec78e8dea301d8b7f27e16528494df261 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 30 May 2025 12:53:08 +0300 Subject: [PATCH 16/18] rotation_matrix -> rotation_angle --- src/Operators/vector_rotation_operators.jl | 40 ++++++----- test/test_vector_rotation_operators.jl | 80 +++++++++++++--------- 2 files changed, 70 insertions(+), 50 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index a43de724df..00f901c408 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/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,7 +48,13 @@ _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) -@inline function rotation_matrix(i, j, grid::OrthogonalSphericalShellGrid) + +""" + rotation_angle(i, j, grid::OrthogonalSphericalShellGrid) + +Return the rotation angle (in degrees) of the `i, j`-th point of the `grid`. +""" +@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()) @@ -58,19 +64,21 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f Δ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) + Δ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 + # θ 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) - return Rcosθ / R, Rsinθ / R + cosθ, sinθ = Rcosθ / R, Rsinθ / R + + return atand(sinθ / cosθ) end # Intrinsic and extrinsic conversion for `OrthogonalSphericalShellGrid`s, @@ -82,14 +90,13 @@ end # 2D vectors @inline function intrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uₑ, vₑ) - cosθ, sinθ = rotation_matrix(i, j, grid) + θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid) + sinθ = sind(θ_degrees) + cosθ = cosd(θ_degrees) 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θ @@ -108,14 +115,13 @@ end # 2D vectors @inline function extrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uᵢ, vᵢ) - cosθ, sinθ = rotation_matrix(i, j, grid) - + θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid) + sinθ = sind(θ_degrees) + cosθ = cosd(θ_degrees) + 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θ diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 23bd83cb5e..2b7dc2f469 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -2,8 +2,8 @@ include("dependencies_for_runtests.jl") using Statistics: mean using Oceananigans.Operators -using Oceananigans.Operators: rotation_matrix -using Oceananigans.Grids: haversine +using Oceananigans.Operators: rotation_angle +using Distances: haversine # To be used in the test below as `KernelFunctionOperation`s @inline intrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) = @@ -78,40 +78,54 @@ function test_vector_rotation(grid) end @testset "Vector rotation" begin - if arch isa CPU - - # Build a custom grid that is rotated and test that the rotation matrix is correct. - # We build a grid that is rotated by θᵢ degrees clockwise from the vertical. - grid = OrthogonalSphericalShellGrid(; size=(4, 4, 1), z=(0, 1), radius=1, conformal_mapping=nothing) - - angles = [-22.5, 30, 45, 60] - - for θᵢ in angles - sinθ = sind(θᵢ) - cosθ = cosd(θᵢ) - - for i in 1:6, j in 1:6 - grid.λᶠᶠᵃ[i, j, 1] = ( (i-1) * cosθ + (j-1) * sinθ) * 1e-4 - grid.φᶠᶠᵃ[i, j, 1] = (- (i-1) * sinθ + (j-1) * cosθ) * 1e-4 - end - - λᶠᶠᵃ = grid.λᶠᶠᵃ - φᶠᶠᵃ = grid.φᶠᶠᵃ - - for i in 1:6, j in 1:6 - grid.Δxᶜᶠᵃ[i, j] = haversine((λᶠᶠᵃ[i+1, j], φᶠᶠᵃ[i+1, j]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), 1) - grid.Δyᶠᶜᵃ[i, j] = haversine((λᶠᶠᵃ[i, j+1], φᶠᶠᵃ[i, j+1]), (λᶠᶠᵃ[i, j], φᶠᶠᵃ[i, j]), 1) - end - - for i in 1:3, j in 1:3 - cosθ, sinθ = rotation_matrix(i, j, grid) - θ = atand(sinθ / cosθ) - @test θ ≈ θᵢ + for arch in archs + @testset "Rotation between Intrinsic to Extrinsic reference frame [$(typeof(arch))]" begin + + if arch isa CPU + @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 end end - end - 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" cubed_sphere_grid = ConformalCubedSphereGrid(arch; panel_size=(10, 10, 1), z=(-1, 0)) From f9f99d656adf245a8868cdca81c55808b13836d2 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 30 May 2025 13:41:43 +0300 Subject: [PATCH 17/18] add rotation_angle docs --- src/Operators/vector_rotation_operators.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/Operators/vector_rotation_operators.jl b/src/Operators/vector_rotation_operators.jl index 00f901c408..c4aa1e484f 100644 --- a/src/Operators/vector_rotation_operators.jl +++ b/src/Operators/vector_rotation_operators.jl @@ -53,6 +53,8 @@ _intrinsic_ coordinate systems are equivalent. However, for other grids (e.g., f 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) @@ -78,7 +80,8 @@ Return the rotation angle (in degrees) of the `i, j`-th point of the `grid`. cosθ, sinθ = Rcosθ / R, Rsinθ / R - return atand(sinθ / cosθ) + θ_degrees = atand(sinθ / cosθ) + return θ_degrees end # Intrinsic and extrinsic conversion for `OrthogonalSphericalShellGrid`s, @@ -90,13 +93,13 @@ end # 2D vectors @inline function intrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uₑ, vₑ) + u = getvalue(uₑ, i, j, k, grid) + v = getvalue(vₑ, i, j, k, grid) + θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid) sinθ = sind(θ_degrees) cosθ = cosd(θ_degrees) - u = getvalue(uₑ, i, j, k, grid) - v = getvalue(vₑ, i, j, k, grid) - uᵢ = u * cosθ - v * sinθ vᵢ = u * sinθ + v * cosθ @@ -115,13 +118,13 @@ end # 2D vectors @inline function extrinsic_vector(i, j, k, grid::OrthogonalSphericalShellGrid, uᵢ, vᵢ) + u = getvalue(uᵢ, i, j, k, grid) + v = getvalue(vᵢ, i, j, k, grid) + θ_degrees = rotation_angle(i, j, grid::OrthogonalSphericalShellGrid) sinθ = sind(θ_degrees) cosθ = cosd(θ_degrees) - u = getvalue(uᵢ, i, j, k, grid) - v = getvalue(vᵢ, i, j, k, grid) - uₑ = + u * cosθ + v * sinθ vₑ = - u * sinθ + v * cosθ From 300c7afe1faa5d7485309d61a160fe844cd5eaaa Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 30 May 2025 13:41:43 +0300 Subject: [PATCH 18/18] cleanup --- test/test_vector_rotation_operators.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl index 2b7dc2f469..c4ad2307ed 100644 --- a/test/test_vector_rotation_operators.jl +++ b/test/test_vector_rotation_operators.jl @@ -84,7 +84,7 @@ end if arch isa CPU @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 + # 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 @@ -130,8 +130,10 @@ end @info " Testing the conversion of a vector between the Intrinsic and Extrinsic reference frame" 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) + + for grid in (cubed_sphere_grid, tripolar_grid) + test_vector_rotation(grid) + end end end end