diff --git a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_vortex_shedding.jl b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_vortex_shedding.jl new file mode 100644 index 00000000..8e53297d --- /dev/null +++ b/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) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index faede1a6..68a21948 100644 --- a/src/TrixiAtmo.jl +++ b/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 diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 17ca7777..8ab7d4c4 100644 --- a/src/meshes/meshes.jl +++ b/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") diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl new file mode 100644 index 00000000..367c808d --- /dev/null +++ b/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 diff --git a/test/test_3d_potential_temperature.jl b/test/test_3d_potential_temperature.jl index 7b4af7e2..1b4ba8f5 100644 --- a/test/test_3d_potential_temperature.jl +++ b/test/test_3d_potential_temperature.jl @@ -210,4 +210,58 @@ end @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) end +@trixi_testset "elixir_potential_temperature_vortex_shedding" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "global_circulation", + "elixir_potential_temperature_vortex_shedding.jl"), + l2=[ + 0.00011127272731957073, + 0.036136807950192446, + 0.036132932940041926, + 0.047765679229021134, + 0.03287203491711795, + 5.112500661735035 + ], + linf=[ + 0.0011895976218974091, + 0.16161469307155263, + 0.15524551919165414, + 0.25696936007842397, + 0.34356951757109755, + 42.337304272491956 + ], + rtol=1e-9, + tspan=(0.0, 0.0001 * SECONDS_PER_DAY), + trees_per_cube_face=(2, 2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) +end + +@trixi_testset "elixir_potential_temperature_vortex_shedding with Sleve" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "global_circulation", + "elixir_potential_temperature_vortex_shedding.jl"), + l2=[ + 0.00011127562057264662, + 0.03613675415332188, + 0.03613293640311338, + 0.04776567584661323, + 0.03287272064361882, + 5.114726513677777 + ], + linf=[ + 0.001189598854849594, + 0.1613051140325089, + 0.15524552078160023, + 0.25697003550490627, + 0.3435698663628841, + 42.337304272491956 + ], + rtol=1e-9, + tspan=(0.0, 0.0001 * SECONDS_PER_DAY), + trees_per_cube_face=(2, 2), adapt_vertical_grid=Sleve(0.7, 0.8)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) +end + end