modified: test/utilities.jl

	modified:   test/runtests.jl

test/runtests.jl

@@ -9,6 +9,7 @@ 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 "Coupler Compatibility" begin @time include("coupler_compatibility.jl") end

test/utilities.jl

@@ -3,6 +3,112 @@ 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, 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
+
+
 @testset "sort_hdf5_files" begin
     day_sec(t) =
         (floor(Int, t / (60 * 60 * 24)), floor(Int, t % (60 * 60 * 24)))
@@ -16,17 +122,217 @@ import ClimaAtmos as CA
 end
 
 @testset "gaussian_smooth" begin
-
     # No smooth on constant
     @test CA.gaussian_smooth(3.0 * ones(132, 157)) ≈ 3.0 * ones(132, 157)
-
     randy = rand(123, 145)
-
     smoothed = CA.gaussian_smooth(randy)
-
     # min
     @test extrema(randy)[1] <= smoothed[1]
-
     # max
     @test extrema(randy)[2] >= smoothed[2]
 end
+
+@testset "kinetic_energy (c.f. analytical function)" begin
+    (; cent_space, face_space) = get_cartesian_spaces()
+    ccoords, fcoords = get_coords(cent_space, face_space)
+    uₕ, uᵥ = get_test_functions(cent_space, face_space)
+    (; x, y, z) = ccoords
+    # Type helpers
+    C1 = Geometry.Covariant1Vector
+    C2 = Geometry.Covariant2Vector
+    C3 = Geometry.Covariant3Vector
+    C12 = Geometry.Covariant12Vector
+    CT123 = Geometry.Contravariant123Vector
+    # Exercise function
+    κ = zeros(cent_space)
+    CA.compute_kinetic!(κ, uₕ, uᵥ)
+    ᶜκ_exact = @. 1 // 2 *
+       cos(z)^2 *
+       ((sin(x)^2) * (cos(y)^2) + (cos(x)^2) * (sin(y)^2))
+    # Test upto machine precision approximation
+    @test ᶜκ_exact ≈ κ
+end
+
+@testset "compute_strain_rate (c.f. analytical function)" begin
+    # Test compute_strain_rate_face
+    (; helem, cent_space, face_space) = get_cartesian_spaces()
+    ccoords, fcoords = get_coords(cent_space, face_space)
+    UVW = Geometry.UVWVector
+    C123 = Geometry.Covariant123Vector
+    # Alloc scratch space
+    ᶜϵ =
+        ᶜtemp_UVWxUVW = Fields.Field(
+            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
+            cent_space,
+        ) # ᶜstrain_rate
+    ᶠϵ =
+        ᶠtemp_UVWxUVW = Fields.Field(
+            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
+            face_space,
+        )
+    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))
+    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)
+    @test ᶜϵ.components.data.:5 == ᶜϵ.components.data.:7 .* FT(0)
+    @test ᶜϵ.components.data.:3 == ᶜϵ.components.data.:7
+    @test ᶜϵ.components.data.:2 == ᶜϵ.components.data.:4
+    @test ᶜϵ.components.data.:6 == ᶜϵ.components.data.:8
+
+    ᶜϵ₁₃ = ᶜϵ.components.data.:3
+    ᶜϵ₂₃ = ᶜϵ.components.data.:6
+    c₁₃ = @. -1 // 2 * sin(x) * cos(y) * sin(z)
+    c₂₃ = @. 1 // 2 * cos(x) * sin(y) * sin(z)
+    maximum(abs.((ᶜϵ₁₃ .- c₁₃) ./ (c₁₃ .+ eps(Float32)) .* 100)) < FT(0.5)
+    maximum(abs.((ᶜϵ₂₃ .- c₂₃) ./ (c₂₃ .+ eps(Float32)) .* 100)) < FT(0.5)
+
+    # Face valued strain-rate
+    (; x, y, z) = fcoords
+    ᶠϵ₁₃ = ᶠϵ.components.data.:3
+    ᶠϵ₂₃ = ᶠϵ.components.data.:6
+    f₁₃ = @. -1 // 2 * sin(x) * cos(y) * sin(z)
+    f₂₃ = @. 1 // 2 * cos(x) * sin(y) * sin(z)
+    # Check boundary conditions (see src/utils/utilities)
+    # `slab` works per element
+    for elem_id in 1:helem
+        @test maximum(
+            abs.(
+                Fields.field_values(
+                    Fields.slab(f₁₃, 1, elem_id) .-
+                    Fields.slab(ᶠϵ₁₃, 1, elem_id),
+                )
+            ),
+        ) < eps(FT) # bottom face
+        @test maximum(
+            abs.(
+                Fields.field_values(
+                    Fields.slab(f₁₃, 11, elem_id) .-
+                    Fields.slab(ᶠϵ₁₃, 11, elem_id),
+                )
+            ),
+        ) < eps(FT) # top face
+    end
+end
+
+@testset "horizontal integral at boundary" begin
+    # Test both `horizontal_integral_at_boundary` methods
+    (; cent_space, face_space) = get_cartesian_spaces()
+    _, fcoords = get_coords(cent_space, face_space)
+    ᶠu, ᶠv, ᶠw = taylor_green_ic(fcoords)
+    FT = eltype(ᶠu)
+    halflevel = ClimaCore.Utilities.half
+    y₁ = CA.horizontal_integral_at_boundary(ᶠu, halflevel)
+    y₂ = CA.horizontal_integral_at_boundary(ᶠv, halflevel)
+    y₃ = CA.horizontal_integral_at_boundary(ᶠw, halflevel)
+    @test y₁ <= sqrt(eps(FT))
+    @test y₂ <= sqrt(eps(FT))
+    @test y₃ == y₃
+    ᶠuₛ = Fields.level(ᶠu, halflevel)
+    ᶠvₛ = Fields.level(ᶠv, halflevel)
+    ᶠwₛ = Fields.level(ᶠw, halflevel)
+    y₁ = CA.horizontal_integral_at_boundary(ᶠuₛ)
+    y₂ = CA.horizontal_integral_at_boundary(ᶠvₛ)
+    y₃ = CA.horizontal_integral_at_boundary(ᶠwₛ)
+    @test y₁ <= sqrt(eps(FT))
+    @test y₂ <= sqrt(eps(FT))
+    @test y₃ == y₃
+end
+
+@testset "get mesh metrics" begin
+    # We have already constructed cent_space and face_space (3d)
+    # > These contain local geometry properties. 
+    # If grid properties change the updates will need to be caught in this 
+    # g³³_field function
+
+    # This just tests getter functions
+    # Correctness is checked in ClimaCore.jl
+    (; cent_space, face_space) = get_cartesian_spaces()
+    lg_gⁱʲ = cent_space.grid.center_local_geometry.gⁱʲ
+    lg_g³³ = lg_gⁱʲ.components.data.:9
+    @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)
+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