Add arbitrary topography on P4estMeshCubedSphere

examples/euler/dry_air/global_circulation/elixir_potential_temperature_vortex_shedding.jl

@@ -0,0 +1,142 @@
+# An idealized mountain-triggered mesoscale flow test case (DCMIP-2025 Test Case 2).
+# This setup examines the interaction between a stratified rotating flow and isolated
+# topography, leading to vortex shedding in the lee of the mountain.
+#
+# The test-case uses a reduced-radius Earth configuration (R = Rₑ/20), scaled rotation rate,
+# and a Gaussian mountain with peak height 2000 m.
+#
+# References:
+# - DCMIP-2025 Organizing Committee (2025)
+#   Test Case 2: Mountain-Triggered Meso-scale Flow Phenomena
+#   https://sites.google.com/umich.edu/dcmip-2025/dcmip-2025-test-cases/test-case-2-mountain-triggered-meso-scale-flow-phenomena
+
+using OrdinaryDiffEqSSPRK
+using Trixi, TrixiAtmo
+using LinearAlgebra: norm
+
+function cartesian_to_sphere(x)
+    r = norm(x)
+    lambda = atan(x[2], x[1])
+    if lambda < 0
+        lambda += 2 * pi
+    end
+    phi = asin(x[3] / r)
+
+    return lambda, phi, r
+end
+
+function initial_condition_vortex_shedding(x, t,
+                                           equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)
+    lon, lat, r = cartesian_to_sphere(x)
+    radius_earth = EARTH_RADIUS / 20
+    R = equations.c_p - equations.c_v
+    k = R / equations.c_p
+    # Convert spherical velocity to Cartesian
+    u0 = 20
+    v1 = -u0 * cos(lat) * sin(lon)
+    v2 = u0 * cos(lat) * cos(lon)
+    v3 = 0
+    g = EARTH_GRAVITATIONAL_ACCELERATION  # g
+    T = 288
+    angular_velocity = EARTH_ROTATION_RATE * 20  # Ω
+    r = norm(x)
+    # Make sure that r is not smaller than radius_earth
+    z = max(r - radius_earth, 0.0)
+    phi = g * z
+    N = 0.0182 # Brunt–Väisälä frequency
+    p = equations.p_0 * exp(-radius_earth * N^2 * u0 / (2 * g^2 * k) *
+            (u0 / radius_earth + 2 * angular_velocity) * (sin(lat)^2 - 1) -
+            N^2 / (g^2 * k) * phi)
+    rho = p / (R * T)
+    return prim2cons(SVector(rho, v1, v2, v3, p, phi), equations)
+end
+
+@inline function source_terms_coriolis(u, x, t,
+                                       equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)
+    radius_earth = EARTH_RADIUS / 20  # a
+    angular_velocity = EARTH_ROTATION_RATE * 20  # Ω
+
+    du1 = zero(eltype(u))
+    du2 = zero(eltype(u))
+    du3 = zero(eltype(u))
+    du4 = zero(eltype(u))
+    du5 = zero(eltype(u))
+    # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4]
+    du2 -= -2 * angular_velocity * u[3]
+    du3 -= 2 * angular_velocity * u[2]
+
+    return SVector(du1, du2, du3, du4, du5, zero(eltype(u)))
+end
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
+
+initial_condition = initial_condition_vortex_shedding
+
+boundary_conditions = (; inside = boundary_condition_slip_wall,
+                       outside = boundary_condition_slip_wall)
+
+polydeg = 5
+surface_flux = (FluxLMARS(340), flux_zero)
+volume_flux = (flux_tec, flux_nonconservative_souza_etal)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+function initial_topography_gaussian(lat, lon)
+    h_0 = 2000
+    h0 = 10000
+    d = 12500
+    lambda_c = pi
+    phi_c = 20 * pi / 180
+
+    r_earth = EARTH_RADIUS / 20
+    dist = r_earth * acos(clamp(sin(phi_c) * sin(lat) +
+                      cos(phi_c) * cos(lat) * cos(lon - lambda_c), -1, 1))
+
+    return h_0 * exp(-(dist / d)^2)
+end
+
+trees_per_cube_face = (16, 8)
+
+mesh = P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS / 20,
+                                      20000,
+                                      polydeg = polydeg,
+                                      initial_refinement_level = 0,
+                                      initial_topography = initial_topography_gaussian,
+                                      adapt_vertical_grid = GalChen())
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_coriolis,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+T = 1.0 # 20 small earth days
+tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 1 standard earth day
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 5760, save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution)
+
+###############################################################################
+# Use a Runge-Kutta method with automatic (error based) time step size control
+# Enable threading of the RK method for better performance on multiple threads
+tol = 1e-6
+sol = solve(ode,
+            SSPRK43(thread = Trixi.True());
+            abstol = tol, reltol = tol, ode_default_options()...,
+            callback = callbacks)

src/TrixiAtmo.jl

@@ -98,6 +98,8 @@ export bottom_topography_isolated_mountain, bottom_topography_unsteady_solid_bod
 
 export AtmosphereLayers, AtmosphereLayersRainyBubble
 
+export P4estMeshCubedSphereTopography, GalChen, Sleve
+
 export examples_dir
 
 end # module TrixiAtmo

src/meshes/meshes.jl

@@ -1,3 +1,4 @@
 include("p4est_cubed_sphere.jl")
 include("p4est_icosahedron_quads.jl")
 include("dgmulti_icosahedron_tri.jl")
+include("p4est_cubed_sphere_topography.jl")

src/meshes/p4est_cubed_sphere_topography.jl

@@ -0,0 +1,213 @@
+using Trixi: connectivity_cubed_sphere, new_p4est
+
+abstract type SmoothVerticalCoordinate end
+
+"""
+    Sleve(eta_H, s)
+
+SLEVE (Smooth LEvel VErtical) coordinate transformation for terrain-following
+vertical grids in atmospheric models. Ensures that the topographic influence
+decays smoothly with height, returning to flat levels well below the model top.
+
+The transformed height is given by:
+- For `η ≤ η_H`:
+  `z = η·H + z_s · sinh((η_H - η) / s / η_H) / sinh(1/s)`
+- For `η > η_H`:
+  `z = η·H`
+
+where `η = z_reference / H` is the normalized reference height.
+
+# Arguments
+- `eta_H`: fraction of the domain height `H` above which the
+           topography influence vanishes completely.
+           Must be in `(0, 1)`. Typical value: `0.7`.
+- `s`: controls the rate of decay of the topographic influence with height. 
+       Smaller values produce faster decay. Typical value: `0.8`.
+
+# References
+- Schär, C., Leuenberger, D., Fuhrer, O., Lüthi, D., & Frei, C. (2002)
+  A new terrain-following vertical coordinate formulation for atmospheric
+  prediction models
+  Monthly Weather Review, 130(10), 2459–2480
+  [DOI: 10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2)
+"""
+struct Sleve{RealT} <: SmoothVerticalCoordinate
+    eta_H::RealT  # fraction of H above which topography influence vanishes
+    s::RealT           # rate of decay of the topography influence with height
+end
+
+@inline function (adapt_sleve::Sleve)(z_reference, z_topography, H)
+    @unpack s, eta_H = adapt_sleve
+    eta = z_reference / H
+    if eta <= eta_H
+        z = eta * H +
+            z_topography * sinh((eta_H - eta) / s / eta_H) / sinh(one(z_reference) / s)
+    else
+        z = eta * H
+    end
+    return z
+end
+
+"""
+    GalChen()
+
+Gal-Chen & Somerville terrain-following vertical coordinate transformation.
+The topographic influence decays linearly with height, reaching zero at the
+top of the domain `H`.
+
+The transformed height is given by:
+  `z = z_reference + (H - z_reference) · z_s / H`
+
+which ensures:
+- At the surface (`z_reference = 0`): `z = z_s` (grid follows topography exactly)
+- At the top (`z_reference = H`): `z = H` (grid is flat, independent of topography)
+
+# References
+- Gal-Chen, T. & Somerville, R. C. J. (1975)
+  On the use of a coordinate transformation for the solution of the
+  Navier-Stokes equations
+  Journal of Computational Physics, 17(2), 209–228
+  [DOI: 10.1016/0021-9991(75)90037-6](https://doi.org/10.1016/0021-9991(75)90037-6)
+"""
+struct GalChen <: SmoothVerticalCoordinate end
+
+@inline function (adapt_galchen::GalChen)(z_reference, z_topography, H)
+    z = z_reference + (H - z_reference) * z_topography / H
+    return z
+end
+
+"""
+    P4estMeshCubedSphereTopography(trees_per_face_dimension, layers, inner_radius, thickness;
+                         polydeg, RealT=Float64,
+                         initial_refinement_level=0, unsaved_changes=true,
+                         p4est_partition_allow_for_coarsening=true, initial_topography, adapt_vertical_grid)
+
+Build a "Cubed Sphere" mesh as `P4estMesh` with
+`6 * trees_per_face_dimension^2 * layers` trees and a given topography.
+
+The mesh will have two boundaries, `:inside` and `:outside`.
+
+# Arguments
+- `trees_per_face_dimension::Integer`: the number of trees in the first two local dimensions of
+                                       each face.
+- `layers::Integer`: the number of trees in the third local dimension of each face, i.e., the number
+                     of layers of the sphere.
+- `inner_radius::RealT`: the inner radius of the sphere.
+- `thickness::RealT`: the thickness of the spherical shell. The outer radius will be `inner_radius + thickness`.
+- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh.
+                      The mapping will be approximated by an interpolation polynomial
+                      of the specified degree for each tree.
+- `RealT::Type`: the type that should be used for coordinates.
+- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts.
+- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
+- `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity
+                                                independent of domain partitioning. Should be `false` for static meshes
+                                                to permit more fine-grained partitioning.
+- `initial_topography`: initial topography as a function of x, y, z coordinates, that modifies the surface spherical layer.
+- `adaptive_vertical_grid`: smoothing of the vertical element size in the radial direction, to gradually restore the spherical shape. Two options are available: GalChen() or Sleve(etaH, s).        
+"""
+function P4estMeshCubedSphereTopography(trees_per_face_dimension, layers, inner_radius,
+                                        thickness;
+                                        polydeg, RealT = Float64,
+                                        initial_refinement_level = 0,
+                                        unsaved_changes = true,
+                                        p4est_partition_allow_for_coarsening = true,
+                                        initial_topography, adapt_vertical_grid = GalChen())
+    connectivity = connectivity_cubed_sphere(trees_per_face_dimension, layers)
+
+    n_trees = 6 * trees_per_face_dimension^2 * layers
+
+    basis = LobattoLegendreBasis(RealT, polydeg)
+    nodes = basis.nodes
+
+    tree_node_coordinates = Array{RealT, 5}(undef, 3,
+                                            ntuple(_ -> length(nodes), 3)...,
+                                            n_trees)
+    calc_tree_node_coordinates!(tree_node_coordinates, nodes,
+                                trees_per_face_dimension, layers,
+                                inner_radius, thickness, initial_topography,
+                                adapt_vertical_grid)
+
+    p4est = new_p4est(connectivity, initial_refinement_level)
+
+    boundary_names = fill(Symbol("---"), 2 * 3, n_trees)
+    boundary_names[5, :] .= Symbol("inside")
+    boundary_names[6, :] .= Symbol("outside")
+
+    return P4estMesh{3}(p4est, tree_node_coordinates, nodes,
+                        boundary_names, "", unsaved_changes,
+                        p4est_partition_allow_for_coarsening)
+end
+
+function cubed_sphere_mapping_topography(xi, eta, zeta, inner_radius, thickness, direction,
+                                         initial_topography, adapt_vertical_grid)
+    alpha = xi * pi / 4
+    beta = eta * pi / 4
+
+    x = tan(alpha)
+    y = tan(beta)
+
+    cube_coordinates = (SVector(-1, -x, y),
+                        SVector(1, x, y),
+                        SVector(x, -1, y),
+                        SVector(-x, 1, y),
+                        SVector(-x, y, -1),
+                        SVector(x, y, 1))
+
+    r = sqrt(1 + x^2 + y^2)
+
+    unit_pt = cube_coordinates[direction] / r
+    lat = asin(unit_pt[3] / r)
+    lon = atan(unit_pt[2], unit_pt[1])
+    z_topography = initial_topography(lat, lon)
+
+    z_reference = thickness * 0.5 * (zeta + 1)
+    z = adapt_vertical_grid(z_reference, z_topography, thickness)
+    return (inner_radius + z) / r * cube_coordinates[direction]
+end
+
+function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5},
+                                     nodes, trees_per_face_dimension, layers,
+                                     inner_radius, thickness, initial_topography,
+                                     adapt_vertical_grid)
+    n_cells_x = n_cells_y = trees_per_face_dimension
+    n_cells_z = layers
+
+    linear_indices = LinearIndices((n_cells_x, n_cells_y, n_cells_z, 6))
+
+    dx = 2 / n_cells_x
+    dy = 2 / n_cells_y
+    dz = 2 / n_cells_z
+
+    for direction in 1:6
+        for cell_z in 1:n_cells_z, cell_y in 1:n_cells_y, cell_x in 1:n_cells_x
+            tree = linear_indices[cell_x, cell_y, cell_z, direction]
+
+            x_offset = -1 + (cell_x - 1) * dx + dx / 2
+            y_offset = -1 + (cell_y - 1) * dy + dy / 2
+            z_offset = -1 + (cell_z - 1) * dz + dz / 2
+
+            for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes)
+                node_coordinates[:, i, j, k, tree] .= cubed_sphere_mapping_topography(x_offset +
+                                                                                      dx /
+                                                                                      2 *
+                                                                                      nodes[i],
+                                                                                      y_offset +
+                                                                                      dy /
+                                                                                      2 *
+                                                                                      nodes[j],
+                                                                                      z_offset +
+                                                                                      dz /
+                                                                                      2 *
+                                                                                      nodes[k],
+                                                                                      inner_radius,
+                                                                                      thickness,
+                                                                                      direction,
+                                                                                      initial_topography,
+                                                                                      adapt_vertical_grid)
+            end
+        end
+    end
+
+    return nothing
+end