From 590b1c9995e9b91074e63001c05c4ed3f4fc1f7c Mon Sep 17 00:00:00 2001 From: Dennis Yatunin Date: Thu, 16 Jan 2025 17:03:25 -0800 Subject: [PATCH] Clarify DSS weights and fix DSS buffer bug --- ext/cuda/topologies_dss.jl | 12 +++++----- src/Fields/Fields.jl | 7 ++---- src/Grids/Grids.jl | 2 +- src/Grids/extruded.jl | 7 ++++-- src/Grids/level.jl | 5 ++++- src/Grids/spectralelement.jl | 40 ++++++++++++++++++--------------- src/Quadratures/Quadratures.jl | 35 +++++++++++------------------ src/Spaces/Spaces.jl | 3 ++- src/Spaces/dss.jl | 24 +++++++++----------- src/Spaces/extruded.jl | 4 ---- src/Spaces/spectralelement.jl | 2 -- src/Topologies/dss.jl | 40 ++++++++++++++++----------------- src/Topologies/dss_transform.jl | 4 ++-- test/Spaces/unit_dss.jl | 32 ++++++++++++++++++++++++++ test/Spaces/unit_spaces.jl | 2 +- 15 files changed, 120 insertions(+), 99 deletions(-) diff --git a/ext/cuda/topologies_dss.jl b/ext/cuda/topologies_dss.jl index 5662e68d5e..e3fa7df9d7 100644 --- a/ext/cuda/topologies_dss.jl +++ b/ext/cuda/topologies_dss.jl @@ -2,7 +2,7 @@ import ClimaCore: DataLayouts, Topologies, Spaces, Fields import ClimaCore.DataLayouts: CartesianFieldIndex using CUDA import ClimaCore.Topologies -import ClimaCore.Topologies: DSSDataTypes, DSSPerimeterTypes, DSSWeightTypes +import ClimaCore.Topologies: DSSDataTypes, DSSPerimeterTypes import ClimaCore.Topologies: perimeter_vertex_node_index _max_threads_cuda() = 256 @@ -198,7 +198,7 @@ function Topologies.dss_transform!( data::DSSDataTypes, perimeter::Topologies.Perimeter2D, local_geometry::DSSDataTypes, - weight::DSSWeightTypes, + dss_weights::DSSDataTypes, localelems::AbstractVector{Int}, ) nlocalelems = length(localelems) @@ -214,7 +214,7 @@ function Topologies.dss_transform!( data, perimeter, local_geometry, - weight, + dss_weights, localelems, Val(nlocalelems), ) @@ -231,7 +231,7 @@ function Topologies.dss_transform!( data, perimeter, local_geometry, - weight, + dss_weights, localelems, ) end @@ -243,7 +243,7 @@ function dss_transform_kernel!( data::DSSDataTypes, perimeter::Topologies.Perimeter2D, local_geometry::DSSDataTypes, - weight::DSSWeightTypes, + dss_weights::DSSDataTypes, localelems::AbstractVector{Int}, ::Val{nlocalelems}, ) where {nlocalelems} @@ -260,7 +260,7 @@ function dss_transform_kernel!( src = Topologies.dss_transform( data[loc], local_geometry[loc], - weight[loc], + dss_weights[loc], ) perimeter_data[CI(p, 1, 1, level, elem)] = Topologies.drop_vert_dim(eltype(perimeter_data), src) diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 9e687f1bd8..5f3f62bddf 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -474,11 +474,8 @@ end Create a buffer for communicating neighbour information of `field`. """ -function Spaces.create_dss_buffer(field::Field) - space = axes(field) - hspace = Spaces.horizontal_space(space) - Spaces.create_dss_buffer(field_values(field), hspace) -end +Spaces.create_dss_buffer(field::Field) = + Spaces.create_dss_buffer(field_values(field), axes(field)) Base.@propagate_inbounds function level( field::Union{ diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl index 39f0f235d2..400aa269fb 100644 --- a/src/Grids/Grids.jl +++ b/src/Grids/Grids.jl @@ -60,7 +60,7 @@ function local_geometry_type end # Fallback, but this requires user error-handling local_geometry_type(::Type{T}) where {T} = Union{} -function local_dss_weights end +function dss_weights end function quadrature_style end function vertical_topology end diff --git a/src/Grids/extruded.jl b/src/Grids/extruded.jl index 105cd887ca..6c583cad47 100644 --- a/src/Grids/extruded.jl +++ b/src/Grids/extruded.jl @@ -129,8 +129,11 @@ topology(grid::ExtrudedFiniteDifferenceGrid) = topology(grid.horizontal_grid) vertical_topology(grid::ExtrudedFiniteDifferenceGrid) = topology(grid.vertical_grid) -local_dss_weights(grid::ExtrudedFiniteDifferenceGrid) = - local_dss_weights(grid.horizontal_grid) +# Since ∂z/∂ξ₃ and r are continuous across element boundaries, we can reuse +# the horizontal weights instead of calling compute_dss_weights on the +# extruded local geometry. +dss_weights(grid::AbstractExtrudedFiniteDifferenceGrid, ::Staggering) = + dss_weights(grid.horizontal_grid, nothing) local_geometry_data(grid::AbstractExtrudedFiniteDifferenceGrid, ::CellCenter) = grid.center_local_geometry diff --git a/src/Grids/level.jl b/src/Grids/level.jl index 7211124a4a..bb18b0d0ea 100644 --- a/src/Grids/level.jl +++ b/src/Grids/level.jl @@ -16,7 +16,10 @@ level( topology(levelgrid::LevelGrid) = topology(levelgrid.full_grid) -local_dss_weights(grid::LevelGrid) = local_dss_weights(grid.full_grid) +dss_weights(levelgrid::LevelGrid{<:Any, Int}, ::Nothing) = + level(dss_weights(levelgrid.full_grid, CellCenter()), levelgrid.level) +dss_weights(levelgrid::LevelGrid{<:Any, PlusHalf{Int}}, ::Nothing) = + level(dss_weights(levelgrid.full_grid, CellFace()), levelgrid.level + half) local_geometry_type(::Type{LevelGrid{G, L}}) where {G, L} = local_geometry_type(G) diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl index 7e20872fd5..583567b739 100644 --- a/src/Grids/spectralelement.jl +++ b/src/Grids/spectralelement.jl @@ -102,17 +102,13 @@ function _SpectralElementGrid1D( ) end end - dss_weights = copy(local_geometry.J) - dss_weights .= one(FT) - Topologies.dss_1d!(topology, dss_weights) - dss_weights = one(FT) ./ dss_weights return SpectralElementGrid1D( topology, quadrature_style, global_geometry, local_geometry, - dss_weights, + compute_dss_weights(local_geometry, topology, quadrature_style), ) end @@ -136,7 +132,7 @@ mutable struct SpectralElementGrid2D{ quadrature_style::Q global_geometry::GG local_geometry::LG - local_dss_weights::D + dss_weights::D internal_surface_geometry::IS boundary_surface_geometries::BS enable_bubble::Bool @@ -418,14 +414,6 @@ function _SpectralElementGrid2D( end end - # dss_weights = J ./ dss(J) - J = DataLayouts.rebuild(local_geometry.J, DA) - dss_local_weights = copy(J) - if quadrature_style isa Quadratures.GLL - Topologies.dss!(dss_local_weights, topology) - end - dss_local_weights .= J ./ dss_local_weights - SG = Geometry.SurfaceGeometry{ FT, Geometry.AxisVector{FT, Geometry.LocalAxis{AIdx}, SVector{2, FT}}, @@ -497,12 +485,14 @@ function _SpectralElementGrid2D( internal_surface_geometry = nothing boundary_surface_geometries = nothing end + + device_local_geometry = DataLayouts.rebuild(local_geometry, DA) return SpectralElementGrid2D( topology, quadrature_style, global_geometry, - DataLayouts.rebuild(local_geometry, DA), - dss_local_weights, + device_local_geometry, + compute_dss_weights(device_local_geometry, topology, quadrature_style), internal_surface_geometry, boundary_surface_geometries, enable_bubble, @@ -578,6 +568,21 @@ function compute_surface_geometry( return Geometry.SurfaceGeometry(sWJ, n) end +function compute_dss_weights(local_geometry, topology, quadrature_style) + is_dss_required = + Quadratures.unique_degrees_of_freedom(quadrature_style) < + Quadratures.degrees_of_freedom(quadrature_style) + # Although the weights are defined as WJ / Σ collocated WJ, we can use J + # instead of WJ if the weights are symmetric across element boundaries. + dss_weights = copy(local_geometry.J) + if topology isa Topologies.IntervalTopology + is_dss_required && Topologies.dss_1d!(topology, dss_weights) + else + is_dss_required && Topologies.dss!(dss_weights, topology) + end + @. dss_weights = local_geometry.J / dss_weights + return dss_weights +end # accessors @@ -588,8 +593,7 @@ local_geometry_data(grid::AbstractSpectralElementGrid, ::Nothing) = global_geometry(grid::AbstractSpectralElementGrid) = grid.global_geometry quadrature_style(grid::AbstractSpectralElementGrid) = grid.quadrature_style -local_dss_weights(grid::SpectralElementGrid1D) = grid.dss_weights -local_dss_weights(grid::SpectralElementGrid2D) = grid.local_dss_weights +dss_weights(grid::AbstractSpectralElementGrid, ::Nothing) = grid.dss_weights ## GPU compatibility struct DeviceSpectralElementGrid2D{Q, GG, LG} <: AbstractSpectralElementGrid diff --git a/src/Quadratures/Quadratures.jl b/src/Quadratures/Quadratures.jl index 3584d8a745..9cb2e053ea 100644 --- a/src/Quadratures/Quadratures.jl +++ b/src/Quadratures/Quadratures.jl @@ -9,21 +9,21 @@ export QuadratureStyle, GLL, GL, polynomial_degree, degrees_of_freedom, quadrature_points """ - QuadratureStyle + QuadratureStyle{Nq} Quadrature style supertype. See sub-types: - [`GLL`](@ref) - [`GL`](@ref) - [`Uniform`](@ref) """ -abstract type QuadratureStyle end +abstract type QuadratureStyle{Nq} end """ polynomial_degree(QuadratureStyle) -> Int Returns the polynomial degree of the `QuadratureStyle` concrete type """ -function polynomial_degree end +@inline polynomial_degree(::QuadratureStyle{Nq}) where {Nq} = Nq - 1 """ @@ -31,7 +31,7 @@ function polynomial_degree end Returns the degrees_of_freedom of the `QuadratureStyle` concrete type """ -function degrees_of_freedom end +@inline degrees_of_freedom(::QuadratureStyle{Nq}) where {Nq} = Nq """ points, weights = quadrature_points(::Type{FT}, quadrature_style) @@ -46,15 +46,12 @@ function quadrature_points end Gauss-Legendre-Lobatto quadrature using `Nq` quadrature points. """ -struct GLL{Nq} <: QuadratureStyle end +struct GLL{Nq} <: QuadratureStyle{Nq} end Base.show(io::IO, ::GLL{Nq}) where {Nq} = print(io, Nq, "-point Gauss-Legendre-Lobatto quadrature") -@inline polynomial_degree(::GLL{Nq}) where {Nq} = Int(Nq - 1) -@inline degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq) -unique_degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq - 1) - +unique_degrees_of_freedom(::GLL{Nq}) where {Nq} = Nq - 1 @generated function quadrature_points(::Type{FT}, ::GLL{Nq}) where {FT, Nq} points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.both) :($(SVector{Nq}(points)), $(SVector{Nq}(weights))) @@ -65,14 +62,12 @@ end Gauss-Legendre quadrature using `Nq` quadrature points. """ -struct GL{Nq} <: QuadratureStyle end +struct GL{Nq} <: QuadratureStyle{Nq} end + Base.show(io::IO, ::GL{Nq}) where {Nq} = print(io, Nq, "-point Gauss-Legendre quadrature") -@inline polynomial_degree(::GL{Nq}) where {Nq} = Int(Nq - 1) -@inline degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq) -unique_degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq) - +unique_degrees_of_freedom(::GL{Nq}) where {Nq} = Nq @generated function quadrature_points(::Type{FT}, ::GL{Nq}) where {FT, Nq} points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.neither) :($(SVector{Nq}(points)), $(SVector{Nq}(weights))) @@ -83,11 +78,9 @@ end Uniformly-spaced quadrature. """ -struct Uniform{Nq} <: QuadratureStyle end - -@inline polynomial_degree(::Uniform{Nq}) where {Nq} = Int(Nq - 1) -@inline degrees_of_freedom(::Uniform{Nq}) where {Nq} = Int(Nq) +struct Uniform{Nq} <: QuadratureStyle{Nq} end +unique_degrees_of_freedom(::Uniform{Nq}) where {Nq} = Nq @generated function quadrature_points(::Type{FT}, ::Uniform{Nq}) where {FT, Nq} points = SVector{Nq}(range(-1 + FT(1 / Nq), step = FT(2 / Nq), length = Nq)) weights = SVector{Nq}(ntuple(i -> FT(2 / Nq), Nq)) @@ -99,11 +92,9 @@ end Uniformly-spaced quadrature including boundary. """ -struct ClosedUniform{Nq} <: QuadratureStyle end - -@inline polynomial_degree(::ClosedUniform{Nq}) where {Nq} = Int(Nq - 1) -@inline degrees_of_freedom(::ClosedUniform{Nq}) where {Nq} = Int(Nq) +struct ClosedUniform{Nq} <: QuadratureStyle{Nq} end +unique_degrees_of_freedom(::ClosedUniform{Nq}) where {Nq} = Nq - 1 @generated function quadrature_points( ::Type{FT}, ::ClosedUniform{Nq}, diff --git a/src/Spaces/Spaces.jl b/src/Spaces/Spaces.jl index 9c955b5266..edc59ba2b8 100644 --- a/src/Spaces/Spaces.jl +++ b/src/Spaces/Spaces.jl @@ -37,7 +37,7 @@ import ..Grids: local_geometry_type, local_geometry_data, global_geometry, - local_dss_weights, + dss_weights, quadrature_style import ClimaComms @@ -69,6 +69,7 @@ vertical_topology(space::AbstractSpace) = vertical_topology(grid(space)) local_geometry_data(space::AbstractSpace) = local_geometry_data(grid(space), staggering(space)) +dss_weights(space::AbstractSpace) = dss_weights(grid(space), staggering(space)) function n_elements_per_panel_direction(space::AbstractSpace) hspace = Spaces.horizontal_space(space) diff --git a/src/Spaces/dss.jl b/src/Spaces/dss.jl index 53b99cf9dd..96caa127f5 100644 --- a/src/Spaces/dss.jl +++ b/src/Spaces/dss.jl @@ -13,8 +13,7 @@ import ..Topologies: load_from_recv_buffer!, DSSTypesAll, DSSDataTypes, - DSSPerimeterTypes, - DSSWeightTypes + DSSPerimeterTypes perimeter(space::AbstractSpectralElementSpace) = Topologies.Perimeter2D( @@ -23,10 +22,7 @@ perimeter(space::AbstractSpectralElementSpace) = Topologies.Perimeter2D( """ - create_dss_buffer( - data::Union{DataLayouts.IJFH, DataLayouts.VIJFH}, - hspace::AbstractSpectralElementSpace, - ) + create_dss_buffer(data, space) Creates a [`DSSBuffer`](@ref) for the field data corresponding to `data` """ @@ -37,13 +33,13 @@ function create_dss_buffer( DataLayouts.VIJFH, DataLayouts.VIJHF, }, - hspace::SpectralElementSpace2D, + space, ) create_dss_buffer( data, - topology(hspace), - local_geometry_data(hspace), - local_dss_weights(hspace), + topology(space), + local_geometry_data(space), + dss_weights(space), ) end @@ -54,7 +50,7 @@ function create_dss_buffer( DataLayouts.VIFH, DataLayouts.VIHF, }, - hspace::SpectralElementSpace1D, + space, ) nothing end @@ -122,7 +118,7 @@ function weighted_dss_prepare!( dss_buffer, data, local_geometry_data(space), - local_dss_weights(hspace), + dss_weights(space), Spaces.perimeter(hspace), dss_buffer.perimeter_elems, ) @@ -236,7 +232,7 @@ function weighted_dss_internal!( topology(hspace), data, local_geometry_data(space), - local_dss_weights(space), + dss_weights(space), ) else device = ClimaComms.device(topology(hspace)) @@ -245,7 +241,7 @@ function weighted_dss_internal!( dss_buffer, data, local_geometry_data(space), - local_dss_weights(space), + dss_weights(space), Spaces.perimeter(hspace), dss_buffer.internal_elems, ) diff --git a/src/Spaces/extruded.jl b/src/Spaces/extruded.jl index 427aa6f76a..ced77f28d7 100644 --- a/src/Spaces/extruded.jl +++ b/src/Spaces/extruded.jl @@ -90,10 +90,6 @@ FaceExtrudedFiniteDifferenceSpace(space::ExtrudedFiniteDifferenceSpace) = CenterExtrudedFiniteDifferenceSpace(space::ExtrudedFiniteDifferenceSpace) = ExtrudedFiniteDifferenceSpace(grid(space), CellCenter()) - -local_dss_weights(space::ExtrudedFiniteDifferenceSpace) = - local_dss_weights(grid(space)) - staggering(space::ExtrudedFiniteDifferenceSpace) = getfield(space, :staggering) grid(space::ExtrudedFiniteDifferenceSpace) = getfield(space, :grid) space(space::ExtrudedFiniteDifferenceSpace, staggering::Staggering) = diff --git a/src/Spaces/spectralelement.jl b/src/Spaces/spectralelement.jl index a1b02e424d..a38b211892 100644 --- a/src/Spaces/spectralelement.jl +++ b/src/Spaces/spectralelement.jl @@ -7,8 +7,6 @@ Topologies.nlocalelems(space::AbstractSpectralElementSpace) = quadrature_style(space::AbstractSpectralElementSpace) = quadrature_style(grid(space)) -local_dss_weights(space::AbstractSpectralElementSpace) = - local_dss_weights(grid(space)) horizontal_space(space::AbstractSpectralElementSpace) = space nlevels(space::AbstractSpectralElementSpace) = 1 diff --git a/src/Topologies/dss.jl b/src/Topologies/dss.jl index 75afab7fec..0b2c389d38 100644 --- a/src/Topologies/dss.jl +++ b/src/Topologies/dss.jl @@ -17,7 +17,6 @@ const DSSDataTypes = Union{ DataLayouts.VIJHF, } const DSSPerimeterTypes = Union{DataLayouts.VIFH, DataLayouts.VIHF} -const DSSWeightTypes = Union{DataLayouts.IJFH, DataLayouts.IJHF} """ DSSBuffer{G, D, A, B} @@ -63,7 +62,7 @@ end data::Union{DataLayouts.IJFH{S}, DataLayouts.IJHF{S}, DataLayouts.VIJFH{S}, DataLayouts.VIJHF{S}}, topology::Topology2D, local_geometry::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF, Nothing} = nothing, - local_weights::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF, Nothing} = nothing, + dss_weights::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF, Nothing} = nothing, ) where {S} Creates a [`DSSBuffer`](@ref) for the field data corresponding to `data` @@ -72,13 +71,13 @@ create_dss_buffer( data::DSSDataTypes, topology::Topology2D, local_geometry::Union{DSSDataTypes, Nothing} = nothing, - local_weights::Union{DSSDataTypes, Nothing} = nothing, + dss_weights::Union{DSSDataTypes, Nothing} = nothing, ) = create_dss_buffer( data, topology, DataLayouts.VIFH, local_geometry, - local_weights, + dss_weights, ) function create_dss_buffer( @@ -86,14 +85,14 @@ function create_dss_buffer( topology::Topology2D, ::Type{PerimeterLayout}, local_geometry::Union{DSSDataTypes, Nothing} = nothing, - local_weights::Union{DSSDataTypes, Nothing} = nothing, + dss_weights::Union{DSSDataTypes, Nothing} = nothing, ) where {PerimeterLayout} S = eltype(data) Nij = DataLayouts.get_Nij(data) Nij_lg = isnothing(local_geometry) ? Nij : DataLayouts.get_Nij(local_geometry) Nij_weights = - isnothing(local_weights) ? Nij : DataLayouts.get_Nij(local_weights) + isnothing(dss_weights) ? Nij : DataLayouts.get_Nij(dss_weights) @assert Nij == Nij_lg == Nij_weights perimeter::Perimeter2D = Perimeter2D(Nij) context = ClimaComms.context(topology) @@ -101,7 +100,6 @@ function create_dss_buffer( convert_to_array = DA isa Array ? false : true (_, _, _, Nv, Nh) = Base.size(data) Np = length(perimeter) - Nf = DataLayouts.ncomponents(data) nfacedof = Nij - 2 T = eltype(parent(data)) # Add TS for Covariant123Vector @@ -112,8 +110,9 @@ function create_dss_buffer( elseif eltype(data) <: Geometry.Contravariant123Vector Geometry.UVWVector{T} else - _transformed_type(data, local_geometry, local_weights, DA) # extract transformed type + _transformed_type(data, local_geometry, dss_weights, DA) # extract transformed type end + Nf = DataLayouts.typesize(T, TS) perimeter_data = if !isnothing(local_geometry) fdim = DataLayouts.field_dim(DataLayouts.singleton(local_geometry)) @@ -188,7 +187,7 @@ assert_same_eltype(::DataLayouts.AbstractData, ::Nothing) = nothing dss_buffer::DSSBuffer, data::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, local_geometry::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, - weight::Union{DataLayouts.IJFH, DataLayouts.IJHF}, + dss_weights::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, perimeter::Perimeter2D, localelems::AbstractVector{Int}, ) @@ -202,7 +201,7 @@ Arguments: - `dss_buffer`: [`DSSBuffer`](@ref) generated by `create_dss_buffer` function for field data - `data`: field data - `local_geometry`: local metric information defined at each node -- `weight`: local dss weights for horizontal space +- `dss_weights`: local dss weights for horizontal space - `perimeter`: perimeter iterator - `localelems`: list of local elements to perform transformation operations on @@ -213,7 +212,7 @@ function dss_transform!( dss_buffer::DSSBuffer, data::DSSDataTypes, local_geometry::DSSDataTypes, - weight::DSSWeightTypes, + dss_weights::DSSDataTypes, perimeter::Perimeter2D, localelems::AbstractVector{Int}, ) @@ -224,7 +223,7 @@ function dss_transform!( data, perimeter, local_geometry, - weight, + dss_weights, localelems, ) end @@ -257,8 +256,9 @@ to `UVVector` if `T <: UVVector`. ::ClimaComms.AbstractCPUDevice, perimeter_data::Union{DataLayouts.VIFH, DataLayouts.VIHF}, data::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, - perimeter::AbstractPerimeter, - bc, + perimeter::Perimeter2D, + local_geometry::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, + dss_weights::Union{DataLayouts.IJFH, DataLayouts.IJHF, DataLayouts.VIJFH, DataLayouts.VIJHF}, localelems::Vector{Int}, ) @@ -270,7 +270,8 @@ Arguments: - `perimeter_data`: contains the perimeter field data, represented on the physical axis, corresponding to the full field data in `data` - `data`: field data - `perimeter`: perimeter iterator -- `bc`: A broadcasted object representing the dss transform. +- `local_geometry`: local metric information defined at each node +- `dss_weights`: local dss weights for horizontal space - `localelems`: list of local elements to perform transformation operations on Part of [`ClimaCore.Spaces.weighted_dss!`](@ref). @@ -280,21 +281,20 @@ function dss_transform!( perimeter_data::DSSPerimeterTypes, data::DSSDataTypes, perimeter::Perimeter2D{Nq}, - local_geometry, - weight::DSSWeightTypes, + local_geometry::DSSDataTypes, + dss_weights::DSSDataTypes, localelems::Vector{Int}, ) where {Nq} (_, _, _, nlevels, _) = DataLayouts.universal_size(perimeter_data) CI = CartesianIndex - (; ∂ξ∂x) = local_geometry @inbounds for elem in localelems for (p, (ip, jp)) in enumerate(perimeter) for level in 1:nlevels loc = CI(ip, jp, 1, level, elem) src = dss_transform( data[loc], - local_geometry[CI(ip, jp, 1, level, elem)], - weight[loc], + local_geometry[loc], + dss_weights[loc], ) perimeter_data[CI(p, 1, 1, level, elem)] = drop_vert_dim(eltype(perimeter_data), src) diff --git a/src/Topologies/dss_transform.jl b/src/Topologies/dss_transform.jl index 1f887a6a38..88858b4314 100644 --- a/src/Topologies/dss_transform.jl +++ b/src/Topologies/dss_transform.jl @@ -226,13 +226,13 @@ end _transformed_type( data::DataLayouts.AbstractData, local_geometry::Union{DataLayouts.AbstractData, Nothing}, - local_weights::Union{DataLayouts.AbstractData, Nothing}, + dss_weights::Union{DataLayouts.AbstractData, Nothing}, ::Type{DA}, ) where {DA} = typeof( dss_transform( _representative_slab(data, DA), _representative_slab(local_geometry, DA), - _representative_slab(local_weights, DA), + _representative_slab(dss_weights, DA), CartesianIndex(1, 1, 1, 1, 1), ), ) diff --git a/test/Spaces/unit_dss.jl b/test/Spaces/unit_dss.jl index 32b0b259fb..1dfd8b741a 100644 --- a/test/Spaces/unit_dss.jl +++ b/test/Spaces/unit_dss.jl @@ -8,6 +8,7 @@ using ClimaComms using Random ClimaComms.@import_required_backends +import ClimaCore import ClimaCore: Domains, Fields, @@ -19,6 +20,11 @@ import ClimaCore: Topologies, DataLayouts +include( + joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"), +) +import .TestUtilities as TU + function get_space_cs(::Type{FT}; context, R = 300.0) where {FT} domain = Domains.SphereDomain{FT}(300.0) mesh = Meshes.EquiangularCubedSphere(domain, 3) @@ -111,3 +117,29 @@ end n_dss_affected_y123 == n_dss_affected_y3 * 3 end + +function test_dss_conservation(space) + Random.seed!(1) # ensures reproducibility + field = zeros(space) + FT = Spaces.undertype(space) + parent(field) .= rand.(FT) + 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) +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) + center_space = + TU.CenterExtrudedFiniteDifferenceSpace(FT; space_kwargs...) + for space in (center_space, Spaces.face_space(center_space)) + test_dss_conservation(space) + end + end + end +end diff --git a/test/Spaces/unit_spaces.jl b/test/Spaces/unit_spaces.jl index 313d137506..4a691db084 100644 --- a/test/Spaces/unit_spaces.jl +++ b/test/Spaces/unit_spaces.jl @@ -242,7 +242,7 @@ end @test Spaces.local_geometry_type(typeof(space)) <: Geometry.LocalGeometry local_geometry_slab = slab(Spaces.local_geometry_data(space), 1) - dss_weights_slab = slab(Spaces.local_dss_weights(space), 1) + dss_weights_slab = slab(Spaces.dss_weights(space), 1) @static if on_gpu adapted_space = Adapt.adapt(CUDA.KernelAdaptor(), space)