From aba90221b3e527a3fbab136bf0dda26cf1ce97a5 Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Fri, 8 Mar 2024 11:35:28 -0800 Subject: [PATCH] modified: test/callbacks.jl modified: test/rrtmgp_interface.jl modified: test/test_helpers.jl modified: test/topography.jl Reorganise modified: .buildkite/pipeline.yml renamed: test/gravity_wave/gw_plotutils.jl -> test/parameterized_tendencies/gravity_wave/gw_plotutils.jl renamed: test/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl -> test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl renamed: test/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl -> test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl renamed: test/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl -> test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl renamed: test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl -> test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl renamed: test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl -> test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl new file: test/parameters/parameter_tests.jl new file: test/parameters/parameter_tests.toml new file: sponge/rayleigh_sponge.jl modified: parameterized_tendencies/sponge/rayleigh_sponge.jl --- .buildkite/pipeline.yml | 10 +- test/callbacks.jl | 72 +++++++ .../gravity_wave/gw_plotutils.jl | 0 .../nogw_test_3d.jl | 0 .../nogw_test_mima.jl | 0 .../nogw_test_single_column.jl | 0 .../orographic_gravity_wave/ogwd_3d.jl | 0 .../orographic_gravity_wave/ogwd_baseflux.jl | 0 .../sponge/rayleigh_sponge.jl | 38 ++++ test/{ => parameters}/parameter_tests.jl | 0 test/{ => parameters}/parameter_tests.toml | 0 test/rrtmgp_interface.jl | 89 ++++++++ test/runtests.jl | 5 +- test/test_helpers.jl | 142 +++++++++++++ test/topography.jl | 60 ++++++ test/utilities.jl | 196 ++++++++---------- 16 files changed, 498 insertions(+), 114 deletions(-) create mode 100644 test/callbacks.jl rename test/{ => parameterized_tendencies}/gravity_wave/gw_plotutils.jl (100%) rename test/{ => parameterized_tendencies}/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl (100%) rename test/{ => parameterized_tendencies}/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl (100%) rename test/{ => parameterized_tendencies}/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl (100%) rename test/{ => parameterized_tendencies}/gravity_wave/orographic_gravity_wave/ogwd_3d.jl (100%) rename test/{ => parameterized_tendencies}/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl (100%) create mode 100644 test/parameterized_tendencies/sponge/rayleigh_sponge.jl rename test/{ => parameters}/parameter_tests.jl (100%) rename test/{ => parameters}/parameter_tests.toml (100%) create mode 100644 test/rrtmgp_interface.jl create mode 100644 test/test_helpers.jl create mode 100644 test/topography.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 56782a0af00..dcc16a5e5c2 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -96,27 +96,27 @@ steps: steps: - label: ":computer: non-orographic gravity wave parameterization unit test 3d" - command: "julia --color=yes --project=examples test/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl" + command: "julia --color=yes --project=examples test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl" artifact_paths: "nonorographic_gravity_wave_test_3d/*" agents: slurm_mem: 20GB - label: ":computer: non-orographic gravity wave parameterization test with MiMA output" - command: "julia --color=yes --project=examples test/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl" + command: "julia --color=yes --project=examples test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl" artifact_paths: "nonorographic_gravity_wave_test_mima/*" agents: slurm_mem: 20GB - label: ":computer: non-orographic gravity wave parameterization unit test single column" - command: "julia --color=yes --project=examples test/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl" + command: "julia --color=yes --project=examples test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl" artifact_paths: "nonorographic_gravity_wave_test_single_column/*" - label: ":computer: orographic gravity wave parameterization unit test for base flux calculation" - command: "julia --color=yes --project=examples test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl" + command: "julia --color=yes --project=examples test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl" artifact_paths: "orographic_gravity_wave_test_baseflux/*" - label: ":computer: orographic gravity wave parameterization unit test for 3d calculation" - command: "julia --color=yes --project=examples test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl" + command: "julia --color=yes --project=examples test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl" artifact_paths: "orographic_gravity_wave_test_3d/*" - label: ":computer: single column non-orographic gravity wave parameterization" diff --git a/test/callbacks.jl b/test/callbacks.jl new file mode 100644 index 00000000000..f30b3233717 --- /dev/null +++ b/test/callbacks.jl @@ -0,0 +1,72 @@ +using Test +import ClimaAtmos as CA +import SciMLBase as SMB + +testfun!() = π +cb_default = CA.call_every_n_steps(testfun!;) +test_nsteps = 999 +test_dt = 1 +test_tend = 999.0 + +cb_1 = CA.call_every_n_steps( + testfun!, + test_nsteps; + skip_first = false, + call_at_end = false, + condition = nothing, +) +cb_2 = + CA.call_every_dt(testfun!, test_dt; skip_first = false, call_at_end = false) +cb_3 = CA.callback_from_affect(cb_2.affect!) +cb_4 = CA.call_every_n_steps( + testfun!, + 3; + skip_first = false, + call_at_end = false, + condition = nothing, +) +cb_set = SMB.CallbackSet(cb_1, cb_2, cb_4) + +@testset "simple default callback" begin + @test cb_default.condition.n == 1 + @test cb_default.affect!.f!() == π +end + +# per n steps +@testset "every n-steps callback" begin + @test cb_1.initialize.skip_first == false + @test cb_1.condition.n == test_nsteps + @test cb_1.affect!.f!() == π + @test_throws AssertionError CA.call_every_n_steps( + testfun!, + Inf; + skip_first = false, + call_at_end = false, + ) +end + +# per dt interval +@testset "dt interval callback" begin + @test cb_2 isa SMB.DiscreteCallback + @test cb_2.affect!.dt == test_dt + @test cb_2.affect!.cb!.f!() == π + @test_throws AssertionError CA.call_every_dt( + testfun!, + Inf; + skip_first = false, + call_at_end = false, + ) +end + +@testset "atmos callbacks and callback sets" begin + # atmoscallbacks from discrete callbacks + @test cb_3.f!() == π + atmos_cbs = CA.atmos_callbacks(cb_set) + # test against expected callback outcomes + tspan = [0, test_tend] + @test CA.n_expected_calls(cb_set, test_dt, tspan)[2] == test_tend + @test CA.n_steps_per_cycle_per_cb(cb_set, test_dt) == + [test_nsteps; test_dt; 3] +end + +# query functions diff --git a/test/gravity_wave/gw_plotutils.jl b/test/parameterized_tendencies/gravity_wave/gw_plotutils.jl similarity index 100% rename from test/gravity_wave/gw_plotutils.jl rename to test/parameterized_tendencies/gravity_wave/gw_plotutils.jl diff --git a/test/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl similarity index 100% rename from test/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl rename to test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl diff --git a/test/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl similarity index 100% rename from test/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl rename to test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl diff --git a/test/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl similarity index 100% rename from test/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl rename to test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl diff --git a/test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl similarity index 100% rename from test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl rename to test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl diff --git a/test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl similarity index 100% rename from test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl rename to test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl diff --git a/test/parameterized_tendencies/sponge/rayleigh_sponge.jl b/test/parameterized_tendencies/sponge/rayleigh_sponge.jl new file mode 100644 index 00000000000..36f8686524e --- /dev/null +++ b/test/parameterized_tendencies/sponge/rayleigh_sponge.jl @@ -0,0 +1,38 @@ + +import ClimaAtmos as CA +import SurfaceFluxes as SF +import ClimaAtmos.Parameters as CAP +import ClimaCore as CC + +### Common Objects ### +@testset begin "Rayleigh-sponge functions" + ### Boilerplate default integrator objects + config = CA.AtmosConfig(Dict("initial_condition" => "DryBaroclinicWave")); + parsed_args = config.parsed_args + simulation = CA.get_simulation(config); + (;integrator) = simulation + Y = integrator.u; + zmax = maximum(CC.Fields.coordinate_field(Y.f).z) + z = CC.Fields.coordinate_field(Y.c).z + Y.c.uₕ.components.data.:1 .= ones(axes(Y.c)) + Y.c.uₕ.components.data.:2 .= ones(axes(Y.c)) + FT = eltype(Y); + ᶜYₜ = Y .* FT(0) + ### Component test begins here + rs = CA.RayleighSponge(; zd = FT(0), α_uₕ=FT(1), α_w=FT(1)) + (; ᶜβ_rayleigh_uₕ) = CA.rayleigh_sponge_cache(Y, rs) + @test ᶜβ_rayleigh_uₕ == @. sin(FT(π)/2 * z/zmax)^2 + test_cache = (;rayleigh_sponge = (;ᶜβ_rayleigh_uₕ = ᶜβ_rayleigh_uₕ)) + CC.Fields.bycolumn(axes(Y.c)) do colidx + CA.rayleigh_sponge_tendency!(ᶜYₜ,Y,test_cache, FT(0), colidx, rs) + end + # Test that only required tendencies are affected + for (var_name) in filter(x -> (x != :uₕ), propertynames(Y.c)) + @test ᶜYₜ.c.:($var_name) == zeros(axes(Y.c)) + end + for (var_name) in propertynames(Y.f) + @test ᶜYₜ.f.:($var_name) == zeros(axes(Y.f)) + end + @test ᶜYₜ.c.uₕ.components.data.:1 == -1 .* (ᶜβ_rayleigh_uₕ) + @test ᶜYₜ.c.uₕ.components.data.:2 == -1 .* (ᶜβ_rayleigh_uₕ) +end diff --git a/test/parameter_tests.jl b/test/parameters/parameter_tests.jl similarity index 100% rename from test/parameter_tests.jl rename to test/parameters/parameter_tests.jl diff --git a/test/parameter_tests.toml b/test/parameters/parameter_tests.toml similarity index 100% rename from test/parameter_tests.toml rename to test/parameters/parameter_tests.toml diff --git a/test/rrtmgp_interface.jl b/test/rrtmgp_interface.jl new file mode 100644 index 00000000000..64f53c4875c --- /dev/null +++ b/test/rrtmgp_interface.jl @@ -0,0 +1,89 @@ +using Test +using Random +Random.seed!(1234) +import ClimaAtmos as CA +import ClimaAtmos.RRTMGPInterface as RRTMGPI +import ClimaParams as CP +import RRTMGP.Parameters.RRTMGPParameters + +param_dict = CP.create_toml_dict(Float64) +params = RRTMGPParameters(param_dict) + +using Statistics + +include("test_helpers.jl") + +@testset "array <-> field" begin + (; bubble_space) = get_spherical_spaces() + test_field = ones(bubble_space) + idx = size(parent(test_field)) # IJFH layout + npts = idx[1] * idx[2] * idx[4] + @assert test_field isa Fields.Field + test_array = RRTMGPI.field2array(test_field) + @test test_array isa Array + @test length(test_array) == npts + new_field = RRTMGPI.array2field(test_array, bubble_space) + @test new_field == test_field + @test size(parent(new_field)) == idx +end + +@testset "arithmetic ops" begin + FloatType = (Float32, Float64) + for FT in FloatType + _T = [FT(1)] + _Tꜛ = [FT(1)] + _Tꜜ = [FT(2)] + _p = [FT(1 // 2)] + _p₁ = [FT(2)] + _p₂ = [FT(3)] + _Tₛ = [FT(5)] + XT = similar(_T) + XP = similar(_p) + @test RRTMGPI.uniform_z_p.(_T, _p₁, _Tꜜ, _p₁, _Tꜛ) == _p₁ + @test RRTMGPI.uniform_z_p.(_T, _p₁, _Tꜜ, _p₁, _Tꜛ) == _p₁ + + RRTMGPI.interp!(RRTMGPI.ArithmeticMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ) + @test XT[1] == (_Tꜜ[1] + _Tꜛ[1]) / 2 + @test XP[1] == (_p₁[1] + _p₂[1]) / 2 + + RRTMGPI.interp!(RRTMGPI.GeometricMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ) + @test XT[1] == sqrt(_Tꜜ[1] * _Tꜛ[1]) + @test XP[1] == sqrt(_p₁[1] * _p₂[1]) + + RRTMGPI.interp!(RRTMGPI.UniformZ(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ) + @test XT[1] == (_Tꜜ[1] + _Tꜛ[1]) / 2 + @test XP[1] == RRTMGPI.uniform_z_p.(XT, _p₁, _Tꜜ, _p₂, _Tꜛ)[] + + RRTMGPI.interp!(RRTMGPI.UniformP(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ) + @test XP[1] == (_p₁[1] + _p₂[1]) / 2 + # The second result entry has the same function form as uniform_z_p, + # with arguments swapped + @test XT[1] == RRTMGPI.uniform_z_p.(XP, _Tꜜ, _p₁, _Tꜛ, _p₂)[] + + RRTMGPI.extrap!(RRTMGPI.ArithmeticMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + @test XP[1] == @. (3 * _p₁[1] - _p₂[1]) / 2 + @test XT[1] == @. (3 * _Tꜜ[1] - _Tꜛ[1]) / 2 + + RRTMGPI.extrap!(RRTMGPI.GeometricMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + @test XP[1] == @. sqrt(_p₁[1]^3 /_p₂[1]) + @test XT[1] == @. sqrt(_Tꜜ[1]^3 /_Tꜛ[1]) + + RRTMGPI.extrap!(RRTMGPI.ArithmeticMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + _XP = XP + RRTMGPI.extrap!(RRTMGPI.UniformZ(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + @test _XP[1] == XP[1] + # The result entry has the same function form as uniform_z_p, + # with arguments swapped + @test XT[1] == RRTMGPI.uniform_z_p.(XP, _Tꜜ, _p₁, _Tꜛ, _p₂)[] + + RRTMGPI.extrap!(RRTMGPI.ArithmeticMean(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + _XT = XT + RRTMGPI.extrap!(RRTMGPI.UniformP(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + @test XP[1] == RRTMGPI.uniform_z_p.(XT, _p₁, _Tꜜ, _p₂, _Tꜛ)[] + @test _XT[1] == XT[1] + + # RRTMGPI.extrap!(RRTMGPI.UseSurfaceTempAtBottom(), XP, XT, _p₁, _Tꜜ, _p₂, _Tꜛ, _Tₛ, params) + # @test XT[1] == _Tₛ[1] + + end +end diff --git a/test/runtests.jl b/test/runtests.jl index eb5af125b6b..0d9d79011ff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,11 +9,14 @@ using Test #! format: off @safetestset "Aqua" begin @time include("aqua.jl") end +@safetestset "Callbacks" begin @time include("callbacks.jl") end @safetestset "Utilities" begin @time include("utilities.jl") end -@safetestset "Parameter tests" begin @time include("parameter_tests.jl") end +@safetestset "Parameter tests" begin @time include("parameters/parameter_tests.jl") end @safetestset "Coupler Compatibility" begin @time include("coupler_compatibility.jl") end @safetestset "Configuration tests" begin @time include("config.jl") end @safetestset "surface albedo tests" begin @time include("surface_albedo.jl") end +@safetestset "Radiation interface tests" begin @time include("rrtmgp_interface.jl") end +@safetestset "Topography tests" begin @time include("topography.jl") end #! format: on diff --git a/test/test_helpers.jl b/test/test_helpers.jl new file mode 100644 index 00000000000..b193594872a --- /dev/null +++ b/test/test_helpers.jl @@ -0,0 +1,142 @@ + +### BoilerPlate Code +using ClimaComms +using IntervalSets + +import ClimaCore: + ClimaCore, + Domains, + Geometry, + Grids, + Fields, + Operators, + Meshes, + Spaces, + Quadratures, + Topologies, + Hypsography + +function get_spherical_spaces(; FT = Float32) + context = ClimaComms.SingletonCommsContext() + radius = FT(10π) + ne = 4 + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, ne) + topology = Topologies.Topology2D(context, mesh) + quad = Quadratures.GLL{Nq}() + space = Spaces.SpectralElementSpace2D(topology, quad) + enable_bubble = false + no_bubble_space = + Spaces.SpectralElementSpace2D(topology, quad; enable_bubble) + + # Now check constructor with bubble enabled + enable_bubble = true + bubble_space = Spaces.SpectralElementSpace2D(topology, quad; enable_bubble) + + lat = Fields.coordinate_field(bubble_space).lat + long = Fields.coordinate_field(bubble_space).long + coords = Fields.coordinate_field(bubble_space) + return (; + bubble_space = bubble_space, + no_bubble_space = no_bubble_space, + lat = lat, + long = long, + coords = coords, + FT = FT, + ) +end + +function get_cartesian_spaces(; FT = Float32) + xlim = (FT(0), FT(π)) + zlim = (FT(0), FT(π)) + helem = 5 + velem = 10 + npoly = 5 + ndims = 3 + stretch = Meshes.Uniform() + device = ClimaComms.CPUSingleThreaded() + comms_context = ClimaComms.SingletonCommsContext(device) + # Horizontal Grid Construction + quad = Quadratures.GLL{npoly + 1}() + horzdomain = Domains.RectangleDomain( + Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]), + Geometry.YPoint{FT}(xlim[1]) .. Geometry.YPoint{FT}(xlim[2]), + x1periodic = true, + x2periodic = true, + ) + # Assume same number of elems (helem) in (x,y) directions + horzmesh = Meshes.RectilinearMesh(horzdomain, helem, helem) + horz_topology = Topologies.Topology2D( + comms_context, + horzmesh, + Topologies.spacefillingcurve(horzmesh), + ) + h_space = + Spaces.SpectralElementSpace2D(horz_topology, quad, enable_bubble = true) + + horz_grid = Spaces.grid(h_space) + + # Vertical Grid Construction + vertdomain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(zlim[1]), + Geometry.ZPoint{FT}(zlim[2]); + boundary_names = (:bottom, :top), + ) + vertmesh = Meshes.IntervalMesh(vertdomain, stretch, nelems = velem) + vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) + vert_topology = Topologies.IntervalTopology( + ClimaComms.SingletonCommsContext(device), + vertmesh, + ) + vert_grid = Grids.FiniteDifferenceGrid(vert_topology) + ArrayType = ClimaComms.array_type(device) + grid = Grids.ExtrudedFiniteDifferenceGrid(horz_grid, vert_grid) + cent_space = Spaces.CenterExtrudedFiniteDifferenceSpace(grid) + face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(grid) + return ( + helem = helem, + cent_space = cent_space, + face_space = face_space, + xlim = xlim, + zlim = zlim, + velem = velem, + npoly = npoly, + quad = quad, + ) +end + +function get_coords(cent_space, face_space) + ccoords = Fields.coordinate_field(cent_space) + fcoords = Fields.coordinate_field(face_space) + return ccoords, fcoords +end + +function taylor_green_ic(coords) + u = @. sin(coords.x) * cos(coords.y) * cos(coords.z) + v = @. -cos(coords.x) * sin(coords.y) * cos(coords.z) + #TODO: If a w field is introduced include it here. + return u, v, u .* 0 +end + +function get_test_functions(cent_space, face_space) + ccoords, fcoords = get_coords(cent_space, face_space) + FT = eltype(ccoords) + Q = zero.(ccoords.x) + # Exact velocity profiles + u, v, w = taylor_green_ic(ccoords) + ᶠu, ᶠv, ᶠw = taylor_green_ic(fcoords) + (; x, y, z) = ccoords + UVW = Geometry.UVWVector + # Assemble (Cartesian) velocity + ᶜu = @. UVW(Geometry.UVector(u)) + + UVW(Geometry.VVector(v)) + + UVW(Geometry.WVector(w)) + ᶠu = @. UVW(Geometry.UVector(ᶠu)) + + UVW(Geometry.VVector(ᶠv)) + + UVW(Geometry.WVector(ᶠw)) + # Get covariant components + uₕ = @. Geometry.Covariant12Vector(ᶜu) + uᵥ = @. Geometry.Covariant3Vector(ᶠu) + return uₕ, uᵥ +end diff --git a/test/topography.jl b/test/topography.jl new file mode 100644 index 00000000000..c6f43929ee0 --- /dev/null +++ b/test/topography.jl @@ -0,0 +1,60 @@ +using Test +using Random +import ClimaAtmos as CA +Random.seed!(1234) + +### BoilerPlate Code +using ClimaComms +using IntervalSets + +import ClimaCore: + ClimaCore, + Domains, + Geometry, + Grids, + Fields, + Operators, + Meshes, + Spaces, + Quadratures, + Topologies, + Hypsography + +function get_spherical_spaces(; FT = Float32) + context = ClimaComms.SingletonCommsContext() + radius = FT(10π) + ne = 4 + Nq = 4 + domain = Domains.SphereDomain(radius) + mesh = Meshes.EquiangularCubedSphere(domain, ne) + topology = Topologies.Topology2D(context, mesh) + quad = Quadratures.GLL{Nq}() + space = Spaces.SpectralElementSpace2D(topology, quad) + enable_bubble = false + no_bubble_space = + Spaces.SpectralElementSpace2D(topology, quad; enable_bubble) + + # Now check constructor with bubble enabled + enable_bubble = true + bubble_space = Spaces.SpectralElementSpace2D(topology, quad; enable_bubble) + + lat = Fields.coordinate_field(bubble_space).lat + long = Fields.coordinate_field(bubble_space).long + coords = Fields.coordinate_field(bubble_space) + return (; + bubble_space = bubble_space, + no_bubble_space = no_bubble_space, + lat = lat, + long = long, + coords = coords, + FT = FT, + ) +end + +@testset "test topography functions" begin + (; FT, coords) = get_spherical_spaces() + @test extrema(CA.topography_dcmip200(coords)) == (FT(0), FT(2000)) + loc = findmax(parent(CA.topography_dcmip200(coords))) + @test parent(coords.lat)[loc[2]] == FT(0) + @test parent(coords.long)[loc[2]] == FT(-90) +end diff --git a/test/utilities.jl b/test/utilities.jl index 13112a15f83..847c38ba42f 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -3,111 +3,7 @@ using Random Random.seed!(1234) import ClimaAtmos as CA -### BoilerPlate Code -using ClimaComms -using IntervalSets - -import ClimaCore: - ClimaCore, - Domains, - Geometry, - Grids, - Fields, - Operators, - Meshes, - Spaces, - Quadratures, - Topologies, - Hypsography - - -FT = Float32 - -function get_cartesian_spaces(; FT = Float32) - xlim = (FT(0), FT(π)) - zlim = (FT(0), FT(π)) - helem = 5 - velem = 10 - npoly = 5 - ndims = 3 - stretch = Meshes.Uniform() - device = ClimaComms.CPUSingleThreaded() - comms_context = ClimaComms.SingletonCommsContext(device) - # Horizontal Grid Construction - quad = Quadratures.GLL{npoly + 1}() - horzdomain = Domains.RectangleDomain( - Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]), - Geometry.YPoint{FT}(xlim[1]) .. Geometry.YPoint{FT}(xlim[2]), - x1periodic = true, - x2periodic = true, - ) - # Assume same number of elems (helem) in (x,y) directions - horzmesh = Meshes.RectilinearMesh(horzdomain, helem, helem) - horz_topology = Topologies.Topology2D( - comms_context, - horzmesh, - Topologies.spacefillingcurve(horzmesh), - ) - h_space = - Spaces.SpectralElementSpace2D(horz_topology, quad, enable_bubble = true) - - horz_grid = Spaces.grid(h_space) - - # Vertical Grid Construction - vertdomain = Domains.IntervalDomain( - Geometry.ZPoint{FT}(zlim[1]), - Geometry.ZPoint{FT}(zlim[2]); - boundary_names = (:bottom, :top), - ) - vertmesh = Meshes.IntervalMesh(vertdomain, stretch, nelems = velem) - vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) - vert_topology = Topologies.IntervalTopology( - ClimaComms.SingletonCommsContext(device), - vertmesh, - ) - vert_grid = Grids.FiniteDifferenceGrid(vert_topology) - ArrayType = ClimaComms.array_type(device) - grid = Grids.ExtrudedFiniteDifferenceGrid(horz_grid, vert_grid) - cent_space = Spaces.CenterExtrudedFiniteDifferenceSpace(grid) - face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(grid) - return (helem = helem, cent_space = cent_space, face_space = face_space) -end - -function get_coords(cent_space, face_space) - ccoords = Fields.coordinate_field(cent_space) - fcoords = Fields.coordinate_field(face_space) - return ccoords, fcoords -end - -function taylor_green_ic(coords) - u = @. sin(coords.x) * cos(coords.y) * cos(coords.z) - v = @. -cos(coords.x) * sin(coords.y) * cos(coords.z) - #TODO: If a w field is introduced include it here. - return u, v, u .* 0 -end - -function get_test_functions(cent_space, face_space) - ccoords, fcoords = get_coords(cent_space, face_space) - FT = eltype(ccoords) - Q = zero.(ccoords.x) - # Exact velocity profiles - u, v, w = taylor_green_ic(ccoords) - ᶠu, ᶠv, ᶠw = taylor_green_ic(fcoords) - (; x, y, z) = ccoords - UVW = Geometry.UVWVector - # Assemble (Cartesian) velocity - ᶜu = @. UVW(Geometry.UVector(u)) + - UVW(Geometry.VVector(v)) + - UVW(Geometry.WVector(w)) - ᶠu = @. UVW(Geometry.UVector(ᶠu)) + - UVW(Geometry.VVector(ᶠv)) + - UVW(Geometry.WVector(ᶠw)) - # Get covariant components - uₕ = @. Geometry.Covariant12Vector(ᶜu) - uᵥ = @. Geometry.Covariant3Vector(ᶠu) - return uₕ, uᵥ -end - +include("test_helpers.jl") @testset "sort_hdf5_files" begin day_sec(t) = @@ -157,6 +53,7 @@ end # Test compute_strain_rate_face (; helem, cent_space, face_space) = get_cartesian_spaces() ccoords, fcoords = get_coords(cent_space, face_space) + FT = eltype(ccoords.x) UVW = Geometry.UVWVector C123 = Geometry.Covariant123Vector # Alloc scratch space @@ -269,7 +166,90 @@ end @test Fields.field_values( CA.g³³_field(Fields.coordinate_field(cent_space).x), ) == lg_g³³ - @test maximum(abs.(lg_g³³ .- CA.g³³.(lg_gⁱʲ).components.data.:1)) == FT(0) - @test maximum(abs.(CA.g³ʰ.(lg_gⁱʲ).components.data.:1)) == FT(0) - @test maximum(abs.(CA.g³ʰ.(lg_gⁱʲ).components.data.:2)) == FT(0) + @test maximum(abs.(lg_g³³ .- CA.g³³.(lg_gⁱʲ).components.data.:1)) == 0 + @test maximum(abs.(CA.g³ʰ.(lg_gⁱʲ).components.data.:1)) == 0 + @test maximum(abs.(CA.g³ʰ.(lg_gⁱʲ).components.data.:2)) == 0 +end + +@testset "interval domain" begin + # Interval Spaces + (; zlim, velem) = get_cartesian_spaces() + line_mesh = CA.periodic_line_mesh(; x_max = zlim[2], x_elem = velem) + @test line_mesh isa Meshes.IntervalMesh + @test Geometry.XPoint(zlim[1]) == Meshes.domain(line_mesh).coord_min + @test Geometry.XPoint(zlim[2]) == Meshes.domain(line_mesh).coord_max + @test velem == Meshes.nelements(line_mesh) +end + +@testset "periodic rectangle meshes (spectral elements)" begin + # Interval Spaces + (; xlim, zlim, velem, helem, npoly) = get_cartesian_spaces() + rectangle_mesh = CA.periodic_rectangle_mesh(; + x_max = xlim[2], + y_max = xlim[2], + x_elem = helem, + y_elem = helem, + ) + @test rectangle_mesh isa Meshes.RectilinearMesh + @test Meshes.domain(rectangle_mesh) isa Meshes.RectangleDomain + @test Meshes.nelements(rectangle_mesh) == helem^2 + @test Meshes.element_horizontal_length_scale(rectangle_mesh) == + eltype(xlim)(π / npoly) + @test Meshes.elements(rectangle_mesh) == CartesianIndices((helem, helem)) +end + +@testset "make horizontal spaces" begin + + (; xlim, zlim, velem, helem, npoly, quad) = get_cartesian_spaces() + device = ClimaComms.CPUSingleThreaded() + comms_ctx = ClimaComms.context(device) + FT = eltype(xlim) + # 1D Space + line_mesh = CA.periodic_line_mesh(; x_max = zlim[2], x_elem = velem) + @test line_mesh isa Meshes.AbstractMesh1D + horz_plane_space = + CA.make_horizontal_space(line_mesh, quad, comms_ctx, true) + @test Spaces.column(horz_plane_space, 1, 1) isa Spaces.PointSpace + + # 2D Space + rectangle_mesh = CA.periodic_rectangle_mesh(; + x_max = xlim[2], + y_max = xlim[2], + x_elem = helem, + y_elem = helem, + ) + @test rectangle_mesh isa Meshes.AbstractMesh2D + horz_plane_space = + CA.make_horizontal_space(rectangle_mesh, quad, comms_ctx, true) + @test Spaces.nlevels(horz_plane_space) == 1 + @test Spaces.node_horizontal_length_scale(horz_plane_space) == + FT(π / npoly / 5) + @test Spaces.column(horz_plane_space, 1, 1, 1) isa Spaces.PointSpace +end + +@testset "make hybrid spaces" begin + (; cent_space, face_space, xlim, zlim, velem, helem, npoly, quad) = + get_cartesian_spaces() + device = ClimaComms.CPUSingleThreaded() + comms_ctx = ClimaComms.context(device) + z_stretch = Meshes.Uniform() + rectangle_mesh = CA.periodic_rectangle_mesh(; + x_max = xlim[2], + y_max = xlim[2], + x_elem = helem, + y_elem = helem, + ) + horz_plane_space = + CA.make_horizontal_space(rectangle_mesh, quad, comms_ctx, true) + test_cent_space, test_face_space = CA.make_hybrid_spaces( + horz_plane_space, + zlim[2], + velem, + z_stretch; + surface_warp = nothing, + topo_smoothing = false, + deep = false, + ) + @test test_cent_space == cent_space + @test test_face_space == face_space end