From 68f07f4f1165c79d6650e87140a45799e760a424 Mon Sep 17 00:00:00 2001 From: Dennis Yatunin Date: Mon, 27 Jan 2025 15:09:26 -0800 Subject: [PATCH] Fix vertical metric terms, add SEM option for horizontal metric terms --- examples/hybrid/sphere/deformation_flow.jl | 2 +- src/Grids/finitedifference.jl | 94 ++++++-------- src/Grids/spectralelement.jl | 141 +++++++++++---------- src/Hypsography/Hypsography.jl | 9 +- test/Operators/integrals.jl | 5 +- test/Spaces/sphere.jl | 42 +++--- test/Spaces/unit_dss.jl | 7 +- test/TestUtilities/TestUtilities.jl | 9 +- 8 files changed, 156 insertions(+), 153 deletions(-) diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl index ef7c7e3d02..70ef627fa4 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -405,7 +405,7 @@ compare_tracer_props(a, b; buffer = 1) = all( @test compare_tracer_props(tracer_roughnesses(lim_third_upwind_sol), tracer_roughnesses(third_upwind_sol); buffer = 1.0) @test compare_tracer_props(tracer_roughnesses(lim_fct_sol) , tracer_roughnesses(third_upwind_sol); buffer = 0.93) @test compare_tracer_props(tracer_ranges(fct_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0) - @test compare_tracer_props(tracer_ranges(lim_third_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0) + @test compare_tracer_props(tracer_ranges(lim_third_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 1.2) @test compare_tracer_props(tracer_ranges(lim_fct_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0) @test compare_tracer_props(tracer_ranges(lim_first_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 0.6) @test compare_tracer_props(tracer_ranges(lim_centered_sol) , tracer_ranges(third_upwind_sol); buffer = 1.3) diff --git a/src/Grids/finitedifference.jl b/src/Grids/finitedifference.jl index 402400568b..0ba0fac27c 100644 --- a/src/Grids/finitedifference.jl +++ b/src/Grids/finitedifference.jl @@ -56,13 +56,21 @@ function _FiniteDifferenceGrid(topology::Topologies.IntervalTopology) CT = Topologies.coordinate_type(mesh) FT = Geometry.float_type(CT) Nv_face = length(mesh.faces) + Nv_cent = Nv_face - 1 # construct on CPU, adapt to GPU + center_coordinates = DataLayouts.VF{CT, Nv_cent}(Array{FT}, Nv_cent) face_coordinates = DataLayouts.VF{CT, Nv_face}(Array{FT}, Nv_face) + for v in 1:Nv_cent + center_coordinates[vindex(v)] = (mesh.faces[v] + mesh.faces[v + 1]) / 2 + end for v in 1:Nv_face face_coordinates[vindex(v)] = mesh.faces[v] end - center_local_geometry, face_local_geometry = - fd_geometry_data(face_coordinates, Val(Topologies.isperiodic(topology))) + center_local_geometry, face_local_geometry = fd_geometry_data( + center_coordinates, + face_coordinates, + Val(Topologies.isperiodic(topology)), + ) return FiniteDifferenceGrid( topology, @@ -74,84 +82,56 @@ end # called by the FiniteDifferenceGrid constructor, and the ExtrudedFiniteDifferenceGrid constructor with Hypsography function fd_geometry_data( + center_coordinates::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, face_coordinates::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, ::Val{periodic}, ) where {FT, periodic} CT = Geometry.ZPoint{FT} AIdx = (3,) LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}} + ∂x∂ξ_axes = (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()) (Ni, Nj, Nk, Nv, Nh) = size(face_coordinates) Nv_face = Nv - periodic Nv_cent = Nv - 1 - center_local_geometry = similar(face_coordinates, LG, Val(Nv_cent)) + center_local_geometry = similar(center_coordinates, LG, Val(Nv_cent)) face_local_geometry = similar(face_coordinates, LG, Val(Nv_face)) - c1(args...) = + cent_coord(args...) = + Geometry.component(center_coordinates[CartesianIndex(args...)], 1) + face_coord(args...) = Geometry.component(face_coordinates[CartesianIndex(args...)], 1) for h in 1:Nh, k in 1:Nk, j in 1:Nj, i in 1:Ni for v in 1:Nv_cent - # centers - coord⁻ = c1(i, j, k, v, h) - coord⁺ = c1(i, j, k, v + 1, h) - # use a "discrete Jacobian" - coord = (coord⁺ + coord⁻) / 2 - Δcoord = coord⁺ - coord⁻ - J = Δcoord - WJ = Δcoord - ∂x∂ξ = SMatrix{1, 1}(J) + J = face_coord(i, j, k, v + 1, h) - face_coord(i, j, k, v, h) + WJ = J + x = CT(cent_coord(i, j, k, v, h)) + ∂x∂ξ = Geometry.AxisTensor(∂x∂ξ_axes, SMatrix{1, 1}(J)) center_local_geometry[CartesianIndex(i, j, k, v, h)] = - Geometry.LocalGeometry( - CT(coord), - J, - WJ, - Geometry.AxisTensor( - ( - Geometry.LocalAxis{AIdx}(), - Geometry.CovariantAxis{AIdx}(), - ), - ∂x∂ξ, - ), - ) + Geometry.LocalGeometry(x, J, WJ, ∂x∂ξ) end for v in 1:Nv_face - coord = c1(i, j, k, v, h) - if v == 1 + if periodic && v == 1 + # periodic boundary face + J⁺ = face_coord(i, j, k, 2, h) - face_coord(i, j, k, 1, h) + J⁻ = face_coord(i, j, k, Nv, h) - face_coord(i, j, k, Nv - 1, h) + J = (J⁺ + J⁻) / 2 + WJ = J + elseif !periodic && v == 1 # bottom face - if periodic - Δcoord⁺ = c1(i, j, k, 2, h) - c1(i, j, k, 1, h) - Δcoord⁻ = c1(i, j, k, Nv, h) - c1(i, j, k, Nv - 1, h) - J = (Δcoord⁺ + Δcoord⁻) / 2 - WJ = J - else - coord⁺ = c1(i, j, k, 2, h) - J = coord⁺ - coord - WJ = J / 2 - end - elseif v == Nv_cent + 1 - @assert !periodic + J = face_coord(i, j, k, 2, h) - face_coord(i, j, k, 1, h) + WJ = J / 2 + elseif v == Nv # top face - coord⁻ = c1(i, j, k, Nv - 1, h) - J = coord - coord⁻ + @assert !periodic + J = face_coord(i, j, k, Nv, h) - face_coord(i, j, k, Nv - 1, h) WJ = J / 2 else - coord⁺ = c1(i, j, k, v + 1, h) - coord⁻ = c1(i, j, k, v - 1, h) - J = (coord⁺ - coord⁻) / 2 + J = cent_coord(i, j, k, v, h) - cent_coord(i, j, k, v - 1, h) WJ = J end - ∂x∂ξ = SMatrix{1, 1}(J) + x = CT(face_coord(i, j, k, v, h)) + ∂x∂ξ = Geometry.AxisTensor(∂x∂ξ_axes, SMatrix{1, 1}(J)) face_local_geometry[CartesianIndex(i, j, k, v, h)] = - Geometry.LocalGeometry( - CT(coord), - J, - WJ, - Geometry.AxisTensor( - ( - Geometry.LocalAxis{AIdx}(), - Geometry.CovariantAxis{AIdx}(), - ), - ∂x∂ξ, - ), - ) + Geometry.LocalGeometry(x, J, WJ, ∂x∂ξ) end end return (center_local_geometry, face_local_geometry) diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl index 583567b739..24e583b7fc 100644 --- a/src/Grids/spectralelement.jl +++ b/src/Grids/spectralelement.jl @@ -136,6 +136,7 @@ mutable struct SpectralElementGrid2D{ internal_surface_geometry::IS boundary_surface_geometries::BS enable_bubble::Bool + autodiff_metric::Bool end Adapt.@adapt_structure SpectralElementGrid2D @@ -145,15 +146,18 @@ local_geometry_type( ) where {T, Q, GG, LG, D, IS, BS} = eltype(LG) # calls eltype from DataLayouts """ - SpectralElementSpace2D(topology, quadrature_style; enable_bubble, horizontal_layout_type = DataLayouts.IJFH) + SpectralElementSpace2D(topology, quadrature_style; enable_bubble, autodiff_metric, horizontal_layout_type = DataLayouts.IJFH) Construct a `SpectralElementSpace2D` instance given a `topology` and `quadrature`. The flag `enable_bubble` enables the `bubble correction` for more accurate element areas. +The flag `autodiff_metric` enables the use of automatic differentiation instead of the +SEM for computing metric terms. # Input arguments: - topology: Topology2D - quadrature_style: QuadratureStyle - enable_bubble: Bool +- autodiff_metric: Bool - horizontal_layout_type: Type{<:AbstractData} The idea behind the so-called `bubble_correction` is that the numerical area @@ -182,6 +186,7 @@ function SpectralElementGrid2D( quadrature_style::Quadratures.QuadratureStyle; horizontal_layout_type = DataLayouts.IJFH, enable_bubble::Bool = false, + autodiff_metric::Bool = true, ) get!( Cache.OBJECT_CACHE, @@ -190,6 +195,7 @@ function SpectralElementGrid2D( topology, quadrature_style, enable_bubble, + autodiff_metric, horizontal_layout_type, ), ) do @@ -198,6 +204,7 @@ function SpectralElementGrid2D( quadrature_style, horizontal_layout_type; enable_bubble, + autodiff_metric, ) end end @@ -217,12 +224,14 @@ _SpectralElementGrid2D( quadrature_style::Quadratures.QuadratureStyle, horizontal_layout_type = DataLayouts.IJFH; enable_bubble::Bool, + autodiff_metric::Bool, ) = _SpectralElementGrid2D( topology, quadrature_style, Val(Topologies.nlocalelems(topology)), horizontal_layout_type; enable_bubble, + autodiff_metric, ) function _SpectralElementGrid2D( @@ -231,6 +240,7 @@ function _SpectralElementGrid2D( ::Val{Nh}, ::Type{horizontal_layout_type}; enable_bubble::Bool, + autodiff_metric::Bool, ) where {Nh, horizontal_layout_type} @assert horizontal_layout_type <: Union{DataLayouts.IJHF, DataLayouts.IJFH} surface_layout_type = if horizontal_layout_type <: DataLayouts.IJFH @@ -274,9 +284,8 @@ function _SpectralElementGrid2D( local_geometry = horizontal_layout_type{LG, Nq}(Array{FT}, Nh) - quad_points, quad_weights = - Quadratures.quadrature_points(FT, quadrature_style) - high_order_quad_points, high_order_quad_weights = + _, quad_weights = Quadratures.quadrature_points(FT, quadrature_style) + _, high_order_quad_weights = Quadratures.quadrature_points(FT, high_order_quadrature_style) for (lidx, elem) in enumerate(Topologies.localelems(topology)) elem_area = zero(FT) @@ -285,16 +294,18 @@ function _SpectralElementGrid2D( interior_elem_area = zero(FT) rel_interior_elem_area_Δ = zero(FT) local_geometry_slab = slab(local_geometry, lidx) + lg_args = + (global_geometry, topology, quadrature_style, autodiff_metric, elem) + high_order_lg_args = ( + global_geometry, + topology, + high_order_quadrature_style, + autodiff_metric, + elem, + ) # high-order quadrature loop for computing geometric element face area. for i in 1:high_order_Nq, j in 1:high_order_Nq - ξ = SVector(high_order_quad_points[i], high_order_quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = local_geometry_at_nodal_point(high_order_lg_args..., i, j) J_high_order = det(Geometry.components(∂u∂ξ)) WJ_high_order = J_high_order * @@ -304,14 +315,7 @@ function _SpectralElementGrid2D( end # low-order quadrature loop for computing numerical element face area for i in 1:Nq, j in 1:Nq - ξ = SVector(quad_points[i], quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = local_geometry_at_nodal_point(lg_args..., i, j) J = det(Geometry.components(∂u∂ξ)) WJ = J * quad_weights[i] * quad_weights[j] elem_area += WJ @@ -325,14 +329,7 @@ function _SpectralElementGrid2D( if enable_bubble if abs(elem_area - high_order_elem_area) ≤ eps(FT) for i in 1:Nq, j in 1:Nq - ξ = SVector(quad_points[i], quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = local_geometry_at_nodal_point(lg_args..., i, j) J = det(Geometry.components(∂u∂ξ)) WJ = J * quad_weights[i] * quad_weights[j] local_geometry_slab[slab_index(i, j)] = @@ -354,14 +351,8 @@ function _SpectralElementGrid2D( # Use uniform bubble correction if Nq == 2 for i in 1:Nq, j in 1:Nq - ξ = SVector(quad_points[i], quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = + local_geometry_at_nodal_point(lg_args..., i, j) J = det(Geometry.components(∂u∂ξ)) J += Δarea / Nq^2 WJ = J * quad_weights[i] * quad_weights[j] @@ -370,14 +361,8 @@ function _SpectralElementGrid2D( end else # Higher-order elements: Use HOMME bubble correction for the interior nodes for i in 2:(Nq - 1), j in 2:(Nq - 1) - ξ = SVector(quad_points[i], quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = + local_geometry_at_nodal_point(lg_args..., i, j) J = det(Geometry.components(∂u∂ξ)) WJ = J * quad_weights[i] * quad_weights[j] interior_elem_area += WJ @@ -391,14 +376,8 @@ function _SpectralElementGrid2D( rel_interior_elem_area_Δ = Δarea / interior_elem_area for i in 1:Nq, j in 1:Nq - ξ = SVector(quad_points[i], quad_points[j]) - u, ∂u∂ξ = compute_local_geometry( - global_geometry, - topology, - elem, - ξ, - Val(AIdx), - ) + u, ∂u∂ξ = + local_geometry_at_nodal_point(lg_args..., i, j) J = det(Geometry.components(∂u∂ξ)) # Modify J only for interior nodes if i != 1 && j != 1 && i != Nq && j != Nq @@ -496,44 +475,74 @@ function _SpectralElementGrid2D( internal_surface_geometry, boundary_surface_geometries, enable_bubble, + autodiff_metric, ) end -function compute_local_geometry( +function ξ_at_nodal_point(FT, quadrature_style, i, j) + quad_points = Quadratures.quadrature_points(FT, quadrature_style)[1] + return SVector(quad_points[i], quad_points[j]) +end + +function ∂f∂ξ_at_nodal_point(f, FT, quadrature_style, autodiff_metric, i, j) + if autodiff_metric + ξ = ξ_at_nodal_point(FT, quadrature_style, i, j) + return ForwardDiff.jacobian(f, ξ) + end + nodal_indices = SOneTo(Quadratures.degrees_of_freedom(quadrature_style)) + deriv_matrix = Quadratures.differentiation_matrix(FT, quadrature_style) + ∂f∂ξ¹ = sum(nodal_indices) do i′ + deriv_matrix[i, i′] * f(ξ_at_nodal_point(FT, quadrature_style, i′, j)) + end + ∂f∂ξ² = sum(nodal_indices) do j′ + deriv_matrix[j, j′] * f(ξ_at_nodal_point(FT, quadrature_style, i, j′)) + end + return hcat(∂f∂ξ¹, ∂f∂ξ²) +end + +function local_geometry_at_nodal_point( global_geometry::Geometry.SphericalGlobalGeometry, topology, + quadrature_style, + autodiff_metric, elem, - ξ, - ::Val{AIdx}, -) where {AIdx} + i, + j, +) + FT = eltype(Topologies.coordinate_type(topology)) + AIdx = Geometry.coordinate_axis(get_CoordType2D(topology)) + ξ = ξ_at_nodal_point(FT, quadrature_style, i, j) x = Meshes.coordinates(topology.mesh, elem, ξ) - u = Geometry.LatLongPoint(x, global_geometry) ∂x∂ξ = Geometry.AxisTensor( (Geometry.Cartesian123Axis(), Geometry.CovariantAxis{AIdx}()), - ForwardDiff.jacobian(ξ) do ξ + ∂f∂ξ_at_nodal_point(FT, quadrature_style, autodiff_metric, i, j) do ξ Geometry.components(Meshes.coordinates(topology.mesh, elem, ξ)) end, ) + u = Geometry.LatLongPoint(x, global_geometry) G = Geometry.local_to_cartesian(global_geometry, u) ∂u∂ξ = Geometry.project(Geometry.LocalAxis{AIdx}(), G' * ∂x∂ξ) - return u, ∂u∂ξ end -function compute_local_geometry( - global_geometry::Geometry.AbstractGlobalGeometry, +function local_geometry_at_nodal_point( + ::Geometry.AbstractGlobalGeometry, topology, + quadrature_style, + autodiff_metric, elem, - ξ, - ::Val{AIdx}, -) where {AIdx} + i, + j, +) + FT = eltype(Topologies.coordinate_type(topology)) + AIdx = Geometry.coordinate_axis(get_CoordType2D(topology)) + ξ = ξ_at_nodal_point(FT, quadrature_style, i, j) u = Meshes.coordinates(topology.mesh, elem, ξ) ∂u∂ξ = Geometry.AxisTensor( (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), - ForwardDiff.jacobian(ξ) do ξ + ∂f∂ξ_at_nodal_point(FT, quadrature_style, autodiff_metric, i, j) do ξ Geometry.components(Meshes.coordinates(topology.mesh, elem, ξ)) end, ) - return u, ∂u∂ξ end diff --git a/src/Hypsography/Hypsography.jl b/src/Hypsography/Hypsography.jl index b30197fe42..a8a6949bb8 100644 --- a/src/Hypsography/Hypsography.jl +++ b/src/Hypsography/Hypsography.jl @@ -173,11 +173,15 @@ function _ExtrudedFiniteDifferenceGrid( @assert Spaces.grid(axes(adaption.surface)) == horizontal_grid z_surface = Fields.field_values(adaption.surface) + center_z_ref = + Grids.local_geometry_data(vertical_grid, Grids.CellCenter()).coordinates face_z_ref = Grids.local_geometry_data(vertical_grid, Grids.CellFace()).coordinates vertical_domain = Topologies.domain(vertical_grid) z_top = vertical_domain.coord_max + center_z = + ref_z_to_physical_z.(Ref(adaption), center_z_ref, z_surface, Ref(z_top)) face_z = ref_z_to_physical_z.(Ref(adaption), face_z_ref, z_surface, Ref(z_top)) @@ -186,16 +190,18 @@ function _ExtrudedFiniteDifferenceGrid( vertical_grid, adaption, global_geometry, + center_z, face_z, ) end -# generic 5-arg hypsography constructor, uses computed face_z points +# generic hypsography constructor, uses computed center_z and face_z points function _ExtrudedFiniteDifferenceGrid( horizontal_grid::Grids.AbstractGrid, vertical_grid::Grids.FiniteDifferenceGrid, adaption::HypsographyAdaption, global_geometry::Geometry.AbstractGlobalGeometry, + center_z::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, face_z::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, ) where {FT} # construct the "flat" grid @@ -213,6 +219,7 @@ function _ExtrudedFiniteDifferenceGrid( ArrayType = ClimaComms.array_type(horizontal_grid.topology) # currently only works on Arrays (center_z_local_geometry, face_z_local_geometry) = Grids.fd_geometry_data( + Adapt.adapt(Array, center_z), Adapt.adapt(Array, face_z), Val(Topologies.isperiodic(vertical_grid.topology)), ) diff --git a/test/Operators/integrals.jl b/test/Operators/integrals.jl index 44f4a20620..da7f893ba7 100644 --- a/test/Operators/integrals.jl +++ b/test/Operators/integrals.jl @@ -189,8 +189,9 @@ end device = ClimaComms.device() context = ClimaComms.SingletonCommsContext(device) for FT in (Float32, Float64) - for topography in (false, true), deep in (false, true) - space_kwargs = (; context, topography, deep) + bools = (false, true) + for topography in bools, deep in bools, autodiff_metric in bools + space_kwargs = (; context, topography, deep, autodiff_metric) space = TU.CenterExtrudedFiniteDifferenceSpace(FT; space_kwargs...) test_fubinis_theorem(space) end diff --git a/test/Spaces/sphere.jl b/test/Spaces/sphere.jl index 35e8e42599..6a95e853d5 100644 --- a/test/Spaces/sphere.jl +++ b/test/Spaces/sphere.jl @@ -113,25 +113,25 @@ end horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem) horztopology = Topologies.Topology2D(context, horzmesh) quad = Quadratures.GLL{Nq}() - horzgrid = Grids.SpectralElementGrid2D(horztopology, quad) - horzspace = Spaces.SpectralElementSpace2D(horzgrid) - @test Spaces.node_horizontal_length_scale(horzspace) ≈ - sqrt((4 * pi * radius^2) / (helem^2 * 6 * (Nq - 1)^2)) - - # "shallow atmosphere" spherical shell: volume = surface area * height - shallow_grid = Grids.ExtrudedFiniteDifferenceGrid(horzgrid, vertgrid) - - @test sum(ones(Spaces.CenterExtrudedFiniteDifferenceSpace(shallow_grid))) ≈ - 4pi * radius^2 * (zlim[2] - zlim[1]) rtol = 1e-3 - @test sum(ones(Spaces.FaceExtrudedFiniteDifferenceSpace(shallow_grid))) ≈ - 4pi * radius^2 * (zlim[2] - zlim[1]) rtol = 1e-3 - - deep_grid = - Grids.ExtrudedFiniteDifferenceGrid(horzgrid, vertgrid; deep = true) - - @test sum(ones(Spaces.CenterExtrudedFiniteDifferenceSpace(deep_grid))) ≈ - 4pi / 3 * ((zlim[2] + radius)^3 - (zlim[1] + radius)^3) rtol = 1e-3 - @test sum(ones(Spaces.FaceExtrudedFiniteDifferenceSpace(deep_grid))) ≈ - 4pi / 3 * ((zlim[2] + radius)^3 - (zlim[1] + radius)^3) rtol = 1e-3 - + for autodiff_metric in (true, false) + horzgrid = + Grids.SpectralElementGrid2D(horztopology, quad; autodiff_metric) + horzspace = Spaces.SpectralElementSpace2D(horzgrid) + @test Spaces.node_horizontal_length_scale(horzspace) ≈ + sqrt((4 * pi * radius^2) / (helem^2 * 6 * (Nq - 1)^2)) + + # In a shallow atmosphere, the volume is given by surface area * height. + shallow_volume = 4pi * radius^2 * (zlim[2] - zlim[1]) + deep_volume = 4pi / 3 * ((zlim[2] + radius)^3 - (zlim[1] + radius)^3) + + for (deep, ref_volume) in ((false, shallow_volume), (true, deep_volume)) + grid = Grids.ExtrudedFiniteDifferenceGrid(horzgrid, vertgrid; deep) + center_volume = + sum(ones(Spaces.CenterExtrudedFiniteDifferenceSpace(grid))) + face_volume = + sum(ones(Spaces.FaceExtrudedFiniteDifferenceSpace(grid))) + @test center_volume ≈ ref_volume rtol = deep ? 2e-5 : 6e-6 + @test face_volume ≈ ref_volume rtol = deep ? 1e-5 : 6e-6 + end + end end diff --git a/test/Spaces/unit_dss.jl b/test/Spaces/unit_dss.jl index 1dfd8b741a..bdf6f8e4e4 100644 --- a/test/Spaces/unit_dss.jl +++ b/test/Spaces/unit_dss.jl @@ -126,15 +126,16 @@ function test_dss_conservation(space) integral_before_dss = sum(field) Spaces.weighted_dss!(field) integral_after_dss = sum(field) - @test integral_after_dss ≈ integral_before_dss rtol = 14 * eps(FT) + @test integral_after_dss ≈ integral_before_dss rtol = 18 * eps(FT) end @testset "DSS Conservation on Cubed Sphere" begin device = ClimaComms.device() context = ClimaComms.SingletonCommsContext(device) for FT in (Float32, Float64) - for topography in (false, true), deep in (false, true) - space_kwargs = (; context, topography, deep) + bools = (false, true) + for topography in bools, deep in bools, autodiff_metric in bools + space_kwargs = (; context, topography, deep, autodiff_metric) center_space = TU.CenterExtrudedFiniteDifferenceSpace(FT; space_kwargs...) for space in (center_space, Spaces.face_space(center_space)) diff --git a/test/TestUtilities/TestUtilities.jl b/test/TestUtilities/TestUtilities.jl index b415db01ca..220f829b17 100644 --- a/test/TestUtilities/TestUtilities.jl +++ b/test/TestUtilities/TestUtilities.jl @@ -113,6 +113,7 @@ function CenterExtrudedFiniteDifferenceSpace( Nq = 4, deep = false, topography = false, + autodiff_metric = true, horizontal_layout_type = DataLayouts.IJFH, ) where {FT} radius = FT(128) @@ -131,8 +132,12 @@ function CenterExtrudedFiniteDifferenceSpace( hmesh = Meshes.EquiangularCubedSphere(hdomain, helem) htopology = Topologies.Topology2D(context, hmesh) quad = Quadratures.GLL{Nq}() - hspace = - Spaces.SpectralElementSpace2D(htopology, quad; horizontal_layout_type) + hspace = Spaces.SpectralElementSpace2D( + htopology, + quad; + autodiff_metric, + horizontal_layout_type, + ) hypsography = if topography # some non-trivial function of latitude and longitude