diff --git a/Project.toml b/Project.toml index ca9fc6e487..d81b10494b 100644 --- a/Project.toml +++ b/Project.toml @@ -78,10 +78,9 @@ julia = "1.9" [extras] DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" -OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimesDates = "bdfc003b-8df8-5c39-adcd-3a9087f5df4a" [targets] -test = ["DataDeps", "Enzyme", "OrthogonalSphericalShellGrids", "SafeTestsets", "Test", "TimesDates"] +test = ["DataDeps", "Enzyme", "SafeTestsets", "Test", "TimesDates"] diff --git a/test/runtests.jl b/test/runtests.jl index c78eeefbb6..816d8d578a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,6 +47,7 @@ CUDA.allowscalar() do @testset "Unit tests" begin include("test_grids.jl") include("test_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_cubed_spheres.jl b/test/test_cubed_spheres.jl deleted file mode 100644 index d8984cb0a1..0000000000 --- a/test/test_cubed_spheres.jl +++ /dev/null @@ -1,232 +0,0 @@ -include("dependencies_for_runtests.jl") -include("data_dependencies.jl") - -using Statistics: mean -using Oceananigans.Operators -using Oceananigans.CubedSpheres -using Oceananigans.Models.HydrostaticFreeSurfaceModels -using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField - -using OrthogonalSphericalShellGrids - -# To be used in the test below as `KernelFunctionOperation`s -@inline intrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) = - @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[1] - -@inline intrinsic_vector_y_component(i, j, k, grid, uₑ, vₑ) = - @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[2] - -@inline extrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) = - @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[1] - -@inline extrinsic_vector_y_component(i, j, k, grid, uₑ, vₑ) = - @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[2] - -function kinetic_energy(u, v) - ke = Field(0.5 * (u * u + v * v)) - return compute!(ke) -end - -function test_vector_rotation(grid) - u = XFaceField(grid) - v = YFaceField(grid) - - # Purely longitudinal flow in the extrinsic coordinate system - set!(u, 1) - set!(v, 0) - - # Convert it to an "Instrinsic" reference frame - uᵢ = KernelFunctionOperation{Face, Center, Center}(intrinsic_vector_x_component, grid, u, v) - vᵢ = KernelFunctionOperation{Center, Face, 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 maximum(uᵢ) ≈ maximum(vᵢ) - @test minimum(uᵢ) ≈ minimum(vᵢ) - @test mean(uᵢ) ≈ mean(vᵢ) - @test mean(uᵢ) > 0 # The mean value should be positive - - # Kinetic energy should remain the same - KE = kinetic_energy(uᵢ, vᵢ) - @test all(on_architecture(CPU(), interior(KE)) .≈ 0.5) - - # Convert it back to a purely zonal velocity (vₑ == 0) - uₑ = KernelFunctionOperation{Face, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ) - vₑ = KernelFunctionOperation{Center, Face, 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) - @test all(on_architecture(CPU(), interior(vₑ)) .≈ 0) - @test all(on_architecture(CPU(), interior(uₑ)) .≈ 1) - - # Purely meridional flow in the extrinsic coordinate system - set!(u, 0) - set!(v, 1) - - # Convert it to an "Instrinsic" reference frame - uᵢ = KernelFunctionOperation{Face, Center, Center}(intrinsic_vector_x_component, grid, u, v) - vᵢ = KernelFunctionOperation{Center, Face, 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 maximum(uᵢ) ≈ maximum(vᵢ) - @test minimum(uᵢ) ≈ minimum(vᵢ) - @test mean(uᵢ) ≈ mean(vᵢ) - @test mean(vᵢ) > 0 # The mean value should be positive - - # Kinetic energy should remain the same - KE = kinetic_energy(uᵢ, vᵢ) - @test all(on_architecture(CPU(), interior(KE)) .≈ 0.5) - - # Convert it back to a purely zonal velocity (vₑ == 0) - uₑ = KernelFunctionOperation{Face, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ) - vₑ = KernelFunctionOperation{Center, Face, 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) - @test all(on_architecture(CPU(), interior(vₑ)) .≈ 1) - @test all(on_architecture(CPU(), interior(uₑ)) .≈ 0) - - # Mixed zonal and meridional flow. - set!(u, 0.5) - set!(v, 0.5) - - # Convert it to an "Instrinsic" reference frame - uᵢ = KernelFunctionOperation{Face, Center, Center}(intrinsic_vector_x_component, grid, u, v) - vᵢ = KernelFunctionOperation{Center, Face, 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 maximum(uᵢ) ≈ maximum(vᵢ) - @test minimum(uᵢ) ≈ minimum(vᵢ) - @test mean(uᵢ) ≈ mean(vᵢ) - @test mean(vᵢ) > 0 # The mean value should be positive - - # Kinetic energy should remain the same - KE = kinetic_energy(uᵢ, vᵢ) - @test all(on_architecture(CPU(), interior(KE)) .≈ 0.25) - - # Convert it back to a purely zonal velocity (vₑ == 0) - uₑ = KernelFunctionOperation{Face, Center, Center}(extrinsic_vector_x_component, grid, uᵢ, vᵢ) - vₑ = KernelFunctionOperation{Center, Face, 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) - @test all(on_architecture(CPU(), interior(vₑ)) .≈ 0.5) - @test all(on_architecture(CPU(), interior(uₑ)) .≈ 0.5) -end - - -@testset "Cubed spheres" begin - - @testset "Conformal cubed sphere grid" begin - @info " Testing conformal cubed sphere grid..." - - grid = ConformalCubedSphereGrid(panel_size=(10, 10, 1), z=(-1, 0)) - @test try show(grid); println(); true; catch; false; end - end - - for arch in archs - - @info " Constructing a ConformalCubedSphereGrid from file [$(typeof(arch))]..." - - # These tests cause an undefined `Bound Access Error` on GPU's CI with the new CUDA version. - # The error is not reproducible neither on Tartarus nor on Sverdrup. - # These are excised for the moment (PR #2253) as Cubed sphere will be reworked - if !(arch isa GPU) - # Prototype grid and model for subsequent tests - cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2" - grid = OldConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0)) - - @info " Constructing a HydrostaticFreeSurfaceModel on a ConformalCubedSphereGrid [$(typeof(arch))]..." - - free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1) - model = HydrostaticFreeSurfaceModel(; grid, free_surface, - momentum_advection = VectorInvariant(), - coriolis = nothing, - closure = nothing, - tracers = :c, - buoyancy = nothing) - - @testset "Constructing a grid from file [$(typeof(arch))]" begin - @test grid isa ConformalCubedSphereGrid - 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" - trg_grid = TripolarGrid(arch, size = (20, 20, 1), z = (0, 1)) - - test_vector_rotation(grid) - test_vector_rotation(trg_grid) - end - - @testset "CubedSphereData and CubedSphereFields [$(typeof(arch))]" begin - @info " Testing CubedSphereData and CubedSphereFields [$(typeof(arch))]..." - c = model.tracers.c - η = model.free_surface.η - - set!(c, 0) - set!(η, 0) - - CUDA.allowscalar(true) - @test all(all(face_c .== 0) for face_c in faces(c)) - @test all(all(face_η .== 0) for face_η in faces(η)) - CUDA.allowscalar(false) - - @test maximum(abs, c) == 0 - @test minimum(abs, c) == 0 - @test mean(c) == 0 - - @test maximum(abs, η) == 0 - @test minimum(abs, η) == 0 - @test mean(η) == 0 - end - - @testset "Constructing a HydrostaticFreeSurfaceModel [$(typeof(arch))]" begin - @test model isa HydrostaticFreeSurfaceModel - end - - @testset "Time stepping a HydrostaticFreeSurfaceModel [$(typeof(arch))]" begin - @info " Time-stepping HydrostaticFreeSurfaceModel on a ConformalCubedSphereGrid [$(typeof(arch))]..." - time_step!(model, 1) - @test try time_step!(model, 1); true; catch; false; end - end - - @testset "VerticalVorticityField on ConformalCubedSphereGrid [$(typeof(arch))]" begin - @info " Testing VerticalVorticityField on a ConformalCubedSphereGrid [$(typeof(arch))]..." - ζ = VerticalVorticityField(model) - - @test ζ isa Field - - set!(model, u = (x, y, z) -> rand()) - - @test try - compute!(ζ) - true - catch err - println(sprint(showerror, err)) - false - end - @test maximum(abs, ζ) > 0 # fingers crossed - end - end - end -end diff --git a/test/test_vector_rotation_operators.jl b/test/test_vector_rotation_operators.jl new file mode 100644 index 0000000000..eeb61a8bb9 --- /dev/null +++ b/test/test_vector_rotation_operators.jl @@ -0,0 +1,178 @@ +include("dependencies_for_runtests.jl") + +using Statistics: mean +using Oceananigans.Operators + +# To be used in the test below as `KernelFunctionOperation`s +@inline intrinsic_vector_x_component(i, j, k, grid, uₑ, vₑ) = + @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[1] + +@inline intrinsic_vector_y_component(i, j, k, grid, uₑ, vₑ) = + @inbounds intrinsic_vector(i, j, k, grid, uₑ, vₑ)[2] + +@inline extrinsic_vector_x_component(i, j, k, grid, uᵢ, vᵢ) = + @inbounds extrinsic_vector(i, j, k, grid, uᵢ, vᵢ)[1] + +@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 +end + +function test_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 + +@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) + end + end +end