diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index cc8dfa699d..47deff0fc9 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/output_active/*" 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/output_active/*" 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/output_active/*" - 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/output_active/*" - 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/output_active/*" - label: ":computer: single column non-orographic gravity wave parameterization" diff --git a/test/callbacks.jl b/test/callbacks.jl new file mode 100644 index 0000000000..f30b323371 --- /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 99% rename from test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl rename to test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl index 69185cd750..acebaf4310 100644 --- a/test/gravity_wave/orographic_gravity_wave/ogwd_3d.jl +++ b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_3d.jl @@ -11,7 +11,9 @@ using ClimaCoreTempestRemap using Interpolations const FT = Float64 -include("../../../post_processing/remap/remap_helpers.jl") +include( + joinpath(pkgdir(ClimaAtmos), "post_processing/remap", "remap_helpers.jl"), +) include("../gw_plotutils.jl") comms_ctx = ClimaComms.SingletonCommsContext() 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 97% rename from test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl rename to test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl index 064bafb2e7..4523e5352f 100644 --- a/test/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl +++ b/test/parameterized_tendencies/gravity_wave/orographic_gravity_wave/ogwd_baseflux.jl @@ -8,7 +8,9 @@ using ClimaCoreTempestRemap const FT = Float64 include("../gw_plotutils.jl") -include("../../../post_processing/remap/remap_helpers.jl") +include( + joinpath(pkgdir(ClimaAtmos), "post_processing/remap", "remap_helpers.jl"), +) comms_ctx = ClimaComms.SingletonCommsContext() diff --git a/test/parameterized_tendencies/microphysics/precipitation.jl b/test/parameterized_tendencies/microphysics/precipitation.jl new file mode 100644 index 0000000000..2c50da307b --- /dev/null +++ b/test/parameterized_tendencies/microphysics/precipitation.jl @@ -0,0 +1,110 @@ + +import ClimaAtmos as CA +import SurfaceFluxes as SF +import ClimaAtmos.Parameters as CAP +import ClimaCore as CC +import Thermodynamics as TD +import CloudMicrophysics as CM + +include("../../test_helpers.jl") + +### Common Objects ### +@testset begin + "Precipitation tendency functions" + ### Boilerplate default integrator objects + config = CA.AtmosConfig( + Dict( + "initial_condition" => "PrecipitatingColumn", + "moist" => "nonequil", + "precip_model" => "0M", + "config" => "column", + "output_default_diagnostics" => false, + ), + ) + (; Y, p, params) = generate_test_simulation(config) + + FT = eltype(Y) + ᶜYₜ = zero(Y) + ### Component test begins here + + @info "0M Scheme" + ### 0-Moment Scheme + precip_model = CA.Microphysics0Moment() + precip_cache = CA.precipitation_cache(Y, precip_model) + # Test cache to verify expected variables exist in tendency function + test_varnames = ( + :ᶜS_ρq_tot, + :ᶜ3d_rain, + :ᶜ3d_snow, + :col_integrated_rain, + :col_integrated_snow, + ) + for var_name in test_varnames + @test var_name ∈ propertynames(precip_cache) + end + turbconv_model = nothing # Extend to other turbulence convection models + CC.Fields.bycolumn(axes(Y.c)) do colidx + CA.compute_precipitation_cache!( + Y, + p, + colidx, + precip_model, + turbconv_model, + ) + end + @test maximum(abs.(p.precipitation.ᶜS_ρq_tot)) <= sqrt(eps(FT)) + + # Test that tendencies result in correct water-mass budget, + # and that the tendency modification corresponds exactly to the + # cached source term. + CC.Fields.bycolumn(axes(Y.c)) do colidx + CA.precipitation_tendency!( + ᶜYₜ, + Y, + p, + FT(0), + colidx, + precip_model, + turbconv_model, + ) + end + @test ᶜYₜ.c.ρ == ᶜYₜ.c.ρq_tot + @test ᶜYₜ.c.ρ == p.precipitation.ᶜS_ρq_tot + + ### 1-Moment Scheme + @info "1M Scheme" + config = CA.AtmosConfig( + Dict( + "initial_condition" => "PrecipitatingColumn", + "moist" => "nonequil", + "precip_model" => "1M", + "config" => "column", + "output_default_diagnostics" => false, + ), + ) + (; Y, p, params) = generate_test_simulation(config) + precip_model = CA.Microphysics1Moment() + (; turbconv_model) = p.atmos + precip_cache = CA.precipitation_cache(Y, precip_model) + ᶜYₜ = Y .* FT(0) + test_varnames = (:ᶜSqₜᵖ, :ᶜSqᵣᵖ, :ᶜSqₛᵖ, :ᶜSeₜᵖ) + for var_name in test_varnames + @test var_name ∈ propertynames(precip_cache) + end + @test CA.qₚ(FT(10), FT(2)) == FT(5) + @test CA.qₚ(FT(-10), FT(2)) == FT(0) + @test CA.limit(FT(10), FT(2)) == FT(1) + CC.Fields.bycolumn(axes(Y.c)) do colidx + CA.precipitation_tendency!( + ᶜYₜ, + Y, + p, + FT(0), + colidx, + precip_model, + turbconv_model, + ) + end + @test ᶜYₜ.c.ρ == ᶜYₜ.c.ρq_tot + @test ᶜYₜ.c.ρ == Y.c.ρ .* p.precipitation.ᶜSqₜᵖ +end diff --git a/test/parameterized_tendencies/sponge/rayleigh_sponge.jl b/test/parameterized_tendencies/sponge/rayleigh_sponge.jl new file mode 100644 index 0000000000..ca4743ec2f --- /dev/null +++ b/test/parameterized_tendencies/sponge/rayleigh_sponge.jl @@ -0,0 +1,37 @@ + +import ClimaAtmos as CA +import SurfaceFluxes as SF +import ClimaAtmos.Parameters as CAP +import ClimaCore as CC + +include("../../test_helpers.jl") +### Common Objects ### +@testset begin + "Rayleigh-sponge functions" + ### Boilerplate default integrator objects + config = CA.AtmosConfig(Dict("initial_condition" => "DryBaroclinicWave")) + (; Y) = generate_test_simulation(config) + 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ₜ = zero(Y) + ### 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/parameterized_tendencies/sponge/viscous_sponge.jl b/test/parameterized_tendencies/sponge/viscous_sponge.jl new file mode 100644 index 0000000000..97534bb93f --- /dev/null +++ b/test/parameterized_tendencies/sponge/viscous_sponge.jl @@ -0,0 +1,30 @@ + +import ClimaAtmos as CA +import SurfaceFluxes as SF +import ClimaAtmos.Parameters as CAP +import ClimaCore as CC +include("../../test_helpers.jl") + +### Common Objects ### +@testset begin + "Rayleigh-sponge functions" + ### Boilerplate default integrator objects + config = CA.AtmosConfig(Dict("initial_condition" => "DryBaroclinicWave")) + (; Y) = generate_test_simulation(config) + z = CC.Fields.coordinate_field(Y.c).z + zmax = maximum(CC.Fields.coordinate_field(Y.f).z) + Y.c.uₕ.components.data.:1 .= ones(axes(Y.c)) + Y.c.uₕ.components.data.:2 .= ones(axes(Y.c)) + ᶜYₜ = Y .* FT(0) + FT = eltype(Y) + ### Component test begins here + rs = CA.RayleighSponge(; zd = FT(0), α_uₕ = FT(1), α_w = FT(1)) + (; ᶜβ_rayleigh_uₕ) = CA.viscous_sponge_cache(Y, rs) + @test ᶜβ_rayleigh_uₕ == @. sin(FT(π) / 2 * z / zmax)^2 + test_cache = (; viscous_sponge = (; ᶜβ_rayleigh_uₕ = ᶜβ_rayleigh_uₕ)) + CC.Fields.bycolumn(axes(Y.c)) do colidx + CA.viscous_sponge_tendency!(ᶜYₜ, Y, test_cache, FT(0), colidx, rs) + 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/rrtmgp_interface.jl b/test/rrtmgp_interface.jl new file mode 100644 index 0000000000..92681da6a6 --- /dev/null +++ b/test/rrtmgp_interface.jl @@ -0,0 +1,150 @@ +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) + +### RRTMGP Interface Tests +# Includes tests for functions defined within the RRTMGPInterface.jl file. +# Assesses that interp / extrap functions are correctly defined. + +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] + + end +end diff --git a/test/runtests.jl b/test/runtests.jl index eb5af125b6..db0c16c1c2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,11 +9,17 @@ 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 "Sponge interface tests" begin @time include("parameterized_tendencies/sponge/rayleigh_sponge.jl") end +@safetestset "Precipitation interface tests" begin @time include("parameterized_tendencies/microphysics/precipitation.jl") end +@safetestset "Model getters" begin @time include("solver/model_getters.jl") end +@safetestset "Topography tests" begin @time include("topography.jl") end #! format: on diff --git a/test/solver/model_getters.jl b/test/solver/model_getters.jl new file mode 100644 index 0000000000..b0ba8d668a --- /dev/null +++ b/test/solver/model_getters.jl @@ -0,0 +1,43 @@ +import ClimaAtmos as CA + +@testset "Model config" begin + config = CA.AtmosConfig(Dict("config" => "box")) + parsed_args = config.parsed_args + @test CA.get_model_config(parsed_args) isa CA.BoxModel + + config = CA.AtmosConfig(Dict("config" => "plane")) + parsed_args = config.parsed_args + @test CA.get_model_config(parsed_args) isa CA.PlaneModel +end + +@testset "Hyperdiffusion config" begin + @info "CAM_SE (Special case of ClimaHyperdiffusion)" + config = CA.AtmosConfig( + Dict( + "hyperdiff" => "CAM_SE", + "vorticity_hyperdiffusion_coefficient" => 0.2, + "scalar_hyperdiffusion_coefficient" => 99.0, + "divergence_damping_factor" => 99.0, + ), + ) + parsed_args = config.parsed_args + FT = eltype(config) + hyperdiff_model = CA.get_hyperdiffusion_model(parsed_args, FT) + # Coefficients are currently hardcoded in the CAM_SE type!!! + # Regardless of user specified coeffs above, the CAM_SE selection + # should override these coefficients to match values in + # https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2017MS001257 + cam_se_ν₄_vort = FT(0.150 * 1.238) + cam_se_ν₄_scalar = FT(0.751 * 1.238) + cam_se_ν₄_dd = FT(5) + @test hyperdiff_model isa CA.ClimaHyperdiffusion + @test hyperdiff_model.ν₄_vorticity_coeff == cam_se_ν₄_vort + @test hyperdiff_model.ν₄_scalar_coeff == cam_se_ν₄_scalar + @test hyperdiff_model.divergence_damping_factor == cam_se_ν₄_dd + + @info "Test unrecognized Hyperdiffusion scheme" + config = CA.AtmosConfig(Dict("hyperdiff" => "UnknownHyperdiffusion")) + parsed_args = config.parsed_args + FT = eltype(config) + @test_throws ErrorException CA.get_hyperdiffusion_model(parsed_args, FT) +end diff --git a/test/test_helpers.jl b/test/test_helpers.jl new file mode 100644 index 0000000000..bb959c9e90 --- /dev/null +++ b/test/test_helpers.jl @@ -0,0 +1,163 @@ +### BoilerPlate Code +using ClimaComms +using IntervalSets + +import ClimaCore: + ClimaCore, + Domains, + Geometry, + Grids, + Fields, + Operators, + Meshes, + Spaces, + Quadratures, + Topologies, + Hypsography + +### Unit Test Helpers +# If wrappers for general operations are required in unit tests +# (e.g. construct spaces, construct simulations from configs, +# specialised test functions with multiple uses, define them here.) + +function generate_test_simulation(config) + parsed_args = config.parsed_args + simulation = CA.get_simulation(config) + (; integrator) = simulation + Y = integrator.u + p = integrator.p + return (; Y = Y, p = p, params = p.params, simulation = simulation) +end + +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 + +""" + get_test_functions(cent_space, face_space) +Given center and face space objects for the staggered grid +construction, generate velocity profiles defined by simple +trigonometric functions (for checks against analytical solutions). +This will be generalised to accept arbitrary input profiles +(e.g. spherical harmonics) for testing purposes in the future. +""" +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 0000000000..f289a28c67 --- /dev/null +++ b/test/topography.jl @@ -0,0 +1,35 @@ +using Test +using Random +import ClimaAtmos as CA +Random.seed!(1234) + +### BoilerPlate Code +using ClimaComms +using IntervalSets + +### Unit Tests for topography +# Ensures that space construction hooks in ClimaAtmos +# result in correct warped spaces with surface elevation. + +import ClimaCore: + ClimaCore, + Domains, + Geometry, + Grids, + Fields, + Operators, + Meshes, + Spaces, + Quadratures, + Topologies, + Hypsography + +include("test_helpers.jl") + +@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 13112a15f8..b6ef39eb67 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) = @@ -133,6 +29,8 @@ end end @testset "kinetic_energy (c.f. analytical function)" begin + # Test kinetic energy function for staggered grids + # given an analytical expression for the velocity profiles (; cent_space, face_space) = get_cartesian_spaces() ccoords, fcoords = get_coords(cent_space, face_space) uₕ, uᵥ = get_test_functions(cent_space, face_space) @@ -157,6 +55,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 @@ -184,10 +83,6 @@ end UVW(Geometry.WVector(ᶠw)) CA.compute_strain_rate_center!(ᶜϵ, Geometry.Covariant123Vector.(ᶠu)) CA.compute_strain_rate_face!(ᶠϵ, Geometry.Covariant123Vector.(ᶜu)) - # Check by component, verify symmetry. - # Strain rate functions only compute vertical derivatives right now. - # Thus, terms in horizontal derivatives must be zero. - # FIXME: (This needs to be updated if a 3d operator is introduced. # Center valued strain rate @test ᶜϵ.components.data.:1 == ᶜϵ.components.data.:1 .* FT(0) @@ -269,7 +164,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