diff --git a/test/Spaces/unit_dss.jl b/test/Spaces/unit_dss.jl new file mode 100644 index 0000000000..41c47dc699 --- /dev/null +++ b/test/Spaces/unit_dss.jl @@ -0,0 +1,131 @@ +#= +julia --project +ENV["CLIMACOMMS_DEVICE"] = "CPU"; +using Revise; include(joinpath("test", "Spaces", "unit_dss.jl")) +=# +using Test +using ClimaComms +using Random +ClimaComms.@import_required_backends + +import ClimaCore: + Domains, + Fields, + Geometry, + Meshes, + Operators, + Spaces, + Quadratures, + Topologies, + DataLayouts + +function get_space_cs(::Type{FT}; context, R = 300.0) where {FT} + domain = Domains.SphereDomain{FT}(300.0) + mesh = Meshes.EquiangularCubedSphere(domain, 3) + topology = Topologies.Topology2D(context, mesh) + quad = Quadratures.GLL{4}() + space = Spaces.SpectralElementSpace2D(topology, quad) + return space +end + +function one_to_n(a::AbstractArray) + _a = Array(a) + Random.seed!(1234) + for i in 1:length(_a) + _a[i] = rand() + end + return typeof(a)(_a) +end + +function test_dss_count(f::Fields.Field, buff::Topologies.DSSBuffer, nc) + parent(f) .= one_to_n(parent(f)) + @test allunique(parent(f)) + cf = copy(f) + Spaces.weighted_dss!(f => buff) + n_dss_unaffected = count(parent(f) .== parent(cf)) + n_dss_affected = length(parent(f)) - n_dss_unaffected + return (; n_dss_affected) +end + +function get_space_and_buffers(::Type{FT}; context) where {FT} + init_state_covariant12(local_geometry, p) = + Geometry.Covariant12Vector(1.0, -1.0) + init_state_covariant123(local_geometry, p) = + Geometry.Covariant123Vector(1.0, -1.0, 1.0) + init_state_covariant3(local_geometry, p) = Geometry.Covariant3Vector(1.0) + + R = FT(6.371229e6) + npoly = 2 + z_max = FT(30e3) + z_elem = 3 + h_elem = 2 + device = ClimaComms.device(context) + @info "running dss-Covariant123Vector test on $(device)" h_elem z_elem npoly R z_max FT + # horizontal space + domain = Domains.SphereDomain{FT}(R) + horizontal_mesh = Meshes.EquiangularCubedSphere(domain, h_elem) + horizontal_topology = Topologies.Topology2D( + context, + horizontal_mesh, + Topologies.spacefillingcurve(horizontal_mesh), + ) + quad = Quadratures.GLL{npoly + 1}() + h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad) + # vertical space + z_domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(zero(z_max)), + Geometry.ZPoint{FT}(z_max); + boundary_names = (:bottom, :top), + ) + z_mesh = Meshes.IntervalMesh(z_domain, nelems = z_elem) + z_topology = Topologies.IntervalTopology(context, z_mesh) + z_center_space = Spaces.CenterFiniteDifferenceSpace(z_topology) + space = Spaces.ExtrudedFiniteDifferenceSpace(h_space, z_center_space) + args = (Fields.local_geometry_field(space), Ref(nothing)) + y12 = init_state_covariant12.(args...) + y123 = init_state_covariant123.(args...) + y3 = init_state_covariant3.(args...) + dss_buffer = (; + y12 = Spaces.create_dss_buffer(y12), + y123 = Spaces.create_dss_buffer(y123), + y3 = Spaces.create_dss_buffer(y3), + ) + return (; space, y12, y123, y3, dss_buffer) +end + +@testset "DSS of AxisTensors on Cubed Sphere" begin + FT = Float64 + device = ClimaComms.device() + nt = get_space_and_buffers(FT; context = ClimaComms.context(device)) + + # test DSS for a Covariant12Vector + # ensure physical velocity is continuous across SE boundary for initial state + n_dss_affected_y12 = + test_dss_count(nt.y12, nt.dss_buffer.y12, 2).n_dss_affected + n_dss_affected_y123 = + test_dss_count(nt.y123, nt.dss_buffer.y123, 3).n_dss_affected + n_dss_affected_y3 = + test_dss_count(nt.y3, nt.dss_buffer.y3, 1).n_dss_affected + + @test n_dss_affected_y12 * 3 / 2 == + n_dss_affected_y123 == + n_dss_affected_y3 * 3 + + @test nt.dss_buffer.y12.scalarfidx == Int[] + @test nt.dss_buffer.y12.covariant12fidx == Int[1] + @test nt.dss_buffer.y12.contravariant12fidx == Int[] + @test nt.dss_buffer.y12.covariant123fidx == Int[] + @test nt.dss_buffer.y12.contravariant123fidx == Int[] + + @test nt.dss_buffer.y123.scalarfidx == Int[] + @test nt.dss_buffer.y123.covariant12fidx == Int[] + @test nt.dss_buffer.y123.contravariant12fidx == Int[] + @test nt.dss_buffer.y123.covariant123fidx == Int[1] + @test nt.dss_buffer.y123.contravariant123fidx == Int[] + + @test nt.dss_buffer.y3.scalarfidx == Int[1] + @test nt.dss_buffer.y3.covariant12fidx == Int[] + @test nt.dss_buffer.y3.contravariant12fidx == Int[] + @test nt.dss_buffer.y3.covariant123fidx == Int[] + @test nt.dss_buffer.y3.contravariant123fidx == Int[] +end diff --git a/test/Topologies/unit_dss_transform.jl b/test/Topologies/unit_dss_transform.jl new file mode 100644 index 0000000000..57c6031987 --- /dev/null +++ b/test/Topologies/unit_dss_transform.jl @@ -0,0 +1,89 @@ +#= +julia --project +using Revise; include(joinpath("test", "Topologies", "unit_dss_transform.jl")) +=# +using Test +using ClimaComms +using ClimaCore.Topologies: dss_transform, dss_untransform +using Random +ClimaComms.@import_required_backends + +import ClimaCore: + Domains, + Fields, + Geometry, + Meshes, + Operators, + Spaces, + Quadratures, + Topologies, + DataLayouts + +function get_space(::Type{FT}; context) where {FT} + R = FT(6.371229e6) + npoly = 2 + z_max = FT(30e3) + z_elem = 3 + h_elem = 2 + device = ClimaComms.device(context) + @info "running dss-Covariant123Vector test on $(device)" h_elem z_elem npoly R z_max FT + # horizontal space + domain = Domains.SphereDomain{FT}(R) + horizontal_mesh = Meshes.EquiangularCubedSphere(domain, h_elem) + horizontal_topology = Topologies.Topology2D( + context, + horizontal_mesh, + Topologies.spacefillingcurve(horizontal_mesh), + ) + quad = Quadratures.GLL{npoly + 1}() + h_space = Spaces.SpectralElementSpace2D(horizontal_topology, quad) + # vertical space + z_domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(zero(z_max)), + Geometry.ZPoint{FT}(z_max); + boundary_names = (:bottom, :top), + ) + z_mesh = Meshes.IntervalMesh(z_domain, nelems = z_elem) + z_topology = Topologies.IntervalTopology(context, z_mesh) + z_center_space = Spaces.CenterFiniteDifferenceSpace(z_topology) + space = Spaces.ExtrudedFiniteDifferenceSpace(h_space, z_center_space) + return space +end + +@testset "dss_transform" begin + device = ClimaComms.device() + space = get_space(Float64; context = ClimaComms.context(device)) + + local_geometry = Fields.local_geometry_field(space) + map(local_geometry) do lg + FT = Geometry.undertype(typeof(lg)) + (; lat, long, z) = lg.coordinates + # Test that vertical component is treated as a scalar: + + arg = Geometry.Covariant123Vector(FT(lat), FT(long), FT(z)) + weight = 2 + dss_t = dss_transform(arg, lg, weight) + dss_ut = dss_untransform(Geometry.Covariant123Vector{FT}, dss_t, lg) + @test dss_t isa Geometry.UVWVector + @test typeof(arg) == typeof(dss_ut) + @test arg ≈ dss_ut / weight + + arg = Geometry.Covariant12Vector(FT(lat), FT(long)) + weight = 2 + dss_t = dss_transform(arg, lg, weight) + dss_ut = dss_untransform(Geometry.Covariant12Vector{FT}, dss_t, lg) + @test dss_t isa Geometry.UVWVector + @test typeof(arg) == typeof(dss_ut) + @test arg ≈ dss_ut / weight + + arg = Geometry.Covariant3Vector(FT(z)) + weight = 2 + dss_t = dss_transform(arg, lg, weight) + dss_ut = dss_untransform(Geometry.Covariant3Vector{FT}, dss_t, lg) + @test typeof(arg) == typeof(dss_ut) + @test dss_t isa Geometry.Covariant3Vector + @test dss_t === arg * weight + @test arg == dss_ut / weight + FT(1) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 136f5a4dfb..d8734d8ad3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,8 +30,10 @@ UnitTest("Rectangle topology" ,"Topologies/rectangle.jl"), UnitTest("Rectangle surface topology" ,"Topologies/rectangle_sfc.jl"), UnitTest("Cubedsphere topology" ,"Topologies/cubedsphere.jl"), UnitTest("Cubedsphere surface topology" ,"Topologies/cubedsphere_sfc.jl"), +UnitTest("dss_transform" ,"Topologies/unit_dss_transform.jl"), UnitTest("Quadratures" ,"Quadratures/Quadratures.jl"), UnitTest("Spaces" ,"Spaces/unit_spaces.jl"), +UnitTest("dss" ,"Spaces/unit_dss.jl"), UnitTest("Spaces - serial CPU DSS" ,"Spaces/ddss1.jl"), UnitTest("Spaces - DSS cubed sphere" ,"Spaces/ddss1_cs.jl"), UnitTest("Sphere spaces" ,"Spaces/sphere.jl"),