From fb4dad77eceb955ba5e5397da17dee7744662d1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 14:18:26 +0100 Subject: [PATCH 01/19] Added non-conservative 2D Euler equations with gravity potential --- src/TrixiAtmo.jl | 4 +++- src/equations/equations.jl | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 2133d518..b9311b47 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -32,12 +32,14 @@ include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D, - SplitCovariantShallowWaterEquations2D + SplitCovariantShallowWaterEquations2D, CompressibleEulerEquationsWithGravity2D export GlobalCartesianCoordinates, GlobalSphericalCoordinates export flux_chandrashekar, FluxLMARS +export flux_nonconservative_waruszewski + export flux_nonconservative_zeros, flux_nonconservative_ec, source_terms_geometric_coriolis diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 53618a94..7fabd517 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -326,4 +326,5 @@ include("covariant_shallow_water_split.jl") include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") include("reference_data.jl") +include("compressible_euler_gravity_2d.jl") end # @muladd From 78fd1d0879ff5b6a84c59e5301045773c0756c24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 14:29:27 +0100 Subject: [PATCH 02/19] Added equation... --- .../compressible_euler_gravity_2d.jl | 1087 +++++++++++++++++ 1 file changed, 1087 insertions(+) create mode 100644 src/equations/compressible_euler_gravity_2d.jl diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl new file mode 100644 index 00000000..dd5be4b5 --- /dev/null +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -0,0 +1,1087 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent +@doc raw""" + CompressibleEulerEquationsWithGravity2D(gamma) + +The compressible Euler equations with gravity, +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \\ \Phi +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \\ 0 +\end{pmatrix} ++ +\frac{\partial}{\partial y} +\begin{pmatrix} +\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \\ 0 +\end{pmatrix} += +\begin{pmatrix} +0 \\ - \rho \nabla \Phi \\ 0 \\ 0 +\end{pmatrix} +``` +for an ideal gas with ratio of specific heat `gamma` in two space dimensions. + +Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific total energy, +which includes the internal, kinetik, and geopotential energy, `\Phi` is the gravitational +geopotential, and +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) - \rho \Phi \right) +``` +the pressure. + +References: +- Souza, A. N., He, J., Bischoff, T., Waruszewski, M., Novak, L., Barra, V., ... & Schneider, T. (2023). The flux‐differencing discontinuous galerkin method applied to an idealized fully compressible nonhydrostatic dry atmosphere. Journal of Advances in Modeling Earth Systems, 15(4), e2022MS003527. https://doi.org/10.1029/2022MS003527. +- Waruszewski, M., Kozdon, J. E., Wilcox, L. C., Gibson, T. H., & Giraldo, F. X. (2022). Entropy stable discontinuous Galerkin methods for balance laws in non-conservative form: Applications to the Euler equations with gravity. Journal of Computational Physics, 468, 111507. https://doi.org/10.1016/j.jcp.2022.111507. +""" +struct CompressibleEulerEquationsWithGravity2D{RealT <: Real} <: + Trixi.AbstractCompressibleEulerEquations{2, 5} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + function CompressibleEulerEquationsWithGravity2D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravity2D) = True() +Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravity2D) = ("rho", + "rho_v1", + "rho_v2", + "rho_e", + "phi") +Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravity2D) = ("rho", + "v1", + "v2", + "p", + "phi") + +# Set initial conditions at physical location `x` for time `t` +# """ +# initial_condition_constant(x, t, equations::CompressibleEulerEquationsWithGravity2D) + +# A constant initial condition to test free-stream preservation. +# """ +# function initial_condition_constant(x, t, equations::CompressibleEulerEquationsWithGravity2D) +# rho = 1.0 +# rho_v1 = 0.1 +# rho_v2 = -0.2 +# rho_e = 10.0 +# return SVector(rho, rho_v1, rho_v2, rho_e) +# end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + +Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book +- Eleuterio F. Toro (2009) + Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + 3rd edition + [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + +Should be used together with [`UnstructuredMesh2D`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_local = Trixi.rotate_to_x(u_inner, normal, equations) + + # compute the primitive variables + rho_local, v_normal, v_tangent, p_local, _ = cons2prim(u_local, equations) + + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Pratical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0.0 + A = 2 / ((equations.gamma + 1) * rho_local) + B = p_local * (equations.gamma - 1) / (equations.gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return (SVector(zero(eltype(u_inner)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_, + SVector(zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_) +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquationsWithGravity2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(1, 0) + else # orientation == 2 + normal_direction = SVector(0, 1) + end + + # compute and return the flux using `boundary_condition_slip_wall` routine above + return boundary_condition_slip_wall(u_inner, normal_direction, direction, + x, t, surface_flux_function, equations) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquationsWithGravity2D) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + fluxes = boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, equations) + boundary_flux = (-fluxes[1], -fluxes[2]) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux +end + +# Calculate 2D flux for a single point +@inline function Trixi.flux(u, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, phi = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = (rho_e + p) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = (rho_e + p) * v2 + end + return SVector(f1, f2, f3, f4, zero(eltype(u))) +end + +# Calculate 2D flux for a single point in the normal direction +# Note, this directional vector is not normalized +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + rho_e = u[4] + rho, v1, v2, p, _ = cons2prim(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + f4 = (rho_e + p) * v_normal + return SVector(f1, f2, f3, f4, zero(eltype(u))) +end + +""" + flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + +This flux is is a modification of the original kinetic energy preserving two-point flux by +- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) + Kinetic energy and entropy preserving schemes for compressible flows + by split convective forms + [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058) + +The modification is in the energy flux to guarantee pressure equilibrium and was developed by +- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) + Preventing spurious pressure oscillations in split convective form discretizations for + compressible flows + [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060) +""" +@inline function Trixi.flux_shima_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Average each factor of products in flux + rho_avg = 1 / 2 * (rho_ll + rho_rr) + v1_avg = 1 / 2 * (v1_ll + v1_rr) + v2_avg = 1 / 2 * (v2_ll + v2_rr) + p_avg = 1 / 2 * (p_ll + p_rr) + kin_avg = 1 / 2 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + pv1_avg = 1 / 2 * (p_ll * v1_rr + p_rr * v1_ll) + f1 = rho_avg * v1_avg + f2 = f1 * v1_avg + p_avg + f3 = f1 * v2_avg + f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + else + pv2_avg = 1 / 2 * (p_ll * v2_rr + p_rr * v2_ll) + f1 = rho_avg * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + f4 = p_avg * v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Average each factor of products in flux + rho_avg = 1 / 2 * (rho_ll + rho_rr) + v1_avg = 1 / 2 * (v1_ll + v1_rr) + v2_avg = 1 / 2 * (v2_ll + v2_rr) + v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) + p_avg = 1 / 2 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_avg * v_dot_n_avg + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = (f1 * velocity_square_avg + + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + +Kinetic energy preserving two-point flux by +- Kennedy and Gruber (2008) + Reduced aliasing formulations of the convective terms within the + Navier-Stokes equations for a compressible fluid + [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020) +""" +@inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_e_ll = u_ll[4] + rho_e_rr = u_rr[4] + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Average each factor of products in flux + rho_avg = 1 / 2 * (rho_ll + rho_rr) + v1_avg = 1 / 2 * (v1_ll + v1_rr) + v2_avg = 1 / 2 * (v2_ll + v2_rr) + p_avg = 1 / 2 * (p_ll + p_rr) + e_avg = 1 / 2 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = rho_avg * v1_avg + f2 = rho_avg * v1_avg * v1_avg + p_avg + f3 = rho_avg * v1_avg * v2_avg + f4 = (rho_avg * e_avg + p_avg) * v1_avg + else + f1 = rho_avg * v2_avg + f2 = rho_avg * v2_avg * v1_avg + f3 = rho_avg * v2_avg * v2_avg + p_avg + f4 = (rho_avg * e_avg + p_avg) * v2_avg + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_e_ll = u_ll[4] + rho_e_rr = u_rr[4] + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Average each factor of products in flux + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + p_avg = 0.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_avg * v_dot_n_avg + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = f1 * e_avg + p_avg * v_dot_n_avg + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) + +Entropy conserving two-point flux by +- Chandrashekar (2013) + Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes + for Compressible Euler and Navier-Stokes Equations + [DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a) +""" +@inline function Trixi.flux_chandrashekar(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2) + specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2) + + # Compute the necessary mean values + rho_avg = 0.5 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_mean = 0.5 * rho_avg / beta_avg + velocity_square_avg = specific_kin_ll + specific_kin_rr + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + p_mean + f3 = f1 * v2_avg + f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + else + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_mean + f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + +Entropy conserving and kinetic energy preserving two-point flux by +- Hendrik Ranocha (2018) + Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods + for Hyperbolic Balance Laws + [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) +See also +- Hendrik Ranocha (2020) + Entropy Conserving and Kinetic Energy Preserving Numerical Methods for + the Euler Equations Using Summation-by-Parts Operators + [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) +""" +@inline function Trixi.flux_ranocha(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` + # in exact arithmetic since + # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) + # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + p_avg + f3 = f1 * v2_avg + f4 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + else + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + f4 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` + # in exact arithmetic since + # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) + # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +function flux_nonconservative_waruszewski(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, phi_rr = u_rr + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) + # noncons = 0.5 * (rho_ll + rho_rr) * (phi_rr - phi_ll) + + f0 = zero(eltype(u_ll)) + return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], + f0, f0) +end + +function flux_nonconservative_waruszewski(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, phi_rr = u_rr + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) + # noncons = 0.5 * (rho_ll + rho_rr) * (phi_rr - phi_ll) + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, noncons, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, noncons, f0, f0) + end +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" +# The struct is already defined in CompressibleEulerEquations2D + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4, _ = v * u_ll + f4 = f4 + p_ll * v + else + f1, f2, f3, f4, _ = v * u_rr + f4 = f4 + p_rr * v + end + + if orientation == 1 + f2 = f2 + p + else # orientation == 2 + f3 = f3 + p + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5f0 * (rho_ll + rho_rr) + p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4, _ = u_ll * v + f4 = f4 + p_ll * v + else + f1, f2, f3, f4, _ = u_rr * v + f4 = f4 + p_rr * v + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4, zero(eltype(u_ll))) +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + 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 Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + 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 + +# Calculate minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + if orientation == 1 # x-direction + λ_min = v1_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v1_rr + sqrt(equations.gamma * p_rr / rho_rr) + else # y-direction + λ_min = v2_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v2_rr + sqrt(equations.gamma * p_rr / rho_rr) + end + + return λ_min, λ_max +end + +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + # The v_normals are already scaled by the norm + λ_min = v_normal_ll - sqrt(equations.gamma * p_ll / rho_ll) * norm_ + λ_max = v_normal_rr + sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + return λ_min, λ_max +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this rotation of the state vector +@inline function Trixi.rotate_to_x(u, normal_vector, + equations::CompressibleEulerEquationsWithGravity2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4], + u[5]) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function Trixi.rotate_from_x(u, normal_vector, + equations::CompressibleEulerEquationsWithGravity2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 t_1 0; + # 0 n_2 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] - s * u[3], + s * u[2] + c * u[3], + u[4], + u[5]) +end + +""" + flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) + +Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro +[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) +Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140) +""" +function Trixi.flux_hllc(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Calculate primitive variables and speed of sound + rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, phi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, phi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + e_ll = rho_e_ll / rho_ll + p_ll = (equations.gamma - 1) * + (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2) - rho_ll * phi_ll) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + e_rr = rho_e_rr / rho_rr + p_rr = (equations.gamma - 1) * + (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2) - rho_rr * phi_rr) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + # Obtain left and right fluxes + f_ll = flux(u_ll, orientation, equations) + f_rr = flux(u_rr, orientation, equations) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + if orientation == 1 # x-direction + vel_L = v1_ll + vel_R = v1_rr + ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 + elseif orientation == 2 # y-direction + vel_L = v2_ll + vel_R = v2_rr + ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + end + vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho + ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) + H_ll = (rho_e_ll + p_ll) / rho_ll + H_rr = (rho_e_rr + p_rr) / rho_rr + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) + Ssl = min(vel_L - c_ll, vel_roe - c_roe) + Ssr = max(vel_R + c_rr, vel_roe + c_roe) + sMu_L = Ssl - vel_L + sMu_R = Ssr - vel_R + + if Ssl >= 0.0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + elseif Ssr <= 0.0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + else + SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / + (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0.0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar4 = densStar * enerStar + if orientation == 1 # x-direction + UStar2 = densStar * SStar + UStar3 = densStar * v2_ll + elseif orientation == 2 # y-direction + UStar2 = densStar * v1_ll + UStar3 = densStar * SStar + end + f1 = f_ll[1] + Ssl * (UStar1 - rho_ll) + f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) + f3 = f_ll[3] + Ssl * (UStar3 - rho_v2_ll) + f4 = f_ll[4] + Ssl * (UStar4 - rho_e_ll) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar4 = densStar * enerStar + if orientation == 1 # x-direction + UStar2 = densStar * SStar + UStar3 = densStar * v2_rr + elseif orientation == 2 # y-direction + UStar2 = densStar * v1_rr + UStar3 = densStar * SStar + end + f1 = f_rr[1] + Ssr * (UStar1 - rho_rr) + f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) + f3 = f_rr[3] + Ssr * (UStar3 - rho_v2_rr) + f4 = f_rr[4] + Ssr * (UStar4 - rho_e_rr) + end + end + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) + +Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. +Special estimates of the signal velocites and linearization of the Riemann problem developed +by Einfeldt to ensure that the internal energy and density remain positive during the computation +of the numerical flux. + +- Bernd Einfeldt (1988) + On Godunov-type methods for gas dynamics. + [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021) +- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991) + On Godunov-type methods near low densities. + [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) +""" +function Trixi.flux_hlle(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Calculate primitive variables, enthalpy and speed of sound + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # `u_ll[4]` is total energy `rho_e_ll` on the left + H_ll = (u_ll[4] + p_ll) / rho_ll + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + + # `u_rr[4]` is total energy `rho_e_rr` on the right + H_rr = (u_rr[4] + p_rr) / rho_rr + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr) + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho + v_roe_mag = v1_roe^2 + v2_roe^2 + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) + + # Compute convenience constant for positivity preservation, see + # https://doi.org/10.1016/0021-9991(91)90211-3 + beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + + # Estimate the edges of the Riemann fan (with positivity conservation) + if orientation == 1 # x-direction + SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(v1_roe)) + SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(v1_roe)) + elseif orientation == 2 # y-direction + SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, zero(v2_roe)) + SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe)) + end + + if SsL >= 0.0 && SsR > 0.0 + # Positive supersonic speed + f_ll = flux(u_ll, orientation, equations) + + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + elseif SsR <= 0.0 && SsL < 0.0 + # Negative supersonic speed + f_rr = flux(u_rr, orientation, equations) + + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + else + # Subsonic case + # Compute left and right fluxes + f_ll = flux(u_ll, orientation, equations) + f_rr = flux(u_rr, orientation, equations) + + f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) / + (SsR - SsL) + f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) / + (SsR - SsL) + f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) / + (SsR - SsL) + f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) / + (SsR - SsL) + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function Trixi.max_abs_speeds(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, v1, v2, p, _ = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, phi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + + return SVector(rho, v1, v2, p, phi) +end + +# Convert conservative variables to entropy (see, e.g., Waruszewski et al. (2022)) +@inline function Trixi.cons2entropy(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, phi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square - rho * phi) + s = log(p) - equations.gamma * log(rho) + rho_p = rho / p + + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + rho_p * (0.5 * v_square - phi) + w2 = rho_p * v1 + w3 = rho_p * v2 + w4 = -rho_p + + return SVector(w1, w2, w3, w4, phi) +end + +@inline function Trixi.entropy2cons(w, + equations::CompressibleEulerEquationsWithGravity2D) + # See Waruszewski et al. (2022) + @unpack gamma = equations + + # convert to entropy `-rho * s` + # instead of `-rho * s / (gamma - 1)` + V1, V2, V3, V5, phi = w .* (gamma - 1) + + # s = specific entropy + s = gamma - V1 + (V2^2 + V3^2) / (2 * V5) - V5 * phi + + rho_iota = ((gamma - 1) / (-V5)^gamma)^(equations.inv_gamma_minus_one) * + exp(-s * equations.inv_gamma_minus_one) + + rho = -rho_iota * V5 + rho_v1 = rho_iota * V2 + rho_v2 = rho_iota * V3 + rho_e = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5)) + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_e, w[5]) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, + equations::CompressibleEulerEquationsWithGravity2D) + rho, v1, v2, p, phi = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_e = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_e, phi) +end + +@inline function Trixi.density(u, equations::CompressibleEulerEquationsWithGravity2D) + rho = u[1] + return rho +end + +@inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, phi = u + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) + return p +end + +@inline function Trixi.density_pressure(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, phi = u + rho_times_p = (equations.gamma - 1) * + (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) + return rho_times_p +end + +# Calculate thermodynamic entropy for a conservative state `cons` +@inline function Trixi.entropy_thermodynamic(cons, + equations::CompressibleEulerEquationsWithGravity2D) + # Pressure + p = (equations.gamma - 1) * + (cons[4] - 1 / 2 * (cons[2]^2 + cons[3]^2) / cons[1] - cons[1] * cons[5]) + + # Thermodynamic entropy + s = log(p) - equations.gamma * log(cons[1]) + + return s +end + +# Calculate mathematical entropy for a conservative state `cons` +@inline function Trixi.entropy_math(cons, + equations::CompressibleEulerEquationsWithGravity2D) + # Mathematical entropy + S = -entropy_thermodynamic(cons, equations) * cons[1] * + equations.inv_gamma_minus_one + + return S +end + +# Default entropy is the mathematical entropy +@inline Trixi.entropy(cons, equations::CompressibleEulerEquationsWithGravity2D) = entropy_math(cons, + equations) + +# Calculate total energy for a conservative state `cons` +@inline Trixi.energy_total(cons, ::CompressibleEulerEquationsWithGravity2D) = cons[4] + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_e, _ = u + return (rho_v1^2 + rho_v2^2) / (2 * rho) +end + +# Calculate internal energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, + equations::CompressibleEulerEquationsWithGravity2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) - + cons[1] * cons[5] +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) +end +end # @muladd From 57ebc3b9fb3a4776ba1506b3403aed80f4495125 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 14:42:23 +0100 Subject: [PATCH 03/19] Added benchmark tests Co-authored-by: Benedict Geihe --- examples/elixir_euler_gravity_equilibrium.jl | 74 ++++++++++ .../elixir_euler_gravity_schaer_mountains.jl | 135 +++++++++++++++++ examples/elixir_euler_gravity_warmbubble.jl | 136 ++++++++++++++++++ 3 files changed, 345 insertions(+) create mode 100644 examples/elixir_euler_gravity_equilibrium.jl create mode 100644 examples/elixir_euler_gravity_schaer_mountains.jl create mode 100644 examples/elixir_euler_gravity_warmbubble.jl diff --git a/examples/elixir_euler_gravity_equilibrium.jl b/examples/elixir_euler_gravity_equilibrium.jl new file mode 100644 index 00000000..9e1b4b76 --- /dev/null +++ b/examples/elixir_euler_gravity_equilibrium.jl @@ -0,0 +1,74 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravity2D(1.4) + +function initial_condition_constant(x, t, + equations::CompressibleEulerEquationsWithGravity2D) + rho = exp(-x[2]) + v1 = 0.0 + v2 = 0.0 + p = exp(-x[2]) + prim = SVector(rho, v1, v2, p, x[2]) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_constant + +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski) + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10_000, + periodicity = false) + +boundary_conditions = (; x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/examples/elixir_euler_gravity_schaer_mountains.jl b/examples/elixir_euler_gravity_schaer_mountains.jl new file mode 100644 index 00000000..0e3a8e61 --- /dev/null +++ b/examples/elixir_euler_gravity_schaer_mountains.jl @@ -0,0 +1,135 @@ +using OrdinaryDiffEq +using Trixi + +#schaer mountain test case + +# Initial condition +function initial_condition_schaer_mountain(x, t, + equations::CompressibleEulerEquationsWithGravity2D) + g = 9.81 + c_p = 1004.0 + c_v = 717.0 + gamma = c_p / c_v + + # Exner pressure from hydrostatic balance for x[2] + potential_temperature_int = 280.0 #constant of integration + bvfrequency = 0.01 #Brunt-Väisälä frequency + + exner = 1 + + g^2 / (c_p * potential_temperature_int * bvfrequency^2) * + (exp(-bvfrequency^2 / g * x[2]) - 1) + + # mean potential temperature + potential_temperature = potential_temperature_int * exp(bvfrequency^2 / g * x[2]) + + # temperature + T = potential_temperature * exner + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # density + rho = p / (R * T) + + #Geopotential + phi = g * x[2] + + v1 = 10.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + phi + return SVector(rho, rho * v1, rho * v2, rho * E, phi) +end + +############################################################################### + +function mapping(xi_, eta_) + xi = xi_ * 25_000 #xi_ * 10_000.0 + eta = eta_ * 10_500 + 10_500# eta_ * 500.0 + 500.0 + #upper boundary + H = 21_000.0 + + #topography + h_c = 250.0 + lambda_c = 4000.0 + a_c = 5000.0 + + topo = -h_c * exp(-(xi / a_c)^2) * cos(pi * xi / lambda_c)^2 + + x = xi + y = H * (eta - topo) / (H - topo) + return SVector(x, y) +end + +# Create curved mesh with 200 x 100 elements +polydeg = 3 +cells_per_dimension = (60, 30) +mesh = P4estMesh(cells_per_dimension; polydeg = polydeg, mapping = mapping, + periodicity = false) + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravity2D(1004.0 / 717.0) + +initial_condition = initial_condition_schaer_mountain + +boundary_conditions_dirichlet = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_schaer_mountain), + :x_pos => BoundaryConditionDirichlet(initial_condition_schaer_mountain), + :y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall) + +basis = LobattoLegendreBasis(polydeg) + +volume_flux = (flux_ranocha, flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-25_000.0, 0.0) +coordinates_max = (25_000.0, 21_000.0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions_dirichlet) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 60 * 60)# * 10) # 10h = 36000 s + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +solution_variables = cons2prim + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = solution_variables) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/elixir_euler_gravity_warmbubble.jl b/examples/elixir_euler_gravity_warmbubble.jl new file mode 100644 index 00000000..9e5dcdb6 --- /dev/null +++ b/examples/elixir_euler_gravity_warmbubble.jl @@ -0,0 +1,136 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquationsWithGravity2D) + RealT = eltype(x) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000 + center_z = 2000 + # radius of perturbation + radius = 2000 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300 + potential_temperature_perturbation = zero(RealT) + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5f0 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + # T = potential_temperature - g / (c_p) * x[2] + + # density + rho = p / (R * T) + + # Geopotential + phi = g * x[2] + + v1 = 20 + v2 = 0 + E = c_v * T + 0.5f0 * (v1^2 + v2^2) + phi + return SVector(rho, rho * v1, rho * v2, rho * E, phi) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma) + +initial_condition = warm_bubble_setup + +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, -5000.0) +coordinates_max = (20_000.0, 15_000.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000, + periodicity = (true, false)) + +boundary_conditions = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 10.0, #interval = 1, #dt = 10.0, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); From 7a7d194b39c6f78e46e915d1f567aad2b3209cd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 14:43:40 +0100 Subject: [PATCH 04/19] Added test module and equilibrium tests --- examples/elixir_euler_gravity_equilibrium.jl | 2 +- test/runtests.jl | 4 ++ test/test_2d_euler_gravity.jl | 59 ++++++++++++++++++++ 3 files changed, 64 insertions(+), 1 deletion(-) create mode 100644 test/test_2d_euler_gravity.jl diff --git a/examples/elixir_euler_gravity_equilibrium.jl b/examples/elixir_euler_gravity_equilibrium.jl index 9e1b4b76..23d587eb 100644 --- a/examples/elixir_euler_gravity_equilibrium.jl +++ b/examples/elixir_euler_gravity_equilibrium.jl @@ -48,7 +48,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 1 +analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, analysis_polydeg = polydeg) diff --git a/test/runtests.jl b/test/runtests.jl index 8e65de24..2f1ef05e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,10 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_2d_moist_euler.jl") end + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "euler_gravity" + include("test_2d_euler_gravity.jl") + end + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" include("test_spherical_advection.jl") end diff --git a/test/test_2d_euler_gravity.jl b/test/test_2d_euler_gravity.jl new file mode 100644 index 00000000..f12531f3 --- /dev/null +++ b/test/test_2d_euler_gravity.jl @@ -0,0 +1,59 @@ +module TestSphericalAdvection + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "elixir_euler_gravity_equilibrium (Shima + LLF)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixiatmo_testset "elixir_euler_gravity_equilibrium (Kennedy-Gruber)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0], + volume_flux=(flux_kennedy_gruber, flux_nonconservative_waruszewski), + surface_flux=(flux_kennedy_gruber, flux_nonconservative_waruszewski)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixiatmo_testset "elixir_euler_gravity_equilibrium (Ranocha)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0], + volume_flux=(flux_ranocha, flux_nonconservative_waruszewski), + surface_flux=(flux_ranocha, flux_nonconservative_waruszewski)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +end # module From 0028e0843919c33f87e34b40a291f808c303eee9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 14:47:46 +0100 Subject: [PATCH 05/19] Fixed typo --- src/equations/compressible_euler_gravity_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index dd5be4b5..5202ecd8 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -119,7 +119,7 @@ Should be used together with [`UnstructuredMesh2D`](@ref). # Get the solution of the pressure Riemann problem # See Section 6.3.3 of # Eleuterio F. Toro (2009) - # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Pratical Introduction + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) if v_normal <= 0.0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed From 8f35c603f9f4b0b706eea0693f748963e1f0d927 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 18 Mar 2025 15:34:20 +0100 Subject: [PATCH 06/19] Fixed bug in entropy2cons --- src/equations/compressible_euler_gravity_2d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 5202ecd8..977fa7cc 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -986,7 +986,8 @@ end # convert to entropy `-rho * s` # instead of `-rho * s / (gamma - 1)` - V1, V2, V3, V5, phi = w .* (gamma - 1) + V1, V2, V3, V5, _ = w .* (gamma - 1) + phi = w[5] # s = specific entropy s = gamma - V1 + (V2^2 + V3^2) / (2 * V5) - V5 * phi @@ -998,7 +999,7 @@ end rho_v1 = rho_iota * V2 rho_v2 = rho_iota * V3 rho_e = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5)) + rho * phi - return SVector(rho, rho_v1, rho_v2, rho_e, w[5]) + return SVector(rho, rho_v1, rho_v2, rho_e, phi) end # Convert primitive to conservative variables From e521782bd44b215f6024dd3a1b09121747ef9b9e Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 19 Mar 2025 14:33:02 +0100 Subject: [PATCH 07/19] try to make spell check happy --- .typos.toml | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .typos.toml diff --git a/.typos.toml b/.typos.toml new file mode 100644 index 00000000..97e44980 --- /dev/null +++ b/.typos.toml @@ -0,0 +1,2 @@ +[default.extend-words] +claus = "claus" From bfbdb237d2ed40494758981678775318c21e49ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 19 Mar 2025 14:54:35 +0100 Subject: [PATCH 08/19] Schaer mountains test runs with Kennedy-Gruber, not with Ranocha nor Shima --- examples/elixir_euler_gravity_schaer_mountains.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/elixir_euler_gravity_schaer_mountains.jl b/examples/elixir_euler_gravity_schaer_mountains.jl index 0e3a8e61..72f4f3be 100644 --- a/examples/elixir_euler_gravity_schaer_mountains.jl +++ b/examples/elixir_euler_gravity_schaer_mountains.jl @@ -81,7 +81,7 @@ boundary_conditions_dirichlet = Dict(:x_neg => BoundaryConditionDirichlet(initia basis = LobattoLegendreBasis(polydeg) -volume_flux = (flux_ranocha, flux_nonconservative_waruszewski) +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) volume_integral = VolumeIntegralFluxDifferencing(volume_flux) From 4a538d8a21db198eadfb4e7a2632b9e0e713fa00 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 20 Mar 2025 11:39:21 +0100 Subject: [PATCH 09/19] minor stuff --- examples/elixir_euler_gravity_warmbubble.jl | 4 ++-- src/equations/compressible_euler_gravity_2d.jl | 11 ++++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/elixir_euler_gravity_warmbubble.jl b/examples/elixir_euler_gravity_warmbubble.jl index 9e5dcdb6..9d4e82ca 100644 --- a/examples/elixir_euler_gravity_warmbubble.jl +++ b/examples/elixir_euler_gravity_warmbubble.jl @@ -1,6 +1,6 @@ -using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK -using Trixi +using OrdinaryDiffEqLowStorageRK +using Trixi, TrixiAtmo # Warm bubble test case from # - Wicker, L. J., and Skamarock, W. C. (1998) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 977fa7cc..303e3fb7 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -555,8 +555,8 @@ References: c = flux_lmars.speed_of_sound # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) if orientation == 1 v_ll = v1_ll @@ -594,8 +594,8 @@ end c = flux_lmars.speed_of_sound # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] @@ -1076,7 +1076,8 @@ end cons[1] * cons[5] end -# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the +# gravitational potential @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsWithGravity2D) From 300063577c6b419db51456328295dbed82e4d8c8 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Apr 2025 10:32:51 +0200 Subject: [PATCH 10/19] add velocity getter --- src/equations/compressible_euler_gravity_2d.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 303e3fb7..1a6ad8dc 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -1076,6 +1076,13 @@ end cons[1] * cons[5] end +@inline function Trixi.velocity(u, equations::CompressibleEulerEquationsWithGravity2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + return SVector(v1, v2) +end + # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the # gravitational potential @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, From 5e27bd9261853a78218551049364b5e972da4a2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 9 Mar 2026 16:30:09 +0100 Subject: [PATCH 11/19] Remove comment --- src/equations/compressible_euler_gravity_2d.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 1a6ad8dc..9f4296e0 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -511,7 +511,6 @@ function flux_nonconservative_waruszewski(u_ll, u_rr, normal_direction::Abstract # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) - # noncons = 0.5 * (rho_ll + rho_rr) * (phi_rr - phi_ll) f0 = zero(eltype(u_ll)) return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], @@ -525,7 +524,6 @@ function flux_nonconservative_waruszewski(u_ll, u_rr, orientation::Integer, # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) - # noncons = 0.5 * (rho_ll + rho_rr) * (phi_rr - phi_ll) f0 = zero(eltype(u_ll)) if orientation == 1 From 8c065859094fe3bae4a5c6f54a17ef4c719ce9a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 09:53:26 +0100 Subject: [PATCH 12/19] Cleaned up the implementation a bit, renamed total energy variables, and removed unused functions --- .../compressible_euler_gravity_2d.jl | 281 +++--------------- 1 file changed, 38 insertions(+), 243 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 9f4296e0..69a8f956 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -7,21 +7,21 @@ @doc raw""" CompressibleEulerEquationsWithGravity2D(gamma) -The compressible Euler equations with gravity, +The compressible Euler equations with gravity and total energy, ```math \frac{\partial}{\partial t} \begin{pmatrix} -\rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \\ \Phi +\rho \\ \rho v_1 \\ \rho v_2 \\ \rho e_{tot} \\ \Phi \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} - \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \\ 0 + \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e_{tot} +p) v_1 \\ 0 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} -\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \\ 0 + \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e_{tot} +p) v_2 \\ 0 \end{pmatrix} = \begin{pmatrix} @@ -30,7 +30,7 @@ The compressible Euler equations with gravity, ``` for an ideal gas with ratio of specific heat `gamma` in two space dimensions. -Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific total energy, +Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e_{tot}`` the specific total energy, which includes the internal, kinetik, and geopotential energy, `\Phi` is the gravitational geopotential, and ```math @@ -57,7 +57,7 @@ Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravity2D) = Tr Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravity2D) = ("rho", "rho_v1", "rho_v2", - "rho_e", + "rho_etot", "phi") Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravity2D) = ("rho", "v1", @@ -65,20 +65,6 @@ Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravity2D) = "p", "phi") -# Set initial conditions at physical location `x` for time `t` -# """ -# initial_condition_constant(x, t, equations::CompressibleEulerEquationsWithGravity2D) - -# A constant initial condition to test free-stream preservation. -# """ -# function initial_condition_constant(x, t, equations::CompressibleEulerEquationsWithGravity2D) -# rho = 1.0 -# rho_v1 = 0.1 -# rho_v2 = -0.2 -# rho_e = 10.0 -# return SVector(rho, rho_v1, rho_v2, rho_e) -# end - """ boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, equations::CompressibleEulerEquationsWithGravity2D) @@ -199,20 +185,21 @@ end # Calculate 2D flux for a single point @inline function Trixi.flux(u, orientation::Integer, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, phi = u + rho, rho_v1, rho_v2, rho_etot, phi = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + p = (equations.gamma - 1) * + (rho_etot - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) if orientation == 1 f1 = rho_v1 f2 = rho_v1 * v1 + p f3 = rho_v1 * v2 - f4 = (rho_e + p) * v1 + f4 = (rho_etot + p) * v1 else f1 = rho_v2 f2 = rho_v2 * v1 f3 = rho_v2 * v2 + p - f4 = (rho_e + p) * v2 + f4 = (rho_etot + p) * v2 end return SVector(f1, f2, f3, f4, zero(eltype(u))) end @@ -221,7 +208,7 @@ end # Note, this directional vector is not normalized @inline function Trixi.flux(u, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravity2D) - rho_e = u[4] + rho_etot = u[4] rho, v1, v2, p, _ = cons2prim(u, equations) v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] @@ -229,7 +216,7 @@ end f1 = rho_v_normal f2 = rho_v_normal * v1 + p * normal_direction[1] f3 = rho_v_normal * v2 + p * normal_direction[2] - f4 = (rho_e + p) * v_normal + f4 = (rho_etot + p) * v_normal return SVector(f1, f2, f3, f4, zero(eltype(u))) end @@ -320,8 +307,8 @@ Kinetic energy preserving two-point flux by @inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_e_ll = u_ll[4] - rho_e_rr = u_rr[4] + rho_etot_ll = u_ll[4] + rho_etot_rr = u_rr[4] rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) @@ -330,19 +317,19 @@ Kinetic energy preserving two-point flux by v1_avg = 1 / 2 * (v1_ll + v1_rr) v2_avg = 1 / 2 * (v2_ll + v2_rr) p_avg = 1 / 2 * (p_ll + p_rr) - e_avg = 1 / 2 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + etot_avg = 1 / 2 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) # Calculate fluxes depending on orientation if orientation == 1 f1 = rho_avg * v1_avg f2 = rho_avg * v1_avg * v1_avg + p_avg f3 = rho_avg * v1_avg * v2_avg - f4 = (rho_avg * e_avg + p_avg) * v1_avg + f4 = (rho_avg * etot_avg + p_avg) * v1_avg else f1 = rho_avg * v2_avg f2 = rho_avg * v2_avg * v1_avg f3 = rho_avg * v2_avg * v2_avg + p_avg - f4 = (rho_avg * e_avg + p_avg) * v2_avg + f4 = (rho_avg * etot_avg + p_avg) * v2_avg end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) @@ -351,8 +338,8 @@ end @inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_e_ll = u_ll[4] - rho_e_rr = u_rr[4] + rho_etot_ll = u_ll[4] + rho_etot_rr = u_rr[4] rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) @@ -362,13 +349,13 @@ end v2_avg = 0.5 * (v2_ll + v2_rr) v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] p_avg = 0.5 * (p_ll + p_rr) - e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + etot_avg = 0.5 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] - f4 = f1 * e_avg + p_avg * v_dot_n_avg + f4 = f1 * etot_avg + p_avg * v_dot_n_avg return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end @@ -743,200 +730,6 @@ end u[5]) end -""" - flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) - -Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro -[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) -Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140) -""" -function Trixi.flux_hllc(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquationsWithGravity2D) - # Calculate primitive variables and speed of sound - rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, phi_ll = u_ll - rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, phi_rr = u_rr - - v1_ll = rho_v1_ll / rho_ll - v2_ll = rho_v2_ll / rho_ll - e_ll = rho_e_ll / rho_ll - p_ll = (equations.gamma - 1) * - (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2) - rho_ll * phi_ll) - c_ll = sqrt(equations.gamma * p_ll / rho_ll) - - v1_rr = rho_v1_rr / rho_rr - v2_rr = rho_v2_rr / rho_rr - e_rr = rho_e_rr / rho_rr - p_rr = (equations.gamma - 1) * - (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2) - rho_rr * phi_rr) - c_rr = sqrt(equations.gamma * p_rr / rho_rr) - - # Obtain left and right fluxes - f_ll = flux(u_ll, orientation, equations) - f_rr = flux(u_rr, orientation, equations) - - # Compute Roe averages - sqrt_rho_ll = sqrt(rho_ll) - sqrt_rho_rr = sqrt(rho_rr) - sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr - if orientation == 1 # x-direction - vel_L = v1_ll - vel_R = v1_rr - ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 - elseif orientation == 2 # y-direction - vel_L = v2_ll - vel_R = v2_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 - end - vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) - H_ll = (rho_e_ll + p_ll) / rho_ll - H_rr = (rho_e_rr + p_rr) / rho_rr - H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) - Ssl = min(vel_L - c_ll, vel_roe - c_roe) - Ssr = max(vel_R + c_rr, vel_roe + c_roe) - sMu_L = Ssl - vel_L - sMu_R = Ssr - vel_R - - if Ssl >= 0.0 - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - f4 = f_ll[4] - elseif Ssr <= 0.0 - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - f4 = f_rr[4] - else - SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / - (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0.0 <= SStar - densStar = rho_ll * sMu_L / (Ssl - SStar) - enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) - UStar1 = densStar - UStar4 = densStar * enerStar - if orientation == 1 # x-direction - UStar2 = densStar * SStar - UStar3 = densStar * v2_ll - elseif orientation == 2 # y-direction - UStar2 = densStar * v1_ll - UStar3 = densStar * SStar - end - f1 = f_ll[1] + Ssl * (UStar1 - rho_ll) - f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) - f3 = f_ll[3] + Ssl * (UStar3 - rho_v2_ll) - f4 = f_ll[4] + Ssl * (UStar4 - rho_e_ll) - else - densStar = rho_rr * sMu_R / (Ssr - SStar) - enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) - UStar1 = densStar - UStar4 = densStar * enerStar - if orientation == 1 # x-direction - UStar2 = densStar * SStar - UStar3 = densStar * v2_rr - elseif orientation == 2 # y-direction - UStar2 = densStar * v1_rr - UStar3 = densStar * SStar - end - f1 = f_rr[1] + Ssr * (UStar1 - rho_rr) - f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) - f3 = f_rr[3] + Ssr * (UStar3 - rho_v2_rr) - f4 = f_rr[4] + Ssr * (UStar4 - rho_e_rr) - end - end - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - -""" - flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) - -Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. -Special estimates of the signal velocites and linearization of the Riemann problem developed -by Einfeldt to ensure that the internal energy and density remain positive during the computation -of the numerical flux. - -- Bernd Einfeldt (1988) - On Godunov-type methods for gas dynamics. - [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021) -- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991) - On Godunov-type methods near low densities. - [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) -""" -function Trixi.flux_hlle(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquationsWithGravity2D) - # Calculate primitive variables, enthalpy and speed of sound - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) - - # `u_ll[4]` is total energy `rho_e_ll` on the left - H_ll = (u_ll[4] + p_ll) / rho_ll - c_ll = sqrt(equations.gamma * p_ll / rho_ll) - - # `u_rr[4]` is total energy `rho_e_rr` on the right - H_rr = (u_rr[4] + p_rr) / rho_rr - c_rr = sqrt(equations.gamma * p_rr / rho_rr) - - # Compute Roe averages - sqrt_rho_ll = sqrt(rho_ll) - sqrt_rho_rr = sqrt(rho_rr) - inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr) - - v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho - v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho - v_roe_mag = v1_roe^2 + v2_roe^2 - - H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) - - # Compute convenience constant for positivity preservation, see - # https://doi.org/10.1016/0021-9991(91)90211-3 - beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) - - # Estimate the edges of the Riemann fan (with positivity conservation) - if orientation == 1 # x-direction - SsL = min(v1_roe - c_roe, v1_ll - beta * c_ll, zero(v1_roe)) - SsR = max(v1_roe + c_roe, v1_rr + beta * c_rr, zero(v1_roe)) - elseif orientation == 2 # y-direction - SsL = min(v2_roe - c_roe, v2_ll - beta * c_ll, zero(v2_roe)) - SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe)) - end - - if SsL >= 0.0 && SsR > 0.0 - # Positive supersonic speed - f_ll = flux(u_ll, orientation, equations) - - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - f4 = f_ll[4] - elseif SsR <= 0.0 && SsL < 0.0 - # Negative supersonic speed - f_rr = flux(u_rr, orientation, equations) - - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - f4 = f_rr[4] - else - # Subsonic case - # Compute left and right fluxes - f_ll = flux(u_ll, orientation, equations) - f_rr = flux(u_rr, orientation, equations) - - f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) / - (SsR - SsL) - f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) / - (SsR - SsL) - f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) / - (SsR - SsL) - f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) / - (SsR - SsL) - end - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - @inline function Trixi.max_abs_speeds(u, equations::CompressibleEulerEquationsWithGravity2D) rho, v1, v2, p, _ = cons2prim(u, equations) @@ -947,11 +740,12 @@ end # Convert conservative variables to primitive @inline function Trixi.cons2prim(u, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, phi = u + rho, rho_v1, rho_v2, rho_etot, phi = u v1 = rho_v1 / rho v2 = rho_v2 / rho - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + p = (equations.gamma - 1) * + (rho_etot - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) return SVector(rho, v1, v2, p, phi) end @@ -959,12 +753,12 @@ end # Convert conservative variables to entropy (see, e.g., Waruszewski et al. (2022)) @inline function Trixi.cons2entropy(u, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, phi = u + rho, rho_v1, rho_v2, rho_etot, phi = u v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - p = (equations.gamma - 1) * (rho_e - 0.5 * rho * v_square - rho * phi) + p = (equations.gamma - 1) * (rho_etot - 0.5 * rho * v_square - rho * phi) s = log(p) - equations.gamma * log(rho) rho_p = rho / p @@ -996,8 +790,8 @@ end rho = -rho_iota * V5 rho_v1 = rho_iota * V2 rho_v2 = rho_iota * V3 - rho_e = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5)) + rho * phi - return SVector(rho, rho_v1, rho_v2, rho_e, phi) + rho_etot = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5)) + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_etot, phi) end # Convert primitive to conservative variables @@ -1006,9 +800,9 @@ end rho, v1, v2, p, phi = prim rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_e = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + - rho * phi - return SVector(rho, rho_v1, rho_v2, rho_e, phi) + rho_etot = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_etot, phi) end @inline function Trixi.density(u, equations::CompressibleEulerEquationsWithGravity2D) @@ -1017,16 +811,17 @@ end end @inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, phi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) + rho, rho_v1, rho_v2, rho_etot, phi = u + p = (equations.gamma - 1) * + (rho_etot - 0.5 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) return p end @inline function Trixi.density_pressure(u, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, phi = u + rho, rho_v1, rho_v2, rho_etot, phi = u rho_times_p = (equations.gamma - 1) * - (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) + (rho * rho_etot - 0.5 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) return rho_times_p end @@ -1063,7 +858,7 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function Trixi.energy_kinetic(u, equations::CompressibleEulerEquationsWithGravity2D) - rho, rho_v1, rho_v2, rho_e, _ = u + rho, rho_v1, rho_v2, rho_etot, _ = u return (rho_v1^2 + rho_v2^2) / (2 * rho) end From f521c62f79d8e73c3f4bdd9d67f0f82d182fa869 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 10:12:06 +0100 Subject: [PATCH 13/19] Added single precision compatibility --- .../compressible_euler_gravity_2d.jl | 88 ++++++++++--------- 1 file changed, 45 insertions(+), 43 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 69a8f956..0fdb7162 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -107,17 +107,17 @@ Should be used together with [`UnstructuredMesh2D`](@ref). # Eleuterio F. Toro (2009) # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) - if v_normal <= 0.0 + if v_normal <= 0 sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * - equations.gamma * - equations.inv_gamma_minus_one) + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) else # v_normal > 0.0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end @@ -189,7 +189,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho p = (equations.gamma - 1) * - (rho_etot - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + (rho_etot - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) if orientation == 1 f1 = rho_v1 f2 = rho_v1 * v1 + p @@ -281,7 +281,7 @@ end v2_avg = 1 / 2 * (v2_ll + v2_rr) v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) p_avg = 1 / 2 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -289,7 +289,7 @@ end f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one - + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end @@ -344,12 +344,12 @@ end rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 0.5 * (rho_ll + rho_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] - p_avg = 0.5 * (p_ll + p_rr) - etot_avg = 0.5 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + etot_avg = 0.5f0 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -374,19 +374,19 @@ Entropy conserving two-point flux by # Unpack left and right state rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2) - specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2) # Compute the necessary mean values - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_mean = 0.5 * rho_avg / beta_avg + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_mean = 0.5f0 * rho_avg / beta_avg velocity_square_avg = specific_kin_ll + specific_kin_rr # Calculate fluxes depending on orientation @@ -394,13 +394,15 @@ Entropy conserving two-point flux by f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_mean f3 = f1 * v2_avg - f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f4 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg else f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_mean - f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f4 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + f2 * v1_avg + f3 * v2_avg end @@ -435,10 +437,10 @@ See also # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -447,14 +449,14 @@ See also f3 = f1 * v2_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) else f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) @@ -475,18 +477,18 @@ end # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end @@ -745,7 +747,7 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho p = (equations.gamma - 1) * - (rho_etot - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + (rho_etot - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) return SVector(rho, v1, v2, p, phi) end @@ -758,12 +760,12 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v_square = v1^2 + v2^2 - p = (equations.gamma - 1) * (rho_etot - 0.5 * rho * v_square - rho * phi) + p = (equations.gamma - 1) * (rho_etot - 0.5f0 * rho * v_square - rho * phi) s = log(p) - equations.gamma * log(rho) rho_p = rho / p w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - - rho_p * (0.5 * v_square - phi) + rho_p * (0.5f0 * v_square - phi) w2 = rho_p * v1 w3 = rho_p * v2 w4 = -rho_p @@ -800,7 +802,7 @@ end rho, v1, v2, p, phi = prim rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_etot = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + + rho_etot = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) + rho * phi return SVector(rho, rho_v1, rho_v2, rho_etot, phi) end @@ -813,7 +815,7 @@ end @inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravity2D) rho, rho_v1, rho_v2, rho_etot, phi = u p = (equations.gamma - 1) * - (rho_etot - 0.5 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) + (rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) return p end @@ -821,7 +823,7 @@ end equations::CompressibleEulerEquationsWithGravity2D) rho, rho_v1, rho_v2, rho_etot, phi = u rho_times_p = (equations.gamma - 1) * - (rho * rho_etot - 0.5 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) + (rho * rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) return rho_times_p end @@ -883,7 +885,7 @@ end equations::CompressibleEulerEquationsWithGravity2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) + diss = -0.5f0 * λ * (u_rr - u_ll) return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) end end # @muladd From a471becd9c66cfbaa30bd63b497b974ac75a0c85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 11:11:19 +0100 Subject: [PATCH 14/19] Fixed Ranocha and Shima fluxes, and removed unused Chandrashekar flux --- .../compressible_euler_gravity_2d.jl | 84 +++++-------------- 1 file changed, 21 insertions(+), 63 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index 0fdb7162..b5f81df0 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -239,8 +239,8 @@ The modification is in the energy flux to guarantee pressure equilibrium and was @inline function Trixi.flux_shima_etal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) # Average each factor of products in flux rho_avg = 1 / 2 * (rho_ll + rho_rr) @@ -248,6 +248,7 @@ The modification is in the energy flux to guarantee pressure equilibrium and was v2_avg = 1 / 2 * (v2_ll + v2_rr) p_avg = 1 / 2 * (p_ll + p_rr) kin_avg = 1 / 2 * (v1_ll * v1_rr + v2_ll * v2_rr) + phi_avg = 1 / 2 * (phi_ll + phi_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -255,13 +256,15 @@ The modification is in the energy flux to guarantee pressure equilibrium and was f1 = rho_avg * v1_avg f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg - f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + + f1 * phi_avg else pv2_avg = 1 / 2 * (p_ll * v2_rr + p_rr * v2_ll) f1 = rho_avg * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg - f4 = p_avg * v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg + f4 = p_avg * v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg + + f1 * phi_avg end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) @@ -270,8 +273,8 @@ end @inline function Trixi.flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] @@ -282,6 +285,7 @@ end v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) p_avg = 1 / 2 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) + phi_avg = 1 / 2 * (phi_ll + phi_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -289,7 +293,7 @@ end f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * velocity_square_avg + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one - + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + f1 * phi_avg return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end @@ -360,55 +364,6 @@ end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end -""" - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquationsWithGravity2D) - -Entropy conserving two-point flux by -- Chandrashekar (2013) - Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes - for Compressible Euler and Navier-Stokes Equations - [DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a) -""" -@inline function Trixi.flux_chandrashekar(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquationsWithGravity2D) - # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) - beta_ll = 0.5f0 * rho_ll / p_ll - beta_rr = 0.5f0 * rho_rr / p_rr - specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2) - specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2) - - # Compute the necessary mean values - rho_avg = 0.5f0 * (rho_ll + rho_rr) - rho_mean = ln_mean(rho_ll, rho_rr) - beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5f0 * (beta_ll + beta_rr) - v1_avg = 0.5f0 * (v1_ll + v1_rr) - v2_avg = 0.5f0 * (v2_ll + v2_rr) - p_mean = 0.5f0 * rho_avg / beta_avg - velocity_square_avg = specific_kin_ll + specific_kin_rr - - # Calculate fluxes depending on orientation - if orientation == 1 - f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean - f3 = f1 * v2_avg - f4 = f1 * 0.5f0 * - (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + - f2 * v1_avg + f3 * v2_avg - else - f1 = rho_mean * v2_avg - f2 = f1 * v1_avg - f3 = f1 * v2_avg + p_mean - f4 = f1 * 0.5f0 * - (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + - f2 * v1_avg + f3 * v2_avg - end - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - """ flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsWithGravity2D) @@ -427,8 +382,8 @@ See also @inline function Trixi.flux_ranocha(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) # Compute the necessary mean values rho_mean = ln_mean(rho_ll, rho_rr) @@ -441,6 +396,7 @@ See also v2_avg = 0.5f0 * (v2_ll + v2_rr) p_avg = 0.5f0 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) + phi_avg = 1 / 2 * (phi_ll + phi_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -449,14 +405,14 @@ See also f3 = f1 * v2_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + f1 * phi_avg else f1 = rho_mean * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg f4 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + f1 * phi_avg end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) @@ -465,8 +421,8 @@ end @inline function Trixi.flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravity2D) # Unpack left and right state - rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] @@ -481,6 +437,7 @@ end v2_avg = 0.5f0 * (v2_ll + v2_rr) p_avg = 0.5f0 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) + phi_avg = 1 / 2 * (phi_ll + phi_rr) # Calculate fluxes depending on normal_direction f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) @@ -488,7 +445,8 @@ end f3 = f1 * v2_avg + p_avg * normal_direction[2] f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + f1 * phi_avg) return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end From 35efc3e99a90537789644d92119fb970ac9f8aea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 11:50:36 +0100 Subject: [PATCH 15/19] Update tests --- .../elixir_euler_gravity_warmbubble.jl | 2 +- .../elixir_euler_gravity_schaer_mountains.jl | 17 ++++++++--------- .../elixir_euler_gravity_equilibrium.jl | 6 +++--- test/test_2d_euler_gravity.jl | 19 ++++++++----------- 4 files changed, 20 insertions(+), 24 deletions(-) rename examples/{ => euler/dry_air/buoyancy}/elixir_euler_gravity_warmbubble.jl (98%) rename examples/{ => euler/dry_air/mountain_flow}/elixir_euler_gravity_schaer_mountains.jl (87%) rename examples/{ => euler/dry_air/tests}/elixir_euler_gravity_equilibrium.jl (96%) diff --git a/examples/elixir_euler_gravity_warmbubble.jl b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl similarity index 98% rename from examples/elixir_euler_gravity_warmbubble.jl rename to examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl index 9d4e82ca..fba4cb7a 100644 --- a/examples/elixir_euler_gravity_warmbubble.jl +++ b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl @@ -81,7 +81,7 @@ equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma) initial_condition = warm_bubble_setup volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) -surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_zero) polydeg = 3 solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, diff --git a/examples/elixir_euler_gravity_schaer_mountains.jl b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl similarity index 87% rename from examples/elixir_euler_gravity_schaer_mountains.jl rename to examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl index 72f4f3be..e93a7404 100644 --- a/examples/elixir_euler_gravity_schaer_mountains.jl +++ b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl @@ -1,4 +1,5 @@ -using OrdinaryDiffEq +using OrdinaryDiffEqLowStorageRK +using TrixiAtmo using Trixi #schaer mountain test case @@ -74,15 +75,15 @@ equations = CompressibleEulerEquationsWithGravity2D(1004.0 / 717.0) initial_condition = initial_condition_schaer_mountain -boundary_conditions_dirichlet = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_schaer_mountain), - :x_pos => BoundaryConditionDirichlet(initial_condition_schaer_mountain), - :y_neg => boundary_condition_slip_wall, - :y_pos => boundary_condition_slip_wall) +boundary_conditions_dirichlet = (x_neg = BoundaryConditionDirichlet(initial_condition_schaer_mountain), + x_pos = BoundaryConditionDirichlet(initial_condition_schaer_mountain), + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) basis = LobattoLegendreBasis(polydeg) -volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) -surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_zero) volume_integral = VolumeIntegralFluxDifferencing(volume_flux) @@ -131,5 +132,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), maxiters = 1.0e7, dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); - -summary_callback() diff --git a/examples/elixir_euler_gravity_equilibrium.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl similarity index 96% rename from examples/elixir_euler_gravity_equilibrium.jl rename to examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl index 23d587eb..0703beb3 100644 --- a/examples/elixir_euler_gravity_equilibrium.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl @@ -1,5 +1,5 @@ - -using OrdinaryDiffEq +using OrdinaryDiffEqLowStorageRK +using TrixiAtmo using Trixi ############################################################################### @@ -19,7 +19,7 @@ end initial_condition = initial_condition_constant volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski) -surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski) +surface_flux = (flux_lax_friedrichs, flux_zero) polydeg = 3 solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, diff --git a/test/test_2d_euler_gravity.jl b/test/test_2d_euler_gravity.jl index f12531f3..b122a570 100644 --- a/test/test_2d_euler_gravity.jl +++ b/test/test_2d_euler_gravity.jl @@ -1,14 +1,11 @@ -module TestSphericalAdvection - -using Test -using TrixiAtmo +module Test2DEulerGravity include("test_trixiatmo.jl") -EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") +EXAMPLES_DIR = joinpath(EXAMPLES_DIR, "euler/dry_air") -@trixiatmo_testset "elixir_euler_gravity_equilibrium (Shima + LLF)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, +@trixi_testset "elixir_euler_gravity_equilibrium (Shima + LLF)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", "elixir_euler_gravity_equilibrium.jl"), l2=[0.0, 0.0, 0.0, 0.0, 0.0], linf=[0.0, 0.0, 0.0, 0.0, 0.0]) @@ -22,8 +19,8 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") end end -@trixiatmo_testset "elixir_euler_gravity_equilibrium (Kennedy-Gruber)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, +@trixi_testset "elixir_euler_gravity_equilibrium (Kennedy-Gruber)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", "elixir_euler_gravity_equilibrium.jl"), l2=[0.0, 0.0, 0.0, 0.0, 0.0], linf=[0.0, 0.0, 0.0, 0.0, 0.0], @@ -39,8 +36,8 @@ end end end -@trixiatmo_testset "elixir_euler_gravity_equilibrium (Ranocha)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, +@trixi_testset "elixir_euler_gravity_equilibrium (Ranocha)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", "elixir_euler_gravity_equilibrium.jl"), l2=[0.0, 0.0, 0.0, 0.0, 0.0], linf=[0.0, 0.0, 0.0, 0.0, 0.0], From d96b3e45dc03c05342851214a8ff2c2c28453642 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 12:08:27 +0100 Subject: [PATCH 16/19] Fixed outdated TRIXIATMO_TEST call --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7e747066..ccd84e88 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -69,7 +69,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) end end - @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "euler_gravity" + @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_gravity" include("test_2d_euler_gravity.jl") end From 8dda4bd0ae18df3a6ca66fb3ded5641af3d44046 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 17:38:06 +0100 Subject: [PATCH 17/19] Added type tests, fixed type inconsistencies, and renamed flux_nonconservative_waruszewski_etal --- .../elixir_euler_gravity_warmbubble.jl | 2 +- .../elixir_euler_gravity_schaer_mountains.jl | 2 +- .../tests/elixir_euler_gravity_equilibrium.jl | 2 +- .../compressible_euler_gravity_2d.jl | 67 ++++++++-------- test/test_type.jl | 76 +++++++++++++++++++ 5 files changed, 113 insertions(+), 36 deletions(-) diff --git a/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl index fba4cb7a..3e211abb 100644 --- a/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl +++ b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl @@ -80,7 +80,7 @@ equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma) initial_condition = warm_bubble_setup -volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski_etal) surface_flux = (FluxLMARS(340.0), flux_zero) polydeg = 3 diff --git a/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl index e93a7404..38fb7398 100644 --- a/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl +++ b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl @@ -82,7 +82,7 @@ boundary_conditions_dirichlet = (x_neg = BoundaryConditionDirichlet(initial_cond basis = LobattoLegendreBasis(polydeg) -volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski) +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski_etal) surface_flux = (FluxLMARS(340.0), flux_zero) volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl index 0703beb3..491d3c05 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl @@ -18,7 +18,7 @@ end initial_condition = initial_condition_constant -volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski) +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski_etal) surface_flux = (flux_lax_friedrichs, flux_zero) polydeg = 3 diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index b5f81df0..f4194073 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -113,7 +113,7 @@ Should be used together with [`UnstructuredMesh2D`](@ref). (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * equations.gamma * equations.inv_gamma_minus_one) - else # v_normal > 0.0 + else # v_normal > 0 A = 2 / ((equations.gamma + 1) * rho_local) B = p_local * (equations.gamma - 1) / (equations.gamma + 1) p_star = p_local + @@ -146,9 +146,9 @@ Should be used together with [`TreeMesh`](@ref). equations::CompressibleEulerEquationsWithGravity2D) # get the appropriate normal vector from the orientation if orientation == 1 - normal_direction = SVector(1, 0) + normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner))) else # orientation == 2 - normal_direction = SVector(0, 1) + normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner))) end # compute and return the flux using `boundary_condition_slip_wall` routine above @@ -243,23 +243,23 @@ The modification is in the energy flux to guarantee pressure equilibrium and was rho_rr, v1_rr, v2_rr, p_rr, phi_rr = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - kin_avg = 1 / 2 * (v1_ll * v1_rr + v2_ll * v2_rr) - phi_avg = 1 / 2 * (phi_ll + phi_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + kin_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) + phi_avg = 0.5f0 * (phi_ll + phi_rr) # Calculate fluxes depending on orientation if orientation == 1 - pv1_avg = 1 / 2 * (p_ll * v1_rr + p_rr * v1_ll) + pv1_avg = 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) f1 = rho_avg * v1_avg f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + f1 * phi_avg else - pv2_avg = 1 / 2 * (p_ll * v2_rr + p_rr * v2_ll) + pv2_avg = 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) f1 = rho_avg * v2_avg f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg @@ -279,13 +279,13 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - v_dot_n_avg = 1 / 2 * (v_dot_n_ll + v_dot_n_rr) - p_avg = 1 / 2 * (p_ll + p_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + p_avg = 0.5f0 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) - phi_avg = 1 / 2 * (phi_ll + phi_rr) + phi_avg = 0.5f0 * (phi_ll + phi_rr) # Calculate fluxes depending on normal_direction f1 = rho_avg * v_dot_n_avg @@ -317,11 +317,11 @@ Kinetic energy preserving two-point flux by rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) # Average each factor of products in flux - rho_avg = 1 / 2 * (rho_ll + rho_rr) - v1_avg = 1 / 2 * (v1_ll + v1_rr) - v2_avg = 1 / 2 * (v2_ll + v2_rr) - p_avg = 1 / 2 * (p_ll + p_rr) - etot_avg = 1 / 2 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + etot_avg = 0.5f0 * (rho_etot_ll / rho_ll + rho_etot_rr / rho_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -396,7 +396,7 @@ See also v2_avg = 0.5f0 * (v2_ll + v2_rr) p_avg = 0.5f0 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) - phi_avg = 1 / 2 * (phi_ll + phi_rr) + phi_avg = 0.5f0 * (phi_ll + phi_rr) # Calculate fluxes depending on orientation if orientation == 1 @@ -437,7 +437,7 @@ end v2_avg = 0.5f0 * (v2_ll + v2_rr) p_avg = 0.5f0 * (p_ll + p_rr) velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr) - phi_avg = 1 / 2 * (phi_ll + phi_rr) + phi_avg = 0.5f0 * (phi_ll + phi_rr) # Calculate fluxes depending on normal_direction f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) @@ -451,8 +451,9 @@ end return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end -function flux_nonconservative_waruszewski(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerEquationsWithGravity2D) +function flux_nonconservative_waruszewski_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) rho_ll, _, _, _, phi_ll = u_ll rho_rr, _, _, _, phi_rr = u_rr @@ -464,8 +465,8 @@ function flux_nonconservative_waruszewski(u_ll, u_rr, normal_direction::Abstract f0, f0) end -function flux_nonconservative_waruszewski(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquationsWithGravity2D) +function flux_nonconservative_waruszewski_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) rho_ll, _, _, _, phi_ll = u_ll rho_rr, _, _, _, phi_rr = u_rr @@ -786,11 +787,11 @@ end end # Calculate thermodynamic entropy for a conservative state `cons` -@inline function Trixi.entropy_thermodynamic(cons, - equations::CompressibleEulerEquationsWithGravity2D) +@inline function entropy_thermodynamic(cons, + equations::CompressibleEulerEquationsWithGravity2D) # Pressure p = (equations.gamma - 1) * - (cons[4] - 1 / 2 * (cons[2]^2 + cons[3]^2) / cons[1] - cons[1] * cons[5]) + (cons[4] - 0.5f0 * (cons[2]^2 + cons[3]^2) / cons[1] - cons[1] * cons[5]) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -799,8 +800,8 @@ end end # Calculate mathematical entropy for a conservative state `cons` -@inline function Trixi.entropy_math(cons, - equations::CompressibleEulerEquationsWithGravity2D) +@inline function entropy_math(cons, + equations::CompressibleEulerEquationsWithGravity2D) # Mathematical entropy S = -entropy_thermodynamic(cons, equations) * cons[1] * equations.inv_gamma_minus_one diff --git a/test/test_type.jl b/test/test_type.jl index 617ece0d..718a1021 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -198,6 +198,82 @@ end end end +@timed_testset "Compressible Euler Total Energy With Gravity 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquationsWithGravity2D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), zero(RealT), zero(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] + flux_lmars = FluxLMARS(RealT(340.0)) + + 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} + end + 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 + @test eltype(@inferred flux(u_ll, orientation, + equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_waruszewski_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred flux_lmars(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_nonconservative_waruszewski_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred boundary_condition_slip_wall(u, normal_direction, + x, t, + surface_flux_function, + equations)) == SVector{5, 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), From 9a3ba670d89097522cfcaa5f4adfb533c990d414 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 18 Mar 2026 23:06:49 +0100 Subject: [PATCH 18/19] Format --- test/test_2d_euler_gravity.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_2d_euler_gravity.jl b/test/test_2d_euler_gravity.jl index b122a570..8dc9bc53 100644 --- a/test/test_2d_euler_gravity.jl +++ b/test/test_2d_euler_gravity.jl @@ -24,8 +24,10 @@ end "elixir_euler_gravity_equilibrium.jl"), l2=[0.0, 0.0, 0.0, 0.0, 0.0], linf=[0.0, 0.0, 0.0, 0.0, 0.0], - volume_flux=(flux_kennedy_gruber, flux_nonconservative_waruszewski), - surface_flux=(flux_kennedy_gruber, flux_nonconservative_waruszewski)) + volume_flux=(flux_kennedy_gruber, + flux_nonconservative_waruszewski_etal), + surface_flux=(flux_kennedy_gruber, + flux_nonconservative_waruszewski_etal)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -41,8 +43,8 @@ end "elixir_euler_gravity_equilibrium.jl"), l2=[0.0, 0.0, 0.0, 0.0, 0.0], linf=[0.0, 0.0, 0.0, 0.0, 0.0], - volume_flux=(flux_ranocha, flux_nonconservative_waruszewski), - surface_flux=(flux_ranocha, flux_nonconservative_waruszewski)) + volume_flux=(flux_ranocha, flux_nonconservative_waruszewski_etal), + surface_flux=(flux_ranocha, flux_nonconservative_waruszewski_etal)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 2825d7e00bd8af8ee014265a9f3592a440d83c1b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Mar 2026 17:40:13 +0100 Subject: [PATCH 19/19] Fixed docs typos --- src/equations/compressible_euler_gravity_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl index f4194073..dbd4c343 100644 --- a/src/equations/compressible_euler_gravity_2d.jl +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -31,10 +31,10 @@ The compressible Euler equations with gravity and total energy, for an ideal gas with ratio of specific heat `gamma` in two space dimensions. Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e_{tot}`` the specific total energy, -which includes the internal, kinetik, and geopotential energy, `\Phi` is the gravitational +which includes the internal, kinetik, and geopotential energy, ``\Phi`` is the gravitational geopotential, and ```math -p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) - \rho \Phi \right) +p = (\gamma - 1) \left( \rho e_{tot} - \frac{1}{2} \rho (v_1^2+v_2^2) - \rho \Phi \right) ``` the pressure.