From e308e40797ec4c52266a3721f1b80e76e5437308 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 5 May 2026 20:38:08 +0200 Subject: [PATCH 01/11] arbitrary horography --- ...emperature_baroclinic_instability_earth.jl | 266 ++++++++++++++++++ src/meshes/meshes.jl | 1 + src/meshes/p4est_cubed_sphere_horography.jl | 156 ++++++++++ 3 files changed, 423 insertions(+) create mode 100644 examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl create mode 100644 src/meshes/p4est_cubed_sphere_horography.jl diff --git a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl new file mode 100644 index 00000000..a48a8a6e --- /dev/null +++ b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl @@ -0,0 +1,266 @@ +# An idealized baroclinic instability test case +# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. +# +# This elixir takes about 8 hours, using 16 threads of an AMD Ryzen 7 7800X3D. +# +# References: +# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) +# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores +# https://doi.org/10.1002/qj.2241 + +using OrdinaryDiffEqSSPRK +using Trixi, TrixiAtmo +using LinearAlgebra: norm +using NetCDF + +# Unperturbed balanced steady-state. +# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). +# The other velocity components are zero. +function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + radius_earth = EARTH_RADIUS # a + half_width_parameter = 2 # b + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g + k = 3 # k + surface_pressure = 1.0f5 # p₀ + gas_constant = 287 # R + surface_equatorial_temperature = 310 # T₀ᴱ + surface_polar_temperature = 240 # T₀ᴾ + lapse_rate = 0.005 # Γ + angular_velocity = EARTH_ROTATION_RATE # Ω + + # Distance to the center of the Earth + r = z + radius_earth + + # In the paper: T₀ + temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) + # In the paper: A, B, C, H + const_a = 1 / lapse_rate + const_b = (temperature0 - surface_polar_temperature) / + (temperature0 * surface_polar_temperature) + const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / + (surface_equatorial_temperature * surface_polar_temperature) + const_h = gas_constant * temperature0 / gravitational_acceleration + + # In the paper: (r - a) / bH + scaled_z = z / (half_width_parameter * const_h) + + # Temporary variables + temp1 = exp(lapse_rate / temperature0 * z) + temp2 = exp(-scaled_z^2) + + # In the paper: ̃τ₁, ̃τ₂ + tau1 = const_a * lapse_rate / temperature0 * temp1 + + const_b * (1 - 2 * scaled_z^2) * temp2 + tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 + + # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' + inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 + inttau2 = const_c * z * temp2 + + # Temporary variables + temp3 = r / radius_earth * cos(lat) + temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) + + # In the paper: T + temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) + + # In the paper: U, u (zonal wind, first component of spherical velocity) + big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * + (temp3^(k - 1) - temp3^(k + 1)) + temp5 = radius_earth * cos(lat) + u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) + + # Hydrostatic pressure + p = surface_pressure * + exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) + + # Density (via ideal gas law) + rho = p / (gas_constant * temperature) + + return rho, u, p +end + +# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) +function perturbation_stream_function(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + perturbation_radius = 1 / 6 # d₀ / a + perturbed_wind_amplitude = 1 # Vₚ + perturbation_lon = pi / 9 # Longitude of perturbation location + perturbation_lat = 2 * pi / 9 # Latitude of perturbation location + pertz = 15000 # Perturbation height cap + + # Great circle distance (d in the paper) divided by a (radius of the Earth) + # because we never actually need d without dividing by a + great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + + cos(perturbation_lat) * cos(lat) * + cos(lon - perturbation_lon)) + + # In the first case, the vertical taper function is per definition zero. + # In the second case, the stream function is per definition zero. + if z > pertz || great_circle_distance_by_a > perturbation_radius + return 0, 0 + end + + # Vertical tapering of stream function + perttaper = 1 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 + + # sin/cos(pi * d / (2 * d_0)) in the paper + sin_, cos_ = sincos(0.5f0 * pi * great_circle_distance_by_a / perturbation_radius) + + # Common factor for both u and v + factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ + + u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + + cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / + sin(great_circle_distance_by_a) + + v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / + sin(great_circle_distance_by_a) + + return u_perturbation, v_perturbation +end + +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 + +# Initial condition for an idealized baroclinic instability test +# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A +function initial_condition_baroclinic_instability(x, t, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + lon, lat, r = cartesian_to_sphere(x) + radius_earth = EARTH_RADIUS + # Make sure that the r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + + # Unperturbed basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Stream function type perturbation + u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) + + u += u_perturbation + v = v_perturbation + + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u - sin(lat) * cos(lon) * v + v2 = cos(lon) * u - sin(lat) * sin(lon) * v + v3 = cos(lat) * v + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0) + if z > 0 + r = norm(x) + else + r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) + end + r = -norm(x) + phi = radius_earth^2 * gravitational_acceleration / r + + return prim2cons(SVector(rho, v1, v2, v3, p, phi), equations) +end + +@inline function source_terms_baroclinic_instability(u, x, t, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + radius_earth = EARTH_RADIUS # a + angular_velocity = EARTH_ROTATION_RATE # Ω + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0) + r = z + radius_earth + + 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_baroclinic_instability + +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)) +trees_per_cube_face = (8, 4) +filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_bed.nc" +topo = ncread(filename, "z") + +function makepositive(x) + if x >= 0 + return x + else + return zero(typeof(x)) + end +end + +topo = makepositive.(topo) +# Add extra column and row for periodicity +topo = [topo; topo[1:1, :]] +topo = [topo topo[:, 1:1]] + +delta_topo = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians + +mesh = TrixiAtmo.P4estMeshCubedSphereHorography(trees_per_cube_face..., EARTH_RADIUS, 30000, + polydeg = polydeg, initial_refinement_level = 0, horography = topo, delta_horography = delta_topo, adapt_vertical_grid = TrixiAtmo.GalChen()) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_baroclinic_instability, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. +T = 0.01 # 10 days +tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 10 days + +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 = 100, 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/meshes/meshes.jl b/src/meshes/meshes.jl index 17ca7777..198d6858 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_horography.jl") diff --git a/src/meshes/p4est_cubed_sphere_horography.jl b/src/meshes/p4est_cubed_sphere_horography.jl new file mode 100644 index 00000000..fc0e7271 --- /dev/null +++ b/src/meshes/p4est_cubed_sphere_horography.jl @@ -0,0 +1,156 @@ +using Trixi: connectivity_cubed_sphere, new_p4est +abstract type AdaptVerticalGrid end + +struct GalChen <: AdaptVerticalGrid end + +struct Sleve{RealT} <: AdaptVerticalGrid + etaH::RealT + s::RealT +end + +@inline function (adapt_galchen::GalChen)(zRef, zs, H) + z = zRef + (H - zRef) * zs / H + return z + end + + @inline function (adapt_sleve::Sleve)(zRef, zs, H) + (; s, etaH) = Sleve + eta = zRef / H + if eta <= etaH + z = eta * H + zs * sinh((etaH - eta) / s / etaH) / sinh(one(zRef) / s) + else + z = eta * H + end + return z + end + +""" + P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; + polydeg, RealT=Float64, + initial_refinement_level=0, unsaved_changes=true, + p4est_partition_allow_for_coarsening=true) + +Build a "Cubed Sphere" mesh as `P4estMesh` with +`6 * trees_per_face_dimension^2 * layers` trees. + +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. +""" +function P4estMeshCubedSphereHorography(trees_per_face_dimension, layers, inner_radius, thickness; + polydeg, RealT = Float64, + initial_refinement_level = 0, unsaved_changes = true, + p4est_partition_allow_for_coarsening = true, horography = nothing, delta_horography = nothing, 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, horography, delta_horography, 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 horography_height(x_cart, y_cart, z_cart, horography, delta_horography) + RealT = eltype(x_cart) + r = sqrt(x_cart^2 + y_cart^2 + z_cart^2) + lat = asin(z_cart / r) + phi = atan(y_cart, x_cart) + + i = round(Int, (phi + pi) / delta_horography + 1) + j = round(Int, (lat + 0.5 * pi) / delta_horography + 1) + + i = clamp(i, 1, size(horography, 1)) + j = clamp(j, 1, size(horography, 2)) + + return Float64(horography[i, j]) +end + +function cubed_sphere_mapping_topo(xi, eta, zeta, inner_radius, thickness, direction, + horography, delta_horography, 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 + zs = horography_height(unit_pt[1], unit_pt[2], unit_pt[3], horography, delta_horography) + + zRef = thickness * 0.5 * (zeta + 1) + z = adapt_vertical_grid(zRef, zs, 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, horography, delta_horography, 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_topo( + x_offset + dx / 2 * nodes[i], + y_offset + dy / 2 * nodes[j], + z_offset + dz / 2 * nodes[k], + inner_radius, thickness, direction, + horography, delta_horography, adapt_vertical_grid) + end + end + end + + return nothing +end From b81a343010023e57020a28c6fdbf6a363b0abff3 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 5 May 2026 21:28:17 +0200 Subject: [PATCH 02/11] rename to topography, add vortex shedding with gaussian hill --- ...emperature_baroclinic_instability_earth.jl | 28 +- ...r_potential_temperature_vortex_shedding.jl | 269 ++++++++++++++++++ src/meshes/meshes.jl | 2 +- ...hy.jl => p4est_cubed_sphere_topography.jl} | 105 +++---- 4 files changed, 349 insertions(+), 55 deletions(-) create mode 100644 examples/euler/dry_air/global_circulation/elixir_potential_temperature_vortex_shedding.jl rename src/meshes/{p4est_cubed_sphere_horography.jl => p4est_cubed_sphere_topography.jl} (53%) diff --git a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl index a48a8a6e..b488ffd3 100644 --- a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl +++ b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl @@ -208,7 +208,9 @@ volume_flux = (flux_tec, flux_nonconservative_souza_etal) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + trees_per_cube_face = (8, 4) + filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_bed.nc" topo = ncread(filename, "z") @@ -225,10 +227,30 @@ topo = makepositive.(topo) topo = [topo; topo[1:1, :]] topo = [topo topo[:, 1:1]] -delta_topo = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians +function initial_topography_earth(x_cart, y_cart, z_cart, earth_topography) + RealT = eltype(x_cart) + delta_topography = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians + r = sqrt(x_cart^2 + y_cart^2 + z_cart^2) + lat = asin(z_cart / r) + phi = atan(y_cart, x_cart) + + i = round(Int, (phi + pi) / delta_topography + 1) + j = round(Int, (lat + 0.5 * pi) / delta_topography + 1) + + i = clamp(i, 1, size(earth_topography, 1)) + j = clamp(j, 1, size(earth_topography, 2)) + + return RealT(earth_topography[i, j]) +end -mesh = TrixiAtmo.P4estMeshCubedSphereHorography(trees_per_cube_face..., EARTH_RADIUS, 30000, - polydeg = polydeg, initial_refinement_level = 0, horography = topo, delta_horography = delta_topo, adapt_vertical_grid = TrixiAtmo.GalChen()) +mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000, + polydeg = polydeg, + initial_refinement_level = 0, + initial_topography = (x, y, z) -> initial_topography_earth(x, + y, + z, + topo), + adapt_vertical_grid = TrixiAtmo.GalChen()) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_baroclinic_instability, 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..081a1fde --- /dev/null +++ b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_vortex_shedding.jl @@ -0,0 +1,269 @@ +# An idealized baroclinic instability test case +# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. +# +# This elixir takes about 8 hours, using 16 threads of an AMD Ryzen 7 7800X3D. +# +# References: +# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) +# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores +# https://doi.org/10.1002/qj.2241 + +using OrdinaryDiffEqSSPRK +using Trixi, TrixiAtmo +using LinearAlgebra: norm +using NetCDF + +# Unperturbed balanced steady-state. +# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). +# The other velocity components are zero. +function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + radius_earth = EARTH_RADIUS # a + half_width_parameter = 2 # b + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g + k = 3 # k + surface_pressure = 1.0f5 # p₀ + gas_constant = 287 # R + surface_equatorial_temperature = 310 # T₀ᴱ + surface_polar_temperature = 240 # T₀ᴾ + lapse_rate = 0.005 # Γ + angular_velocity = EARTH_ROTATION_RATE # Ω + + # Distance to the center of the Earth + r = z + radius_earth + + # In the paper: T₀ + temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) + # In the paper: A, B, C, H + const_a = 1 / lapse_rate + const_b = (temperature0 - surface_polar_temperature) / + (temperature0 * surface_polar_temperature) + const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / + (surface_equatorial_temperature * surface_polar_temperature) + const_h = gas_constant * temperature0 / gravitational_acceleration + + # In the paper: (r - a) / bH + scaled_z = z / (half_width_parameter * const_h) + + # Temporary variables + temp1 = exp(lapse_rate / temperature0 * z) + temp2 = exp(-scaled_z^2) + + # In the paper: ̃τ₁, ̃τ₂ + tau1 = const_a * lapse_rate / temperature0 * temp1 + + const_b * (1 - 2 * scaled_z^2) * temp2 + tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 + + # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' + inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 + inttau2 = const_c * z * temp2 + + # Temporary variables + temp3 = r / radius_earth * cos(lat) + temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) + + # In the paper: T + temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) + + # In the paper: U, u (zonal wind, first component of spherical velocity) + big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * + (temp3^(k - 1) - temp3^(k + 1)) + temp5 = radius_earth * cos(lat) + u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) + + # Hydrostatic pressure + p = surface_pressure * + exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) + + # Density (via ideal gas law) + rho = p / (gas_constant * temperature) + + return rho, u, p +end + +# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) +function perturbation_stream_function(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + perturbation_radius = 1 / 6 # d₀ / a + perturbed_wind_amplitude = 1 # Vₚ + perturbation_lon = pi / 9 # Longitude of perturbation location + perturbation_lat = 2 * pi / 9 # Latitude of perturbation location + pertz = 15000 # Perturbation height cap + + # Great circle distance (d in the paper) divided by a (radius of the Earth) + # because we never actually need d without dividing by a + great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + + cos(perturbation_lat) * cos(lat) * + cos(lon - perturbation_lon)) + + # In the first case, the vertical taper function is per definition zero. + # In the second case, the stream function is per definition zero. + if z > pertz || great_circle_distance_by_a > perturbation_radius + return 0, 0 + end + + # Vertical tapering of stream function + perttaper = 1 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 + + # sin/cos(pi * d / (2 * d_0)) in the paper + sin_, cos_ = sincos(0.5f0 * pi * great_circle_distance_by_a / perturbation_radius) + + # Common factor for both u and v + factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ + + u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + + cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / + sin(great_circle_distance_by_a) + + v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / + sin(great_circle_distance_by_a) + + return u_perturbation, v_perturbation +end + +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 + +# Initial condition for an idealized baroclinic instability test +# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A +function initial_condition_baroclinic_instability(x, t, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + lon, lat, r = cartesian_to_sphere(x) + radius_earth = EARTH_RADIUS + # Make sure that the r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + + # Unperturbed basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Stream function type perturbation + u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) + + u += u_perturbation + v = v_perturbation + + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u - sin(lat) * cos(lon) * v + v2 = cos(lon) * u - sin(lat) * sin(lon) * v + v3 = cos(lat) * v + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0) + if z > 0 + r = norm(x) + else + r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) + end + r = -norm(x) + phi = radius_earth^2 * gravitational_acceleration / r + + return prim2cons(SVector(rho, v1, v2, v3, p, phi), equations) +end + +@inline function source_terms_baroclinic_instability(u, x, t, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + radius_earth = EARTH_RADIUS # a + angular_velocity = EARTH_ROTATION_RATE # Ω + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0) + r = z + radius_earth + + 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_baroclinic_instability + +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)) + +trees_per_cube_face = (8, 4) + +function initial_topography_gaussian(x, y, z) + h_0 = 2000 + d = 12500 + lambda_c = pi + phi_c = -20 * pi / 180 + + r_earth = sqrt(x^2 + y^2 + z^2) + lat = asin(z / r_earth) + lon = atan(y, x) + + 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 + +mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000, + polydeg = polydeg, + initial_refinement_level = 0, + initial_topography = initial_topography_gaussian, + adapt_vertical_grid = TrixiAtmo.GalChen()) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_baroclinic_instability, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. +T = 0.01 # 10 days +tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 10 days + +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 = 100, 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/meshes/meshes.jl b/src/meshes/meshes.jl index 198d6858..8ab7d4c4 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -1,4 +1,4 @@ include("p4est_cubed_sphere.jl") include("p4est_icosahedron_quads.jl") include("dgmulti_icosahedron_tri.jl") -include("p4est_cubed_sphere_horography.jl") +include("p4est_cubed_sphere_topography.jl") diff --git a/src/meshes/p4est_cubed_sphere_horography.jl b/src/meshes/p4est_cubed_sphere_topography.jl similarity index 53% rename from src/meshes/p4est_cubed_sphere_horography.jl rename to src/meshes/p4est_cubed_sphere_topography.jl index fc0e7271..0b7f0b20 100644 --- a/src/meshes/p4est_cubed_sphere_horography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -1,4 +1,5 @@ using Trixi: connectivity_cubed_sphere, new_p4est + abstract type AdaptVerticalGrid end struct GalChen <: AdaptVerticalGrid end @@ -8,24 +9,25 @@ struct Sleve{RealT} <: AdaptVerticalGrid s::RealT end -@inline function (adapt_galchen::GalChen)(zRef, zs, H) - z = zRef + (H - zRef) * zs / H - return z - end +@inline function (adapt_galchen::GalChen)(z_reference, z_topography, H) + z = z_reference + (H - z_reference) * z_topography / H + return z +end - @inline function (adapt_sleve::Sleve)(zRef, zs, H) +@inline function (adapt_sleve::Sleve)(z_reference, z_topography, H) (; s, etaH) = Sleve - eta = zRef / H - if eta <= etaH - z = eta * H + zs * sinh((etaH - eta) / s / etaH) / sinh(one(zRef) / s) - else - z = eta * H - end - return z + eta = z_reference / H + if eta <= etaH + z = eta * H + + z_topography * sinh((etaH - eta) / s / etaH) / sinh(one(z_reference) / s) + else + z = eta * H end + return z +end """ - P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness; + 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) @@ -52,10 +54,13 @@ The mesh will have two boundaries, `:inside` and `:outside`. independent of domain partitioning. Should be `false` for static meshes to permit more fine-grained partitioning. """ -function P4estMeshCubedSphereHorography(trees_per_face_dimension, layers, inner_radius, thickness; - polydeg, RealT = Float64, - initial_refinement_level = 0, unsaved_changes = true, - p4est_partition_allow_for_coarsening = true, horography = nothing, delta_horography = nothing, adapt_vertical_grid = GalChen()) +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 @@ -68,7 +73,8 @@ function P4estMeshCubedSphereHorography(trees_per_face_dimension, layers, inner_ n_trees) calc_tree_node_coordinates!(tree_node_coordinates, nodes, trees_per_face_dimension, layers, - inner_radius, thickness, horography, delta_horography, adapt_vertical_grid) + inner_radius, thickness, initial_topography, + adapt_vertical_grid) p4est = new_p4est(connectivity, initial_refinement_level) @@ -81,49 +87,35 @@ function P4estMeshCubedSphereHorography(trees_per_face_dimension, layers, inner_ p4est_partition_allow_for_coarsening) end -function horography_height(x_cart, y_cart, z_cart, horography, delta_horography) - RealT = eltype(x_cart) - r = sqrt(x_cart^2 + y_cart^2 + z_cart^2) - lat = asin(z_cart / r) - phi = atan(y_cart, x_cart) - - i = round(Int, (phi + pi) / delta_horography + 1) - j = round(Int, (lat + 0.5 * pi) / delta_horography + 1) - - i = clamp(i, 1, size(horography, 1)) - j = clamp(j, 1, size(horography, 2)) - - return Float64(horography[i, j]) -end - -function cubed_sphere_mapping_topo(xi, eta, zeta, inner_radius, thickness, direction, - horography, delta_horography, adapt_vertical_grid) +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 + 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)) + 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 - zs = horography_height(unit_pt[1], unit_pt[2], unit_pt[3], horography, delta_horography) + z_topography = initial_topography(unit_pt[1], unit_pt[2], unit_pt[3]) - zRef = thickness * 0.5 * (zeta + 1) - z = adapt_vertical_grid(zRef, zs, thickness) + 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, horography, delta_horography, adapt_vertical_grid) + inner_radius, thickness, initial_topography, + adapt_vertical_grid) n_cells_x = n_cells_y = trees_per_face_dimension n_cells_z = layers @@ -142,12 +134,23 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, 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_topo( - x_offset + dx / 2 * nodes[i], - y_offset + dy / 2 * nodes[j], - z_offset + dz / 2 * nodes[k], - inner_radius, thickness, direction, - horography, delta_horography, adapt_vertical_grid) + 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 From 73336a7209563a262acaae84607722b03dcaa1f9 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 5 May 2026 22:02:01 +0200 Subject: [PATCH 03/11] add vortex shedding --- ...r_potential_temperature_vortex_shedding.jl | 169 +++--------------- 1 file changed, 21 insertions(+), 148 deletions(-) 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 index 081a1fde..2dcf3b16 100644 --- 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 @@ -1,127 +1,6 @@ -# An idealized baroclinic instability test case -# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. -# -# This elixir takes about 8 hours, using 16 threads of an AMD Ryzen 7 7800X3D. -# -# References: -# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) -# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores -# https://doi.org/10.1002/qj.2241 - using OrdinaryDiffEqSSPRK using Trixi, TrixiAtmo using LinearAlgebra: norm -using NetCDF - -# Unperturbed balanced steady-state. -# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). -# The other velocity components are zero. -function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) - # Parameters from Table 1 in the paper - # Corresponding names in the paper are commented - radius_earth = EARTH_RADIUS # a - half_width_parameter = 2 # b - gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g - k = 3 # k - surface_pressure = 1.0f5 # p₀ - gas_constant = 287 # R - surface_equatorial_temperature = 310 # T₀ᴱ - surface_polar_temperature = 240 # T₀ᴾ - lapse_rate = 0.005 # Γ - angular_velocity = EARTH_ROTATION_RATE # Ω - - # Distance to the center of the Earth - r = z + radius_earth - - # In the paper: T₀ - temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) - # In the paper: A, B, C, H - const_a = 1 / lapse_rate - const_b = (temperature0 - surface_polar_temperature) / - (temperature0 * surface_polar_temperature) - const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / - (surface_equatorial_temperature * surface_polar_temperature) - const_h = gas_constant * temperature0 / gravitational_acceleration - - # In the paper: (r - a) / bH - scaled_z = z / (half_width_parameter * const_h) - - # Temporary variables - temp1 = exp(lapse_rate / temperature0 * z) - temp2 = exp(-scaled_z^2) - - # In the paper: ̃τ₁, ̃τ₂ - tau1 = const_a * lapse_rate / temperature0 * temp1 + - const_b * (1 - 2 * scaled_z^2) * temp2 - tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 - - # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' - inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 - inttau2 = const_c * z * temp2 - - # Temporary variables - temp3 = r / radius_earth * cos(lat) - temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) - - # In the paper: T - temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) - - # In the paper: U, u (zonal wind, first component of spherical velocity) - big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * - (temp3^(k - 1) - temp3^(k + 1)) - temp5 = radius_earth * cos(lat) - u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) - - # Hydrostatic pressure - p = surface_pressure * - exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) - - # Density (via ideal gas law) - rho = p / (gas_constant * temperature) - - return rho, u, p -end - -# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) -function perturbation_stream_function(lon, lat, z) - # Parameters from Table 1 in the paper - # Corresponding names in the paper are commented - perturbation_radius = 1 / 6 # d₀ / a - perturbed_wind_amplitude = 1 # Vₚ - perturbation_lon = pi / 9 # Longitude of perturbation location - perturbation_lat = 2 * pi / 9 # Latitude of perturbation location - pertz = 15000 # Perturbation height cap - - # Great circle distance (d in the paper) divided by a (radius of the Earth) - # because we never actually need d without dividing by a - great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + - cos(perturbation_lat) * cos(lat) * - cos(lon - perturbation_lon)) - - # In the first case, the vertical taper function is per definition zero. - # In the second case, the stream function is per definition zero. - if z > pertz || great_circle_distance_by_a > perturbation_radius - return 0, 0 - end - - # Vertical tapering of stream function - perttaper = 1 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 - - # sin/cos(pi * d / (2 * d_0)) in the paper - sin_, cos_ = sincos(0.5f0 * pi * great_circle_distance_by_a / perturbation_radius) - - # Common factor for both u and v - factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ - - u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + - cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / - sin(great_circle_distance_by_a) - - v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / - sin(great_circle_distance_by_a) - - return u_perturbation, v_perturbation -end function cartesian_to_sphere(x) r = norm(x) @@ -134,30 +13,22 @@ function cartesian_to_sphere(x) return lambda, phi, r end -# Initial condition for an idealized baroclinic instability test -# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A -function initial_condition_baroclinic_instability(x, t, +function initial_condition_vortex_shedding(x, t, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) lon, lat, r = cartesian_to_sphere(x) - radius_earth = EARTH_RADIUS + radius_earth = EARTH_RADIUS/20 # Make sure that the r is not smaller than radius_earth z = max(r - radius_earth, 0.0) - - # Unperturbed basic state - rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) - - # Stream function type perturbation - u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) - - u += u_perturbation - v = v_perturbation - + R = equations.c_p - equations.c_v + k = R/equations.c_p # Convert spherical velocity to Cartesian - v1 = -sin(lon) * u - sin(lat) * cos(lon) * v - v2 = cos(lon) * u - sin(lat) * sin(lon) * v - v3 = cos(lat) * v - gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g - + 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) @@ -167,15 +38,17 @@ function initial_condition_baroclinic_instability(x, t, r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) end r = -norm(x) - phi = radius_earth^2 * gravitational_acceleration / r - + phi = radius_earth^2 * g/ r + N = 0.0182 + 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_baroclinic_instability(u, x, t, +@inline function source_terms_coriolis(u, x, t, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) - radius_earth = EARTH_RADIUS # a - angular_velocity = EARTH_ROTATION_RATE # Ω + radius_earth = EARTH_RADIUS/20 # a + angular_velocity = EARTH_ROTATION_RATE*20 # Ω r = norm(x) # Make sure that r is not smaller than radius_earth @@ -197,7 +70,7 @@ equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 10 c_v = 717, gravity = EARTH_GRAVITATIONAL_ACCELERATION) -initial_condition = initial_condition_baroclinic_instability +initial_condition = initial_condition_vortex_shedding boundary_conditions = (; inside = boundary_condition_slip_wall, outside = boundary_condition_slip_wall) @@ -227,14 +100,14 @@ function initial_topography_gaussian(x, y, z) return h_0 * exp(-(dist / d)^2) end -mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000, +mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS/20, 20000, polydeg = polydeg, initial_refinement_level = 0, initial_topography = initial_topography_gaussian, adapt_vertical_grid = TrixiAtmo.GalChen()) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms = source_terms_baroclinic_instability, + source_terms = source_terms_coriolis, boundary_conditions = boundary_conditions) ############################################################################### From 0b0828ac39416276d917b8faefe7ec101ea8ae52 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 10:13:49 +0200 Subject: [PATCH 04/11] fix mapping for the topography --- .../elixir_potential_temperature_vortex_shedding.jl | 9 +++++---- src/meshes/p4est_cubed_sphere_topography.jl | 2 ++ 2 files changed, 7 insertions(+), 4 deletions(-) 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 index 2dcf3b16..0f155a37 100644 --- 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 @@ -82,21 +82,22 @@ volume_flux = (flux_tec, flux_nonconservative_souza_etal) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) -trees_per_cube_face = (8, 4) +trees_per_cube_face = (16, 8) function initial_topography_gaussian(x, y, z) h_0 = 2000 + h0 = 10000 d = 12500 lambda_c = pi - phi_c = -20 * pi / 180 + phi_c = 20 * pi / 180 r_earth = sqrt(x^2 + y^2 + z^2) lat = asin(z / r_earth) lon = atan(y, x) + 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 @@ -124,7 +125,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(dt = 100, save_initial_solution = true, +save_solution = SaveSolutionCallback(dt = 5760, save_initial_solution = true, save_final_solution = true) callbacks = CallbackSet(summary_callback, diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 0b7f0b20..7a91c549 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -53,6 +53,8 @@ The mesh will have two boundaries, `:inside` and `:outside`. - `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 r estore the spherical shape. Two options are available: GalChen() or Sleve(etaH, s). """ function P4estMeshCubedSphereTopography(trees_per_face_dimension, layers, inner_radius, thickness; From 7750d3cd852e3b430a1e320b5fdc486ddfadcbbb Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 12:28:29 +0200 Subject: [PATCH 05/11] topography as function of longitude and latitute, small adjustments --- ...emperature_baroclinic_instability_earth.jl | 33 ++++---- ...r_potential_temperature_vortex_shedding.jl | 36 ++++----- src/meshes/p4est_cubed_sphere_topography.jl | 80 +++++++++++++++---- 3 files changed, 102 insertions(+), 47 deletions(-) diff --git a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl index b488ffd3..02278afb 100644 --- a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl +++ b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl @@ -211,8 +211,14 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, trees_per_cube_face = (8, 4) -filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_bed.nc" -topo = ncread(filename, "z") +filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_surface.nc" +#= +mesh_file = Trixi.download( + "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf/ETOPO_2022_v1_60s_N90W180_surface.nc", + joinpath(@__DIR__, "ETOPO_2022_v1_60s_N90W180_surface.nc") +) +=# +earth_topography = ncread(filename, "z") function makepositive(x) if x >= 0 @@ -222,19 +228,16 @@ function makepositive(x) end end -topo = makepositive.(topo) +earth_topography = makepositive.(earth_topography) # Add extra column and row for periodicity -topo = [topo; topo[1:1, :]] -topo = [topo topo[:, 1:1]] +earth_topography = [earth_topography; earth_topography[1:1, :]] +earth_topography = [earth_topography earth_topography[:, 1:1]] -function initial_topography_earth(x_cart, y_cart, z_cart, earth_topography) +function initial_topography_earth(lat, lon, earth_topography) RealT = eltype(x_cart) delta_topography = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians - r = sqrt(x_cart^2 + y_cart^2 + z_cart^2) - lat = asin(z_cart / r) - phi = atan(y_cart, x_cart) - i = round(Int, (phi + pi) / delta_topography + 1) + i = round(Int, (lon + pi) / delta_topography + 1) j = round(Int, (lat + 0.5 * pi) / delta_topography + 1) i = clamp(i, 1, size(earth_topography, 1)) @@ -246,11 +249,11 @@ end mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000, polydeg = polydeg, initial_refinement_level = 0, - initial_topography = (x, y, z) -> initial_topography_earth(x, - y, - z, - topo), - adapt_vertical_grid = TrixiAtmo.GalChen()) + initial_topography = (lat, lon) -> initial_topography_earth(lat, + lon, + earth_topography), + adapt_vertical_grid = TrixiAtmo.Sleve(0.7, + 0.8)) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_baroclinic_instability, 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 index 0f155a37..921352f3 100644 --- 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 @@ -14,21 +14,21 @@ function cartesian_to_sphere(x) end function initial_condition_vortex_shedding(x, t, - equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) lon, lat, r = cartesian_to_sphere(x) - radius_earth = EARTH_RADIUS/20 + radius_earth = EARTH_RADIUS / 20 # Make sure that the r is not smaller than radius_earth z = max(r - radius_earth, 0.0) R = equations.c_p - equations.c_v - k = R/equations.c_p + k = R / equations.c_p # Convert spherical velocity to Cartesian u0 = 20 - v1 = -u0*cos(lat) * sin(lon) + 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 # Ω + 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) @@ -38,17 +38,19 @@ function initial_condition_vortex_shedding(x, t, r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) end r = -norm(x) - phi = radius_earth^2 * g/ r + phi = radius_earth^2 * g / r N = 0.0182 - 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) + 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 # Ω + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + radius_earth = EARTH_RADIUS / 20 # a + angular_velocity = EARTH_ROTATION_RATE * 20 # Ω r = norm(x) # Make sure that r is not smaller than radius_earth @@ -84,24 +86,22 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, trees_per_cube_face = (16, 8) -function initial_topography_gaussian(x, y, z) +function initial_topography_gaussian(lat, lon) h_0 = 2000 h0 = 10000 d = 12500 lambda_c = pi phi_c = 20 * pi / 180 - r_earth = sqrt(x^2 + y^2 + z^2) - lat = asin(z / r_earth) - lon = atan(y, x) - - r_earth = EARTH_RADIUS/20 + 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 -mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS/20, 20000, +mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS / 20, + 20000, polydeg = polydeg, initial_refinement_level = 0, initial_topography = initial_topography_gaussian, diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 7a91c549..df8880f0 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -1,31 +1,81 @@ using Trixi: connectivity_cubed_sphere, new_p4est -abstract type AdaptVerticalGrid end +abstract type SmoothVerticalCoordinate end -struct GalChen <: AdaptVerticalGrid end +""" + Sleve(decay_start_fraction, decay_scale) -struct Sleve{RealT} <: AdaptVerticalGrid - etaH::RealT - s::RealT -end +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. -@inline function (adapt_galchen::GalChen)(z_reference, z_topography, H) - z = z_reference + (H - z_reference) * z_topography / H - return z +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) - (; s, etaH) = Sleve + (; s, eta_H) = Sleve eta = z_reference / H if eta <= etaH z = eta * H + - z_topography * sinh((etaH - eta) / s / etaH) / sinh(one(z_reference) / s) + 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, @@ -53,8 +103,8 @@ The mesh will have two boundaries, `:inside` and `:outside`. - `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 r estore the spherical shape. Two options are available: GalChen() or Sleve(etaH, s). +- `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; @@ -107,7 +157,9 @@ function cubed_sphere_mapping_topography(xi, eta, zeta, inner_radius, thickness, r = sqrt(1 + x^2 + y^2) unit_pt = cube_coordinates[direction] / r - z_topography = initial_topography(unit_pt[1], unit_pt[2], unit_pt[3]) + 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) From f98459ef0088b4aeb77791b7f43960c1fe91d933 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 14:32:54 +0200 Subject: [PATCH 06/11] add tests --- ...r_potential_temperature_vortex_shedding.jl | 4 +-- test/test_3d_potential_temperature.jl | 27 +++++++++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) 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 index 921352f3..1138df60 100644 --- 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 @@ -113,8 +113,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks etc. -T = 0.01 # 10 days -tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 10 days +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) diff --git a/test/test_3d_potential_temperature.jl b/test/test_3d_potential_temperature.jl index 7b4af7e2..7e6f27e5 100644 --- a/test/test_3d_potential_temperature.jl +++ b/test/test_3d_potential_temperature.jl @@ -210,4 +210,31 @@ 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=[ + 2.496983981196748e12, + 8.709311184559809e14, + 8.708527613659252e14, + 1.0364047379154301e15, + 1.5716246388841625e10, + 4.804110544770574 + ], + linf=[ + 2.6402826399432e13, + 3.77696630836204e15, + 3.6352451849496805e15, + 5.624197001818547e15, + 1.6198347751471875e11, + 16.685684057418257 + ], + 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 + end From 61a1d155a68f7b47b476b591bb8921de31994199 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 14:51:43 +0200 Subject: [PATCH 07/11] add tests, fix small bugs --- ...emperature_baroclinic_instability_earth.jl | 291 ------------------ ...r_potential_temperature_vortex_shedding.jl | 47 ++- src/TrixiAtmo.jl | 2 + src/meshes/p4est_cubed_sphere_topography.jl | 6 +- test/test_3d_potential_temperature.jl | 51 ++- 5 files changed, 67 insertions(+), 330 deletions(-) delete mode 100644 examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl diff --git a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl b/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl deleted file mode 100644 index 02278afb..00000000 --- a/examples/euler/dry_air/global_circulation/elixir_potential_temperature_baroclinic_instability_earth.jl +++ /dev/null @@ -1,291 +0,0 @@ -# An idealized baroclinic instability test case -# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. -# -# This elixir takes about 8 hours, using 16 threads of an AMD Ryzen 7 7800X3D. -# -# References: -# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) -# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores -# https://doi.org/10.1002/qj.2241 - -using OrdinaryDiffEqSSPRK -using Trixi, TrixiAtmo -using LinearAlgebra: norm -using NetCDF - -# Unperturbed balanced steady-state. -# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). -# The other velocity components are zero. -function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) - # Parameters from Table 1 in the paper - # Corresponding names in the paper are commented - radius_earth = EARTH_RADIUS # a - half_width_parameter = 2 # b - gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g - k = 3 # k - surface_pressure = 1.0f5 # p₀ - gas_constant = 287 # R - surface_equatorial_temperature = 310 # T₀ᴱ - surface_polar_temperature = 240 # T₀ᴾ - lapse_rate = 0.005 # Γ - angular_velocity = EARTH_ROTATION_RATE # Ω - - # Distance to the center of the Earth - r = z + radius_earth - - # In the paper: T₀ - temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) - # In the paper: A, B, C, H - const_a = 1 / lapse_rate - const_b = (temperature0 - surface_polar_temperature) / - (temperature0 * surface_polar_temperature) - const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / - (surface_equatorial_temperature * surface_polar_temperature) - const_h = gas_constant * temperature0 / gravitational_acceleration - - # In the paper: (r - a) / bH - scaled_z = z / (half_width_parameter * const_h) - - # Temporary variables - temp1 = exp(lapse_rate / temperature0 * z) - temp2 = exp(-scaled_z^2) - - # In the paper: ̃τ₁, ̃τ₂ - tau1 = const_a * lapse_rate / temperature0 * temp1 + - const_b * (1 - 2 * scaled_z^2) * temp2 - tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 - - # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' - inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 - inttau2 = const_c * z * temp2 - - # Temporary variables - temp3 = r / radius_earth * cos(lat) - temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) - - # In the paper: T - temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) - - # In the paper: U, u (zonal wind, first component of spherical velocity) - big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * - (temp3^(k - 1) - temp3^(k + 1)) - temp5 = radius_earth * cos(lat) - u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) - - # Hydrostatic pressure - p = surface_pressure * - exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) - - # Density (via ideal gas law) - rho = p / (gas_constant * temperature) - - return rho, u, p -end - -# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) -function perturbation_stream_function(lon, lat, z) - # Parameters from Table 1 in the paper - # Corresponding names in the paper are commented - perturbation_radius = 1 / 6 # d₀ / a - perturbed_wind_amplitude = 1 # Vₚ - perturbation_lon = pi / 9 # Longitude of perturbation location - perturbation_lat = 2 * pi / 9 # Latitude of perturbation location - pertz = 15000 # Perturbation height cap - - # Great circle distance (d in the paper) divided by a (radius of the Earth) - # because we never actually need d without dividing by a - great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + - cos(perturbation_lat) * cos(lat) * - cos(lon - perturbation_lon)) - - # In the first case, the vertical taper function is per definition zero. - # In the second case, the stream function is per definition zero. - if z > pertz || great_circle_distance_by_a > perturbation_radius - return 0, 0 - end - - # Vertical tapering of stream function - perttaper = 1 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 - - # sin/cos(pi * d / (2 * d_0)) in the paper - sin_, cos_ = sincos(0.5f0 * pi * great_circle_distance_by_a / perturbation_radius) - - # Common factor for both u and v - factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ - - u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + - cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / - sin(great_circle_distance_by_a) - - v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / - sin(great_circle_distance_by_a) - - return u_perturbation, v_perturbation -end - -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 - -# Initial condition for an idealized baroclinic instability test -# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A -function initial_condition_baroclinic_instability(x, t, - equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) - lon, lat, r = cartesian_to_sphere(x) - radius_earth = EARTH_RADIUS - # Make sure that the r is not smaller than radius_earth - z = max(r - radius_earth, 0.0) - - # Unperturbed basic state - rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) - - # Stream function type perturbation - u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) - - u += u_perturbation - v = v_perturbation - - # Convert spherical velocity to Cartesian - v1 = -sin(lon) * u - sin(lat) * cos(lon) * v - v2 = cos(lon) * u - sin(lat) * sin(lon) * v - v3 = cos(lat) * v - gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g - - r = norm(x) - # Make sure that r is not smaller than radius_earth - z = max(r - radius_earth, 0) - if z > 0 - r = norm(x) - else - r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) - end - r = -norm(x) - phi = radius_earth^2 * gravitational_acceleration / r - - return prim2cons(SVector(rho, v1, v2, v3, p, phi), equations) -end - -@inline function source_terms_baroclinic_instability(u, x, t, - equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) - radius_earth = EARTH_RADIUS # a - angular_velocity = EARTH_ROTATION_RATE # Ω - - r = norm(x) - # Make sure that r is not smaller than radius_earth - z = max(r - radius_earth, 0) - r = z + radius_earth - - 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_baroclinic_instability - -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)) - -trees_per_cube_face = (8, 4) - -filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_surface.nc" -#= -mesh_file = Trixi.download( - "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf/ETOPO_2022_v1_60s_N90W180_surface.nc", - joinpath(@__DIR__, "ETOPO_2022_v1_60s_N90W180_surface.nc") -) -=# -earth_topography = ncread(filename, "z") - -function makepositive(x) - if x >= 0 - return x - else - return zero(typeof(x)) - end -end - -earth_topography = makepositive.(earth_topography) -# Add extra column and row for periodicity -earth_topography = [earth_topography; earth_topography[1:1, :]] -earth_topography = [earth_topography earth_topography[:, 1:1]] - -function initial_topography_earth(lat, lon, earth_topography) - RealT = eltype(x_cart) - delta_topography = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians - - i = round(Int, (lon + pi) / delta_topography + 1) - j = round(Int, (lat + 0.5 * pi) / delta_topography + 1) - - i = clamp(i, 1, size(earth_topography, 1)) - j = clamp(j, 1, size(earth_topography, 2)) - - return RealT(earth_topography[i, j]) -end - -mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000, - polydeg = polydeg, - initial_refinement_level = 0, - initial_topography = (lat, lon) -> initial_topography_earth(lat, - lon, - earth_topography), - adapt_vertical_grid = TrixiAtmo.Sleve(0.7, - 0.8)) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms = source_terms_baroclinic_instability, - boundary_conditions = boundary_conditions) - -############################################################################### -# ODE solvers, callbacks etc. -T = 0.01 # 10 days -tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 10 days - -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 = 100, 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/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 index 1138df60..8e53297d 100644 --- 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 @@ -1,3 +1,15 @@ +# 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 @@ -17,8 +29,6 @@ function initial_condition_vortex_shedding(x, t, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) lon, lat, r = cartesian_to_sphere(x) radius_earth = EARTH_RADIUS / 20 - # Make sure that the r is not smaller than radius_earth - z = max(r - radius_earth, 0.0) R = equations.c_p - equations.c_v k = R / equations.c_p # Convert spherical velocity to Cartesian @@ -31,15 +41,9 @@ function initial_condition_vortex_shedding(x, t, 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) - if z > 0 - r = norm(x) - else - r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) - end - r = -norm(x) - phi = radius_earth^2 * g / r - N = 0.0182 + 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) @@ -52,11 +56,6 @@ end radius_earth = EARTH_RADIUS / 20 # a 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) - r = z + radius_earth - du1 = zero(eltype(u)) du2 = zero(eltype(u)) du3 = zero(eltype(u)) @@ -84,8 +83,6 @@ volume_flux = (flux_tec, flux_nonconservative_souza_etal) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) -trees_per_cube_face = (16, 8) - function initial_topography_gaussian(lat, lon) h_0 = 2000 h0 = 10000 @@ -100,12 +97,14 @@ function initial_topography_gaussian(lat, lon) return h_0 * exp(-(dist / d)^2) end -mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS / 20, - 20000, - polydeg = polydeg, - initial_refinement_level = 0, - initial_topography = initial_topography_gaussian, - adapt_vertical_grid = TrixiAtmo.GalChen()) +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, 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/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index df8880f0..4158a9c1 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -37,9 +37,9 @@ struct Sleve{RealT} <: SmoothVerticalCoordinate end @inline function (adapt_sleve::Sleve)(z_reference, z_topography, H) - (; s, eta_H) = Sleve + @unpack s, eta_H = adapt_sleve eta = z_reference / H - if eta <= etaH + if eta <= eta_H z = eta * H + z_topography * sinh((eta_H - eta) / s / eta_H) / sinh(one(z_reference) / s) else @@ -80,7 +80,7 @@ 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) + 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. diff --git a/test/test_3d_potential_temperature.jl b/test/test_3d_potential_temperature.jl index 7e6f27e5..1b4ba8f5 100644 --- a/test/test_3d_potential_temperature.jl +++ b/test/test_3d_potential_temperature.jl @@ -214,20 +214,20 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "global_circulation", "elixir_potential_temperature_vortex_shedding.jl"), l2=[ - 2.496983981196748e12, - 8.709311184559809e14, - 8.708527613659252e14, - 1.0364047379154301e15, - 1.5716246388841625e10, - 4.804110544770574 + 0.00011127272731957073, + 0.036136807950192446, + 0.036132932940041926, + 0.047765679229021134, + 0.03287203491711795, + 5.112500661735035 ], linf=[ - 2.6402826399432e13, - 3.77696630836204e15, - 3.6352451849496805e15, - 5.624197001818547e15, - 1.6198347751471875e11, - 16.685684057418257 + 0.0011895976218974091, + 0.16161469307155263, + 0.15524551919165414, + 0.25696936007842397, + 0.34356951757109755, + 42.337304272491956 ], rtol=1e-9, tspan=(0.0, 0.0001 * SECONDS_PER_DAY), @@ -237,4 +237,31 @@ end @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 From cbe322d81f2f35f7bed00720f2b2e7bba9da0010 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 14:55:03 +0200 Subject: [PATCH 08/11] fix function' signature --- src/meshes/p4est_cubed_sphere_topography.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 4158a9c1..8b685a78 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -3,7 +3,7 @@ using Trixi: connectivity_cubed_sphere, new_p4est abstract type SmoothVerticalCoordinate end """ - Sleve(decay_start_fraction, decay_scale) + Sleve(eta_H, s) SLEVE (Smooth LEvel VErtical) coordinate transformation for terrain-following vertical grids in atmospheric models. Ensures that the topographic influence From d12fe2bd9ec0597f246c3a3e3308404e080269a7 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 17:15:44 +0200 Subject: [PATCH 09/11] fix indentation --- src/meshes/p4est_cubed_sphere_topography.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 8b685a78..69c1274a 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -19,10 +19,10 @@ 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`. + 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) From 1fa227c4ca30c0ed9d91f96f25a7c25c34dd1838 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 6 May 2026 17:25:04 +0200 Subject: [PATCH 10/11] fix indentation again --- src/meshes/p4est_cubed_sphere_topography.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 69c1274a..619eeb9c 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -19,8 +19,8 @@ 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`. + 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`. From 95e36358235d437306d9b861c7ad6a1aa80af452 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Thu, 7 May 2026 07:52:17 +0200 Subject: [PATCH 11/11] add a detail in the docs --- src/meshes/p4est_cubed_sphere_topography.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_cubed_sphere_topography.jl b/src/meshes/p4est_cubed_sphere_topography.jl index 619eeb9c..367c808d 100644 --- a/src/meshes/p4est_cubed_sphere_topography.jl +++ b/src/meshes/p4est_cubed_sphere_topography.jl @@ -83,7 +83,7 @@ end 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. +`6 * trees_per_face_dimension^2 * layers` trees and a given topography. The mesh will have two boundaries, `:inside` and `:outside`.