From aec1fbf2aa581f81aac4e8ea0c5467163bcebf7f Mon Sep 17 00:00:00 2001 From: akshaysridhar Date: Wed, 31 Jul 2024 09:46:53 -0700 Subject: [PATCH] For Covariant123Vector types, transforms the 3rd velocity component instead of treating it as a scalar variable modified: ext/cuda/topologies_dss.jl modified: src/Topologies/dss.jl modified: src/Topologies/dss_transform.jl modified: test/Spaces/terrain_warp.jl modified: test/Spaces/ddss1_cs.jl +format Remove dss-op on surface elevation horizontal gradient modified: src/Hypsography/Hypsography.jl modified: test/Spaces/terrain_warp.jl modified: test/Spaces/terrain_warp.jl --- ext/cuda/topologies_dss.jl | 114 ++++++++++++++++++ src/Hypsography/Hypsography.jl | 2 - src/Topologies/dss.jl | 205 ++++++++++++++++++++++++++++++-- src/Topologies/dss_transform.jl | 37 ++++-- test/Spaces/ddss1_cs.jl | 2 +- test/Spaces/terrain_warp.jl | 107 +++++++++++++++++ 6 files changed, 443 insertions(+), 24 deletions(-) diff --git a/ext/cuda/topologies_dss.jl b/ext/cuda/topologies_dss.jl index 1ce3355c16..e0ba96f04f 100644 --- a/ext/cuda/topologies_dss.jl +++ b/ext/cuda/topologies_dss.jl @@ -188,6 +188,8 @@ function Topologies.dss_transform!( scalarfidx::AbstractVector{Int}, covariant12fidx::AbstractVector{Int}, contravariant12fidx::AbstractVector{Int}, + covariant123fidx::AbstractVector{Int}, + contravariant123fidx::AbstractVector{Int}, localelems::AbstractVector{Int}, ) nlocalelems = length(localelems) @@ -212,6 +214,8 @@ function Topologies.dss_transform!( scalarfidx, covariant12fidx, contravariant12fidx, + covariant123fidx, + contravariant123fidx, localelems, ) auto_launch!( @@ -236,6 +240,8 @@ function dss_transform_kernel!( scalarfidx::AbstractVector{Int}, covariant12fidx::AbstractVector{Int}, contravariant12fidx::AbstractVector{Int}, + covariant123fidx::AbstractVector{Int}, + contravariant123fidx::AbstractVector{Int}, localelems::AbstractVector{Int}, ) where {FT <: AbstractFloat, Nq} gidx = threadIdx().x + (blockIdx().x - 1) * blockDim().x @@ -288,6 +294,70 @@ function dss_transform_kernel!( p∂x∂ξ[idx22] * pdata[data_idx2] ) * weight end + for fidx in covariant123fidx + data_idx1 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + + (idx11, idx12, idx13, idx21, idx22, idx23, idx31, idx32, idx33) = + Topologies._get_idx_metric_3d( + sizet_metric, + (level, ip, jp, elem), + ) + + # Covariant to physical transformation + pperimeter_data[level, p, fidx, elem] = + ( + p∂ξ∂x[idx11] * pdata[data_idx1] + + p∂ξ∂x[idx12] * pdata[data_idx2] + + p∂ξ∂x[idx13] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 1, elem] = + ( + p∂ξ∂x[idx21] * pdata[data_idx1] + + p∂ξ∂x[idx22] * pdata[data_idx2] + + p∂ξ∂x[idx23] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 2, elem] = + ( + p∂ξ∂x[idx31] * pdata[data_idx1] + + p∂ξ∂x[idx32] * pdata[data_idx2] + + p∂ξ∂x[idx33] * pdata[data_idx3] + ) * pw + end + + for fidx in contravariant123fidx + data_idx1 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + Topologies._get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + (idx11, idx12, idx13, idx21, idx22, idx23, idx31, idx32, idx33) = + _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + # Contravariant to physical transformation + pperimeter_data[level, p, fidx, elem] = + ( + p∂x∂ξ[idx11] * pdata[data_idx1] + + p∂x∂ξ[idx21] * pdata[data_idx2] + + p∂x∂ξ[idx31] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 1, elem] = + ( + p∂x∂ξ[idx12] * pdata[data_idx1] + + p∂x∂ξ[idx22] * pdata[data_idx2] + + p∂x∂ξ[idx32] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 2, elem] = + ( + p∂x∂ξ[idx13] * pdata[data_idx1] + + p∂x∂ξ[idx23] * pdata[data_idx2] + + p∂x∂ξ[idx33] * pdata[data_idx3] + ) * pw + end end return nothing end @@ -302,6 +372,8 @@ function Topologies.dss_untransform!( scalarfidx::AbstractVector{Int}, covariant12fidx::AbstractVector{Int}, contravariant12fidx::AbstractVector{Int}, + covariant123fidx::AbstractVector{Int}, + contravariant123fidx::AbstractVector{Int}, localelems::AbstractVector{Int}, ) where {Nq} nlocalelems = length(localelems) @@ -324,6 +396,8 @@ function Topologies.dss_untransform!( scalarfidx, covariant12fidx, contravariant12fidx, + covariant123fidx, + contravariant123fidx, localelems, ) auto_launch!( @@ -347,6 +421,8 @@ function dss_untransform_kernel!( scalarfidx::AbstractVector{Int}, covariant12fidx::AbstractVector{Int}, contravariant12fidx::AbstractVector{Int}, + covariant123fidx::AbstractVector{Int}, + contravariant123fidx::AbstractVector{Int}, localelems::AbstractVector{Int}, ) where {FT <: AbstractFloat, Nq} gidx = threadIdx().x + (blockIdx().x - 1) * blockDim().x @@ -389,6 +465,44 @@ function dss_untransform_kernel!( p∂ξ∂x[idx12] * pperimeter_data[level, p, fidx, elem] + p∂ξ∂x[idx22] * pperimeter_data[level, p, fidx + 1, elem] end + for fidx in covariant123fidx + data_idx1 = _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + (idx11, idx12, idx13, idx21, idx22, idx23, idx31, idx32, idx33) = + _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + pdata[data_idx1] = + p∂x∂ξ[idx11] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx12] * pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx13] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx2] = + p∂x∂ξ[idx21] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx22] * pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx23] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx3] = + p∂x∂ξ[idx31] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx32] * pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx33] * pperimeter_data[level, p, fidx + 2, elem] + end + for fidx in contravariant123fidx + data_idx1 = _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + (idx11, idx12, idx13, idx21, idx22, idx23, idx31, idx32, idx33) = + _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + pdata[data_idx1] = + p∂ξ∂x[idx11] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx21] * pperimeter_data[level, p, fidx + 1, elem] + + p∂ξ∂x[idx31] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx2] = + p∂ξ∂x[idx12] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx22] * pperimeter_data[level, p, fidx + 1, elem] + p∂ξ∂x[idx32] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx3] = + p∂ξ∂x[idx13] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx23] * pperimeter_data[level, p, fidx + 1, elem] + + p∂ξ∂x[idx33] * pperimeter_data[level, p, fidx + 2, elem] + end end return nothing end diff --git a/src/Hypsography/Hypsography.jl b/src/Hypsography/Hypsography.jl index 75ede63f9c..c1ea038a5c 100644 --- a/src/Hypsography/Hypsography.jl +++ b/src/Hypsography/Hypsography.jl @@ -239,8 +239,6 @@ function _ExtrudedFiniteDifferenceGrid( f = Spaces.create_dss_buffer(face_∇Z_field), ) - Spaces.weighted_dss!(center_∇Z_field => buffer.c, face_∇Z_field => buffer.f) - # construct full local geometry center_local_geometry = Geometry.product_geometry.( diff --git a/src/Topologies/dss.jl b/src/Topologies/dss.jl index dce465412d..49b985118c 100644 --- a/src/Topologies/dss.jl +++ b/src/Topologies/dss.jl @@ -29,6 +29,10 @@ struct DSSBuffer{S, G, D, A, B, VI} "field id for all contravariant12vector fields stored in the `data` array" contravariant12fidx::VI "internal local elements (lidx)" + covariant123fidx::VI + "field id for all contravariant123vector fields stored in the `data` array" + contravariant123fidx::VI + "internal local elements (lidx)" internal_elems::VI "local elements (lidx) located on process boundary" perimeter_elems::VI @@ -107,7 +111,11 @@ function create_dss_buffer( internal_elems = DA(topology.internal_elems) perimeter_elems = DA(topology.perimeter_elems) end - scalarfidx, covariant12fidx, contravariant12fidx = Int[], Int[], Int[] + scalarfidx, + covariant12fidx, + contravariant12fidx, + covariant123fidx, + contravariant123fidx = Int[], Int[], Int[], Int[], Int[] supportedvectortypes = Union{ Geometry.UVector, Geometry.VVector, @@ -119,6 +127,7 @@ function create_dss_buffer( Geometry.Covariant12Vector, Geometry.Covariant3Vector, Geometry.Covariant123Vector, + Geometry.Contravariant123Vector, Geometry.Contravariant12Vector, Geometry.Contravariant3Vector, } @@ -136,10 +145,11 @@ function create_dss_buffer( if fieldtype <: Geometry.Covariant12Vector push!(covariant12fidx, offset + 1) elseif fieldtype <: Geometry.Covariant123Vector - push!(covariant12fidx, offset + 1) - push!(scalarfidx, offset + 3) + push!(covariant123fidx, offset + 1) elseif fieldtype <: Geometry.Contravariant12Vector push!(contravariant12fidx, offset + 1) + elseif fieldtype <: Geometry.Contravariant123Vector + push!(contravariant123fidx, offset + 1) else append!( scalarfidx, @@ -163,10 +173,11 @@ function create_dss_buffer( if S <: Geometry.Covariant12Vector push!(covariant12fidx, 1) elseif S <: Geometry.Covariant123Vector - push!(covariant12fidx, 1) - push!(scalarfidx, 3) + push!(covariant123fidx, 1) elseif S <: Geometry.Contravariant12Vector push!(contravariant12fidx, 1) + elseif S <: Geometry.Contravariant123Vector + push!(contravariant123fidx, 1) else append!(scalarfidx, Vector(1:ncomponents)) end @@ -178,7 +189,9 @@ function create_dss_buffer( end scalarfidx = DA(scalarfidx) covariant12fidx = DA(covariant12fidx) + covariant123fidx = DA(covariant123fidx) contravariant12fidx = DA(contravariant12fidx) + contravariant123fidx = DA(contravariant123fidx) G = typeof(graph_context) D = typeof(perimeter_data) A = typeof(send_data) @@ -194,6 +207,8 @@ function create_dss_buffer( scalarfidx, covariant12fidx, contravariant12fidx, + covariant123fidx, + contravariant123fidx, internal_elems, perimeter_elems, ) @@ -243,8 +258,14 @@ function dss_transform!( localelems::AbstractVector{Int}, ) if !isempty(localelems) - (; scalarfidx, covariant12fidx, contravariant12fidx, perimeter_data) = - dss_buffer + (; + scalarfidx, + covariant12fidx, + contravariant12fidx, + covariant123fidx, + contravariant123fidx, + perimeter_data, + ) = dss_buffer (; ∂ξ∂x, ∂x∂ξ) = local_geometry dss_transform!( device, @@ -257,6 +278,8 @@ function dss_transform!( scalarfidx, covariant12fidx, contravariant12fidx, + covariant123fidx, + contravariant123fidx, localelems, ) end @@ -293,8 +316,14 @@ function dss_untransform!( perimeter::Perimeter2D, localelems::AbstractVector{Int}, ) - (; scalarfidx, covariant12fidx, contravariant12fidx, perimeter_data) = - dss_buffer + (; + scalarfidx, + covariant12fidx, + contravariant12fidx, + covariant123fidx, + contravariant123fidx, + perimeter_data, + ) = dss_buffer (; ∂ξ∂x, ∂x∂ξ) = local_geometry dss_untransform!( device, @@ -306,6 +335,8 @@ function dss_untransform!( scalarfidx, covariant12fidx, contravariant12fidx, + covariant123fidx, + contravariant123fidx, localelems, ) return nothing @@ -323,6 +354,8 @@ end scalarfidx::Vector{Int}, covariant12fidx::Vector{Int}, contravariant12fidx::Vector{Int}, + covariant123fidx::Vector{Int}, + contravariant123fidx::Vector{Int}, localelems::Vector{Int}, ) @@ -338,6 +371,7 @@ Arguments: - `perimeter`: perimeter iterator - `scalarfidx`: field index for scalar fields in the data layout - `covariant12fidx`: field index for Covariant12 vector fields in the data layout +- `covariant123fidx`: field index for Covariant123 vector fields in the data layout - `localelems`: list of local elements to perform transformation operations on Part of [`ClimaCore.Spaces.weighted_dss!`](@ref). @@ -353,6 +387,8 @@ function dss_transform!( scalarfidx::Vector{Int}, covariant12fidx::Vector{Int}, contravariant12fidx::Vector{Int}, + covariant123fidx::Vector{Int}, + contravariant123fidx::Vector{Int}, localelems::Vector{Int}, ) where {Nq} pdata = parent(data) @@ -410,6 +446,84 @@ function dss_transform!( p∂x∂ξ[idx22] * pdata[data_idx2] ) * pw end + + for fidx in covariant123fidx, level in 1:nlevels + data_idx1 = _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + + ( + idx11, + idx12, + idx13, + idx21, + idx22, + idx23, + idx31, + idx32, + idx33, + ) = _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + + # Covariant to physical transformation + pperimeter_data[level, p, fidx, elem] = + ( + p∂ξ∂x[idx11] * pdata[data_idx1] + + p∂ξ∂x[idx12] * pdata[data_idx2] + + p∂ξ∂x[idx13] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 1, elem] = + ( + p∂ξ∂x[idx21] * pdata[data_idx1] + + p∂ξ∂x[idx22] * pdata[data_idx2] + + p∂ξ∂x[idx23] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 2, elem] = + ( + p∂ξ∂x[idx31] * pdata[data_idx1] + + p∂ξ∂x[idx32] * pdata[data_idx2] + + p∂ξ∂x[idx33] * pdata[data_idx3] + ) * pw + end + + for fidx in contravariant123fidx, level in 1:nlevels + data_idx1 = _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + ( + idx11, + idx12, + idx13, + idx21, + idx22, + idx23, + idx31, + idx32, + idx33, + ) = _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + # Contravariant to physical transformation + pperimeter_data[level, p, fidx, elem] = + ( + p∂x∂ξ[idx11] * pdata[data_idx1] + + p∂x∂ξ[idx21] * pdata[data_idx2] + + p∂x∂ξ[idx31] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 1, elem] = + ( + p∂x∂ξ[idx12] * pdata[data_idx1] + + p∂x∂ξ[idx22] * pdata[data_idx2] + + p∂x∂ξ[idx32] * pdata[data_idx3] + ) * pw + pperimeter_data[level, p, fidx + 2, elem] = + ( + p∂x∂ξ[idx13] * pdata[data_idx1] + + p∂x∂ξ[idx23] * pdata[data_idx2] + + p∂x∂ξ[idx33] * pdata[data_idx3] + ) * pw + end end end return nothing @@ -453,6 +567,8 @@ function dss_untransform!( scalarfidx::Vector{Int}, covariant12fidx::Vector{Int}, contravariant12fidx::Vector{Int}, + covariant123fidx::Vector{Int}, + contravariant123fidx::Vector{Int}, localelems::Vector{Int}, ) where {Nq} pdata = parent(data) @@ -505,6 +621,77 @@ function dss_untransform!( p∂ξ∂x[idx22] * pperimeter_data[level, p, fidx + 1, elem] end end + for fidx in covariant123fidx + for level in 1:nlevels + data_idx1 = + _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + ( + idx11, + idx12, + idx13, + idx21, + idx22, + idx23, + idx31, + idx32, + idx33, + ) = _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + pdata[data_idx1] = + p∂x∂ξ[idx11] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx12] * + pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx13] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx2] = + p∂x∂ξ[idx21] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx22] * + pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx23] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx3] = + p∂x∂ξ[idx31] * pperimeter_data[level, p, fidx, elem] + + p∂x∂ξ[idx32] * + pperimeter_data[level, p, fidx + 1, elem] + + p∂x∂ξ[idx33] * pperimeter_data[level, p, fidx + 2, elem] + end + end + for fidx in contravariant123fidx + for level in 1:nlevels + data_idx1 = + _get_idx(sizet_data, (level, ip, jp, fidx, elem)) + data_idx2 = + _get_idx(sizet_data, (level, ip, jp, fidx + 1, elem)) + data_idx3 = + _get_idx(sizet_data, (level, ip, jp, fidx + 2, elem)) + + ( + idx11, + idx12, + idx13, + idx21, + idx22, + idx23, + idx31, + idx32, + idx33, + ) = _get_idx_metric_3d(sizet_metric, (level, ip, jp, elem)) + pdata[data_idx1] = + p∂ξ∂x[idx11] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx21] * pperimeter_data[level, p, fidx + 1, elem] + p∂ξ∂x[idx31] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx2] = + p∂ξ∂x[idx12] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx22] * pperimeter_data[level, p, fidx + 1, elem] + p∂ξ∂x[idx32] * pperimeter_data[level, p, fidx + 2, elem] + pdata[data_idx3] = + p∂ξ∂x[idx13] * pperimeter_data[level, p, fidx, elem] + + p∂ξ∂x[idx23] * + pperimeter_data[level, p, fidx + 1, elem] + + p∂ξ∂x[idx33] * pperimeter_data[level, p, fidx + 2, elem] + end + end end end return nothing diff --git a/src/Topologies/dss_transform.jl b/src/Topologies/dss_transform.jl index a03ab3c79c..282ae28a88 100644 --- a/src/Topologies/dss_transform.jl +++ b/src/Topologies/dss_transform.jl @@ -233,18 +233,31 @@ end # helper functions for DSS2 function _get_idx_metric(sizet::NTuple{5, Int}, loc::NTuple{4, Int}) - @inbounds begin - nmetric = sizet[4] - (i11, i12, i21, i22) = nmetric == 4 ? (1, 2, 3, 4) : (1, 2, 4, 5) - (level, i, j, elem) = loc - inds = ( - linear_ind(sizet, (level, i, j, i11, elem)), - linear_ind(sizet, (level, i, j, i12, elem)), - linear_ind(sizet, (level, i, j, i21, elem)), - linear_ind(sizet, (level, i, j, i22, elem)), - ) - return inds - end + nmetric = sizet[4] + (i11, i12, i21, i22) = nmetric == 4 ? (1, 2, 3, 4) : (1, 2, 4, 5) + (level, i, j, elem) = loc + return ( + _get_idx(sizet, (level, i, j, i11, elem)), + _get_idx(sizet, (level, i, j, i12, elem)), + _get_idx(sizet, (level, i, j, i21, elem)), + _get_idx(sizet, (level, i, j, i22, elem)), + ) +end + +function _get_idx_metric_3d(sizet::NTuple{5, Int}, loc::NTuple{4, Int}) + nmetric = sizet[4] + (level, i, j, elem) = loc + return ( + _get_idx(sizet, (level, i, j, 1, elem)), + _get_idx(sizet, (level, i, j, 2, elem)), + _get_idx(sizet, (level, i, j, 3, elem)), + _get_idx(sizet, (level, i, j, 4, elem)), + _get_idx(sizet, (level, i, j, 5, elem)), + _get_idx(sizet, (level, i, j, 6, elem)), + _get_idx(sizet, (level, i, j, 7, elem)), + _get_idx(sizet, (level, i, j, 8, elem)), + _get_idx(sizet, (level, i, j, 9, elem)), + ) end function _representative_slab( diff --git a/test/Spaces/ddss1_cs.jl b/test/Spaces/ddss1_cs.jl index 89462bdcf6..cf1f6b8024 100644 --- a/test/Spaces/ddss1_cs.jl +++ b/test/Spaces/ddss1_cs.jl @@ -58,7 +58,7 @@ import ClimaCore: ) end -@testset "DSS of Covarinat12Vector & Covariant123Vector on extruded Cubed Sphere mesh (ne = 3, serial run)" begin +@testset "DSS of Covariant12Vector & Covariant123Vector on extruded Cubed Sphere mesh (ne = 3, serial run)" begin FT = Float64 context = ClimaComms.SingletonCommsContext(ClimaComms.CUDADevice()) context_cpu = diff --git a/test/Spaces/terrain_warp.jl b/test/Spaces/terrain_warp.jl index f1ac451fe2..d4baefacba 100644 --- a/test/Spaces/terrain_warp.jl +++ b/test/Spaces/terrain_warp.jl @@ -369,6 +369,113 @@ end end end +@testset "3D Extruded Terrain Warped Space: DSS" begin + context_gen = ClimaComms.SingletonCommsContext(ClimaComms.device()) + context_cpu = + ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()) # CPU context for comparison + for FT in (Float32, Float64) + # Extruded FD-Spectral Hybrid + xmin, xmax = FT(0), FT(π) + ymin, ymax = FT(0), FT(π) + zmin, zmax = FT(0), FT(1) + levels = (5, 10) + polynom = (2, 5, 10) + horzelem = (2, 5, 10) + for nl in levels, np in polynom, nh in horzelem + hv_center_space, hv_face_space = warpedspace_3D( + FT, + (xmin, xmax), + (ymin, ymax), + (zmin, zmax), + nh, + nl, + np; + context = context_gen, + ) + hv_center_space_cpu, hv_face_space_cpu = warpedspace_3D( + FT, + (xmin, xmax), + (ymin, ymax), + (zmin, zmax), + nh, + nl, + np; + context = context_cpu, + ) + ᶜcoords = Fields.coordinate_field(hv_center_space) + ᶠcoords = Fields.coordinate_field(hv_face_space) + ᶜcoords_cpu = Fields.coordinate_field(hv_center_space_cpu) + ᶠcoords_cpu = Fields.coordinate_field(hv_face_space_cpu) + + y1 = one.(Fields.coordinate_field(hv_center_space).z) + y2 = -1 .* y1 + y3 = y1 + y12 = @. Geometry.Covariant12Vector(y1, y2) + + y1_cpu = one.(Fields.coordinate_field(hv_center_space_cpu).z) + y2_cpu = -1 .* y1_cpu + y3_cpu = @. y1_cpu + y12_cpu = @. Geometry.Covariant12Vector(y1_cpu, y2_cpu) + + dss_buffer12 = Spaces.create_dss_buffer(y12) + dss_buffer12_cpu = Spaces.create_dss_buffer(y12_cpu) + + Spaces.weighted_dss!(y12 => dss_buffer12) + Spaces.weighted_dss!(y12_cpu => dss_buffer12_cpu) + + yinit12 = copy(y12) + yinit12_cpu = copy(y12_cpu) + + @test yinit12 ≈ y12 + @test yinit12_cpu ≈ y12_cpu + @test parent(y12_cpu) ≈ Array(parent(y12)) + + # test DSS for a Covariant123Vector + y123 = @. Geometry.Covariant123Vector(y1, y2, y3) + y123_cpu = @. Geometry.Covariant123Vector(y1_cpu, y2_cpu, y3_cpu) + + dss_buffer123 = Spaces.create_dss_buffer(y123) + dss_buffer123_cpu = Spaces.create_dss_buffer(y123_cpu) + + # ensure physical velocity is continous across SE boundary for initial state + Spaces.weighted_dss!(y123 => dss_buffer123) + Spaces.weighted_dss!(y123_cpu => dss_buffer123_cpu) + + yinit123 = copy(y123) + yinit123_cpu = copy(y123_cpu) + + Spaces.weighted_dss!(y123, dss_buffer123) + Spaces.weighted_dss!(y123_cpu, dss_buffer123_cpu) + + @test yinit123 ≈ y123 + @test yinit123_cpu ≈ y123_cpu + @test parent(y123_cpu) ≈ Array(parent(y123)) + + # test DSS for a Contravariant123Vector + y123 = @. Geometry.Contravariant123Vector(y1, y2, y3) + y123_cpu = + @. Geometry.Contravariant123Vector(y1_cpu, y2_cpu, y3_cpu) + + dss_buffer123 = Spaces.create_dss_buffer(y123) + dss_buffer123_cpu = Spaces.create_dss_buffer(y123_cpu) + + # ensure physical velocity is continous across SE boundary for initial state + Spaces.weighted_dss!(y123 => dss_buffer123) + Spaces.weighted_dss!(y123_cpu => dss_buffer123_cpu) + + yinit123 = copy(y123) + yinit123_cpu = copy(y123_cpu) + + Spaces.weighted_dss!(y123, dss_buffer123) + Spaces.weighted_dss!(y123_cpu, dss_buffer123_cpu) + + @test yinit123 ≈ y123 + @test yinit123_cpu ≈ y123_cpu + @test parent(y123_cpu) ≈ Array(parent(y123)) + end + end +end + @testset "3D Extruded Terrain Laplacian Smoothing" begin # Test smoothing for known parameters device = ClimaComms.device()