Add turbo version of baroclinic instability

examples/elixir_euler_energy_baroclinic_instability_turbo.jl

@@ -0,0 +1,280 @@
+# An idealized baroclinic instability test case
+# For optimal results consider increasing the resolution to 16x16x8 trees per cube face, i.e.,
+# set `trees_per_cube_face = (16, 8)` below.
+#
+# 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 OrdinaryDiffEqLowStorageRK
+using Trixi, TrixiAtmo
+using LinearAlgebra: norm
+
+@inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_kennedy_gruber_souza_etal_turbo),
+equations::CompressibleEulerEnergyEquationsWithGravity3D) = Trixi.True()
+
+@inline function flux_lmars_turbo(u_ll, u_rr, normal_direction, equations)
+    flux = FluxLMARS(340)(u_ll, u_rr, normal_direction, equations)
+    return flux, flux
+end
+
+@inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_lmars_turbo), equations::CompressibleEulerEnergyEquationsWithGravity3D) = Trixi.True()
+
+@inline function boundary_condition_slip_wall(u_inner,
+                                              normal_direction::AbstractVector,
+                                              x, t,
+                                              surface_flux_function::typeof(flux_lmars_turbo),
+                                              equations::CompressibleEulerEnergyEquationsWithGravity3D)
+    # normalize the outward pointing direction
+    normal = normal_direction / norm(normal_direction)
+    # compute the normal momentum
+    u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3] + normal[3] * u_inner[4]
+
+    # create the "external" boundary solution state
+    u_boundary = SVector(u_inner[1],
+                         u_inner[2] - 2 * u_normal * normal[1],
+                         u_inner[3] - 2 * u_normal * normal[2],
+                         u_inner[4] - 2 * u_normal * normal[3],
+                         u_inner[5], u_inner[6])
+
+    # calculate the boundary flux
+    flux, _ = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    return flux
+end
+
+# 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 = 6.371229e6  # a
+    half_width_parameter = 2           # b
+    gravitational_acceleration = 9.81     # 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 = 7.29212e-5  # Ω
+
+    # 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::CompressibleEulerEnergyEquationsWithGravity3D)
+    lon, lat, r = cartesian_to_sphere(x)
+    radius_earth = 6.371229e6
+    # 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
+    radius_earth = 6.371229e6  # a
+    gravitational_acceleration = 9.81    # 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::CompressibleEulerEnergyEquationsWithGravity3D)
+    radius_earth = 6.371229e6  # a
+    angular_velocity = 7.29212e-5  # Ω
+
+    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 = CompressibleEulerEnergyEquationsWithGravity3D(c_p = 1004,
+                                                          c_v = 717,
+                                                          gravity = 9.81)
+
+initial_condition = initial_condition_baroclinic_instability
+
+boundary_conditions = Dict(:inside => boundary_condition_slip_wall,
+                           :outside => boundary_condition_slip_wall)
+
+polydeg = 5
+surface_flux = flux_lmars_turbo
+volume_flux = flux_kennedy_gruber_souza_etal_turbo
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+trees_per_cube_face = (8, 4)
+mesh = P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000,
+                            polydeg = polydeg, initial_refinement_level = 0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_baroclinic_instability,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+T = 10 # 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 = 3840, 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
+sol = solve(ode,
+            RDPK3SpFSAL49(thread = Trixi.True());
+            abstol = 1.0e-5, reltol = 1.0e-5, ode_default_options()...,
+            callback = callbacks);

src/TrixiAtmo.jl

@@ -65,6 +65,8 @@ export flux_nonconservative_zeros, flux_nonconservative_ec,
        flux_nonconservative_artiano_etal,
        flux_nonconservative_waruszewski_etal, flux_zero
 
+export flux_kennedy_gruber_souza_etal_turbo
+
 export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier!
 
 export cons2prim_and_vorticity, contravariant2global

src/equations/compressible_euler_energy_with_gravity_3d.jl

@@ -137,12 +137,12 @@ The well balanced on curvilinear coordinates was proven by
                                                        equations::CompressibleEulerEnergyEquationsWithGravity3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _, phi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _, phi_rr = u_rr
-    v1_ll = rho_v1_ll
-    v2_ll = rho_v2_ll
-    v3_ll = rho_v3_ll
-    v1_rr = rho_v1_rr
-    v2_rr = rho_v2_rr
-    v3_rr = rho_v3_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v3_rr = rho_v3_rr / rho_rr
     v1_avg = 0.5f0 * (v1_ll + v1_rr)
     v2_avg = 0.5f0 * (v2_ll + v2_rr)
     v3_avg = 0.5f0 * (v3_ll + v3_rr)
@@ -174,12 +174,12 @@ Well-balanced gravity term for constant potential temperature background state b
                                                    equations::CompressibleEulerEnergyEquationsWithGravity3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _, phi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _, phi_rr = u_rr
-    v1_ll = rho_v1_ll
-    v2_ll = rho_v2_ll
-    v3_ll = rho_v3_ll
-    v1_rr = rho_v1_rr
-    v2_rr = rho_v2_rr
-    v3_rr = rho_v3_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v3_rr = rho_v3_rr / rho_rr
     v1_avg = 0.5f0 * (v1_ll + v1_rr)
     v2_avg = 0.5f0 * (v2_ll + v2_rr)
     v3_avg = 0.5f0 * (v3_ll + v3_rr)
@@ -208,12 +208,12 @@ end
                                                  equations::CompressibleEulerEnergyEquationsWithGravity3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _, phi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _, phi_rr = u_rr
-    v1_ll = rho_v1_ll
-    v2_ll = rho_v2_ll
-    v3_ll = rho_v3_ll
-    v1_rr = rho_v1_rr
-    v2_rr = rho_v2_rr
-    v3_rr = rho_v3_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v3_rr = rho_v3_rr / rho_rr
     v1_avg = 0.5f0 * (v1_ll + v1_rr)
     v2_avg = 0.5f0 * (v2_ll + v2_rr)
     v3_avg = 0.5f0 * (v3_ll + v3_rr)
@@ -289,8 +289,8 @@ Kinetic energy preserving two-point flux by
 @inline function flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector,
                                      equations::CompressibleEulerEnergyEquationsWithGravity3D)
     # Unpack left and right state
-    rho_e_ll = last(u_ll)
-    rho_e_rr = last(u_rr)
+    rho_e_ll = u_ll[5]
+    rho_e_rr = u_rr[5]
     rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations)
 

src/solvers/dgsem_p4est/dg.jl

@@ -50,3 +50,4 @@ include("containers_2d_manifold_in_3d_covariant.jl")
 
 include("dg_2d_manifold_in_3d_cartesian.jl")
 include("dg_2d_manifold_in_3d_covariant.jl")
+include("dg_3d_compressible_euler.jl")