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
CliMA · Mar 15, 2024 · aba9022 · aba9022
1 parent b14862c
commit aba9022
Show file tree

Hide file tree

Showing 16 changed files with 498 additions and 114 deletions.
diff --git 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
@@ -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 → ...d_tendencies/gravity_wave/gw_plotutils.jl b/test/gravity_wave/gw_plotutils.jl → ...d_tendencies/gravity_wave/gw_plotutils.jl
diff --git a/...n_orographic_gravity_wave/nogw_test_3d.jl → ...n_orographic_gravity_wave/nogw_test_3d.jl b/...n_orographic_gravity_wave/nogw_test_3d.jl → ...n_orographic_gravity_wave/nogw_test_3d.jl
diff --git a/...orographic_gravity_wave/nogw_test_mima.jl → ...orographic_gravity_wave/nogw_test_mima.jl b/...orographic_gravity_wave/nogw_test_mima.jl → ...orographic_gravity_wave/nogw_test_mima.jl
diff --git a/...c_gravity_wave/nogw_test_single_column.jl → ...c_gravity_wave/nogw_test_single_column.jl b/...c_gravity_wave/nogw_test_single_column.jl → ...c_gravity_wave/nogw_test_single_column.jl
diff --git a/...y_wave/orographic_gravity_wave/ogwd_3d.jl → ...y_wave/orographic_gravity_wave/ogwd_3d.jl b/...y_wave/orographic_gravity_wave/ogwd_3d.jl → ...y_wave/orographic_gravity_wave/ogwd_3d.jl
diff --git a/.../orographic_gravity_wave/ogwd_baseflux.jl → .../orographic_gravity_wave/ogwd_baseflux.jl b/.../orographic_gravity_wave/ogwd_baseflux.jl → .../orographic_gravity_wave/ogwd_baseflux.jl
diff --git a/test/parameterized_tendencies/sponge/rayleigh_sponge.jl 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 → test/parameters/parameter_tests.jl b/test/parameter_tests.jl → test/parameters/parameter_tests.jl
diff --git a/test/parameter_tests.toml → test/parameters/parameter_tests.toml b/test/parameter_tests.toml → test/parameters/parameter_tests.toml
diff --git a/test/rrtmgp_interface.jl 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
@@ -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
@@ -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