|
| 1 | +# An idealized baroclinic instability test case |
| 2 | +# For optimal results consider increasing the resolution to 16x16x8 trees per cube face. |
| 3 | +# |
| 4 | +# This elixir takes about 8 hours, using 16 threads of an AMD Ryzen 7 7800X3D. |
| 5 | +# |
| 6 | +# References: |
| 7 | +# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013) |
| 8 | +# A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores |
| 9 | +# https://doi.org/10.1002/qj.2241 |
| 10 | + |
| 11 | +using OrdinaryDiffEqSSPRK |
| 12 | +using Trixi, TrixiAtmo |
| 13 | +using LinearAlgebra: norm |
| 14 | + |
| 15 | +# Unperturbed balanced steady-state. |
| 16 | +# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). |
| 17 | +# The other velocity components are zero. |
| 18 | +function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) |
| 19 | + # Parameters from Table 1 in the paper |
| 20 | + # Corresponding names in the paper are commented |
| 21 | + radius_earth = 6.371229e6 # a |
| 22 | + half_width_parameter = 2 # b |
| 23 | + gravitational_acceleration = 9.81 # g |
| 24 | + k = 3 # k |
| 25 | + surface_pressure = 1.0f5 # p₀ |
| 26 | + gas_constant = 287 # R |
| 27 | + surface_equatorial_temperature = 310 # T₀ᴱ |
| 28 | + surface_polar_temperature = 240 # T₀ᴾ |
| 29 | + lapse_rate = 0.005 # Γ |
| 30 | + angular_velocity = 7.29212e-5 # Ω |
| 31 | + |
| 32 | + # Distance to the center of the Earth |
| 33 | + r = z + radius_earth |
| 34 | + |
| 35 | + # In the paper: T₀ |
| 36 | + temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) |
| 37 | + # In the paper: A, B, C, H |
| 38 | + const_a = 1 / lapse_rate |
| 39 | + const_b = (temperature0 - surface_polar_temperature) / |
| 40 | + (temperature0 * surface_polar_temperature) |
| 41 | + const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / |
| 42 | + (surface_equatorial_temperature * surface_polar_temperature) |
| 43 | + const_h = gas_constant * temperature0 / gravitational_acceleration |
| 44 | + |
| 45 | + # In the paper: (r - a) / bH |
| 46 | + scaled_z = z / (half_width_parameter * const_h) |
| 47 | + |
| 48 | + # Temporary variables |
| 49 | + temp1 = exp(lapse_rate / temperature0 * z) |
| 50 | + temp2 = exp(-scaled_z^2) |
| 51 | + |
| 52 | + # In the paper: ̃τ₁, ̃τ₂ |
| 53 | + tau1 = const_a * lapse_rate / temperature0 * temp1 + |
| 54 | + const_b * (1 - 2 * scaled_z^2) * temp2 |
| 55 | + tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 |
| 56 | + |
| 57 | + # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' |
| 58 | + inttau1 = const_a * (temp1 - 1) + const_b * z * temp2 |
| 59 | + inttau2 = const_c * z * temp2 |
| 60 | + |
| 61 | + # Temporary variables |
| 62 | + temp3 = r / radius_earth * cos(lat) |
| 63 | + temp4 = temp3^k - k / (k + 2) * temp3^(k + 2) |
| 64 | + |
| 65 | + # In the paper: T |
| 66 | + temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4)) |
| 67 | + |
| 68 | + # In the paper: U, u (zonal wind, first component of spherical velocity) |
| 69 | + big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 * |
| 70 | + (temp3^(k - 1) - temp3^(k + 1)) |
| 71 | + temp5 = radius_earth * cos(lat) |
| 72 | + u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) |
| 73 | + |
| 74 | + # Hydrostatic pressure |
| 75 | + p = surface_pressure * |
| 76 | + exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4)) |
| 77 | + |
| 78 | + # Density (via ideal gas law) |
| 79 | + rho = p / (gas_constant * temperature) |
| 80 | + |
| 81 | + return rho, u, p |
| 82 | +end |
| 83 | + |
| 84 | +# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) |
| 85 | +function perturbation_stream_function(lon, lat, z) |
| 86 | + # Parameters from Table 1 in the paper |
| 87 | + # Corresponding names in the paper are commented |
| 88 | + perturbation_radius = 1 / 6 # d₀ / a |
| 89 | + perturbed_wind_amplitude = 1 # Vₚ |
| 90 | + perturbation_lon = pi / 9 # Longitude of perturbation location |
| 91 | + perturbation_lat = 2 * pi / 9 # Latitude of perturbation location |
| 92 | + pertz = 15000 # Perturbation height cap |
| 93 | + |
| 94 | + # Great circle distance (d in the paper) divided by a (radius of the Earth) |
| 95 | + # because we never actually need d without dividing by a |
| 96 | + great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + |
| 97 | + cos(perturbation_lat) * cos(lat) * |
| 98 | + cos(lon - perturbation_lon)) |
| 99 | + |
| 100 | + # In the first case, the vertical taper function is per definition zero. |
| 101 | + # In the second case, the stream function is per definition zero. |
| 102 | + if z > pertz || great_circle_distance_by_a > perturbation_radius |
| 103 | + return 0, 0 |
| 104 | + end |
| 105 | + |
| 106 | + # Vertical tapering of stream function |
| 107 | + perttaper = 1 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 |
| 108 | + |
| 109 | + # sin/cos(pi * d / (2 * d_0)) in the paper |
| 110 | + sin_, cos_ = sincos(0.5f0 * pi * great_circle_distance_by_a / perturbation_radius) |
| 111 | + |
| 112 | + # Common factor for both u and v |
| 113 | + factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ |
| 114 | + |
| 115 | + u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + |
| 116 | + cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) / |
| 117 | + sin(great_circle_distance_by_a) |
| 118 | + |
| 119 | + v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / |
| 120 | + sin(great_circle_distance_by_a) |
| 121 | + |
| 122 | + return u_perturbation, v_perturbation |
| 123 | +end |
| 124 | + |
| 125 | +function cartesian_to_sphere(x) |
| 126 | + r = norm(x) |
| 127 | + lambda = atan(x[2], x[1]) |
| 128 | + if lambda < 0 |
| 129 | + lambda += 2 * pi |
| 130 | + end |
| 131 | + phi = asin(x[3] / r) |
| 132 | + |
| 133 | + return lambda, phi, r |
| 134 | +end |
| 135 | + |
| 136 | +# Initial condition for an idealized baroclinic instability test |
| 137 | +# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A |
| 138 | +function initial_condition_baroclinic_instability(x, t, |
| 139 | + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) |
| 140 | + lon, lat, r = cartesian_to_sphere(x) |
| 141 | + radius_earth = 6.371229e6 |
| 142 | + # Make sure that the r is not smaller than radius_earth |
| 143 | + z = max(r - radius_earth, 0.0) |
| 144 | + |
| 145 | + # Unperturbed basic state |
| 146 | + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) |
| 147 | + |
| 148 | + # Stream function type perturbation |
| 149 | + u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) |
| 150 | + |
| 151 | + u += u_perturbation |
| 152 | + v = v_perturbation |
| 153 | + |
| 154 | + # Convert spherical velocity to Cartesian |
| 155 | + v1 = -sin(lon) * u - sin(lat) * cos(lon) * v |
| 156 | + v2 = cos(lon) * u - sin(lat) * sin(lon) * v |
| 157 | + v3 = cos(lat) * v |
| 158 | + radius_earth = 6.371229e6 # a |
| 159 | + gravitational_acceleration = 9.81 # g |
| 160 | + |
| 161 | + r = norm(x) |
| 162 | + # Make sure that r is not smaller than radius_earth |
| 163 | + z = max(r - radius_earth, 0) |
| 164 | + if z > 0 |
| 165 | + r = norm(x) |
| 166 | + else |
| 167 | + r = -(2 * radius_earth^3) / (x[1]^2 + x[2]^2 + x[3]^2) |
| 168 | + end |
| 169 | + r = -norm(x) |
| 170 | + phi = radius_earth^2 * gravitational_acceleration / r |
| 171 | + |
| 172 | + return prim2cons(SVector(rho, v1, v2, v3, p, phi), equations) |
| 173 | +end |
| 174 | + |
| 175 | +@inline function source_terms_baroclinic_instability(u, x, t, |
| 176 | + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) |
| 177 | + radius_earth = 6.371229e6 # a |
| 178 | + angular_velocity = 7.29212e-5 # Ω |
| 179 | + |
| 180 | + r = norm(x) |
| 181 | + # Make sure that r is not smaller than radius_earth |
| 182 | + z = max(r - radius_earth, 0) |
| 183 | + r = z + radius_earth |
| 184 | + |
| 185 | + du1 = zero(eltype(u)) |
| 186 | + du2 = zero(eltype(u)) |
| 187 | + du3 = zero(eltype(u)) |
| 188 | + du4 = zero(eltype(u)) |
| 189 | + du5 = zero(eltype(u)) |
| 190 | + # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4] |
| 191 | + du2 -= -2 * angular_velocity * u[3] |
| 192 | + du3 -= 2 * angular_velocity * u[2] |
| 193 | + |
| 194 | + return SVector(du1, du2, du3, du4, du5, zero(eltype(u))) |
| 195 | +end |
| 196 | +equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D() |
| 197 | + |
| 198 | +initial_condition = initial_condition_baroclinic_instability |
| 199 | + |
| 200 | +boundary_conditions = Dict(:inside => boundary_condition_slip_wall, |
| 201 | + :outside => boundary_condition_slip_wall) |
| 202 | + |
| 203 | +polydeg = 5 |
| 204 | +surface_flux = (FluxLMARS(340), flux_zero) |
| 205 | +volume_flux = (flux_tec, flux_nonconservative_souza_etal) |
| 206 | + |
| 207 | +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, |
| 208 | + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) |
| 209 | +trees_per_cube_face = (8, 4) |
| 210 | +mesh = P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000, |
| 211 | + polydeg = polydeg, initial_refinement_level = 0) |
| 212 | + |
| 213 | +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, |
| 214 | + source_terms = source_terms_baroclinic_instability, |
| 215 | + boundary_conditions = boundary_conditions) |
| 216 | + |
| 217 | +############################################################################### |
| 218 | +# ODE solvers, callbacks etc. |
| 219 | +T = 10 # 10 days |
| 220 | +tspan = (0.0, T * SECONDS_PER_DAY) # time in seconds for 10 days |
| 221 | + |
| 222 | +ode = semidiscretize(semi, tspan) |
| 223 | + |
| 224 | +summary_callback = SummaryCallback() |
| 225 | + |
| 226 | +analysis_interval = 1000 |
| 227 | +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) |
| 228 | + |
| 229 | +alive_callback = AliveCallback(analysis_interval = analysis_interval) |
| 230 | + |
| 231 | +save_solution = SaveSolutionCallback(dt = 100, save_initial_solution = true, |
| 232 | + save_final_solution = true) |
| 233 | + |
| 234 | +callbacks = CallbackSet(summary_callback, |
| 235 | + analysis_callback, |
| 236 | + alive_callback, |
| 237 | + save_solution) |
| 238 | + |
| 239 | +############################################################################### |
| 240 | +# Use a Runge-Kutta method with automatic (error based) time step size control |
| 241 | +# Enable threading of the RK method for better performance on multiple threads |
| 242 | +tol = 1e-6 |
| 243 | +sol = solve(ode, |
| 244 | + SSPRK43(thread = Trixi.True()); |
| 245 | + abstol = tol, reltol = tol, ode_default_options()..., |
| 246 | + callback = callbacks) |
0 commit comments