diff --git a/examples/elixir_euler_potential_temperature_baroclinic_instability.jl b/examples/elixir_euler_potential_temperature_baroclinic_instability.jl index 984ab80a..c04bdc18 100644 --- 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 index 7a273408..eaf79016 100644 --- 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 index 50ba70bb..aa365ab8 100644 --- 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 index f322729e..5859aadd 100644 --- 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 index e2e8d9f3..9ed7e74b 100644 --- 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 index 617bbe04..229c5590 100644 --- 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 index d793eba1..9c166c70 100644 --- 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 index bdc3054c..9a26fbe8 100644 --- 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 index b3fa61c7..3d8cb16a 100644 --- 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 index 1f1795f1..da234415 100644 --- 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 index b58d3903..cdbea70e 100644 --- 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 index ec6ee0ae..4d61c1b4 100644 --- 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 index 0753be27..30b566d0 100644 --- 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 index 8ba54fbe..90da724d 100644 --- 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 index df5fa2f9..73b68533 100644 --- 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 index 36ea7ecb..7187cce1 100644 --- 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 index 25fd649b..1bc1af6e 100644 --- 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 index 05bde4c2..2692e1e0 100644 --- 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 diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index 17a0fbd8..16187541 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -10,21 +10,19 @@ struct CompressibleEulerPotentialTemperatureEquations2D{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 CompressibleEulerPotentialTemperatureEquations2D(; 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 CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, R, - gamma, - inv_gamma_minus_one, - K, stolarsky_factor) + function CompressibleEulerPotentialTemperatureEquations2D(; 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), @@ -38,6 +36,46 @@ varnames(::typeof(cons2prim), "v2", "p1") +# 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 + p = pressure(u, equations) + + 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 + +# Calculate 1D flux for a single point +@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 = pressure(u, equations) + 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 + @inline function boundary_condition_slip_wall(u_inner, orientation, direction, x, t, surface_flux_function, @@ -287,4 +325,73 @@ end λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll) + c_ll, abs(v_rr) + c_rr) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speeds + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + norm_ = norm(normal_direction) + return max(abs(v_ll) + c_ll * norm_, + abs(v_rr) + c_rr * norm_) +end +@inline function pressure(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) + p = equations.K * exp(equations.gamma * log(cons[4])) + return p +end end # @muladd diff --git a/src/equations/compressible_euler_potential_temperature_3d.jl b/src/equations/compressible_euler_potential_temperature_3d.jl index 549c851b..ed492d76 100644 --- a/src/equations/compressible_euler_potential_temperature_3d.jl +++ b/src/equations/compressible_euler_potential_temperature_3d.jl @@ -10,21 +10,19 @@ struct CompressibleEulerPotentialTemperatureEquations3D{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 CompressibleEulerPotentialTemperatureEquations3D(; 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 CompressibleEulerPotentialTemperatureEquations3D{RealT}(p_0, c_p, c_v, R, - gamma, - inv_gamma_minus_one, - K, stolarsky_factor) + function CompressibleEulerPotentialTemperatureEquations3D(; 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), @@ -39,6 +37,27 @@ varnames(::typeof(cons2prim), "v3", "p1") +# 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::CompressibleEulerPotentialTemperatureEquations3D) + rho, rho_v1, rho_v2, rho_v3, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + + v3 * normal_direction[3] + rho_v_normal = rho * v_normal + p = pressure(u, equations) + + f1 = rho_v_normal + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + f4 = rho_v_normal * v3 + p * normal_direction[3] + f5 = rho_theta * v_normal + return SVector(f1, f2, f3, f4, f5) +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 @@ -240,4 +259,127 @@ end S = -s * cons[1] / (equations.gamma - 1) return S end + +@inline function pressure(cons, + equations::CompressibleEulerPotentialTemperatureEquations3D) + p = equations.K * exp(equations.gamma * log(cons[5])) + return p +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function energy_kinetic(u, + equations::CompressibleEulerPotentialTemperatureEquations3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + return 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho +end + +@inline function energy_total(cons, + equations::CompressibleEulerPotentialTemperatureEquations3D) + # Mathematical entropy + p = equations.p_0 * (equations.R * cons[5] / equations.p_0)^equations.gamma + + U = (p / (equations.gamma - 1) + + 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / (cons[1])) + + return U +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquations3D) + 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) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations3D) + 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) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquations3D) + 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) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll) + c_ll, abs(v_rr) + c_rr) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations3D) + 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) + + # Calculate normal velocities and sound speeds + # left + v_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + norm_ = norm(normal_direction) + return max(abs(v_ll) + c_ll * norm_, abs(v_rr) + c_rr * norm_) +end + +@inline function max_abs_speeds(u, + equations::CompressibleEulerPotentialTemperatureEquations3D) + rho, v1, v2, v3, p = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c, abs(v3) + c +end end # @muladd diff --git a/src/equations/compressible_euler_potential_temperature_gravity_1d.jl b/src/equations/compressible_euler_potential_temperature_gravity_1d.jl index 1b35941c..e4aebdb7 100644 --- a/src/equations/compressible_euler_potential_temperature_gravity_1d.jl +++ b/src/equations/compressible_euler_potential_temperature_gravity_1d.jl @@ -12,24 +12,20 @@ struct CompressibleEulerPotentialTemperatureEquationsWithGravity1D{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 CompressibleEulerPotentialTemperatureEquationsWithGravity1D(; g = 9.81, - 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 CompressibleEulerPotentialTemperatureEquationsWithGravity1D{RealT}(p_0, c_p, - c_v, g, R, - gamma, - inv_gamma_minus_one, - K, - stolarsky_factor) + function CompressibleEulerPotentialTemperatureEquationsWithGravity1D(; c_p, c_v, + gravity) + c_p, c_v, g = promote(c_p, c_v, gravity) + 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(c_p)}(p_0, c_p, c_v, g, R, + gamma, + inv_gamma_minus_one, + K, stolarsky_factor) + end end function varnames(::typeof(cons2cons), @@ -65,6 +61,19 @@ have_nonconservative_terms(::CompressibleEulerPotentialTemperatureEquationsWithG return flux, noncons_flux end +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) + 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, 0) +end + """ flux_nonconservative_waruzewski_etal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) @@ -284,7 +293,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]) + cons[1] * cons[4] + U = p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2) / (cons[1]) + cons[1] * cons[4] return U end @@ -303,11 +312,25 @@ end return S end -@inline function Trixi.energy_kinetic(cons, - equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) +@inline function energy_kinetic(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) return 0.5f0 * (cons[2]^2) / (cons[1]) end +@inline function max_abs_speed(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) + rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations) + + # Calculate primitive variables and speed of sound + v_mag_ll = abs(v1_ll) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + v_mag_rr = abs(v1_rr) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(v_mag_ll + c_ll, v_mag_rr + c_rr) +end + @inline function max_abs_speeds(u, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) rho, v1, p = cons2prim(u, equations) @@ -330,9 +353,9 @@ end λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) end -@inline function Trixi.pressure(cons, - equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) - _, _, p = cons2prim(cons, equations) +@inline function pressure(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D) + p = equations.K * exp(equations.gamma * log(cons[3])) return p end end # @muladd diff --git a/src/equations/compressible_euler_potential_temperature_gravity_2d.jl b/src/equations/compressible_euler_potential_temperature_gravity_2d.jl index 69c69014..fade6c07 100644 --- a/src/equations/compressible_euler_potential_temperature_gravity_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_gravity_2d.jl @@ -11,24 +11,20 @@ struct CompressibleEulerPotentialTemperatureEquationsWithGravity2D{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 CompressibleEulerPotentialTemperatureEquationsWithGravity2D(; g = 9.81, - 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 CompressibleEulerPotentialTemperatureEquationsWithGravity2D{RealT}(p_0, c_p, - c_v, g, R, - gamma, - inv_gamma_minus_one, - K, - stolarsky_factor) + function CompressibleEulerPotentialTemperatureEquationsWithGravity2D(; c_p, c_v, + gravity) + c_p, c_v, g = promote(c_p, c_v, gravity) + 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(c_p)}(p_0, c_p, c_v, g, R, + gamma, + inv_gamma_minus_one, + K, stolarsky_factor) + end end function varnames(::typeof(cons2cons), @@ -86,7 +82,7 @@ end u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4], u_inner[5]) end - surface_flux_function = surface_flux_functions[1] + # Calculate boundary flux if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary flux = surface_flux_function(u_inner, u_boundary, orientation, equations) @@ -101,6 +97,46 @@ end return flux, noncons_flux end +# 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::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + 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 + p = pressure(u, equations) + + 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, zero(eltype(u))) +end + +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = pressure(u, equations) + 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, zero(eltype(u))) +end + """ flux_nonconservative_waruzewski_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) @@ -337,15 +373,21 @@ end rho, rho_v1, rho_v2, rho_theta = u w1 = log(equations.K * (rho_theta / rho)^equations.gamma) - equations.gamma - w2 = 0.0 - w3 = 0.0 w4 = rho / rho_theta * equations.gamma - return SVector(w1, w2, w3, w4, zero(eltype(u))) + return SVector(w1, zero(eltype(u)), zero(eltype(u)), w4, zero(eltype(u))) end @inline function entropy_thermodynamic(cons, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + p = equations.K * (cons[4])^equations.gamma + s = log(p) - equations.gamma * log(cons[1]) + S = -s * cons[1] / (equations.gamma - 1) + return S +end + +@inline function energy_total(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) # Mathematical entropy p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma @@ -360,6 +402,17 @@ end entropy_thermodynamic(cons, equations) end +@inline function pressure(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + p = equations.K * exp(equations.gamma * log(cons[4])) + return p +end + +@inline function energy_kinetic(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + return 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1]) +end + @inline function max_abs_speeds(u, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) rho, v1, v2, p = cons2prim(u, equations) @@ -367,4 +420,88 @@ end return abs(v1) + c, abs(v2) + c end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll) + c_ll, abs(v_rr) + c_rr) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speeds + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + norm_ = norm(normal_direction) + return max(abs(v_ll) + c_ll * norm_, + abs(v_rr) + c_rr * norm_) +end end # @muladd diff --git a/src/equations/compressible_euler_potential_temperature_gravity_3d.jl b/src/equations/compressible_euler_potential_temperature_gravity_3d.jl index a2db2846..22b887c1 100644 --- a/src/equations/compressible_euler_potential_temperature_gravity_3d.jl +++ b/src/equations/compressible_euler_potential_temperature_gravity_3d.jl @@ -11,24 +11,20 @@ struct CompressibleEulerPotentialTemperatureEquationsWithGravity3D{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 CompressibleEulerPotentialTemperatureEquationsWithGravity3D(; g = 9.81, - 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 CompressibleEulerPotentialTemperatureEquationsWithGravity3D{RealT}(p_0, c_p, - c_v, g, R, - gamma, - inv_gamma_minus_one, - K, - stolarsky_factor) + function CompressibleEulerPotentialTemperatureEquationsWithGravity3D(; c_p, c_v, + gravity) + c_p, c_v, g = promote(c_p, c_v, gravity) + 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(c_p)}(p_0, c_p, c_v, g, R, + gamma, + inv_gamma_minus_one, + K, stolarsky_factor) + end end function varnames(::typeof(cons2cons), @@ -70,6 +66,27 @@ have_nonconservative_terms(::CompressibleEulerPotentialTemperatureEquationsWithG return flux, noncons_flux end +# 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::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + + v3 * normal_direction[3] + rho_v_normal = rho * v_normal + p = pressure(u, equations) + + f1 = rho_v_normal + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + f4 = rho_v_normal * v3 + p * normal_direction[3] + f5 = rho_theta * v_normal + return SVector(f1, f2, f3, f4, f5, zero(eltype(u))) +end + """ flux_nonconservative_waruzewski_etal(u_ll, u_rr, normal_direction::AbstractVector, @@ -346,4 +363,127 @@ end S = -s * cons[1] / (equations.gamma - 1) return S end + +@inline function pressure(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + p = equations.K * exp(equations.gamma * log(cons[5])) + return p +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function energy_kinetic(u, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + return 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho +end + +@inline function energy_total(cons, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + # Mathematical entropy + p = equations.p_0 * (equations.R * cons[5] / equations.p_0)^equations.gamma + + U = (p / (equations.gamma - 1) + + 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / (cons[1])) + + return U +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + 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) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + 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) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + 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) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll) + c_ll, abs(v_rr) + c_rr) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + 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) + + # Calculate normal velocities and sound speeds + # left + v_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + norm_ = norm(normal_direction) + return max(abs(v_ll) + c_ll * norm_, abs(v_rr) + c_rr * norm_) +end + +@inline function max_abs_speeds(u, + equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) + rho, v1, v2, v3, p = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c, abs(v3) + c +end end # @muladd diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 9ba7c951..32a50a7e 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -16,22 +16,18 @@ 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 -end - -function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) - p_0 = 100000.0 - c_pd = 1004.0 - c_vd = 717.0 - R_d = c_pd - c_vd - c_pv = 1885.0 - c_vv = 1424.0 - R_v = c_pv - c_vv - c_pl = 4186.0 - gamma = c_pd / c_vd # = 1/(1 - kappa) - kappa = 1 - inv(gamma) - L_00 = 3147620 - return CompressibleMoistEulerEquations2D{RealT}(p_0, c_pd, c_vd, R_d, c_pv, c_vv, - R_v, c_pl, g, kappa, gamma, L_00) + function CompressibleMoistEulerEquations2D(; c_pd, c_vd, c_pv, c_vv, gravity) + c_pd, c_vd, c_pv, c_vv, g = promote(c_pd, c_vd, c_pv, c_vv, gravity) + p_0 = 100_000 + R_d = c_pd - c_vd + R_v = c_pv - c_vv + c_pl = 4186 + gamma = c_pd / c_vd # = 1/(1 - kappa) + kappa = 1 - inv(gamma) + L_00 = 3147620 + return new{typeof(gamma)}(p_0, c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, g, kappa, + gamma, L_00) + end end # Calculate 1D flux for a single point. @@ -523,18 +519,18 @@ end # Workaround if an individual density is zero # Thermodynamic entropy - s_d = 0 - s_v = 0 - s_l = 0 + s_d = zero(eltype(u)) + s_v = zero(eltype(u)) + s_l = zero(eltype(u)) # Thermodynamic entropy - if (rho_qd > 0.0) + if (rho_qd > 0) s_d = c_pd * log(T) - R_d * log(rho_qd * R_d * T) end - if (rho_qv > 0.0) + if (rho_qv > 0) s_v = c_pv * log(T) - R_v * log(rho_qv * R_v * T) end - if (rho_ql > 0.0) + if (rho_ql > 0) s_l = c_pl * log(T) end @@ -559,7 +555,7 @@ end rho_v2 = rho * v2 rho_qv = rho * qv rho_ql = rho * ql - rho_E = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) + rho_E = p / (equations.gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end diff --git a/test/runtests.jl b/test/runtests.jl index 509630e3..c1447abc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,6 +34,14 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_trixi_consistency.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "unit_fluxes" + include("test_unit.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "type_stable_tests" + include("test_type.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "moist_euler" include("test_2d_moist_euler.jl") end diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl index 947f874f..f360cc17 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -56,4 +56,5 @@ isdir(outdir) && rm(outdir, recursive = true) @test isapprox(error_trixi, error_atmo, rtol = 1.1e-10) end end + end diff --git a/test/test_type.jl b/test/test_type.jl new file mode 100644 index 00000000..90402c69 --- /dev/null +++ b/test/test_type.jl @@ -0,0 +1,372 @@ +module TestTypeStable + +using TrixiAtmo +using Trixi +using Test + +include("test_trixiatmo.jl") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "Test Type Stability" begin + @timed_testset "Compressible Euler Potential Temperature 1D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquations1D(c_p = RealT(1004), + c_v = RealT(717)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT)) + orientation = 1 + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, orientation, equations)) == + RealT + + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + end + end + + @timed_testset "Compressible Euler Potential Temperature With Gravity 1D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = RealT(1004), + c_v = RealT(717), + gravity = RealT(9.81)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) + orientation = 1 + directions = [1, 2] + + surface_flux_function = (flux_lax_friedrichs, flux_zero) + + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, + direction, + x, t, + surface_flux_function, + equations)) == + SVector{4, RealT} + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_nonconservative_waruzewski_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_artiano_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_souza_etal(u_ll, u_rr, orientation, + equations)) == + RealT + flux_lmars = FluxLMARS(RealT(340)) + @test eltype(@inferred flux_lmars(u_ll, u_rr, orientation, equations)) == + RealT + + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + end + end + + @timed_testset "Compressible Euler Potential Temperature 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquations2D(c_p = RealT(1004), + c_v = RealT(717)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT)) + surface_flux_function = flux_lax_friedrichs + orientations = [1, 2] + directions = [1, 2] + + for direction in directions, orientation in orientations + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, + direction, + x, t, + surface_flux_function, + equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction, + equations)) == RealT + end + end + + @timed_testset "Compressible Euler Potential Temperature With Gravity 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = RealT(1004), + c_v = RealT(717), + gravity = RealT(9.81)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT)) + surface_flux_function = (flux_lax_friedrichs, flux_zero) + orientations = [1, 2] + directions = [1, 2] + + for direction in directions, orientation in orientations + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, + direction, + x, t, + surface_flux_function, + equations)) == + SVector{5, RealT} + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, orientation, + equations)) == RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction, + equations)) == RealT + end + end + + @timed_testset "Compressible Euler Potential Temperature 3D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquations3D(c_p = RealT(1004), + c_v = RealT(717)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) + + surface_flux_function = flux_lax_friedrichs + orientations = [1, 2, 3] + + for orientation in orientations + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, orientation, + equations)) == RealT + end + normal_direction = SVector(one(RealT), one(RealT), one(RealT)) + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction, + equations)) == RealT + end + end + @timed_testset "Compressible Euler Potential Temperature With Gravity 3D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = RealT(1004), + c_v = RealT(717), + gravity = RealT(9.81)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT), one(RealT)) + + orientations = [1, 2, 3] + + for orientation in orientations + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, orientation, + equations)) == RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_tec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_etec(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction, + equations)) == RealT + end + end + + @timed_testset "Compressible Euler Moist Euler 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleMoistEulerEquations2D(c_pd = RealT(2.1), + c_vd = RealT(2), + c_pv = RealT(2.1), + c_vv = RealT(2), + gravity = RealT(1)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT)) + surface_flux_function = FluxLMARS(RealT(340)) + directions = [1, 2] + + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + direction, + x, t, + surface_flux_function, + equations)) == RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, + equations)) == RealT + + @test eltype(@inferred TrixiAtmo.cons2temp(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.cons2drypot(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.cons2moistpot(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.moist_pottemp_thermodynamic(u, equations)) == + RealT + @test eltype(@inferred TrixiAtmo.dry_pottemp_thermodynamic(u, equations)) == + RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.density(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.density_dry(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.density_vapor(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.temperature(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.density_liquid(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.ratio_liquid(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.ratio_vapor(u, equations)) == RealT + @test eltype(@inferred TrixiAtmo.density_pressure(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test eltype(@inferred energy_internal(u, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + end + end + + @timed_testset "Shallow Water 3D" begin + for RealT in (Float32, Float64) + equations = @inferred ShallowWaterEquations3D(gravity = RealT(1), + rotation_rate = RealT(1), + H0 = RealT(1)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT), one(RealT)) + surface_flux_function = FluxLMARS(RealT(340)) + directions = [1, 2] + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction, + equations)) == RealT + end + end +end +end diff --git a/test/test_unit.jl b/test/test_unit.jl new file mode 100644 index 00000000..8cc426ad --- /dev/null +++ b/test/test_unit.jl @@ -0,0 +1,526 @@ +module TestUnit + +using TrixiAtmo +using Trixi + +include("test_trixiatmo.jl") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@timed_testset "Consistency check for EC flux with Potential Temperature: CEPTE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquations2D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 330.0) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_ec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4]) + u_2d = SVector(u[1], u[2], 0, u[4]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquations1D(c_p = 1004.0, + c_v = 717.0) + equations_2d = equations + flux_1d = normal_1d[1] * flux_ec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_ec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3]) + flux_1d = normal_1d[1] * flux_ec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_ec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquations3D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_ec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_ec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_ec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5]] +end + +@timed_testset "Consistency check for TEC flux with Potential Temperature: CEPTE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquations2D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 330.0) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_tec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4]) + u_2d = SVector(u[1], u[2], 0, u[4]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquations1D(c_p = 1004.0, + c_v = 717.0) + equations_2d = equations + flux_1d = normal_1d[1] * flux_tec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_tec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3]) + flux_1d = normal_1d[1] * flux_tec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_tec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquations3D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_tec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_tec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_tec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5]] +end + +@timed_testset "Consistency check for ETEC flux with Potential Temperature: CEPTE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquations2D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 330.0) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_etec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4]) + u_2d = SVector(u[1], u[2], 0, u[4]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquations1D(c_p = 1004.0, + c_v = 717.0) + equations_2d = equations + flux_1d = normal_1d[1] * flux_etec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_etec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3]) + flux_1d = normal_1d[1] * flux_etec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_etec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquations3D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_etec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_etec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_etec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5]] +end + +@timed_testset "Consistency check for LMARS flux with Potential Temperature: CEPTE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquations2D(c_p = 1004.0, c_v = 717.0) + flux_lmars = FluxLMARS(340) + u = SVector(1.1, -0.5, 2.34, 330.0) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4]) + u_2d = SVector(u[1], u[2], 0, u[4]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquations1D(c_p = 1004.0, + c_v = 717.0) + equations_2d = equations + flux_1d = normal_1d[1] * flux_lmars(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_lmars(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquations3D(c_p = 1004.0, c_v = 717.0) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], -u[2], 0.0, 0.0, u[5]) + u_1d = SVector(u[1], -u[2], u[5]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_lmars(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_lmars(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5]] +end + +@timed_testset "Consistency check for EC flux with Potential Temperature with gravity: CEPTEWG" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_ec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4], u[5]) + u_2d = SVector(u[1], u[2], 0, u[4], u[5]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + equations_2d = equations + flux_1d = normal_1d[1] * flux_ec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_ec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5, 1700) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3], u_rr_1d[4]) + flux_1d = normal_1d[1] * flux_ec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_ec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_ec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5], u[6]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_ec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_ec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5, 6]] +end + +@timed_testset "Consistency check for TEC flux with Potential Temperature with gravity: CEPTEWG" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_tec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4], u[5]) + u_2d = SVector(u[1], u[2], 0, u[4], u[5]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + equations_2d = equations + flux_1d = normal_1d[1] * flux_tec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_tec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5, 1700) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3], u_rr_1d[4]) + flux_1d = normal_1d[1] * flux_tec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_tec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_tec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5], u[6]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_tec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_tec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5, 6]] +end + +@timed_testset "Consistency check for ETEC flux with Potential Temperature with gravity: CEPTEWG" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_etec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4], u[5]) + u_2d = SVector(u[1], u[2], 0, u[4], u[5]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + equations_2d = equations + flux_1d = normal_1d[1] * flux_etec(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_etec(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # test when u_ll is not the same as u_rr + u_rr_1d = SVector(2.1, 0.3, 280.5, 1700) + u_rr_2d = SVector(u_rr_1d[1], u_rr_1d[2], 0.0, u_rr_1d[3], u_rr_1d[4]) + flux_1d = normal_1d[1] * flux_etec(u_1d, u_rr_1d, 1, equations_1d) + flux_2d = flux_etec(u_2d, u_rr_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_etec(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5], u[6]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_etec(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_etec(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5, 6]] +end + +@timed_testset "Consistency check for LMARS flux with Potential Temperature with gravity: CEPTEWG" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + flux_lmars = FluxLMARS(340) + u = SVector(1.1, -0.5, 2.34, 330.0, 1700) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + #check consistency between 1D and 2D EC fluxes + u_1d = SVector(u[1], u[2], u[4], u[5]) + u_2d = SVector(u[1], u[2], 0, u[4], u[5]) + normal_1d = SVector(-0.3) + normal_2d = SVector(normal_1d[1], 0.0) + equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + equations_2d = equations + flux_1d = normal_1d[1] * flux_lmars(u_1d, u_1d, 1, equations_1d) + flux_2d = flux_lmars(u_2d, u_2d, normal_2d, equations_2d) + @test flux_1d ≈ flux(u_1d, normal_1d, equations_1d) + @test flux_1d ≈ flux_2d[[1, 2, 4, 5]] + + # check consistency for 3D EC flux + equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, + c_v = 717.0, + gravity = 9.81) + u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1700) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + # check consistency between 1D and 3D EC fluxes + u_3d = SVector(u[1], u[2], 0.0, 0.0, u[5], u[6]) + normal_3d = SVector(normal_1d[1], 0.0, 0.0) + equations_3d = equations + flux_3d = flux_lmars(u_3d, u_3d, normal_3d, equations_3d) + flux_1d = normal_1d[1] * flux_lmars(u_1d, u_1d, 1, equations_1d) + @test flux_1d ≈ flux_3d[[1, 2, 5, 6]] +end + +@timed_testset "Consistency check for 3D shallow water fluxes: SWE" begin + # Set up equations and dummy conservative variables state + equations = ShallowWaterEquations3D(gravity = 1.0) + u = SVector(1.1, -0.5, 2.34, -3.5, 120.0) + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_wintermeyer_etal(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + for normal_direction in normal_directions + @test flux_fjordholm_etal(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + +@timed_testset "Consistency check for split covariant shallow water fluxes: SWE" begin + # Set up equations and dummy conservative variables state + equations = SplitCovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, + EARTH_ROTATION_RATE) + u = SVector(1.1, -0.5, 2.34) + aux_vars = SVector{26}(ones(26)) + orientation = 1 + @test flux_ec(u, u, aux_vars, aux_vars, orientation, equations) ≈ + flux(u, aux_vars, orientation, equations) +end + +end