From bfa6098b5ebe4447a1dbdb7b0f60768bf7daed18 Mon Sep 17 00:00:00 2001 From: Akshay Sridhar Date: Thu, 15 Jun 2023 13:02:57 -0700 Subject: [PATCH] Add SLEVE coordinate, update API docs. modified: docs/refs.bib modified: src/Hypsography/Hypsography.jl Restructure tests --- docs/refs.bib | 13 ++ docs/src/api.md | 2 + src/Hypsography/Hypsography.jl | 39 ++++- test/Spaces/terrain_warp.jl | 270 +++++++++++++++++++++------------ 4 files changed, 230 insertions(+), 94 deletions(-) diff --git a/docs/refs.bib b/docs/refs.bib index ba39434f8e..354bdd5689 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -146,6 +146,19 @@ @article{Ronchi1996 doi = {10.1006/jcph.1996.0047} } +@article{Schar2002, + author = {Christoph Schär and Daniel Leuenberger and Oliver Fuhrer and Daniel Lüthi and Claude Girard"}, + title = {A New Terrain-Following Vertical Coordinate Formulation for Atmospheric Prediction Models"}, + journal = {Monthly Weather Review}, + year = {2002}, + publisher = {American Meteorological Society}, + address = {Boston MA, USA}, + volume = {130}, + number = {10}, + doi = {https://doi.org/10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2}, + pages= {2459 - 2480}, +} + @article{Taylor2010, title={A compatible and conservative spectral element method on unstructured grids}, author={Taylor, Mark A and Fournier, Aim\'{e}}, diff --git a/docs/src/api.md b/docs/src/api.md index 0be6bf0080..51555c17d3 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -250,6 +250,8 @@ Fields.Δz_field ## Hypsography ```@docs +Hypsography.LinearAdaption +Hypsography.SLEVEAdaption Hypsography.diffuse_surface_elevation! ``` diff --git a/src/Hypsography/Hypsography.jl b/src/Hypsography/Hypsography.jl index 8c87c9cabf..f8c61d84f4 100644 --- a/src/Hypsography/Hypsography.jl +++ b/src/Hypsography/Hypsography.jl @@ -18,6 +18,23 @@ struct LinearAdaption{F <: Union{Fields.Field, Nothing}} <: HypsographyAdaption surface::F end +""" + SLEVEAdaption(surface::Field, ηₕ::FT, s::FT) + +Locate vertical levels using an exponential function between the surface field and the top +of the domain, using the method of [Schar2002](@cite). This method is modified +such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper +bounds represent the domain bottom and top respectively. `s` governs the decay rate. +If the decay-scale is poorly specified (i.e., `s * zₜ` is lower than the maximum +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 +end + # deprecated in 0.10.12 LinearAdaption() = LinearAdaption(nothing) @@ -66,11 +83,31 @@ function ExtrudedFiniteDifferenceSpace( grad = Operators.Gradient() wdiv = Operators.WeakDivergence() + z_surface = Fields.field_values(adaption.surface) + + FT = eltype(z_surface) + # TODO: Function dispatch if adaption isa LinearAdaption - z_surface = Fields.field_values(adaption.surface) fZ_data = @. z_ref + (1 - z_ref / z_top) * z_surface fZ = Fields.Field(fZ_data, face_space) + elseif adaption isa SLEVEAdaption + ηₕ = adaption.ηₕ + s = adaption.s + @assert FT(0) <= ηₕ <= FT(1) + @assert s >= FT(0) + η = @. z_ref ./ z_top + 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))). + \n Returning (5/3)*max(zₛ)" + s = @. maximum(z_surface) / z_top * eltype(z_surface) * (5 // 3) + end + fZ_data = @. ifelse( + η <= ηₕ, + η * z_top + z_surface * (sinh((ηₕ - η) / s / ηₕ)) / (sinh(1 / s)), + η * z_top, + ) + fZ = Fields.Field(fZ_data, face_space) end # Take the horizontal gradient for the Z surface field diff --git a/test/Spaces/terrain_warp.jl b/test/Spaces/terrain_warp.jl index 817264724f..fb2efd7670 100644 --- a/test/Spaces/terrain_warp.jl +++ b/test/Spaces/terrain_warp.jl @@ -1,5 +1,6 @@ using Test using ClimaComms +using IntervalSets import ClimaCore: ClimaCore, @@ -36,17 +37,17 @@ function warp_sinsq_3d(coord) y = Geometry.component(coord, 2) eltype(x)(0.5) * sin(x)^2 * sin(y)^2 end -function warpedspace_2D( - FT = Float64, - xlim = (0, π), - zlim = (0, 1), - helem = 2, - velem = 10, - npoly = 5; - stretch = Meshes.Uniform(), - warp_fn = warp_sin_2d, - test_smoothing = false, +function generate_base_spaces( + xlim, + zlim, + helem, + velem, + npoly, + stretch = Meshes.Uniform(); + ndims = 3, ) + comms_context = ClimaComms.SingletonCommsContext(ClimaComms.CPUDevice()) + FT = eltype(xlim) vertdomain = Domains.IntervalDomain( Geometry.ZPoint{FT}(zlim[1]), Geometry.ZPoint{FT}(zlim[2]); @@ -56,22 +57,43 @@ function warpedspace_2D( vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) # Generate Horizontal Space - horzdomain = Domains.IntervalDomain( - Geometry.XPoint{FT}(xlim[1]), - Geometry.XPoint{FT}(xlim[2]); - periodic = true, - ) - horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem) - horztopology = Topologies.IntervalTopology(horzmesh) quad = Spaces.Quadratures.GLL{npoly + 1}() - hspace = Spaces.SpectralElementSpace1D(horztopology, quad) - + if ndims == 2 + horzdomain = Domains.IntervalDomain( + Geometry.XPoint{FT}(xlim[1]), + Geometry.XPoint{FT}(xlim[2]); + periodic = true, + ) + horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem) + horztopology = Topologies.IntervalTopology(comms_context, horzmesh) + hspace = Spaces.SpectralElementSpace1D(horztopology, quad) + elseif ndims == 3 + horzdomain = Domains.RectangleDomain( + Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]), + Geometry.YPoint{FT}(xlim[1]) .. Geometry.YPoint{FT}(xlim[2]), + x1periodic = true, + x2periodic = true, + ) + # Assume same number of elems (helem) in (x,y) directions + horzmesh = Meshes.RectilinearMesh(horzdomain, helem, helem) + horztopology = Topologies.Topology2D(comms_context, horzmesh) + hspace = Spaces.SpectralElementSpace2D(horztopology, quad) + end + return vert_face_space, hspace +end +function generate_smoothed_orography( + hspace, + warp_fn::Function, + helem; + test_smoothing::Bool = false, +) # Extrusion z_surface = 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) dx = x_array[2] - x_array[1] + FT = eltype(x_array) κ = FT(1 / helem) test_smoothing ? Hypsography.diffuse_surface_elevation!( @@ -80,16 +102,47 @@ function warpedspace_2D( maxiter = 10^5, dt = FT(dx / 100), ) : nothing + return z_surface +end + +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), + ) + end +end + +function warpedspace_2D( + FT = Float64, + xlim = (0, π), + zlim = (0, 1), + helem = 2, + velem = 10, + npoly = 5, + stretch = Meshes.Uniform(); + warp_fn = warp_sin_2d, + test_smoothing = false, + adaption = Hypsography.LinearAdaption, +) + vert_face_space, hspace = + generate_base_spaces(xlim, zlim, helem, velem, npoly, ndims = 2) + z_surface = + generate_smoothed_orography(hspace, warp_fn, helem; test_smoothing) + mesh_adapt = get_adaptation(adaption, z_surface) f_space = Spaces.ExtrudedFiniteDifferenceSpace( hspace, vert_face_space, - Hypsography.LinearAdaption(z_surface), + mesh_adapt, ) c_space = Spaces.CenterExtrudedFiniteDifferenceSpace(f_space) return (c_space, f_space) end - function hybridspace_2D( FT = Float64, xlim = (0, π), @@ -99,87 +152,38 @@ function hybridspace_2D( npoly = 5; stretch = Meshes.Uniform(), ) - vertdomain = Domains.IntervalDomain( - Geometry.ZPoint{FT}(zlim[1]), - Geometry.ZPoint{FT}(zlim[2]); - boundary_tags = (:bottom, :top), - ) - vertmesh = Meshes.IntervalMesh(vertdomain, stretch, nelems = velem) - vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) - - # Generate Horizontal Space - horzdomain = Domains.IntervalDomain( - Geometry.XPoint{FT}(xlim[1]), - Geometry.XPoint{FT}(xlim[2]); - periodic = true, - ) - horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem) - horztopology = Topologies.IntervalTopology(horzmesh) - quad = Spaces.Quadratures.GLL{npoly + 1}() - hspace = Spaces.SpectralElementSpace1D(horztopology, quad) - + vert_face_space, hspace = + generate_base_spaces(xlim, zlim, helem, velem, npoly, ndims = 2) # Extrusion f_space = Spaces.ExtrudedFiniteDifferenceSpace(hspace, vert_face_space) c_space = Spaces.CenterExtrudedFiniteDifferenceSpace(f_space) return (c_space, f_space) end - function warpedspace_3D( FT = Float64, xlim = (0, π), ylim = (0, π), zlim = (0, 1), - xelem = 2, - yelem = 2, - zelem = 10, + helem = 2, + velem = 10, npoly = 5; stretch = Meshes.Uniform(), warp_fn = warp_sincos_3d, test_smoothing = false, + adaption = Hypsography.LinearAdaption, ) - vertdomain = Domains.IntervalDomain( - Geometry.ZPoint{FT}(zlim[1]), - Geometry.ZPoint{FT}(zlim[2]); - boundary_tags = (:bottom, :top), - ) - vertmesh = Meshes.IntervalMesh(vertdomain, stretch, nelems = zelem) - vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) - - # Generate Horizontal Space - xdomain = Domains.IntervalDomain( - Geometry.XPoint{FT}(xlim[1]), - Geometry.XPoint{FT}(xlim[2]), - periodic = true, - ) - ydomain = Domains.IntervalDomain( - Geometry.YPoint{FT}(ylim[1]), - Geometry.YPoint{FT}(ylim[2]), - periodic = true, - ) - horzdomain = Domains.RectangleDomain(xdomain, ydomain) - horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem) - horztopology = - Topologies.Topology2D(ClimaComms.SingletonCommsContext(), horzmesh) - quad = Spaces.Quadratures.GLL{npoly + 1}() - hspace = Spaces.SpectralElementSpace2D(horztopology, quad) + vert_face_space, hspace = + generate_base_spaces(xlim, zlim, helem, velem, npoly) # Extrusion - z_surface = warp_fn.(Fields.coordinate_field(hspace)) - x_array = parent(Fields.coordinate_field(hspace).x) - dx = x_array[2] - x_array[1] - κ = FT(1 / max(xelem, yelem)) - test_smoothing ? - Hypsography.diffuse_surface_elevation!( - z_surface; - κ, - maxiter = 10^5, - dt = FT(dx / 100), - ) : nothing + z_surface = + generate_smoothed_orography(hspace, warp_fn, helem; test_smoothing) + mesh_adapt = get_adaptation(adaption, z_surface) f_space = Spaces.ExtrudedFiniteDifferenceSpace( hspace, vert_face_space, - Hypsography.LinearAdaption(z_surface), + mesh_adapt, ) c_space = Spaces.CenterExtrudedFiniteDifferenceSpace(f_space) @@ -199,17 +203,6 @@ end for nl in levels, np in polynom, nh in horzelem ʷhv_center_space, ʷhv_face_space = warpedspace_2D(FT, (xmin, xmax), (zmin, zmax), nh, nl, np;) - hv_center_space, hv_face_space = - hybridspace_2D(FT, (xmin, xmax), (zmin, zmax), nh, nl, np) - ⁿhv_center_space, ⁿhv_face_space = warpedspace_2D( - FT, - (xmin, xmax), - (zmin, zmax), - nh, - nl, - np; - warp_fn = warp_sin_2d, - ) ʷᶜcoords = Fields.coordinate_field(ʷhv_center_space) ʷᶠcoords = Fields.coordinate_field(ʷhv_face_space) z₀ = ClimaCore.Fields.level(ʷᶜcoords.z, 1) @@ -321,7 +314,6 @@ end (ymin, ymax), (zmin, zmax), nh, - nh, nl, np; ) @@ -381,3 +373,95 @@ end end end end + +@testset "Interior Mesh `Adaption` ηₕ Test" begin + # Test interior mesh in different adaptation types + for meshadapt in (Hypsography.SLEVEAdaption,) + for FT in (Float32, Float64) + xmin, xmax = FT(0), FT(π) + zmin, zmax = FT(0), FT(1) + nl = 10 + np = 3 + nh = 4 + ʷhv_center_space, ʷhv_face_space = warpedspace_2D( + FT, + (xmin, xmax), + (zmin, zmax), + nh, + nl, + np; + warp_fn = warp_sin_2d, + adaption = meshadapt, + ) + hv_center_space, hv_face_space = + hybridspace_2D(FT, (xmin, xmax), (zmin, zmax), nh, nl, np) + ʷᶜcoords = Fields.coordinate_field(ʷhv_center_space) + ʷᶠcoords = Fields.coordinate_field(ʷhv_face_space) + ᶜcoords = Fields.coordinate_field(hv_center_space) + ᶠcoords = Fields.coordinate_field(hv_face_space) + # Check ηₛ = 0.75 is correctly applied. + # Expectation: ≈zero difference between unwarped and warped coordinates for η >= ηₕ, where η = z / zₜ + r1 = + ( + parent(ʷᶜcoords)[8:10, :, 2, :] .- + parent(ᶜcoords)[8:10, :, 2, :] + ) ./ parent(ᶜcoords)[8:10, :, 2, :] + @test maximum(r1) <= FT(0.015) + end + end +end + +@testset "Interior Mesh `Adaption` (ηₕ=1, s=1) Test" begin + # Test interior mesh in different adaptation types + for meshadapt in (Hypsography.SLEVEAdaption,) + for FT in (Float32, Float64) + xlim = (FT(0), FT(π)) + zlim = (FT(0), FT(1)) + nl = 10 + np = 3 + nh = 4 + vertdomain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(zlim[1]), + Geometry.ZPoint{FT}(zlim[2]); + boundary_names = (:bottom, :top), + ) + vertmesh = Meshes.IntervalMesh(vertdomain, nelems = nl) + vert_face_space = Spaces.FaceFiniteDifferenceSpace(vertmesh) + + horzdomain = Domains.IntervalDomain( + Geometry.XPoint{FT}(xlim[1]), + Geometry.XPoint{FT}(xlim[2]); + periodic = true, + ) + horzmesh = Meshes.IntervalMesh(horzdomain, nelems = nh) + horztopology = Topologies.IntervalTopology(horzmesh) + + quad = Spaces.Quadratures.GLL{np + 1}() + hspace = Spaces.SpectralElementSpace1D(horztopology, quad) + + # Generate surface elevation profile + z_surface = warp_sin_2d.(Fields.coordinate_field(hspace)) + # Generate space with known mesh-warp parameters ηₕ = 1; s = 1 + fspace = Spaces.ExtrudedFiniteDifferenceSpace( + hspace, + vert_face_space, + Hypsography.SLEVEAdaption(z_surface, FT(1), FT(1)), + ) + for i in 1:(nl + 1) + z_extracted = Fields.Field( + Fields.level(fspace.face_local_geometry.coordinates.z, i), + fspace, + ) + η = FT((i - 1) / 10) + z_surface_known = + @. FT(η) + z_surface * FT(sinh(1 - η) / sinh(1)) + @test maximum( + abs.( + Fields.field_values(z_extracted) .- + Fields.field_values(z_surface_known) + ), + ) <= FT(1e-6) + end + end + end +end