add lengthscale functions to mesh and spaces

.buildkite/pipeline.yml

@@ -5,7 +5,7 @@ agents:
 
 env:
   JULIA_LOAD_PATH: "${JULIA_LOAD_PATH}:${BUILDKITE_BUILD_CHECKOUT_PATH}/.buildkite"
-  JULIA_DEPOT_PATH: "${BUILDKITE_BUILD_PATH}/${BUILDKITE_PIPELINE_SLUG}/depot/default"
+  # JULIA_DEPOT_PATH: "${BUILDKITE_BUILD_PATH}/${BUILDKITE_PIPELINE_SLUG}/depot/default"
   JULIA_MAX_NUM_PRECOMPILE_FILES: 100
   JULIA_CPU_TARGET: 'broadwell;skylake'
   JULIA_NVTX_CALLBACKS: gc

docs/src/api.md

@@ -119,6 +119,7 @@ Meshes.SharedVertices
 Meshes.face_connectivity_matrix
 Meshes.vertex_connectivity_matrix
 Meshes.linearindices
+Meshes.element_horizontal_length_scale
 ```
 
 ## Topologies
@@ -195,6 +196,10 @@ Spaces.SpectralElementSpace2D
 Spaces.SpectralElementSpaceSlab
 ```
 
+```@docs
+Spaces.node_horizontal_length_scale
+```
+
 ### Quadratures
 
 

examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl

@@ -39,7 +39,6 @@ const parameters = (
     ρ₀ = 1.0, # reference density
     c = 2,
     g = 10,
-    D₄ = 1e-4, # hyperdiffusion coefficient
 )
 
 domain = Domains.RectangleDomain(
@@ -105,7 +104,10 @@ function energy(state, p, local_geometry)
 end
 
 function rhs!(dydt, y, _, t)
-    D₄ = parameters.D₄
+    space = axes(y)
+    c = sqrt(parameters.g * parameters.ρ₀)
+    D₄ = 0.0015 * c * Spaces.node_horizontal_length_scale(space)^3 # hyperdiffusion coefficient
+
     g = parameters.g
 
     sdiv = Operators.Divergence()

examples/sphere/shallow_water.jl

@@ -42,7 +42,7 @@ Base.@kwdef struct PhysicalParameters{FT} # rename to PhysicalParameters
     "Gravitational constant"
     g::FT = FT(9.80616)
     "Hyperdiffusion coefficient"
-    D₄::FT = FT(1.0e16)
+    ν₄::FT = FT(0.25)
 end
 #This example solves the shallow-water equations on a cubed-sphere manifold.
 #This file contains five test cases:
@@ -256,7 +256,7 @@ end
 # Set initial condition
 function set_initial_condition(space, test::RossbyHaurwitzTest)
     (; a, h0, ω, K, α, params) = test
-    (; R, Ω, g, D₄) = params
+    (; R, Ω, g) = params
     Y = map(Fields.local_geometry_field(space)) do local_geometry
         coord = local_geometry.coordinates
         ϕ = coord.lat
@@ -296,7 +296,7 @@ end
 
 function set_initial_condition(space, test::SteadyStateCompactTest)
     (; u0, h0, ϕᵦ, ϕₑ, xₑ, α, params) = test
-    (; R, Ω, g, D₄) = params
+    (; R, Ω, g) = params
     Y = map(Fields.local_geometry_field(space)) do local_geometry
         coord = local_geometry.coordinates
 
@@ -369,7 +369,7 @@ function set_initial_condition(space, test::SteadyStateCompactTest)
 end
 function set_initial_condition(space, test::BarotropicInstabilityTest)
     (; u_max, αₚ, βₚ, h0, h_hat, ϕ₀, ϕ₁, ϕ₂, eₙ, α, params) = test
-    (; R, Ω, g, D₄) = params
+    (; R, Ω, g) = params
 
     Y = map(Fields.local_geometry_field(space)) do local_geometry
         coord = local_geometry.coordinates
@@ -428,7 +428,7 @@ function set_initial_condition(
     test::T,
 ) where {T <: Union{MountainTest, SteadyStateTest}}
     (; u0, h0, α, params) = test
-    (; R, Ω, g, D₄) = params
+    (; R, Ω, g) = params
     Y = map(Fields.local_geometry_field(space)) do local_geometry
         coord = local_geometry.coordinates
 
@@ -455,7 +455,11 @@ end
 function rhs!(dYdt, y, parameters, t)
     @nvtx "rhs!" color = colorant"red" begin
         (; f, h_s, ghost_buffer, params) = parameters
-        (; D₄, g) = params
+        (; ν₄, g) = params
+
+        space = axes(y)
+        D₄ = ν₄ * Spaces.node_horizontal_length_scale(space)^3 # hyperdiffusion coefficient
+
 
         div = Operators.Divergence()
         wdiv = Operators.WeakDivergence()
@@ -551,6 +555,8 @@ function shallow_water_driver(ARGS, ::Type{FT}) where {FT}
         global_space =
             space = Spaces.SpectralElementSpace2D(grid_topology, quad)
     end
+    @show Spaces.node_horizontal_length_scale(space)^3
+
     coords = Fields.coordinate_field(space)
     f = set_coriolis_parameter(space, test)
     h_s = surface_topography(space, test)

src/Meshes/Meshes.jl

@@ -107,7 +107,12 @@ The element and face (`opelem`, `opface`) that oppose face `face` of element `el
 """
 function opposing_face end
 
+"""
+    Meshes.element_horizontal_length_scale(mesh::AbstractMesh)
 
+The approximate length scale (in units of distance) of the elements of the mesh.
+"""
+function element_horizontal_length_scale end
 
 include("common.jl")
 include("interval.jl")

src/Meshes/cubedsphere.jl

@@ -89,6 +89,12 @@ end
 
 domain(mesh::AbstractCubedSphere) = mesh.domain
 elements(mesh::AbstractCubedSphere) = CartesianIndices((mesh.ne, mesh.ne, 6))
+nelements(mesh::AbstractCubedSphere) = mesh.ne * mesh.ne * 6
+
+function element_horizontal_length_scale(mesh::AbstractCubedSphere)
+    FT = typeof(mesh.domain.radius)
+    return FT(sqrt(4 * pi / 6)) * mesh.domain.radius / mesh.ne
+end
 
 is_boundary_face(mesh::AbstractCubedSphere, elem, face) = false
 boundary_face_name(mesh::AbstractCubedSphere, elem, face) = nothing

src/Meshes/interval.jl

@@ -34,6 +34,11 @@ function Base.show(io::IO, mesh::IntervalMesh)
     print(io, " of ")
     print(io, mesh.domain)
 end
+function element_horizontal_length_scale(mesh::IntervalMesh)
+    cmax = Geometry.component(mesh.domain.coord_max, 1)
+    cmin = Geometry.component(mesh.domain.coord_min, 1)
+    return (cmax - cmin) / nelements(mesh)
+end
 
 coordinates(mesh::IntervalMesh, elem::Integer, vert::Integer) =
     mesh.faces[elem + vert - 1]

src/Meshes/rectangle.jl

@@ -35,6 +35,11 @@ domain(mesh::RectilinearMesh) =
     RectangleDomain(domain(mesh.intervalmesh1), domain(mesh.intervalmesh2))
 nelements(mesh::RectilinearMesh) =
     nelements(mesh.intervalmesh1) * nelements(mesh.intervalmesh2)
+
+element_horizontal_length_scale(mesh::RectilinearMesh) = sqrt(
+    element_horizontal_length_scale(mesh.intervalmesh1) *
+    element_horizontal_length_scale(mesh.intervalmesh2),
+)
 function elements(mesh::RectilinearMesh)
     # we use the Base Julia CartesianIndices object to index elements in the mesh
     CartesianIndices((

src/Quadratures/Quadratures.jl

@@ -53,6 +53,7 @@ Base.show(io::IO, ::GLL{Nq}) where {Nq} =
 
 @inline polynomial_degree(::GLL{Nq}) where {Nq} = Int(Nq - 1)
 @inline degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq)
+unique_degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq - 1)
 
 @generated function quadrature_points(::Type{FT}, ::GLL{Nq}) where {FT, Nq}
     points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.both)
@@ -70,6 +71,7 @@ Base.show(io::IO, ::GL{Nq}) where {Nq} =
 
 @inline polynomial_degree(::GL{Nq}) where {Nq} = Int(Nq - 1)
 @inline degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq)
+unique_degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq)
 
 @generated function quadrature_points(::Type{FT}, ::GL{Nq}) where {FT, Nq}
     points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.neither)

src/Spaces/spectralelement.jl

@@ -216,6 +216,19 @@ nlevels(space::SpectralElementSpaceSlab1D) = 1
 nlevels(space::SpectralElementSpaceSlab2D) = 1
 
 
+"""
+    Spaces.node_horizontal_length_scale(space::AbstractSpectralElementSpace)
+
+The approximate length scale of the distance between nodes. This is defined as the
+length scale of the mesh (see [`Meshes.element_horizontal_length_scale`](@ref)), divided by the
+number of unique quadrature points along each dimension.
+"""
+function node_horizontal_length_scale(space::AbstractSpectralElementSpace)
+    quad = quadrature_style(space)
+    Nu = Quadratures.unique_degrees_of_freedom(quad)
+    return Meshes.element_horizontal_length_scale(space.topology.mesh) / Nu
+end
+
 
 
 

test/Meshes/cubedsphere.jl

@@ -80,6 +80,8 @@ end
                     end
                 end
             end
+            @test Meshes.element_horizontal_length_scale(mesh) ≈
+                  sqrt((4pi * 5^2) / Meshes.nelements(mesh))
         end
     end
 end

test/Meshes/interval.jl

@@ -48,6 +48,7 @@ end
         @test length(mesh.faces) == nelems + 1 # number of faces is +1 elements
         @test Meshes.nelements(mesh) == nelems
         @test Meshes.elements(mesh) == UnitRange(1, nelems)
+        @test Meshes.element_horizontal_length_scale(mesh) ≈ 1 / nelems
     end
     dom, mesh = unit_intervalmesh(nelems = 2)
     @test Meshes.domain(mesh) isa Domains.IntervalDomain

test/Meshes/rectangle.jl

@@ -51,6 +51,13 @@ meshes = [
           CartesianIndices((2, 3))
 end
 
+@testset "element_horizontal_length_scale" begin
+    for mesh in meshes
+        @test Meshes.element_horizontal_length_scale(mesh) ≈
+              sqrt(1 / Meshes.nelements(mesh))
+    end
+end
+
 @testset "opposing face" begin
     for mesh in meshes
         for elem in Meshes.elements(mesh)

test/Spaces/sphere.jl

@@ -107,6 +107,9 @@ end
     quad = Quadratures.GLL{Nq}()
     horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
 
+    @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)