trixi-framework · benegee · Oct 9, 2025 · Oct 2, 2025 · Oct 2, 2025 · Oct 2, 2025
diff --git a/examples/elixir_euler_potential_temperature_baroclinic_instability.jl b/examples/elixir_euler_potential_temperature_baroclinic_instability.jl
@@ -193,7 +193,9 @@ end
 
     return SVector(du1, du2, du3, du4, du5, zero(eltype(u)))
 end
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 
 initial_condition = initial_condition_baroclinic_instability
 

diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl
@@ -9,7 +9,7 @@ function initial_condition_density_wave(x, t,
     return prim2cons(SVector(rho, v1, p), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquations1D()
+equations = CompressibleEulerPotentialTemperatureEquations1D(c_p = 1004, c_v = 717)
 
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)

diff --git a/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl b/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl
@@ -36,7 +36,9 @@ function initial_condition_gravity_waves(x, t,
     return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 
 # We have an isothermal background state with T0 = 250 K. 
 # The reference speed of sound can be computed as:

diff --git a/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl b/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl
@@ -95,7 +95,9 @@ end
     return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 alpha = 0.035
 xr_B = 60000.0
 linear_hydrostatic_setup = HydrostaticSetup(alpha, xr_B, equations)

diff --git a/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl b/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl
@@ -91,7 +91,9 @@ end
     return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 alpha = 0.03
 xr_B = 40000.0
 

diff --git a/examples/elixir_euler_potential_temperature_robert_bubble.jl b/examples/elixir_euler_potential_temperature_robert_bubble.jl
@@ -51,7 +51,7 @@ end
     return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, zero(eltype(u)))
 end
 
-equations = CompressibleEulerPotentialTemperatureEquations2D()
+equations = CompressibleEulerPotentialTemperatureEquations2D(c_p = 1004, c_v = 717)
 
 surface_flux = FluxLMARS(340)
 volume_flux = flux_tec

diff --git a/examples/elixir_euler_potential_temperature_schaer_mountain.jl b/examples/elixir_euler_potential_temperature_schaer_mountain.jl
@@ -91,7 +91,9 @@ end
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 alpha = 0.03
 xr_B = 20000
 schär_setup = SchärSetup(alpha, xr_B)

diff --git a/examples/elixir_euler_potential_temperature_taylor_green_vortex.jl b/examples/elixir_euler_potential_temperature_taylor_green_vortex.jl
@@ -14,7 +14,7 @@ function initial_condition_taylor_green_vortex(x, t,
     return prim2cons(SVector(rho, v1, v2, v3, p), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquations3D()
+equations = CompressibleEulerPotentialTemperatureEquations3D(c_p = 1004, c_v = 717)
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 surface_flux = flux_etec

diff --git a/examples/elixir_euler_potential_temperature_well_balanced.jl b/examples/elixir_euler_potential_temperature_well_balanced.jl
@@ -32,7 +32,9 @@ function initial_condition_adiabatic(x, t,
     return prim2cons(SVector(rho, v1, p, g * x[1]), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity1D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 polydeg = 2
 basis = LobattoLegendreBasis(polydeg)
 cs = 340.0

diff --git a/examples/elixir_euler_potential_temperature_well_balanced_curvilinear.jl b/examples/elixir_euler_potential_temperature_well_balanced_curvilinear.jl
@@ -33,7 +33,9 @@ function initial_condition_adiabatic(x, t,
     return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
 end
 
-equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D()
+equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
+                                                                        c_v = 717,
+                                                                        gravity = 9.81)
 polydeg = 2
 basis = LobattoLegendreBasis(polydeg)
 cs = 340

diff --git a/examples/elixir_moist_euler_EC_bubble.jl b/examples/elixir_moist_euler_EC_bubble.jl
@@ -6,7 +6,12 @@ using NLsolve: nlsolve
 ###############################################################################
 # semidiscretization of the compressible moist Euler equations
 
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 function moist_state(y, dz, y0, r_t0, theta_e0,
                      equations::CompressibleMoistEulerEquations2D)

diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl
@@ -5,8 +5,12 @@ using TrixiAtmo: source_terms_geopotential, cons2drypot
 
 ###############################################################################
 # semidiscretization of the compressible moist Euler equations
-
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 # Warm bubble test from paper:
 # Wicker, L. J., and W. C. Skamarock, 1998: A time-splitting scheme

diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl
@@ -6,7 +6,12 @@ using NLsolve: nlsolve
 ###############################################################################
 # semidiscretization of the compressible moist Euler equations
 
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 function moist_state(y, dz, y0, r_t0, theta_e0,
                      equations::CompressibleMoistEulerEquations2D)

diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl
@@ -35,7 +35,12 @@ function initial_condition_nonhydrostatic_gravity_wave(x, t,
     return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql)
 end
 
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 initial_condition = initial_condition_nonhydrostatic_gravity_wave
 

diff --git a/examples/elixir_moist_euler_source_terms_dry.jl b/examples/elixir_moist_euler_source_terms_dry.jl
@@ -5,7 +5,12 @@ using TrixiAtmo: source_terms_convergence_test_dry, initial_condition_convergenc
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
-equations = CompressibleMoistEulerEquations2D(; g = 0.0)
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 0.0)
 
 initial_condition = initial_condition_convergence_test_dry
 

diff --git a/examples/elixir_moist_euler_source_terms_moist.jl b/examples/elixir_moist_euler_source_terms_moist.jl
@@ -6,7 +6,12 @@ using TrixiAtmo: source_terms_convergence_test_moist,
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 initial_condition = initial_condition_convergence_test_moist
 

diff --git a/examples/elixir_moist_euler_source_terms_split_moist.jl b/examples/elixir_moist_euler_source_terms_split_moist.jl
@@ -6,7 +6,12 @@ using TrixiAtmo: source_terms_convergence_test_moist,
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
-equations = CompressibleMoistEulerEquations2D()
+c_pd = 1004 # specific heat at constant pressure for dry air
+c_vd = 717  # specific heat at constant volume for dry air
+c_pv = 1885 # specific heat at constant pressure for moist air
+c_vv = 1424 # specific heat at constant volume for moist air
+equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
+                                              c_vv = c_vv, gravity = 9.81)
 
 initial_condition = initial_condition_convergence_test_moist
 

diff --git a/src/equations/compressible_euler_potential_temperature_1d.jl b/src/equations/compressible_euler_potential_temperature_1d.jl
@@ -11,21 +11,19 @@ struct CompressibleEulerPotentialTemperatureEquations1D{RealT <: Real} <:
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
     K::RealT # = p_0 * (R / p_0)^gamma; scaling factor between pressure and weighted potential temperature
     stolarsky_factor::RealT # = (gamma - 1) / gamma; used in the stolarsky mean
-end
-
-function CompressibleEulerPotentialTemperatureEquations1D(; RealT = Float64)
-    p_0 = 100_000
-    c_p = 1004
-    c_v = 717
-    R = c_p - c_v
-    gamma = c_p / c_v
-    inv_gamma_minus_one = inv(gamma - 1)
-    K = p_0 * (R / p_0)^gamma
-    stolarsky_factor = (gamma - 1) / gamma
-    return CompressibleEulerPotentialTemperatureEquations1D{RealT}(p_0, c_p, c_v, R,
-                                                                   gamma,
-                                                                   inv_gamma_minus_one,
-                                                                   K, stolarsky_factor)
+    function CompressibleEulerPotentialTemperatureEquations1D(; c_p, c_v)
+        c_p, c_v = promote(c_p, c_v)
+        p_0 = 100_000
+        R = c_p - c_v
+        gamma = c_p / c_v
+        inv_gamma_minus_one = inv(gamma - 1)
+        K = p_0 * (R / p_0)^gamma
+        stolarsky_factor = (gamma - 1) / gamma
+        return new{typeof(gamma)}(p_0, c_p, c_v, R,
+                                  gamma,
+                                  inv_gamma_minus_one,
+                                  K, stolarsky_factor)
+    end
 end
 
 function varnames(::typeof(cons2cons),
@@ -36,6 +34,47 @@ end
 varnames(::typeof(cons2prim),
 ::CompressibleEulerPotentialTemperatureEquations1D) = ("rho", "v1", "p1")
 
+# Calculate 1D flux for a single point
+@inline function flux(u, orientation::Integer,
+                      equations::CompressibleEulerPotentialTemperatureEquations1D)
+    rho, rho_v1, rho_theta = u
+    v1 = rho_v1 / rho
+    p = pressure(u, equations)
+    f1 = rho_v1
+    f2 = rho_v1 * v1 + p
+    f3 = rho_theta * v1
+
+    return SVector(f1, f2, f3)
+end
+
+# Low Mach number approximate Riemann solver (LMARS) from
+# X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S.
+# 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::FluxLMARS)(u_ll, u_rr, orientation::Integer,
+                                         equations::CompressibleEulerPotentialTemperatureEquations1D)
+    a = flux_lmars.speed_of_sound
+    rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
+
+    rho = 0.5f0 * (rho_ll + rho_rr)
+
+    p_interface = 0.5f0 * (p_ll + p_rr) - 0.5f0 * a * rho * (v1_rr - v1_ll)
+    v_interface = 0.5f0 * (v1_ll + v1_rr) - 1 / (2 * a * rho) * (p_rr - p_ll)
+
+    if (v_interface > 0)
+        f1, f2, f3 = u_ll * v_interface
+    else
+        f1, f2, f3 = u_rr * v_interface
+    end
+
+    return SVector(f1,
+                   f2 + p_interface,
+                   f3)
+end
+
 """
 	flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature1D)
 
@@ -166,7 +205,7 @@ end
     # Mathematical entropy
     p = equations.p_0 * (equations.R * cons[3] / equations.p_0)^equations.gamma
 
-    U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2) / (cons[1]))
+    U = (p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2) / (cons[1]))
 
     return U
 end
@@ -192,7 +231,7 @@ end
 
 @inline function pressure(cons,
                           equations::CompressibleEulerPotentialTemperatureEquations1D)
-    _, _, p = cons2prim(cons, equations)
+    p = equations.K * exp(equations.gamma * log(cons[3]))
     return p
 end
 end # @muladd