trixi-framework · benegee · Jun 24, 2025 · Jun 25, 2025 · Jun 25, 2025 · Jun 27, 2025
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -26,7 +26,7 @@ concurrency:
 
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -35,16 +35,14 @@ jobs:
           - '1.10'
         os:
           - ubuntu-latest
-          - macOS-latest
+          - macos-latest
           - windows-latest
-        arch:
-          - x64
     steps:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
+          arch: default
           show-versioninfo: true
       - uses: julia-actions/cache@v2
 

diff --git a/examples/elixir_euler_potential_temperature_baroclinic_instability.jl b/examples/elixir_euler_potential_temperature_baroclinic_instability.jl
@@ -18,16 +18,16 @@ using LinearAlgebra: norm
 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  # Ω
+    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
@@ -138,7 +138,7 @@ end
 function initial_condition_baroclinic_instability(x, t,
                                                   equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)
     lon, lat, r = cartesian_to_sphere(x)
-    radius_earth = 6.371229e6
+    radius_earth = EARTH_RADIUS
     # Make sure that the r is not smaller than radius_earth
     z = max(r - radius_earth, 0.0)
 
@@ -155,8 +155,7 @@ function initial_condition_baroclinic_instability(x, t,
     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
+    gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION  # g
 
     r = norm(x)
     # Make sure that r is not smaller than radius_earth
@@ -174,8 +173,8 @@ end
 
 @inline function source_terms_baroclinic_instability(u, x, t,
                                                      equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)
-    radius_earth = 6.371229e6  # a
-    angular_velocity = 7.29212e-5  # Ω
+    radius_earth = EARTH_RADIUS  # a
+    angular_velocity = EARTH_ROTATION_RATE  # Ω
 
     r = norm(x)
     # Make sure that r is not smaller than radius_earth
@@ -195,7 +194,7 @@ end
 end
 equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004,
                                                                         c_v = 717,
-                                                                        gravity = 9.81)
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 
 initial_condition = initial_condition_baroclinic_instability
 
@@ -209,7 +208,7 @@ 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)
-mesh = P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000,
+mesh = P4estMeshCubedSphere(trees_per_cube_face..., EARTH_RADIUS, 30000,
                             polydeg = polydeg, initial_refinement_level = 0)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

diff --git a/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl b/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl
@@ -38,7 +38,7 @@ end
 
 equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
                                                                         c_v = 717,
-                                                                        gravity = 9.81)
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 
 # 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
@@ -97,7 +97,7 @@ end
 
 equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
                                                                         c_v = 717,
-                                                                        gravity = 9.81)
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 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
@@ -93,7 +93,7 @@ end
 
 equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
                                                                         c_v = 717,
-                                                                        gravity = 9.81)
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 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
@@ -12,7 +12,7 @@ function initial_condition_robert_bubble(x, t,
     # center of perturbation
     center_x = 500.0
     center_z = 260.0
-    g = 9.81
+    g = EARTH_GRAVITATIONAL_ACCELERATION
     # radius of perturbation
     radius = 250.0
     # distance of current x to center of perturbation
@@ -47,7 +47,7 @@ end
 @inline function source_terms_gravity(u, x, t,
                                       equations::CompressibleEulerPotentialTemperatureEquations2D)
     rho, _, _, _ = u
-    g = 9.81
+    g = EARTH_GRAVITATIONAL_ACCELERATION
     return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, zero(eltype(u)))
 end
 

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

diff --git a/examples/elixir_euler_potential_temperature_well_balanced.jl b/examples/elixir_euler_potential_temperature_well_balanced.jl
@@ -34,7 +34,7 @@ end
 
 equations = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004,
                                                                         c_v = 717,
-                                                                        gravity = 9.81)
+                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 polydeg = 2
 basis = LobattoLegendreBasis(polydeg)
 cs = 340.0

diff --git a/examples/elixir_moist_euler_EC_bubble.jl b/examples/elixir_moist_euler_EC_bubble.jl
@@ -11,104 +11,8 @@ 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)
-    @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations
-    (p, rho, T, r_t, r_v, rho_qv, theta_e) = y
-    p0 = y0[1]
-
-    F = zeros(7, 1)
-    rho_d = rho / (1 + r_t)
-    p_d = R_d * rho_d * T
-    T_C = T - 273.15
-    p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C))
-    L = L_00 - (c_pl - c_pv) * T
-
-    F[1] = (p - p0) / dz + g * rho
-    F[2] = p - (R_d * rho_d + R_v * rho_qv) * T
-    # H = 1 is assumed
-    F[3] = (theta_e -
-            T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) *
-            exp(L * r_v / ((c_pd + c_pl * r_t) * T)))
-    F[4] = r_t - r_t0
-    F[5] = rho_qv - rho_d * r_v
-    F[6] = theta_e - theta_e0
-    a = p_vs / (R_v * T) - rho_qv
-    b = rho - rho_qv - rho_d
-    # H=1 => phi=0
-    F[7] = a + b - sqrt(a * a + b * b)
-
-    return F
-end
-
-function generate_function_of_y(dz, y0, r_t0, theta_e0,
-                                equations::CompressibleMoistEulerEquations2D)
-    function function_of_y(y)
-        return moist_state(y, dz, y0, r_t0, theta_e0, equations)
-    end
-end
-
-#Create Initial atmosphere by generating a layer data set
-struct AtmosphereLayers{RealT <: Real}
-    equations::CompressibleMoistEulerEquations2D
-    # structure:  1--> i-layer (z = total_height/precision *(i-1)),  2--> rho, rho_theta, rho_qv, rho_ql
-    layer_data::Matrix{RealT} # Contains the layer data for each height
-    total_height::RealT # Total height of the atmosphere
-    preciseness::Int # Space between each layer data (dz)
-    layers::Int # Amount of layers (total height / dz)
-    ground_state::NTuple{2, RealT} # (rho_0, p_tilde_0) to define the initial values at height z=0
-    equivalent_potential_temperature::RealT  # Value for theta_e since we have a constant temperature theta_e0=theta_e.
-    mixing_ratios::NTuple{2, RealT} # Constant mixing ratio. Also defines initial guess for rho_qv0 = r_v0 * rho_0.
-end
-
-function AtmosphereLayers(equations; total_height = 10010.0, preciseness = 10,
-                          ground_state = (1.4, 100000.0),
-                          equivalent_potential_temperature = 320,
-                          mixing_ratios = (0.02, 0.02), RealT = Float64)
-    @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl = equations
-    rho0, p0 = ground_state
-    r_t0, r_v0 = mixing_ratios
-    theta_e0 = equivalent_potential_temperature
-
-    rho_qv0 = rho0 * r_v0
-    T0 = theta_e0
-    y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0]
-
-    n = convert(Int, total_height / preciseness)
-    dz = 0.01
-    layer_data = zeros(RealT, n + 1, 4)
-
-    F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
-    sol = nlsolve(F, y0)
-    p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero
-
-    rho_d = rho / (1 + r_t)
-    rho_ql = rho - rho_d - rho_qv
-    kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
-    rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t)
-
-    layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql]
-    for i in (1:n)
-        y0 = deepcopy(sol.zero)
-        dz = preciseness
-        F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
-        sol = nlsolve(F, y0)
-        p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero
-
-        rho_d = rho / (1 + r_t)
-        rho_ql = rho - rho_d - rho_qv
-        kappa_M = (R_d * rho_d + R_v * rho_qv) /
-                  (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
-        rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t)
-
-        layer_data[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql]
-    end
-
-    return AtmosphereLayers{RealT}(equations, layer_data, total_height, dz, n, ground_state,
-                                   theta_e0, mixing_ratios)
-end
+                                              c_vv = c_vv,
+                                              gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 
 # Generate background state from the Layer data set by linearely interpolating the layers
 function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D,
@@ -137,8 +41,8 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE
     rho, rho_e, rho_qv, rho_ql = perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql,
                                                         equations::CompressibleMoistEulerEquations2D)
 
-    v1 = 60.0
-    v2 = 60.0
+    v1 = 60
+    v2 = 60
     rho_v1 = rho * v1
     rho_v2 = rho * v2
     rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2)
@@ -148,7 +52,7 @@ end
 
 # Add perturbation to the profile
 function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql,
-                                equations::CompressibleMoistEulerEquations2D)
+                                equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT}
     @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl, L_00 = equations
     xc = 2000
     zc = 2000
@@ -166,16 +70,17 @@ function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql,
     # Assume pressure stays constant
     if (r < rc && Δθ > 0)
         θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa)
-        θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300)
+        θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5f0 * r / rc)^2 / 300)
         rt = (rho_qv + rho_ql) / rho_d
         rv = rho_qv / rho_d
         θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv)
         if rt > 0
             while true
                 T_loc = θ_loc * (p_loc / p_0)^kappa
-                T_C = T_loc - 273.15
+                T_C = T_loc - convert(RealT, 273.15)
                 # SaturVapor
-                pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C))
+                pvs = convert(RealT, 611.2) *
+                      exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C))
                 rho_d_new = (p_loc - pvs) / (R_d * T_loc)
                 rvs = pvs / (R_v * rho_d_new * T_loc)
                 θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs)
@@ -223,14 +128,13 @@ source_term = source_terms_geopotential
 # Get the DG approximation space
 
 polydeg = 4
-basis = LobattoLegendreBasis(polydeg)
 
 surface_flux = flux_chandrashekar
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
-solver = DGSEM(basis, surface_flux, volume_integral)
+solver = DGSEM(polydeg, surface_flux, volume_integral)
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (4000.0, 4000.0)

diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl
@@ -10,23 +10,24 @@ 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)
+                                              c_vv = c_vv,
+                                              gravity = EARTH_GRAVITATIONAL_ACCELERATION)
 
 # Warm bubble test from paper:
 # Wicker, L. J., and W. C. Skamarock, 1998: A time-splitting scheme
 # for the elastic equations incorporating second-order Runge–Kutta
 # time differencing. Mon. Wea. Rev., 126, 1992–1999.
 function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEquations2D)
     @unpack p_0, kappa, g, c_pd, c_vd, R_d, R_v = equations
-    xc = 10000.0
-    zc = 2000.0
+    xc = 10000
+    zc = 2000
     r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2)
-    rc = 2000.0
-    θ_ref = 300.0
-    Δθ = 0.0
+    rc = 2000
+    θ_ref = 300
+    Δθ = 0
 
     if r <= rc
-        Δθ = 2 * cospi(0.5 * r / rc)^2
+        Δθ = 2 * cospi(0.5f0 * r / rc)^2
     end
 
     # Perturbed state:
@@ -41,9 +42,9 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEq
     rho = p / ((p / p_0)^kappa * R_d * θ)
     T = p / (R_d * rho)
 
-    v1 = 20.0
-    # v1 = 0.0
-    v2 = 0.0
+    v1 = 20
+    # v1 = 0
+    v2 = 0
     rho_v1 = rho * v1
     rho_v2 = rho * v2
     rho_E = rho * c_vd * T + 1 / 2 * rho * (v1^2 + v2^2)
@@ -61,7 +62,6 @@ source_term = source_terms_geopotential
 ###############################################################################
 # Get the DG approximation space
 polydeg = 4
-basis = LobattoLegendreBasis(polydeg)
 
 surface_flux = FluxLMARS(360.0)
 volume_flux = flux_chandrashekar
@@ -70,7 +70,7 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
 # Create DG solver with polynomial degree = 4 and LMARS flux as surface flux
 # and the EC flux (chandrashekar) as volume flux
-solver = DGSEM(basis, surface_flux, volume_integral)
+solver = DGSEM(polydeg, surface_flux, volume_integral)
 
 coordinates_min = (0.0, -5000.0)
 coordinates_max = (20000.0, 15000.0)