From b3780961fab739dfada1a05fa700cfb32d84f127 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Thu, 11 Jan 2024 15:55:33 -0800 Subject: [PATCH] Add option for deep atmosphere Co-authored-by: Akshay Sridhar Co-authored-by: Gabriele Bozzola --- examples/hybrid/plane/topo_agnesi_nh.jl | 5 +- examples/hybrid/plane/topo_schar_nh.jl | 5 +- examples/hybrid/sphere/deformation_flow.jl | 11 +- src/DataLayouts/broadcast.jl | 10 +- src/Geometry/globalgeometry.jl | 145 ++++++-- src/Geometry/localgeometry.jl | 130 +++++-- src/Grids/extruded.jl | 36 +- src/Grids/finitedifference.jl | 220 ++++++++---- src/Hypsography/Hypsography.jl | 395 ++++++++++----------- src/InputOutput/readers.jl | 7 + src/InputOutput/writers.jl | 16 + src/Remapping/distributed_remapping.jl | 2 +- src/Remapping/interpolate_array.jl | 16 +- test/Geometry/geometry.jl | 36 +- test/Hypsography/2d.jl | 30 ++ test/Operators/hybrid/3d.jl | 44 ++- test/Spaces/sphere.jl | 30 +- test/Spaces/terrain_warp.jl | 23 +- 18 files changed, 767 insertions(+), 394 deletions(-) diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl index 709eff23a9..283469fc31 100644 --- a/examples/hybrid/plane/topo_agnesi_nh.jl +++ b/examples/hybrid/plane/topo_agnesi_nh.jl @@ -69,12 +69,11 @@ function hvspace_2D( quad = Quadratures.GLL{npoly + 1}() horzspace = Spaces.SpectralElementSpace1D(horztopology, quad) - z_surface = warp_fn.(Fields.coordinate_field(horzspace)) + z_surface = Geometry.ZPoint.(warp_fn.(Fields.coordinate_field(horzspace))) hv_face_space = Spaces.ExtrudedFiniteDifferenceSpace( horzspace, vert_face_space, - Hypsography.LinearAdaption(), - z_surface, + Hypsography.LinearAdaption(z_surface), ) hv_center_space = Spaces.CenterExtrudedFiniteDifferenceSpace(hv_face_space) return (hv_center_space, hv_face_space) diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl index bbd9bd17d2..5c0225f851 100644 --- a/examples/hybrid/plane/topo_schar_nh.jl +++ b/examples/hybrid/plane/topo_schar_nh.jl @@ -92,12 +92,11 @@ function hvspace_2D( quad = Quadratures.GLL{npoly + 1}() horzspace = Spaces.SpectralElementSpace1D(horztopology, quad) - z_surface = warp_fn.(Fields.coordinate_field(horzspace)) + z_surface = Geometry.ZPoint.(warp_fn.(Fields.coordinate_field(horzspace))) hv_face_space = Spaces.ExtrudedFiniteDifferenceSpace( horzspace, vert_face_space, - Hypsography.LinearAdaption(), - z_surface, + Hypsography.LinearAdaption(z_surface), ) hv_center_space = Spaces.CenterExtrudedFiniteDifferenceSpace(hv_face_space) return (hv_center_space, hv_face_space) diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl index 34ebb2a117..6f102e8d27 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -221,12 +221,15 @@ function run_deformation_flow(use_limiter, fct_op) zd = z - z_c centers = ( - Geometry.LatLongZPoint(ϕ_c, λ_c1, FT(0)), - Geometry.LatLongZPoint(ϕ_c, λ_c2, FT(0)), + Geometry.LatLongZPoint(ϕ_c, λ_c1, z), + Geometry.LatLongZPoint(ϕ_c, λ_c2, z), ) - horz_geometry = Spaces.global_geometry(horz_space) rds = map(centers) do center - Geometry.great_circle_distance(coord, center, horz_geometry) + Geometry.great_circle_distance( + coord, + center, + Spaces.global_geometry(cent_space), + ) end ds = @. min(1, (rds / R_t)^2 + (zd / Z_t)^2) # scaled distance functions diff --git a/src/DataLayouts/broadcast.jl b/src/DataLayouts/broadcast.jl index a6cdac52f5..4a8b722fd2 100644 --- a/src/DataLayouts/broadcast.jl +++ b/src/DataLayouts/broadcast.jl @@ -234,8 +234,8 @@ end function Base.similar( bc::Union{IJFH{<:Any, Nij, A}, Broadcast.Broadcasted{IJFHStyle{Nij, A}}}, ::Type{Eltype}, + (_, _, _, _, Nh) = size(bc), ) where {Nij, A, Eltype} - _, _, _, _, Nh = size(bc) PA = parent_array_type(A) array = similar(PA, (Nij, Nij, typesize(eltype(A), Eltype), Nh)) return IJFH{Eltype, Nij}(array) @@ -244,8 +244,8 @@ end function Base.similar( bc::Union{IFH{<:Any, Ni, A}, Broadcast.Broadcasted{IFHStyle{Ni, A}}}, ::Type{Eltype}, + (_, _, _, _, Nh) = size(bc), ) where {Ni, A, Eltype} - _, _, _, _, Nh = size(bc) PA = parent_array_type(A) array = similar(PA, (Ni, typesize(eltype(A), Eltype), Nh)) return IFH{Eltype, Ni}(array) @@ -272,8 +272,8 @@ end function Base.similar( bc::Union{VF{<:Any, A}, Broadcast.Broadcasted{VFStyle{A}}}, ::Type{Eltype}, + (_, _, _, Nv, _) = size(bc), ) where {A, Eltype} - _, _, _, Nv, _ = size(bc) PA = parent_array_type(A) array = similar(PA, (Nv, typesize(eltype(A), Eltype))) return VF{Eltype}(array) @@ -282,8 +282,8 @@ end function Base.similar( bc::Union{VIFH{<:Any, Ni, A}, Broadcast.Broadcasted{VIFHStyle{Ni, A}}}, ::Type{Eltype}, + (_, _, _, Nv, Nh) = size(bc), ) where {Ni, A, Eltype} - _, _, _, Nv, Nh = size(bc) PA = parent_array_type(A) array = similar(PA, (Nv, Ni, typesize(eltype(A), Eltype), Nh)) return VIFH{Eltype, Ni}(array) @@ -292,8 +292,8 @@ end function Base.similar( bc::Union{VIJFH{<:Any, Nij, A}, Broadcast.Broadcasted{VIJFHStyle{Nij, A}}}, ::Type{Eltype}, + (_, _, _, Nv, Nh) = size(bc), ) where {Nij, A, Eltype} - _, _, _, Nv, Nh = size(bc) PA = parent_array_type(A) array = similar(PA, (Nv, Nij, Nij, typesize(eltype(A), Eltype), Nh)) return VIJFH{Eltype, Nij}(array) diff --git a/src/Geometry/globalgeometry.jl b/src/Geometry/globalgeometry.jl index 0c5875fb33..7673084d03 100644 --- a/src/Geometry/globalgeometry.jl +++ b/src/Geometry/globalgeometry.jl @@ -61,6 +61,7 @@ LocalVector(u::CartesianVector{T,I}, ::CartesianGlobalGeometry) where {T,I} = AxisVector(LocalAxis{I}(), components(u)) =# +abstract type AbstractSphericalGlobalGeometry <: AbstractGlobalGeometry end """ SphericalGlobalGeometry(radius) @@ -79,19 +80,44 @@ along the zero longitude line: direction, `v` being aligned with the `x1` direction, and `w` being aligned with the negative `x3` direction. """ -struct SphericalGlobalGeometry{FT} <: AbstractGlobalGeometry +struct SphericalGlobalGeometry{FT} <: AbstractSphericalGlobalGeometry + radius::FT +end + + +""" + ShallowSphericalGlobalGeometry(radius) + +Similar to [`SphericalGlobalGeometry`](@ref), but for extruded spheres. In this +case, it uses the "shallow-atmosphere" assumption that circumference is the same +at all `z`. +""" +struct ShallowSphericalGlobalGeometry{FT} <: AbstractSphericalGlobalGeometry + radius::FT +end + +""" + DeepSphericalGlobalGeometry(radius) + +Similar to [`SphericalGlobalGeometry`](@ref), but for extruded spheres. In this +case, it uses the "deep-atmosphere" assumption that circumference increases with `z`. +""" +struct DeepSphericalGlobalGeometry{FT} <: AbstractSphericalGlobalGeometry radius::FT end # coordinates -function CartesianPoint(pt::LatLongPoint, global_geom::SphericalGlobalGeometry) +function CartesianPoint( + pt::LatLongPoint, + global_geom::AbstractSphericalGlobalGeometry, +) r = global_geom.radius x1 = r * cosd(pt.long) * cosd(pt.lat) x2 = r * sind(pt.long) * cosd(pt.lat) x3 = r * sind(pt.lat) Cartesian123Point(x1, x2, x3) end -function LatLongPoint(pt::Cartesian123Point, ::SphericalGlobalGeometry) +function LatLongPoint(pt::Cartesian123Point, ::AbstractSphericalGlobalGeometry) ϕ = atand(pt.x3, hypot(pt.x2, pt.x1)) # IEEE754 spec states that atand(±0.0, −0.0) == ±180, however to make the UV # orienation consistent, we define the longitude to be zero at the poles @@ -104,7 +130,10 @@ function LatLongPoint(pt::Cartesian123Point, ::SphericalGlobalGeometry) end -function CartesianPoint(pt::LatLongZPoint, global_geom::SphericalGlobalGeometry) +function CartesianPoint( + pt::LatLongZPoint, + global_geom::AbstractSphericalGlobalGeometry, +) r = global_geom.radius z = pt.z x1 = (r + z) * cosd(pt.long) * cosd(pt.lat) @@ -114,30 +143,21 @@ function CartesianPoint(pt::LatLongZPoint, global_geom::SphericalGlobalGeometry) end function LatLongZPoint( pt::Cartesian123Point, - global_geom::SphericalGlobalGeometry, + global_geom::AbstractSphericalGlobalGeometry, ) llpt = LatLongPoint(pt, global_geom) z = hypot(pt.x1, pt.x2, pt.x3) - global_geom.radius LatLongZPoint(llpt.lat, llpt.long, z) end -""" - great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry) -Compute the great circle (spherical geodesic) distance between `pt1` and `pt2`. -""" -function great_circle_distance( - pt1::LatLongPoint, - pt2::LatLongPoint, - global_geom::SphericalGlobalGeometry, -) - r = global_geom.radius +function unit_great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint) ϕ1 = pt1.lat λ1 = pt1.long ϕ2 = pt2.lat λ2 = pt2.long Δλ = λ1 - λ2 - return r * atan( + return atan( hypot( cosd(ϕ2) * sind(Δλ), cosd(ϕ1) * sind(ϕ2) - sind(ϕ1) * cosd(ϕ2) * cosd(Δλ), @@ -146,6 +166,21 @@ function great_circle_distance( ) end + +""" + great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry) + +Compute the great circle (spherical geodesic) distance between `pt1` and `pt2`. +""" +function great_circle_distance( + pt1::LatLongPoint, + pt2::LatLongPoint, + global_geom::AbstractSphericalGlobalGeometry, +) + r = global_geom.radius + return r * unit_great_circle_distance(pt1, pt2) +end + """ great_circle_distance(pt1::LatLongZPoint, pt2::LatLongZPoint, global_geometry::SphericalGlobalGeometry) @@ -154,15 +189,26 @@ Compute the great circle (spherical geodesic) distance between `pt1` and `pt2`. function great_circle_distance( pt1::LatLongZPoint, pt2::LatLongZPoint, - global_geom::SphericalGlobalGeometry, + global_geom::ShallowSphericalGlobalGeometry, ) - ϕ1 = pt1.lat - λ1 = pt1.long - ϕ2 = pt2.lat - λ2 = pt2.long - latlong_pt1 = LatLongPoint(ϕ1, λ1) - latlong_pt2 = LatLongPoint(ϕ2, λ2) - return great_circle_distance(latlong_pt1, latlong_pt2, global_geom) + r = global_geom.radius + return r * unit_great_circle_distance( + LatLongPoint(pt1.lat, pt1.long), + LatLongPoint(pt2.lat, pt2.long), + ) +end + +function great_circle_distance( + pt1::LatLongZPoint, + pt2::LatLongZPoint, + global_geom::DeepSphericalGlobalGeometry, +) + r = global_geom.radius + R = r + (pt1.z + pt2.z) / 2 + return R * unit_great_circle_distance( + LatLongPoint(pt1.lat, pt1.long), + LatLongPoint(pt2.lat, pt2.long), + ) end """ @@ -180,7 +226,7 @@ end # vectors CartesianVector( u::CartesianVector, - ::SphericalGlobalGeometry, + ::AbstractSphericalGlobalGeometry, ::LocalGeometry, ) = u CartesianVector( @@ -193,7 +239,7 @@ CartesianVector( local_geometry.coordinates, ) function local_to_cartesian( - ::SphericalGlobalGeometry, + ::AbstractSphericalGlobalGeometry, coord::Union{LatLongPoint, LatLongZPoint}, ) ϕ = coord.lat @@ -213,7 +259,7 @@ end function CartesianVector( u::UVWVector, - geom::SphericalGlobalGeometry, + geom::AbstractSphericalGlobalGeometry, coord::Union{LatLongPoint, LatLongZPoint}, ) G = local_to_cartesian(geom, coord) @@ -221,9 +267,52 @@ function CartesianVector( end function LocalVector( u::Cartesian123Vector, - geom::SphericalGlobalGeometry, + geom::AbstractSphericalGlobalGeometry, coord::Union{LatLongPoint, LatLongZPoint}, ) G = local_to_cartesian(geom, coord) G' * u end + +function product_geometry( + horizontal_local_geometry::Geometry.LocalGeometry, + vertical_local_geometry::Geometry.LocalGeometry, + global_geometry::AbstractGlobalGeometry, + ∇z = nothing, +) + coordinates = Geometry.product_coordinates( + horizontal_local_geometry.coordinates, + vertical_local_geometry.coordinates, + ) + J = horizontal_local_geometry.J * vertical_local_geometry.J + WJ = horizontal_local_geometry.WJ * vertical_local_geometry.WJ + ∂x∂ξ = blockmat( + horizontal_local_geometry.∂x∂ξ, + vertical_local_geometry.∂x∂ξ, + ∇z, + ) + return Geometry.LocalGeometry(coordinates, J, WJ, ∂x∂ξ) +end +function product_geometry( + horizontal_local_geometry::Geometry.LocalGeometry, + vertical_local_geometry::Geometry.LocalGeometry, + global_geometry::DeepSphericalGlobalGeometry, + ∇z = nothing, +) + r = global_geometry.radius + z = vertical_local_geometry.coordinates.z + scale = ((r + z) / r) + + coordinates = Geometry.product_coordinates( + horizontal_local_geometry.coordinates, + vertical_local_geometry.coordinates, + ) + J = scale^2 * horizontal_local_geometry.J * vertical_local_geometry.J + WJ = scale^2 * horizontal_local_geometry.WJ * vertical_local_geometry.WJ + ∂x∂ξ = blockmat( + scale * horizontal_local_geometry.∂x∂ξ, + vertical_local_geometry.∂x∂ξ, + ∇z, + ) + return Geometry.LocalGeometry(coordinates, J, WJ, ∂x∂ξ) +end diff --git a/src/Geometry/localgeometry.jl b/src/Geometry/localgeometry.jl index 098ef94f48..31169b27c6 100644 --- a/src/Geometry/localgeometry.jl +++ b/src/Geometry/localgeometry.jl @@ -52,7 +52,34 @@ end undertype(::Type{LocalGeometry{I, C, FT, S}}) where {I, C, FT, S} = FT undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT +""" + blockmat(m11, m22[, m12]) +Construct an `Axis2Tensor` from sub-blocks +""" +function blockmat( + a::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.UAxis, Geometry.Covariant1Axis}, + SMatrix{1, 1, FT, 1}, + }, + b::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, + SMatrix{1, 1, FT, 1}, + }, + c::Nothing = nothing, +) where {FT} + A = Geometry.components(a) + B = Geometry.components(b) + Geometry.AxisTensor( + (Geometry.UWAxis(), Geometry.Covariant13Axis()), + @SMatrix [ + A[1, 1] zero(FT) + zero(FT) B[1, 1] + ] + ) +end function blockmat( a::Geometry.Axis2Tensor{ FT, @@ -64,12 +91,21 @@ function blockmat( Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, SMatrix{1, 1, FT, 1}, }, + c::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant1Axis}, + SMatrix{1, 1, FT, 1}, + }, ) where {FT} A = Geometry.components(a) B = Geometry.components(b) + C = Geometry.components(c) Geometry.AxisTensor( (Geometry.UWAxis(), Geometry.Covariant13Axis()), - SMatrix{2, 2}(A[1, 1], zero(FT), zero(FT), B[1, 1]), + @SMatrix [ + A[1, 1] zero(FT) + C[1, 1] B[1, 1] + ] ) end @@ -84,15 +120,46 @@ function blockmat( Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, SMatrix{1, 1, FT, 1}, }, + c::Nothing = nothing, ) where {FT} A = Geometry.components(a) B = Geometry.components(b) Geometry.AxisTensor( (Geometry.VWAxis(), Geometry.Covariant23Axis()), - SMatrix{2, 2}(A[1, 1], zero(FT), zero(FT), B[1, 1]), + @SMatrix [ + A[1, 1] zero(FT) + zero(FT) B[1, 1] + ] + ) +end +function blockmat( + a::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.VAxis, Geometry.Covariant2Axis}, + SMatrix{1, 1, FT, 1}, + }, + b::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, + SMatrix{1, 1, FT, 1}, + }, + c::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant2Axis}, + SMatrix{1, 1, FT, 1}, + }, +) where {FT} + A = Geometry.components(a) + B = Geometry.components(b) + C = Geometry.components(c) + Geometry.AxisTensor( + (Geometry.VWAxis(), Geometry.Covariant23Axis()), + @SMatrix [ + A[1, 1] zero(FT) + C[1, 1] B[1, 1] + ] ) end - function blockmat( a::Geometry.Axis2Tensor{ FT, @@ -104,36 +171,45 @@ function blockmat( Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, SMatrix{1, 1, FT, 1}, }, + c::Nothing = nothing, ) where {FT} A = Geometry.components(a) B = Geometry.components(b) Geometry.AxisTensor( (Geometry.UVWAxis(), Geometry.Covariant123Axis()), - SMatrix{3, 3}( - A[1, 1], - A[2, 1], - zero(FT), - A[1, 2], - A[2, 2], - zero(FT), - zero(FT), - zero(FT), - B[1, 1], - ), + @SMatrix [ + A[1, 1] A[1, 2] zero(FT) + A[2, 1] A[2, 2] zero(FT) + zero(FT) zero(FT) B[1, 1] + ] ) end - -function product_geometry( - horizontal_local_geometry::Geometry.LocalGeometry, - vertical_local_geometry::Geometry.LocalGeometry, -) - coordinates = Geometry.product_coordinates( - horizontal_local_geometry.coordinates, - vertical_local_geometry.coordinates, +function blockmat( + a::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.UVAxis, Geometry.Covariant12Axis}, + SMatrix{2, 2, FT, 4}, + }, + b::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant3Axis}, + SMatrix{1, 1, FT, 1}, + }, + c::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.WAxis, Geometry.Covariant12Axis}, + SMatrix{1, 2, FT, 2}, + }, +) where {FT} + A = Geometry.components(a) + B = Geometry.components(b) + C = Geometry.components(c) + Geometry.AxisTensor( + (Geometry.UVWAxis(), Geometry.Covariant123Axis()), + @SMatrix [ + A[1, 1] A[1, 2] zero(FT) + A[2, 1] A[2, 2] zero(FT) + C[1, 1] C[1, 2] B[1, 1] + ] ) - J = horizontal_local_geometry.J * vertical_local_geometry.J - WJ = horizontal_local_geometry.WJ * vertical_local_geometry.WJ - ∂x∂ξ = - blockmat(horizontal_local_geometry.∂x∂ξ, vertical_local_geometry.∂x∂ξ) - return Geometry.LocalGeometry(coordinates, J, WJ, ∂x∂ξ) end diff --git a/src/Grids/extruded.jl b/src/Grids/extruded.jl index c7d0cda82a..25ba51c6d2 100644 --- a/src/Grids/extruded.jl +++ b/src/Grids/extruded.jl @@ -37,12 +37,36 @@ mutable struct ExtrudedFiniteDifferenceGrid{ face_local_geometry::LG end -# non-view grids are cached based on their input arguments -# this means that if data is saved in two different files, reloading will give fields which live on the same grid function ExtrudedFiniteDifferenceGrid( horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D}, vertical_grid::FiniteDifferenceGrid, hypsography::HypsographyAdaption = Flat(); + deep = false, +) + if horizontal_grid.global_geometry isa Geometry.SphericalGlobalGeometry + radius = horizontal_grid.global_geometry.radius + if deep + global_geometry = Geometry.DeepSphericalGlobalGeometry(radius) + else + global_geometry = Geometry.ShallowSphericalGlobalGeometry(radius) + end + else + global_geometry = horizontal_grid.global_geometry + end + ExtrudedFiniteDifferenceGrid( + horizontal_grid, + vertical_grid, + hypsography, + global_geometry, + ) +end + +# memoized constructor +function ExtrudedFiniteDifferenceGrid( + horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D}, + vertical_grid::FiniteDifferenceGrid, + hypsography::HypsographyAdaption, + global_geometry::Geometry.AbstractGlobalGeometry, ) get!( Cache.OBJECT_CACHE, @@ -51,32 +75,36 @@ function ExtrudedFiniteDifferenceGrid( horizontal_grid, vertical_grid, hypsography, + global_geometry, ), ) do _ExtrudedFiniteDifferenceGrid( horizontal_grid, vertical_grid, hypsography, + global_geometry, ) end end -# Non-memoized constructor. Should not generally be called. +# Non-memoized constructor. Should not generally be called, but can be defined for other Hypsography types function _ExtrudedFiniteDifferenceGrid( horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D}, vertical_grid::FiniteDifferenceGrid, hypsography::Flat, + global_geometry::Geometry.AbstractGlobalGeometry, ) - global_geometry = horizontal_grid.global_geometry center_local_geometry = Geometry.product_geometry.( horizontal_grid.local_geometry, vertical_grid.center_local_geometry, + Ref(global_geometry), ) face_local_geometry = Geometry.product_geometry.( horizontal_grid.local_geometry, vertical_grid.face_local_geometry, + Ref(global_geometry), ) return ExtrudedFiniteDifferenceGrid( diff --git a/src/Grids/finitedifference.jl b/src/Grids/finitedifference.jl index 9e1b20fc3e..9daab70f5a 100644 --- a/src/Grids/finitedifference.jl +++ b/src/Grids/finitedifference.jl @@ -49,84 +49,22 @@ end function _FiniteDifferenceGrid(topology::Topologies.IntervalTopology) global_geometry = Geometry.CartesianGlobalGeometry() - mesh = topology.mesh - CT = Meshes.coordinate_type(mesh) - AIdx = Geometry.coordinate_axis(CT) - # TODO: FD operators hardcoded to work over the 3-axis, need to generalize - # similar to spectral operators - @assert AIdx == (3,) "FiniteDifference operations only work over the 3-axis (ZPoint) domain" - FT = eltype(CT) ArrayType = ClimaComms.array_type(topology) - face_coordinates = collect(mesh.faces) - LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}} - nface = length(face_coordinates) - Topologies.isperiodic(topology) - ncent = length(face_coordinates) - 1 - # contstruct on CPU, copy to device at end - center_local_geometry = DataLayouts.VF{LG}(Array{FT}, ncent) - face_local_geometry = DataLayouts.VF{LG}(Array{FT}, nface) - for i in 1:ncent - # centers - coord⁻ = Geometry.component(face_coordinates[i], 1) - coord⁺ = Geometry.component(face_coordinates[i + 1], 1) - # at the moment we use a "discrete Jacobian" - # ideally we should use the continuous quantity via the derivative of the warp function - # could we just define this then as deriv on the mesh element coordinates? - coord = (coord⁺ + coord⁻) / 2 - Δcoord = coord⁺ - coord⁻ - J = Δcoord - WJ = Δcoord - ∂x∂ξ = SMatrix{1, 1}(J) - center_local_geometry[i] = Geometry.LocalGeometry( - CT(coord), - J, - WJ, - Geometry.AxisTensor( - (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), - ∂x∂ξ, - ), - ) - end - for i in 1:nface - coord = Geometry.component(face_coordinates[i], 1) - if i == 1 - # bottom face - if Topologies.isperiodic(topology) - Δcoord⁺ = - Geometry.component(face_coordinates[2], 1) - - Geometry.component(face_coordinates[1], 1) - Δcoord⁻ = - Geometry.component(face_coordinates[end], 1) - - Geometry.component(face_coordinates[end - 1], 1) - J = (Δcoord⁺ + Δcoord⁻) / 2 - WJ = J - else - coord⁺ = Geometry.component(face_coordinates[2], 1) - J = coord⁺ - coord - WJ = J / 2 - end - elseif !Topologies.isperiodic(topology) && i == nface - # top face - coord⁻ = Geometry.component(face_coordinates[i - 1], 1) - J = coord - coord⁻ - WJ = J / 2 - else - coord⁺ = Geometry.component(face_coordinates[i + 1], 1) - coord⁻ = Geometry.component(face_coordinates[i - 1], 1) - J = (coord⁺ - coord⁻) / 2 - WJ = J - end - ∂x∂ξ = SMatrix{1, 1}(J) - ∂ξ∂x = SMatrix{1, 1}(inv(J)) - face_local_geometry[i] = Geometry.LocalGeometry( - CT(coord), - J, - WJ, - Geometry.AxisTensor( - (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), - ∂x∂ξ, - ), - ) + + mesh = Topologies.mesh(topology) + CT = Topologies.coordinate_type(mesh) + FT = Geometry.float_type(CT) + Nv_face = length(mesh.faces) + # construct on CPU, adapt to GPU + face_coordinates = DataLayouts.VF{CT}(Array{FT}, Nv_face) + for v in 1:Nv_face + face_coordinates[v] = mesh.faces[v] end + center_local_geometry, face_local_geometry = fd_geometry_data( + face_coordinates; + periodic = Topologies.isperiodic(topology), + ) + return FiniteDifferenceGrid( topology, global_geometry, @@ -135,6 +73,136 @@ function _FiniteDifferenceGrid(topology::Topologies.IntervalTopology) ) end +# called by the FiniteDifferenceGrid constructor, and the ExtrudedFiniteDifferenceGrid constructor with Hypsography +function fd_geometry_data( + face_coordinates::DataLayouts.AbstractData{Geometry.ZPoint{FT}}; + periodic, +) where {FT} + CT = Geometry.ZPoint{FT} + AIdx = (3,) + LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}} + (Ni, Nj, Nk, Nv, Nh) = size(face_coordinates) + Nv_face = Nv - periodic + Nv_cent = Nv - 1 + + center_local_geometry = + similar(face_coordinates, LG, (Ni, Nj, Nk, Nv_cent, Nh)) + face_local_geometry = + similar(face_coordinates, LG, (Ni, Nj, Nk, Nv_face, Nh)) + + 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⁻ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, v, h)], + 1, + ) + coord⁺ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, v + 1, h)], + 1, + ) + # use a "discrete Jacobian" + coord = (coord⁺ + coord⁻) / 2 + Δcoord = coord⁺ - coord⁻ + J = Δcoord + WJ = Δcoord + ∂x∂ξ = 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∂ξ, + ), + ) + end + + for v in 1:Nv_face + coord = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, v, h)], + 1, + ) + if v == 1 + # bottom face + if periodic + Δcoord⁺ = + Geometry.component( + face_coordinates[CartesianIndex(i, j, k, 2, h)], + 1, + ) - Geometry.component( + face_coordinates[CartesianIndex(i, j, k, 1, h)], + 1, + ) + Δcoord⁻ = + Geometry.component( + face_coordinates[CartesianIndex(i, j, k, Nv, h)], + 1, + ) - Geometry.component( + face_coordinates[CartesianIndex( + i, + j, + k, + Nv - 1, + h, + )], + 1, + ) + J = (Δcoord⁺ + Δcoord⁻) / 2 + WJ = J + else + coord⁺ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, 2, h)], + 1, + ) + J = coord⁺ - coord + WJ = J / 2 + end + elseif v == Nv_cent + 1 + @assert !periodic + # top face + coord⁻ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, Nv - 1, h)], + 1, + ) + J = coord - coord⁻ + WJ = J / 2 + else + coord⁺ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, v + 1, h)], + 1, + ) + coord⁻ = Geometry.component( + face_coordinates[CartesianIndex(i, j, k, v - 1, h)], + 1, + ) + J = (coord⁺ - coord⁻) / 2 + WJ = J + end + ∂x∂ξ = 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∂ξ, + ), + ) + end + end + return (center_local_geometry, face_local_geometry) +end + FiniteDifferenceGrid(mesh::Meshes.IntervalMesh) = FiniteDifferenceGrid(Topologies.IntervalTopology(mesh)) diff --git a/src/Hypsography/Hypsography.jl b/src/Hypsography/Hypsography.jl index ef34ddd996..c156507bb0 100644 --- a/src/Hypsography/Hypsography.jl +++ b/src/Hypsography/Hypsography.jl @@ -1,9 +1,18 @@ module Hypsography +import ClimaComms, Adapt + import ..slab, ..column import ..Geometry, - ..Domains, ..Topologies, ..Grids, ..Spaces, ..Fields, ..Operators + ..DataLayouts, + ..Domains, + ..Topologies, + ..Grids, + ..Spaces, + ..Fields, + ..Operators import ..Spaces: ExtrudedFiniteDifferenceSpace +import ClimaCore.Utilities: half import ..Grids: _ExtrudedFiniteDifferenceGrid, @@ -15,28 +24,52 @@ using StaticArrays, LinearAlgebra """ - ref_z_to_physical_z(adaption::HypsographyAdaption, z_ref, z_surface, z_top) + ref_z_to_physical_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint Convert reference `z`s to physical `z`s as prescribed by the given adaption. + +This function has to be the inverse of `physical_z_to_ref_z`. """ function ref_z_to_physical_z( adaption::HypsographyAdaption, - z_ref, - z_surface, - z_top, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, ) end -# For no adaption, z_physical is z_ref. -function ref_z_to_physical_z(adaption::Nothing, z_ref, z_surface, z_top) - return z_ref -end - """ - _check_adaption_constraints(adaption, z_ref, z_surface, z_top) + physical_z_to_ref_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint -Ensures that the `adaption` is well defined for the given parameters. +Convert physical `z`s to reference `z`s as prescribed by the given adaption. + +This function has to be the inverse of `ref_z_to_physical_z`. """ -function _check_adaption_constraints(adaption, z_ref, z_surface, z_top) end +function physical_z_to_ref_z( + adaption::HypsographyAdaption, + z_phys::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, +) end + +# Flat, z_ref = z_physical + +function ref_z_to_physical_z( + ::Flat, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, +) + return z_ref +end + +function physical_z_to_ref_z( + ::Flat, + z_physical::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, +) + return z_physical +end """ LinearAdaption(surface::Field) @@ -44,13 +77,35 @@ function _check_adaption_constraints(adaption, z_ref, z_surface, z_top) end Locate the levels by linear interpolation between the surface field and the top of the domain, using the method of [GalChen1975](@cite). """ -struct LinearAdaption{F <: Union{Fields.Field, Nothing}} <: HypsographyAdaption - # Union can be removed once deprecation removed. +struct LinearAdaption{F <: Fields.Field} <: HypsographyAdaption surface::F + function LinearAdaption(surface::Fields.Field) + if eltype(surface) <: Real + @warn "`LinearAdaptation`: `surface` argument scalar field has been deprecated. Use a field `ZPoint`s." + surface = Geometry.ZPoint.(surface) + end + new{typeof(surface)}(surface) + end end -function ref_z_to_physical_z(adaption::LinearAdaption, z_ref, z_surface, z_top) - return z_ref + (1 - z_ref / z_top) * z_surface +# This method is invoked by the ExtrudedFiniteDifferenceGrid constructor +function ref_z_to_physical_z( + ::LinearAdaption, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, +) + Geometry.ZPoint(z_ref.z + (1 - z_ref.z / z_top.z) * z_surface.z) +end + +# This method is used for remapping +function physical_z_to_ref_z( + ::LinearAdaption, + z_physical::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, +) + Geometry.ZPoint((z_physical.z - z_surface.z) / (1 - z_surface.z / z_top.z)) end """ @@ -64,180 +119,154 @@ If the decay-scale is poorly specified (i.e., `s * zₜ` is lower than the maxim surface elevation), a warning is thrown and `s` is adjusted such that it `szₜ > maximum(z_surface)`. """ struct SLEVEAdaption{F <: Fields.Field, FT <: Real} <: HypsographyAdaption - # Union can be removed once deprecation removed. surface::F ηₕ::FT s::FT + function SLEVEAdaption( + surface::Fields.Field, + ηₕ::FT, + s::FT, + ) where {FT <: Real} + @assert 0 <= ηₕ <= 1 + @assert s >= 0 + if eltype(surface) <: Real + @warn "`SLEVEAdaption`: `surface` argument scalar field has been deprecated. Use a field `ZPoint`s." + surface = Geometry.ZPoint.(surface) + end + new{typeof(surface), FT}(surface, ηₕ, s) + end end -function _check_adaption_constraints( + +function ref_z_to_physical_z( adaption::SLEVEAdaption, - z_ref, - z_surface, - z_top, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, ) - ηₕ = adaption.ηₕ - s = adaption.s - FT = eltype(s) - @assert FT(0) <= ηₕ <= FT(1) - @assert s >= FT(0) - if s * z_top <= maximum(z_surface) - @warn "Decay scale (s*z_top = $(s*z_top)) must be higher than max surface elevation (max(z_surface) = $(maximum(z_surface))). Returning s = FT(0.8). Scale height is therefore s=$(0.8 * z_top) m." + (; ηₕ, s) = adaption + if s * z_top.z <= z_surface.z + error("Decay scale (s*z_top) must be higher than max surface elevation") end -end -function ref_z_to_physical_z(adaption::SLEVEAdaption, z_ref, z_surface, z_top) - ηₕ = adaption.ηₕ - s = adaption.s - η = z_ref ./ z_top - FT = eltype(s) - # Fix scale height if needed - s = ifelse(s * z_top <= maximum(z_surface), FT(8 / 10), s) - return ifelse( - η <= ηₕ, - η * z_top + z_surface * (sinh((ηₕ - η) / s / ηₕ)) / (sinh(1 / s)), - η * z_top, - ) + η = z_ref.z / z_top.z + if η <= ηₕ + return Geometry.ZPoint( + η * z_top.z + + z_surface.z * (sinh((ηₕ - η) / s / ηₕ)) / (sinh(1 / s)), + ) + else + return Geometry.ZPoint(η * z_top.z) + end end -# deprecated in 0.10.12 -LinearAdaption() = LinearAdaption(nothing) - -@deprecate( - ExtrudedFiniteDifferenceSpace( - horizontal_space::Spaces.AbstractSpace, - vertical_space::Spaces.FiniteDifferenceSpace, - adaption::LinearAdaption{Nothing}, - z_surface::Fields.Field, - ), - ExtrudedFiniteDifferenceSpace( - horizontal_space, - vertical_space, - LinearAdaption(z_surface), - ), - false +function physical_z_to_ref_z( + adaption::SLEVEAdaption, + z_physical::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, ) + error("This method is not implemented") +end -# linear coordinates +# can redefine this constructor for e.g. multi-arg SLEVE function _ExtrudedFiniteDifferenceGrid( horizontal_grid::Grids.AbstractGrid, vertical_grid::Grids.FiniteDifferenceGrid, adaption::HypsographyAdaption, + global_geometry::Geometry.AbstractGlobalGeometry, ) - if adaption isa LinearAdaption - if isnothing(adaption.surface) - error("LinearAdaption requires a Field argument") - end - if axes(adaption.surface).grid !== horizontal_grid - error("Terrain must be defined on the horizontal space") - end - end + @assert Spaces.grid(axes(adaption.surface)) == horizontal_grid + z_surface = Fields.field_values(adaption.surface) + + face_z_ref = + Grids.local_geometry_data(vertical_grid, Grids.CellFace()).coordinates + vertical_domain = Topologies.domain(vertical_grid) + z_top = vertical_domain.coord_max - # construct initial flat space, then mutate - ref_grid = Grids.ExtrudedFiniteDifferenceGrid( + face_z = + ref_z_to_physical_z.(Ref(adaption), face_z_ref, z_surface, Ref(z_top)) + + return _ExtrudedFiniteDifferenceGrid( horizontal_grid, vertical_grid, - Flat(), + adaption, + global_geometry, + face_z, ) - face_ref_space = Spaces.FaceExtrudedFiniteDifferenceSpace(ref_grid) - face_ref_coords = Spaces.coordinates_data(face_ref_space) - coord_type = eltype(face_ref_coords) - z_ref = face_ref_coords.z - - vertical_domain = Topologies.domain(vertical_grid) - z_top = vertical_domain.coord_max.z - - grad = Operators.Gradient() - z_surface = Fields.field_values(adaption.surface) +end - FT = eltype(z_surface) +# generic 5-arg hypsography constructor, uses computed face_z points +function _ExtrudedFiniteDifferenceGrid( + horizontal_grid::Grids.AbstractGrid, + vertical_grid::Grids.FiniteDifferenceGrid, + adaption::HypsographyAdaption, + global_geometry::Geometry.AbstractGlobalGeometry, + face_z::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, +) where {FT} + # construct the "flat" grid + # avoid cached constructor so that it gets cleaned up automatically + flat_grid = _ExtrudedFiniteDifferenceGrid( + horizontal_grid, + vertical_grid, + Flat(), + global_geometry, + ) + center_flat_space = Spaces.space(flat_grid, Grids.CellCenter()) + face_flat_space = Spaces.space(flat_grid, Grids.CellFace()) + + # compute the "z-only local geometry" based on face z coords + 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, face_z); + periodic = Topologies.isperiodic(vertical_grid.topology), + ) - _check_adaption_constraints(adaption, z_ref, z_surface, z_top) - fZ_data = ref_z_to_physical_z.(Ref(adaption), z_ref, z_surface, z_top) - fZ = Fields.Field(fZ_data, face_ref_space) + center_z_local_geometry = Adapt.adapt(ArrayType, center_z_local_geometry) + face_z_local_geometry = Adapt.adapt(ArrayType, face_z_local_geometry) - # Take the horizontal gradient for the Z surface field - # for computing updated ∂x∂ξ₃₁, ∂x∂ξ₃₂ terms + # compute ∇Z at face and centers grad = Operators.Gradient() - If2c = Operators.InterpolateF2C() - - # DSS the horizontal gradient of Z surface field to force - # deriv continuity along horizontal element boundaries - f∇Z = grad.(fZ) - Spaces.weighted_dss!(f∇Z) - - # Interpolate horizontal gradient surface field to centers - # used to compute ∂x∂ξ₃₃ (Δz) metric term - cZ = If2c.(fZ) - - # DSS the interpolated horizontal gradients as well - c∇Z = If2c.(f∇Z) - Spaces.weighted_dss!(c∇Z) - - face_local_geometry = copy(ref_grid.face_local_geometry) - center_local_geometry = copy(ref_grid.center_local_geometry) - - Ni, Nj, _, Nv, Nh = size(center_local_geometry) - for h in 1:Nh, j in 1:Nj, i in 1:Ni - face_column = column(face_local_geometry, i, j, h) - fZ_column = column(Fields.field_values(fZ), i, j, h) - f∇Z_column = column(Fields.field_values(f∇Z), i, j, h) - center_column = column(center_local_geometry, i, j, h) - cZ_column = column(Fields.field_values(cZ), i, j, h) - c∇Z_column = column(Fields.field_values(c∇Z), i, j, h) - - # update face metrics - for v in 1:(Nv + 1) - local_geom = face_column[v] - coord = if coord_type <: Geometry.Abstract2DPoint - c1 = Geometry.components(local_geom.coordinates)[1] - coord_type(c1, fZ_column[v]) - elseif coord_type <: Geometry.Abstract3DPoint - c1 = Geometry.components(local_geom.coordinates)[1] - c2 = Geometry.components(local_geom.coordinates)[2] - coord_type(c1, c2, fZ_column[v]) - end - Δz = if v == 1 - # if this is the domain min face level compute the metric - # extrapolating from the bottom face level of the domain - 2 * (cZ_column[v] - fZ_column[v]) - elseif v == Nv + 1 - # if this is the domain max face level compute the metric - # extrapolating from the top face level of the domain - 2 * (fZ_column[v] - cZ_column[v - 1]) - else - cZ_column[v] - cZ_column[v - 1] - end - ∂x∂ξ = reconstruct_metric(local_geom.∂x∂ξ, f∇Z_column[v], Δz) - W = local_geom.WJ / local_geom.J - J = det(Geometry.components(∂x∂ξ)) - face_column[v] = Geometry.LocalGeometry(coord, J, W * J, ∂x∂ξ) - end - # update center metrics - for v in 1:Nv - local_geom = center_column[v] - coord = if coord_type <: Geometry.Abstract2DPoint - c1 = Geometry.components(local_geom.coordinates)[1] - coord_type(c1, cZ_column[v]) - elseif coord_type <: Geometry.Abstract3DPoint - c1 = Geometry.components(local_geom.coordinates)[1] - c2 = Geometry.components(local_geom.coordinates)[2] - coord_type(c1, c2, cZ_column[v]) - end - Δz = fZ_column[v + 1] - fZ_column[v] - ∂x∂ξ = reconstruct_metric(local_geom.∂x∂ξ, c∇Z_column[v], Δz) - W = local_geom.WJ / local_geom.J - J = det(Geometry.components(∂x∂ξ)) - center_column[v] = Geometry.LocalGeometry(coord, J, W * J, ∂x∂ξ) - end - end + center_∇Z_field = + grad.( + Fields.Field( + center_z_local_geometry, + center_flat_space, + ).coordinates.z + ) + Spaces.weighted_dss!(center_∇Z_field) + + face_∇Z_field = + grad.( + Fields.Field(face_z_local_geometry, face_flat_space).coordinates.z + ) + Spaces.weighted_dss!(face_∇Z_field) + + # construct full local geometry + center_local_geometry = + Geometry.product_geometry.( + horizontal_grid.local_geometry, + center_z_local_geometry, + Ref(global_geometry), + Ref(Geometry.WVector(1)) .* + adjoint.(Fields.field_values(center_∇Z_field)), + ) + face_local_geometry = + Geometry.product_geometry.( + horizontal_grid.local_geometry, + face_z_local_geometry, + Ref(global_geometry), + Ref(Geometry.WVector(1)) .* + adjoint.(Fields.field_values(face_∇Z_field)), + ) return ExtrudedFiniteDifferenceGrid( horizontal_grid, vertical_grid, adaption, - ref_grid.global_geometry, + global_geometry, center_local_geometry, face_local_geometry, ) @@ -265,57 +294,25 @@ function diffuse_surface_elevation!( maxiter::Int = 100, dt::T = 1e-1, ) where {T} + if eltype(f) <: Real + f_z = f + elseif eltype(f) <: Geometry.ZPoint + f_z = f.z + end # Define required ops wdiv = Operators.WeakDivergence() grad = Operators.Gradient() - FT = eltype(f) # Create dss buffer - ghost_buffer = (bf = Spaces.create_dss_buffer(f),) + ghost_buffer = (bf = Spaces.create_dss_buffer(f_z),) # Apply smoothing for iter in 1:maxiter # Euler steps - χf = @. wdiv(grad(f)) + χf = @. wdiv(grad(f_z)) Spaces.weighted_dss!(χf, ghost_buffer.bf) - @. f += κ * dt * χf + @. f_z += κ * dt * χf end # Return mutated surface elevation profile return f end -function reconstruct_metric( - ∂x∂ξ::Geometry.Axis2Tensor{ - T, - Tuple{Geometry.UWAxis, Geometry.Covariant13Axis}, - }, - ∇z::Geometry.Covariant1Vector, - Δz::Real, -) where {T} - v∂x∂ξ = Geometry.components(∂x∂ξ) - v∇z = Geometry.components(∇z) - Geometry.AxisTensor(axes(∂x∂ξ), @SMatrix [ - v∂x∂ξ[1, 1] 0 - v∇z[1] Δz - ]) -end - -function reconstruct_metric( - ∂x∂ξ::Geometry.Axis2Tensor{ - T, - Tuple{Geometry.UVWAxis, Geometry.Covariant123Axis}, - }, - ∇z::Geometry.Covariant12Vector, - Δz::Real, -) where {T} - v∂x∂ξ = Geometry.components(∂x∂ξ) - v∇z = Geometry.components(∇z) - Geometry.AxisTensor( - axes(∂x∂ξ), - @SMatrix [ - v∂x∂ξ[1, 1] v∂x∂ξ[1, 2] 0 - v∂x∂ξ[2, 1] v∂x∂ξ[2, 2] 0 - v∇z[1] v∇z[2] Δz - ] - ) -end - end diff --git a/src/InputOutput/readers.jl b/src/InputOutput/readers.jl index ab0d9b866e..cb810ebf73 100644 --- a/src/InputOutput/readers.jl +++ b/src/InputOutput/readers.jl @@ -346,6 +346,13 @@ function read_grid_new(reader, name) hypsography = Hypsography.LinearAdaption( read_field(reader, attrs(group)["hypsography_surface"]), ) + elseif hypsography_type == "SLEVEAdaption" + # Store hyps object for general use ? + hypsography = Hypsography.SLEVEAdaption( + read_field(reader, attrs(group)["hypsography_surface"]), + get(attrs(group), "hypsography_ηₕ", "hypsography_surface_ηₕ"), + get(attrs(group), "hypsography_s", "hypsography_surface_s"), + ) else error("Unsupported hypsography type $hypsography_type") end diff --git a/src/InputOutput/writers.jl b/src/InputOutput/writers.jl index 64a3207104..f888f74082 100644 --- a/src/InputOutput/writers.jl +++ b/src/InputOutput/writers.jl @@ -363,6 +363,12 @@ function write_new!( return name end + +""" + write_new!(writer, domain, name) + +Writes an object of type 'Hypsography' and name 'name' to the HDF5 file. +""" function write_new!( writer::HDF5Writer, space::Grids.ExtrudedFiniteDifferenceGrid, @@ -378,6 +384,15 @@ function write_new!( write_attribute(group, "vertical_grid", write!(writer, space.vertical_grid)) if space.hypsography isa Hypsography.LinearAdaption write_attribute(group, "hypsography_type", "LinearAdaption") + write_attribute( + group, + "hypsography_surface", + write!(writer, space.hypsography.surface, "_z_surface/$name"), # Change to save "space.hyps" + ) + elseif space.hypsography isa Hypsography.SLEVEAdaption + write_attribute(group, "hypsography_type", "SLEVEAdaption") + write_attribute(group, "hypsography_ηₕ", space.hypsography.ηₕ) + write_attribute(group, "hypsography_s", space.hypsography.s) write_attribute( group, "hypsography_surface", @@ -469,6 +484,7 @@ function write!( return name end + """ write!(writer::HDF5Writer, name => value...) diff --git a/src/Remapping/distributed_remapping.jl b/src/Remapping/distributed_remapping.jl index 76a5cb867b..14a7358eb4 100644 --- a/src/Remapping/distributed_remapping.jl +++ b/src/Remapping/distributed_remapping.jl @@ -235,7 +235,7 @@ function interpolate( axes(field) == remapper.space || error("Field is defined on a different space than remapper") - FT = eltype(field) + FT = Spaces.undertype(axes(field)) if length(remapper.target_zcoords) == 0 out_local_array = zeros(FT, size(remapper.local_target_hcoords_bitmask)) diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl index 725de231ab..4cf2342268 100644 --- a/src/Remapping/interpolate_array.jl +++ b/src/Remapping/interpolate_array.jl @@ -590,10 +590,6 @@ function interpolate_column( # we have to make sure that we are passing around views of the same array (as opposed to # copied of it). if physical_z - # We are hardcoding the transformation from Hypsography.LinearAdaption - space.hypsography isa Hypsography.LinearAdaption || - error("Cannot interpolate $(space.hypsography) hypsography") - FT = Spaces.undertype(axes(field)) # interpolate_slab! takes a vector @@ -601,7 +597,7 @@ function interpolate_column( interpolate_slab!( z_surface, - space.hypsography.surface, + space.hypsography.surface.z, [Fields.SlabIndex(nothing, gidx)], [Is], ) @@ -609,11 +605,15 @@ function interpolate_column( z_top = Spaces.vertical_topology(space).mesh.domain.coord_max.z zpts_ref = [ - Geometry.ZPoint((z.z - z_surface) / (1 - z_surface / z_top)) for - z in zpts if z.z > z_surface + Hypsography.physical_z_to_ref_z( + Spaces.grid(space).hypsography, + z_ref, + Geometry.ZPoint(z_surface), + Geometry.ZPoint(z_top), + ) for z_ref in zpts ] - # When zpts = zpts_ref, all the points are above the surface + # Edge case check: when zpts = zpts_ref, all the points are above the surface num_points_below_surface = length(zpts) - length(zpts_ref) fill!( diff --git a/test/Geometry/geometry.jl b/test/Geometry/geometry.jl index 5a6f1457f2..97724b998f 100644 --- a/test/Geometry/geometry.jl +++ b/test/Geometry/geometry.jl @@ -589,14 +589,44 @@ end global_geom, ) ≈ 2 * deg2rad(180.0) rtol = 2eps() +end + + +@testset "shallow great circle distance" begin + global_geom = Geometry.ShallowSphericalGlobalGeometry(2.0) + # test between two LatLongZPoints @test Geometry.great_circle_distance( - Geometry.LatLongZPoint(22.0, 32.0, 2.0), - Geometry.LatLongZPoint(22.0, 32.0, 2.0), + Geometry.LatLongZPoint(0.0, 30.0, 0.0), + Geometry.LatLongZPoint(0.0, 40.0, 0.0), global_geom, - ) == 0.0 + ) ≈ 2 * deg2rad(10.0) rtol = 2eps() + + @test Geometry.great_circle_distance( + Geometry.LatLongZPoint(0.0, 30.0, 10.0), + Geometry.LatLongZPoint(0.0, 40.0, 10.0), + global_geom, + ) ≈ 2 * deg2rad(10.0) rtol = 2eps() + end +@testset "deep great circle distance" begin + global_geom = Geometry.DeepSphericalGlobalGeometry(2.0) + + # test between two LatLongZPoints + @test Geometry.great_circle_distance( + Geometry.LatLongZPoint(0.0, 30.0, 0.0), + Geometry.LatLongZPoint(0.0, 40.0, 0.0), + global_geom, + ) ≈ 2 * deg2rad(10.0) rtol = 2eps() + + @test Geometry.great_circle_distance( + Geometry.LatLongZPoint(0.0, 30.0, 10.0), + Geometry.LatLongZPoint(0.0, 40.0, 10.0), + global_geom, + ) ≈ (10 + 2) * deg2rad(10.0) rtol = 2eps() + +end @testset "1D XPoint Euclidean distance" begin for FT in (Float32, Float64, BigFloat) pt_1 = Geometry.XPoint(one(FT)) diff --git a/test/Hypsography/2d.jl b/test/Hypsography/2d.jl index 30ee941228..6e178cb192 100644 --- a/test/Hypsography/2d.jl +++ b/test/Hypsography/2d.jl @@ -9,6 +9,7 @@ import ClimaCore: Topologies, Spaces, Quadratures, + Grids, Fields, Operators, InputOutput, @@ -65,3 +66,32 @@ face_x = Fields.coordinate_field(hv_face_space).x @test map(u -> u.u, ∇x) ≈ ones(hv_center_space) @test maximum(u -> abs(u.w), ∇x) < 1e-10 + +# Check transformations back and forth +z_ref1 = Geometry.ZPoint{FT}(40.0) +z_top1 = Geometry.ZPoint{FT}(400.0) +z_surface1 = Geometry.ZPoint{FT}(10.0) + +point_space = Spaces.PointSpace(z_surface1) +z_surface_point = Fields.coordinate_field(point_space) + +@test Hypsography.physical_z_to_ref_z( + Grids.Flat(), + Hypsography.ref_z_to_physical_z(Grids.Flat(), z_ref1, z_surface1, z_top1), + z_surface1, + z_top1, +).z ≈ z_ref1.z + +# Linear adaption +linear_adaption = Hypsography.LinearAdaption(z_surface_point) +@test Hypsography.physical_z_to_ref_z( + linear_adaption, + Hypsography.ref_z_to_physical_z( + linear_adaption, + z_ref1, + z_surface1, + z_top1, + ), + z_surface1, + z_top1, +).z ≈ z_ref1.z diff --git a/test/Operators/hybrid/3d.jl b/test/Operators/hybrid/3d.jl index d802316951..bf8bc512d0 100644 --- a/test/Operators/hybrid/3d.jl +++ b/test/Operators/hybrid/3d.jl @@ -9,6 +9,7 @@ import ClimaCore: Meshes, Geometry, Topologies, + Grids, Spaces, Quadratures, Fields, @@ -32,24 +33,53 @@ device = ClimaComms.device() ClimaComms.SingletonCommsContext(device), vertmesh, ) - vert_center_space = Spaces.CenterFiniteDifferenceSpace(verttopology) + vertgrid = Grids.FiniteDifferenceGrid(verttopology) - horzdomain = Domains.SphereDomain(30.0) + radius = 30.0 + horzdomain = Domains.SphereDomain(radius) horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 4) horztopology = Topologies.Topology2D( ClimaComms.SingletonCommsContext(device), horzmesh, ) quad = Quadratures.GLL{3 + 1}() - horzspace = Spaces.SpectralElementSpace2D(horztopology, quad) + horzgrid = Grids.SpectralElementGrid2D(horztopology, quad) - hv_center_space = - Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space) + # shallow + shallow_grid = Grids.ExtrudedFiniteDifferenceGrid(horzgrid, vertgrid) + + coords = Fields.coordinate_field( + Spaces.CenterExtrudedFiniteDifferenceSpace(shallow_grid), + ) + x = Geometry.UVWVector.(cosd.(coords.lat), 0.0, 0.0) + div = Operators.Divergence() + @test div.(x) ≈ zero(coords.z) atol = 1e-4 + + fcoords = Fields.coordinate_field( + Spaces.FaceExtrudedFiniteDifferenceSpace(shallow_grid), + ) + y = map(coord -> Geometry.WVector(0.7), fcoords) + divf2c = Operators.DivergenceF2C() + @test divf2c.(y) ≈ zero(coords.z) atol = 100 * eps(FT) + + # deep + deep_grid = + Grids.ExtrudedFiniteDifferenceGrid(horzgrid, vertgrid; deep = true) - coords = Fields.coordinate_field(hv_center_space) + coords = Fields.coordinate_field( + Spaces.CenterExtrudedFiniteDifferenceSpace(deep_grid), + ) x = Geometry.UVWVector.(cosd.(coords.lat), 0.0, 0.0) div = Operators.Divergence() - @test norm(div.(x)) < 2e-2 + @test div.(x) ≈ zero(coords.z) atol = 1e-4 + + # divergence of a constant outward vector field = 2 w / (r + z) + fcoords = Fields.coordinate_field( + Spaces.FaceExtrudedFiniteDifferenceSpace(deep_grid), + ) + y = map(coord -> Geometry.WVector(0.7), fcoords) + divf2c = Operators.DivergenceF2C() + @test divf2c.(y) ≈ (2 * 0.7) ./ (radius .+ coords.z) atol = 100 * eps(FT) end function hvspace_3D( diff --git a/test/Spaces/sphere.jl b/test/Spaces/sphere.jl index 34331baf4e..a2dedb922a 100644 --- a/test/Spaces/sphere.jl +++ b/test/Spaces/sphere.jl @@ -1,7 +1,7 @@ using LinearAlgebra, IntervalSets using ClimaComms import ClimaCore: - Domains, Topologies, Meshes, Spaces, Geometry, column, Quadratures + Domains, Topologies, Meshes, Grids, Spaces, Geometry, column, Quadratures using Test @@ -87,7 +87,7 @@ end @testset "Volume of a spherical shell" begin FT = Float64 context = ClimaComms.SingletonCommsContext() - radius = FT(128) + radius = FT(10) zlim = (0, 1) helem = 4 zelem = 10 @@ -99,21 +99,31 @@ end boundary_tags = (:bottom, :top), ) vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem) - vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh) + vertgrid = Grids.FiniteDifferenceGrid(vertmesh) horzdomain = Domains.SphereDomain(radius) horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem) horztopology = Topologies.Topology2D(context, horzmesh) quad = Quadratures.GLL{Nq}() - horzspace = Spaces.SpectralElementSpace2D(horztopology, quad) - + 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)) - hv_center_space = - Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space) - # "shallow atmosphere" spherical shell: volume = surface area * height - @test sum(ones(hv_center_space)) ≈ 4pi * radius^2 * (zlim[2] - zlim[1]) rtol = - 1e-3 + 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 + end diff --git a/test/Spaces/terrain_warp.jl b/test/Spaces/terrain_warp.jl index 21e46205d4..0f5fb52cb1 100644 --- a/test/Spaces/terrain_warp.jl +++ b/test/Spaces/terrain_warp.jl @@ -90,7 +90,7 @@ function generate_smoothed_orography( test_smoothing::Bool = false, ) # Extrusion - z_surface = warp_fn.(Fields.coordinate_field(hspace)) + z_surface = Geometry.ZPoint.(warp_fn.(Fields.coordinate_field(hspace))) # An Euler step defines the diffusion coefficient # (See e.g. cfl condition for diffusive terms). x_array = parent(Fields.coordinate_field(hspace).x) @@ -111,11 +111,8 @@ function get_adaptation(adaption, z_surface::Fields.Field) if adaption <: Hypsography.LinearAdaption return adaption(z_surface) elseif adaption <: Hypsography.SLEVEAdaption - return adaption( - z_surface, - eltype(z_surface)(0.75), - eltype(z_surface)(0.60), - ) + FT = eltype(eltype(z_surface)) + return adaption(z_surface, FT(0.75), FT(0.60)) end end @@ -499,16 +496,10 @@ end z_surface = warp_sin_2d.(Fields.coordinate_field(hspace)) # Generate space with known mesh-warp parameters ηₕ = 1; s = 0.1 # Scale height is poorly specified, so code should throw warning. - @test_logs ( - :warn, - "Decay scale (s*z_top = 0.1) must be higher than max surface elevation (max(z_surface) = 0.5). Returning s = FT(0.8). Scale height is therefore s=0.8 m.", - ) - ( - fspace = Spaces.ExtrudedFiniteDifferenceSpace( - hspace, - vert_face_space, - Hypsography.SLEVEAdaption(z_surface, FT(1), FT(0.1)), - ) + @test_throws ErrorException Spaces.ExtrudedFiniteDifferenceSpace( + hspace, + vert_face_space, + Hypsography.SLEVEAdaption(z_surface, FT(1), FT(0.1)), ) end end