Add 2D Potential Temperature Formulation for the Euler Eqs.

examples/elixir_euler_potential_temperature_dry_bubble.jl

@@ -54,7 +54,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
-surface_flux = flux_LMARS
+surface_flux = flux_LMARS(340.0)
 
 volume_flux = flux_theta
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

examples/elixir_euler_potential_temperature_ec.jl

@@ -0,0 +1,77 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+using TrixiAtmo: flux_theta
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerPotentialTemperatureEquations2D()
+
+function initial_condition_weak_blast_wave(x, t,
+                                           equations::CompressibleEulerPotentialTemperatureEquations2D)
+    # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
+    # Set up polar coordinates
+    inicenter = SVector(0, 0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+    phi = atan(y_norm, x_norm)
+    sin_phi, cos_phi = sincos(phi)
+
+    # Calculate primitive variables
+    RealT = eltype(x)
+    rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
+    v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi
+    p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245)
+
+    return TrixiAtmo.prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+volume_flux = flux_theta
+solver = DGSEM(polydeg = 3, surface_flux = flux_theta,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+
+cells_per_dimension = (32, 32)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = (true, true))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/elixir_moist_euler_dry_bubble.jl

@@ -61,7 +61,7 @@ source_term = source_terms_geopotential
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
 
-surface_flux = flux_LMARS
+surface_flux = flux_LMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

examples/elixir_moist_euler_moist_bubble.jl

@@ -242,7 +242,7 @@ source_term = source_terms_moist_bubble
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
 
-surface_flux = flux_LMARS
+surface_flux = flux_LMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl

@@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
 
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
-surface_flux = flux_LMARS
+surface_flux = flux_LMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

src/equations/compressible_euler_potential_temperature_2d.jl

@@ -12,7 +12,6 @@
     g::RealT
     R::RealT
     gamma::RealT
-    a::RealT
 end
 
 function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Float64)
@@ -21,9 +20,8 @@
     c_v = 717.0
     R = c_p - c_v
     gamma = c_p / c_v
-    a = 340.0
     return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R,
-                                                                   gamma, a)
+                                                                   gamma)
 end
 
 function varnames(::typeof(cons2cons),
@@ -38,40 +36,40 @@
 
 # Calculate 1D flux for a single point in the normal direction.
 # Note, this directional vector is not normalized.
 @inline function flux(u, normal_direction::AbstractVector,
                       equations::CompressibleEulerPotentialTemperatureEquations2D)
     rho, rho_v1, rho_v2, rho_theta = u
     v1 = rho_v1 / rho
     v2 = rho_v2 / rho
     v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
     rho_v_normal = rho * v_normal
     f1 = rho_v_normal
     f2 = (rho_v_normal) * v1 + p * normal_direction[1]
     f3 = (rho_v_normal) * v2 + p * normal_direction[2]
     f4 = (rho_theta) * v_normal
     return SVector(f1, f2, f3, f4)
 end
 
 @inline function flux(u, orientation::Integer,
                       equations::CompressibleEulerPotentialTemperatureEquations2D)
     rho, rho_v1, rho_v2, rho_theta = u
     v1 = rho_v1 / rho
     v2 = rho_v2 / rho
     p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma
 
     if orientation == 1
         f1 = rho_v1
         f2 = rho_v1 * v1 + p
         f3 = rho_v1 * v2
         f4 = (rho_theta) * v1
     else
         f1 = rho_v2
         f2 = rho_v2 * v1
         f3 = rho_v2 * v2 + p
         f4 = (rho_theta) * v2
     end
 
     return SVector(f1, f2, f3, f4)
 end
 
 # Slip-wall boundary condition
@@ -99,13 +97,13 @@
     if v_normal <= 0.0
         sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed
         p_star = p_local *
-                 (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma *
-                                                                     inv(gamma - 1))
+                 (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma *
+                                                                       inv(gamma - 1))
     else # v_normal > 0.0
         A = 2.0 / ((gamma + 1) * rho_local)
         B = p_local * (gamma - 1) / (gamma + 1)
         p_star = p_local +
-                 0.5 * v_normal / A *
+                 0.5f0 * v_normal / A *
                  (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B)))
     end
 
@@ -169,23 +167,23 @@
 # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian
 # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013,
 # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml.
-@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector,
-                            equations::CompressibleEulerPotentialTemperatureEquations2D)
-    @unpack a = equations
+
+@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector,
+                                          equations::CompressibleEulerPotentialTemperatureEquations2D)
+    a = flux_lmars.speed_of_sound
     # Unpack left and right state
     rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
 
     v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
     v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
 
-    # Compute the necessary interface flux components
     norm_ = norm(normal_direction)
 
-    rho = 0.5 * (rho_ll + rho_rr)
+    rho = 0.5f0 * (rho_ll + rho_rr)
 
-    p_interface = 0.5 * (p_ll + p_rr) - 0.5 * a * rho * (v_rr - v_ll) / norm_
-    v_interface = 0.5 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_
+    p_interface = 0.5f0 * (p_ll + p_rr) - 0.5f0 * a * rho * (v_rr - v_ll) / norm_
+    v_interface = 0.5f0 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_
 
     if (v_interface > 0)
         f1, f2, f3, f4 = u_ll * v_interface
@@ -200,6 +198,7 @@
 end
 
 ## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation
+
 @inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector,
                             equations::CompressibleEulerPotentialTemperatureEquations2D)
     # Unpack left and right state
@@ -211,20 +210,18 @@
     _, _, _, rho_theta_rr = u_rr
     # Compute the necessary mean values
     rho_mean = ln_mean(rho_ll, rho_rr)
+
     gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma)
-    # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
-    # in exact arithmetic since
-    #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
-    #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.5 * (p_ll + p_rr)
+
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
 
     # Calculate fluxes depending on normal_direction
-    f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
+    f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
     f2 = f1 * v1_avg + p_avg * normal_direction[1]
     f3 = f1 * v2_avg + p_avg * normal_direction[2]
-    f4 = gammamean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
+    f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
     return SVector(f1, f2, f3, f4)
 end
 
@@ -257,7 +254,7 @@
     rho, rho_v1, rho_v2, rho_theta = u
 
     k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma
-    w1 = -0.5 * rho_v1^2 / (rho)^2 - 0.5 * rho_v2^2 / (rho)^2
+    w1 = -0.5f0 * rho_v1^2 / (rho)^2 - 0.5f0 * rho_v2^2 / (rho)^2
     w2 = rho_v1 / rho
     w3 = rho_v2 / rho
     w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1)
@@ -265,30 +262,30 @@
     return SVector(w1, w2, w3, w4)
 end
 
 @inline function entropy_math(cons,
                               equations::CompressibleEulerPotentialTemperatureEquations2D)
     # Mathematical entropy
     p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma
 
     U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1]))
 
     return U
 end
 
 # Default entropy is the mathematical entropy
 @inline function entropy(cons,
                          equations::CompressibleEulerPotentialTemperatureEquations2D)
     entropy_math(cons, equations)
 end
 
 @inline function energy_total(cons,
                               equations::CompressibleEulerPotentialTemperatureEquations2D)
     entropy(cons, equations)
 end
 
 @inline function energy_kinetic(cons,
                                 equations::CompressibleEulerPotentialTemperatureEquations2D)
-    return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1])
+    return 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1])
 end
 
 @inline function max_abs_speeds(u,

src/equations/compressible_moist_euler_2d_lucas.jl

@@ -23,7 +23,6 @@ struct CompressibleMoistEulerEquations2D{RealT <: Real} <:
     kappa::RealT # ratio of the gas constant R_d
     gamma::RealT # = inv(kappa- 1); can be used to write slow divisions as fast multiplications
     L_00::RealT # latent heat of evaporation  at 0 K
-    a::RealT
 end
 
 function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64)
@@ -38,9 +37,8 @@ function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64)
     gamma = c_pd / c_vd # = 1/(1 - kappa)
     kappa = 1 - inv(gamma)
     L_00 = 3147620.0
-    a = 360.0
     return CompressibleMoistEulerEquations2D{RealT}(p_0, c_pd, c_vd, R_d, c_pv, c_vv,
-                                                    R_v, c_pl, g, kappa, gamma, L_00, a)
+                                                    R_v, c_pl, g, kappa, gamma, L_00)
 end
 
 # Calculate 1D flux for a single point.
@@ -452,9 +450,10 @@ end
 # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian
 # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013,
 # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml.
-@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector,
-                            equations::CompressibleMoistEulerEquations2D)
-    @unpack a = equations
+
+@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector,
+                                          equations::CompressibleMoistEulerEquations2D)
+    a = flux_lmars.speed_of_sound
     # Unpack left and right state
     rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, rho_qv_rr, rho_ql_rr = u_rr

src/equations/equations.jl

@@ -4,5 +4,11 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
 abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
+
+struct flux_LMARS{SpeedOfSound}
+    # Estimate for the speed of sound
+    speed_of_sound::SpeedOfSound
+end
+
 include("compressible_moist_euler_2d_lucas.jl")
 include("compressible_euler_potential_temperature_2d.jl")

test/test_2d_euler_potential_temperature.jl

@@ -34,4 +34,29 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory
     end
 end
 
+@trixiatmo_testset "elixir_euler_potential_temperature_ec" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_potential_temperature_ec.jl"),
+                        l2=[
+                            0.06174365254988615,
+                            0.05018008812040643,
+                            0.05018774429827882,
+                            0.0057820937341387935,
+                        ],
+                        linf=[
+                            0.2932352942992503,
+                            0.3108954213591686,
+                            0.31082193857647905,
+                            0.027465217019157606,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100
+    end
+end
+
 end # module

test/test_trixi_consistency.jl

@@ -49,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true)
     trixi_include(trixi_elixir,
                   equations = equations_moist,
                   volume_flux = flux_chandrashekar,
-                  surface_flux = flux_LMARS,
+                  surface_flux = flux_LMARS(360.0),
                   maxiters = maxiters)
 
     errors_atmo = Main.analysis_callback(Main.sol)