Add option for deep atmosphere

Co-authored-by: Akshay Sridhar <asridhar@caltech.edu>
Co-authored-by: Gabriele Bozzola <gbozzola@caltech.edu>

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)

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)

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
 

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)

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,17 +259,60 @@ end
 
 function CartesianVector(
     u::UVWVector,
-    geom::SphericalGlobalGeometry,
+    geom::AbstractSphericalGlobalGeometry,
     coord::Union{LatLongPoint, LatLongZPoint},
 )
     G = local_to_cartesian(geom, coord)
     G * u
 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