From 2c1b45c954436b2f219220d5c9f4af8c2976b5ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Mar 2026 18:39:42 +0100 Subject: [PATCH 01/18] Added Euler equations with total energy and gravity when peeling of pressure from flux --- ...ompressible_euler_gravity_nopressure_2d.jl | 967 ++++++++++++++++++ 1 file changed, 967 insertions(+) create mode 100644 src/equations/compressible_euler_gravity_nopressure_2d.jl diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl new file mode 100644 index 00000000..a1cd8215 --- /dev/null +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -0,0 +1,967 @@ +# 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""" + CompressibleEulerEquationsWithGravityNoPressure2D(gamma) + +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_{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_{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_{tot} +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_{tot}`` 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_{tot} - \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 CompressibleEulerEquationsWithGravityNoPressure2D{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 CompressibleEulerEquationsWithGravityNoPressure2D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravityNoPressure2D) = True() +Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravityNoPressure2D) = ("rho", + "rho_v1", + "rho_v2", + "rho_etot", + "phi") +Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravityNoPressure2D) = ("rho", + "v1", + "v2", + "p", + "phi") + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + +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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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 Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0 + sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0 + A = 2 / ((equations.gamma + 1) * rho_local) + B = p_local * (equations.gamma - 1) / (equations.gamma + 1) + p_star = p_local + + 0.5f0 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # Compute non-conservative term: Since the geopotential is a continuous function, it does not act at the BC + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + # TODO: compute with surface_flux_function + # surface_flux, surface_noncons = surface_flux_function + # flux_nonconservative_waruszewski_etal(u_ll, u_rr, + # normal_direction::AbstractVector, + # equations::CompressibleEulerEquationsWithGravityNoPressure2D) + noncons = (p_star - p_local) + + # For the slip wall we directly set the flux as the normal velocity is zero + return (SVector(zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_, + SVector(zero(eltype(u_inner)), + noncons, + noncons, + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_) +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquationsWithGravityNoPressure2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner))) + else # orientation == 2 + normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner))) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + +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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, rho_v1, rho_v2, rho_etot, phi = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * + (rho_etot - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + f3 = rho_v1 * v2 + f4 = (rho_etot + p) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + f4 = (rho_etot + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_etot = 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 + f3 = rho_v_normal * v2 + f4 = (rho_etot + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + +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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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) + + # Average each factor of products in flux + 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 = 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + f1 = rho_avg * v1_avg + f2 = f1 * v1_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 = 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + f1 = rho_avg * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_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))) +end + +@inline function Trixi.flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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_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 = 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 = 0.5f0 * (phi_ll + phi_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_avg * v_dot_n_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + 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)) + f1 * phi_avg + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + +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::CompressibleEulerEquationsWithGravityNoPressure2D) + # Unpack left and right state + 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) + + # Average each factor of products in flux + 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 + f1 = rho_avg * v1_avg + f2 = rho_avg * v1_avg * v1_avg + f3 = rho_avg * v1_avg * v2_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 + f4 = (rho_avg * etot_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::CompressibleEulerEquationsWithGravityNoPressure2D) + # Unpack left and right state + 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) + + # Average each factor of products in flux + 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.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 + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * etot_avg + p_avg * v_dot_n_avg + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + +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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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) + + # 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.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_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 + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + 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) + f1 * phi_avg + else + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_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) + f1 * phi_avg + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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_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.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_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) + f2 = f1 * v1_avg + f3 = f1 * v2_avg + 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) + + f1 * phi_avg) + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +function flux_nonconservative_waruszewski_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + p_jump = (p_rr - p_ll) + noncons = p_jump + ln_mean(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_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_rr, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + p_jump = (p_rr - p_ll) + noncons = p_jump + ln_mean(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 + +# For `VolumeIntegralSubcellLimiting` the nonconservative flux is created as a callable struct to +# enable dispatch on the type of the nonconservative term (local / jump). +struct FluxNonConservativeChandrashekarIsothermal <: + Trixi.FluxNonConservative{Trixi.NonConservativeJump()} +end + +Trixi.n_nonconservative_terms(::FluxNonConservativeChandrashekarIsothermal) = 2 + +const flux_nonconservative_chandrashekar_isothermal = FluxNonConservativeChandrashekarIsothermal() + +@inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, + orientation::Integer, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) + _, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + + # TODO: we assume R*T constant... Read from aux vars + RT = 1.0 # p_ll / rho_ll + + e_ll = exp(phi_ll / RT) + e_rr = exp(phi_rr / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = -rho_ll * RT * e_ll * (1 / e_rr - 1 / e_ll) + (p_rr - p_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 + +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, + orientation::Integer, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + nonconservative_type::Trixi.NonConservativeLocal, + nonconservative_term::Integer) + # TODO: we assume R*T constant... Read from aux vars + RT = 1.0 + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + if noncons == 1 + rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + e_ll = exp(phi_ll / RT) + + flux = -rho_ll * RT * e_ll + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, flux, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, flux, f0, f0) + end + else #if noncons == 2 + flux = 1 + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, flux, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, flux, f0, f0) + end + end +end + +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, u_rr, + orientation::Integer, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + nonconservative_type::Trixi.NonConservativeJump, + nonconservative_term::Integer) + 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) + + # TODO: we assume R*T constant... Read from aux vars + RT = 1.0 + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + if noncons == 1 + inv_e_ll = exp(-phi_ll / RT) + inv_e_rr = exp(-phi_rr / RT) + + flux = (inv_e_rr - inv_e_ll) + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, flux, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, flux, f0, f0) + end + else #if noncons == 2 + flux = p_rr - p_ll + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, flux, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, flux, f0, f0) + end + end +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + +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 + +# We add the "upwinding" pressure terms here, but not the average pressure terms as in CompressibleEulerEquationsWithGravity2D +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + 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 * 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 + +# We add the "upwinding" pressure terms here, but not the average pressure terms as in CompressibleEulerEquationsWithGravity2D +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # 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 * 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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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 + +@inline function Trixi.max_abs_speeds(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + 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::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, rho_v1, rho_v2, rho_etot, phi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * + (rho_etot - 0.5f0 * (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::CompressibleEulerEquationsWithGravityNoPressure2D) + 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_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.5f0 * 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::CompressibleEulerEquationsWithGravityNoPressure2D) + # See Waruszewski et al. (2022) + @unpack gamma = equations + + # convert to entropy `-rho * s` + # instead of `-rho * s / (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 + + 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_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 +@inline function Trixi.prim2cons(prim, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, v1, v2, p, phi = prim + rho_v1 = rho * v1 + rho_v2 = rho * 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 + +@inline function Trixi.density(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho = u[1] + return rho +end + +@inline function Trixi.pressure(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, rho_v1, rho_v2, rho_etot, phi = u + p = (equations.gamma - 1) * + (rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) + return p +end + +@inline function Trixi.density_pressure(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, rho_v1, rho_v2, rho_etot, phi = u + rho_times_p = (equations.gamma - 1) * + (rho * rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) + return rho_times_p +end + +# Calculate thermodynamic entropy for a conservative state `cons` +@inline function entropy_thermodynamic(cons, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # Pressure + p = (equations.gamma - 1) * + (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]) + + return s +end + +# Calculate mathematical entropy for a conservative state `cons` +@inline function entropy_math(cons, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # 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::CompressibleEulerEquationsWithGravityNoPressure2D) = entropy_math(cons, + equations) + +# Calculate total energy for a conservative state `cons` +@inline Trixi.energy_total(cons, ::CompressibleEulerEquationsWithGravityNoPressure2D) = cons[4] + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho, rho_v1, rho_v2, rho_etot, _ = 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::CompressibleEulerEquationsWithGravityNoPressure2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) - + cons[1] * cons[5] +end + +@inline function Trixi.velocity(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + 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, + orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + 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 62a1a63bfea768d250230b3d0b97d5192802a94f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Mar 2026 18:52:14 +0100 Subject: [PATCH 02/18] Include and export equation --- src/TrixiAtmo.jl | 3 ++- src/equations/equations.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index f9e4dff6..abaf3bbd 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -59,7 +59,8 @@ export CompressibleMoistEulerEquations2D, CompressibleEulerPotentialTemperatureEquationsWithGravity3D, CompressibleEulerEnergyEquationsWithGravity2D, CompressibleEulerEnergyEquationsWithGravity3D, - CompressibleEulerEquationsWithGravity2D + CompressibleEulerEquationsWithGravity2D, + CompressibleEulerEquationsWithGravityNoPressure2D export GlobalCartesianCoordinates, GlobalSphericalCoordinates diff --git a/src/equations/equations.jl b/src/equations/equations.jl index f36a3422..71423c51 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -353,4 +353,5 @@ include("compressible_euler_energy_with_gravity_3d.jl") include("shallow_water_3d.jl") include("reference_data.jl") include("compressible_euler_gravity_2d.jl") +include("compressible_euler_gravity_nopressure_2d.jl") end # @muladd From 2d5bd12e34b1dcb70df0bb4135674d5ba6526feb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Mar 2026 19:22:52 +0100 Subject: [PATCH 03/18] Fixed bugs and added equilibrium test --- ...lixir_euler_gravity_equilibrium_subcell.jl | 86 +++++++++++++++++++ src/TrixiAtmo.jl | 2 +- ...ompressible_euler_gravity_nopressure_2d.jl | 14 +-- 3 files changed, 94 insertions(+), 8 deletions(-) create mode 100644 examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl new file mode 100644 index 00000000..cd398711 --- /dev/null +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -0,0 +1,86 @@ + +using OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiAtmo + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravityNoPressure2D(1.4) + +function initial_condition_constant(x, t, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + 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_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) +surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"],) +# volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; +# volume_flux_dg = volume_flux, +# volume_flux_fv = surface_flux) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + 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, 100.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg, + save_analysis = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) #extra_node_variables = (:limiting_coefficient,) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = () # SubcellLimiterIDPCorrection(), + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., + callback = callbacks); diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index abaf3bbd..a86818c0 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -68,7 +68,7 @@ export NonlinearSolveDG export flux_chandrashekar, FluxLMARS -export flux_nonconservative_waruszewski +export flux_nonconservative_waruszewski, flux_nonconservative_chandrashekar_isothermal export flux_nonconservative_zeros, flux_nonconservative_ec, flux_nonconservative_surface_simplified, source_terms_geometric_coriolis, diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index a1cd8215..d6e939b6 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -125,7 +125,7 @@ Should be used together with [`UnstructuredMesh2D`](@ref). # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 # TODO: compute with surface_flux_function # surface_flux, surface_noncons = surface_flux_function - # flux_nonconservative_waruszewski_etal(u_ll, u_rr, + # surface_noncons(u_ll, u_rr, # normal_direction::AbstractVector, # equations::CompressibleEulerEquationsWithGravityNoPressure2D) noncons = (p_star - p_local) @@ -137,8 +137,8 @@ Should be used together with [`UnstructuredMesh2D`](@ref). zero(eltype(u_inner)), zero(eltype(u_inner))) * norm_, SVector(zero(eltype(u_inner)), - noncons, - noncons, + noncons * normal[1], + noncons * normal[2], zero(eltype(u_inner)), zero(eltype(u_inner))) * norm_) end @@ -532,7 +532,7 @@ end RT = 1.0 # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 - if noncons == 1 + if nonconservative_term == 1 rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) e_ll = exp(phi_ll / RT) @@ -544,7 +544,7 @@ end else #if orientation == 2 return SVector(f0, f0, flux, f0, f0) end - else #if noncons == 2 + else #if nonconservative_term == 2 flux = 1 f0 = zero(eltype(u_ll)) @@ -568,7 +568,7 @@ end RT = 1.0 # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 - if noncons == 1 + if nonconservative_term == 1 inv_e_ll = exp(-phi_ll / RT) inv_e_rr = exp(-phi_rr / RT) @@ -580,7 +580,7 @@ end else #if orientation == 2 return SVector(f0, f0, flux, f0, f0) end - else #if noncons == 2 + else #if nonconservative_term == 2 flux = p_rr - p_ll f0 = zero(eltype(u_ll)) From 5456cb9656e5a8e58a6de0390f9b1444a21c43f5 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 19 Mar 2026 19:38:54 +0100 Subject: [PATCH 04/18] add hydrostatic reconstruction for euler gravity (no pressure if this does not work) --- ...lixir_euler_gravity_equilibrium_subcell.jl | 13 +-- hydrostatic_reconstruction.jl | 78 +++++++++++++++++ ...ompressible_euler_gravity_nopressure_2d.jl | 83 +++++++++++++++++++ 3 files changed, 169 insertions(+), 5 deletions(-) create mode 100644 hydrostatic_reconstruction.jl diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index cd398711..d95330bf 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -22,14 +22,17 @@ initial_condition = initial_condition_constant volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) +surface_flux = (FluxHydrostaticReconstruction(flux_kennedy_gruber, hydrostatic_reconstruction), + FluxHydrostaticreconstruction(flux_nonconservative_chandrashekar_isothermal, hydrostatic_reconstruction)) + polydeg = 3 basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"],) -# volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; -# volume_flux_dg = volume_flux, -# volume_flux_fv = surface_flux) -volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +#volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) @@ -78,7 +81,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = () # SubcellLimiterIDPCorrection(), +stage_callbacks = (SubcellLimiterIDPCorrection(),) sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback diff --git a/hydrostatic_reconstruction.jl b/hydrostatic_reconstruction.jl new file mode 100644 index 00000000..b476c1a6 --- /dev/null +++ b/hydrostatic_reconstruction.jl @@ -0,0 +1,78 @@ +struct FluxHydrostaticReconstruction{NumericalFlux, HydrostaticReconstruction} + numerical_flux::NumericalFlux + hydrostatic_reconstruction::HydrostaticReconstruction +end + +@inline function (numflux::FluxHydrostaticReconstruction)(u_ll, u_rr, + orientation_or_normal_direction, + equations::Trixi.AbstractEquations) + @unpack numerical_flux, hydrostatic_reconstruction = numflux + + # Create the reconstructed left/right solution states in conservative form + u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) + + # Use the reconstructed states to compute the numerical surface flux + return numerical_flux(u_ll_star, u_rr_star, orientation_or_normal_direction, + equations) +end + +# Hydrostatic reconstruction from the paper: +# Ziming Chen, Yingjuan Zhang, Gang Li, Shouguo Qian (2022) +# "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal +# hydrostatic state under gravitational field" +# [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) +@inline function hydrostatic_reconstruction(u_ll, u_rr, equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # Unpack left and right states + rho_ll, rho_v_ll, rho_e_ll, phi_ll = u_ll + rho_rr, rho_v_rr, rho_e_rr, phi_rr = u_rr + + # Compute equilibrium potential (# TODO: For general case we need to add phi_num - phi_exact)) + psi_eq_ll = psi_ll # phi_ll - phi_exact_ll + psi_eq_rr = psi_rr # phi_rr - phi_exact_rr + + # Set RT0 to 1 for now. + RT0 = 1 + + # Compute equlibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0*rho_eq_ll + p_eq_rr = RT0*rho_eq_rr + rho_v_eq_ll = 0.0 + rho_v_eq_rr = 0.0 + rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + + u_eq_ll = (rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) + u_eq_rr = (rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + + # Compute residual contribution + u_res_ll = u_ll - u_eq_ll + u_res_rr = u_rr - u_eq_rr + + # Reconstruct phi + phi_star = max(phi_ll, phi_rr) + + # Compute reconstructed equilibrium potential + psi_eq_ll = psi_ll + phi_ll - phi_star + psi_eq_rr = psi_rr + phi_rr - phi_star + + # Compute reconstructed equilibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0*rho_eq_ll + p_eq_rr = RT0*rho_eq_rr + rho_v_eq_ll = 0.0 + rho_v_eq_rr = 0.0 + rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + + u_eq_ll = SVector(rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) + u_eq_rr = SVector(rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + + # Compute reconstructed state + u_star_ll = u_eq_ll + u_res_ll + u_star_rr = u_eq_rr + u_res_rr + + return u_star_ll, u_star_rr +end \ No newline at end of file diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index d6e939b6..5a466686 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -964,4 +964,87 @@ end diss = -0.5f0 * λ * (u_rr - u_ll) return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) end + +####################################################################### +# Hydrostatic reconstruction + +struct FluxHydrostaticReconstruction{NumericalFlux, HydrostaticReconstruction} + numerical_flux::NumericalFlux + hydrostatic_reconstruction::HydrostaticReconstruction +end + +@inline function (numflux::FluxHydrostaticReconstruction)(u_ll, u_rr, + orientation_or_normal_direction, + equations::Trixi.AbstractEquations) + @unpack numerical_flux, hydrostatic_reconstruction = numflux + + # Create the reconstructed left/right solution states in conservative form + u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) + + # Use the reconstructed states to compute the numerical surface flux + return numerical_flux(u_ll_star, u_rr_star, orientation_or_normal_direction, + equations) +end + +# Hydrostatic reconstruction from the paper: +# Ziming Chen, Yingjuan Zhang, Gang Li, Shouguo Qian (2022) +# "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal +# hydrostatic state under gravitational field" +# [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) +@inline function hydrostatic_reconstruction(u_ll, u_rr, equations::CompressibleEulerEquationsWithGravityNoPressure2D) + # Unpack left and right states + rho_ll, rho_v_ll, rho_e_ll, phi_ll = u_ll + rho_rr, rho_v_rr, rho_e_rr, phi_rr = u_rr + + # Compute equilibrium potential (# TODO: For general case we need to add phi_num - phi_exact)) + psi_eq_ll = psi_ll # phi_ll - phi_exact_ll + psi_eq_rr = psi_rr # phi_rr - phi_exact_rr + + # Set RT0 to 1 for now. + RT0 = 1 + + # Compute equlibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0*rho_eq_ll + p_eq_rr = RT0*rho_eq_rr + rho_v_eq_ll = 0.0 + rho_v_eq_rr = 0.0 + rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + + u_eq_ll = (rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) + u_eq_rr = (rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + + # Compute residual contribution + u_res_ll = u_ll - u_eq_ll + u_res_rr = u_rr - u_eq_rr + + # Reconstruct phi + phi_star = max(phi_ll, phi_rr) + + # Compute reconstructed equilibrium potential + psi_eq_ll = psi_ll + phi_ll - phi_star + psi_eq_rr = psi_rr + phi_rr - phi_star + + # Compute reconstructed equilibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0*rho_eq_ll + p_eq_rr = RT0*rho_eq_rr + rho_v_eq_ll = 0.0 + rho_v_eq_rr = 0.0 + rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + + u_eq_ll = SVector(rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) + u_eq_rr = SVector(rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + + # Compute reconstructed state + u_star_ll = u_eq_ll + u_res_ll + u_star_rr = u_eq_rr + u_res_rr + + return u_star_ll, u_star_rr +end + end # @muladd From afd6b09def6797790534bcd5b8820d648a8826a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 19 Mar 2026 19:49:19 +0100 Subject: [PATCH 05/18] Fixed some bugs --- ...lixir_euler_gravity_equilibrium_subcell.jl | 6 ++- src/TrixiAtmo.jl | 3 +- ...ompressible_euler_gravity_nopressure_2d.jl | 42 ++++++++++--------- 3 files changed, 28 insertions(+), 23 deletions(-) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index d95330bf..695b5b01 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -22,8 +22,10 @@ initial_condition = initial_condition_constant volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) -surface_flux = (FluxHydrostaticReconstruction(flux_kennedy_gruber, hydrostatic_reconstruction), - FluxHydrostaticreconstruction(flux_nonconservative_chandrashekar_isothermal, hydrostatic_reconstruction)) +surface_flux = (FluxHydrostaticReconstruction(flux_kennedy_gruber, + hydrostatic_reconstruction), + FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, + hydrostatic_reconstruction)) polydeg = 3 basis = LobattoLegendreBasis(polydeg) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index a86818c0..f4c831b3 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -66,7 +66,8 @@ export GlobalCartesianCoordinates, GlobalSphericalCoordinates export NonlinearSolveDG -export flux_chandrashekar, FluxLMARS +export flux_chandrashekar, FluxLMARS, FluxHydrostaticReconstruction, + hydrostatic_reconstruction export flux_nonconservative_waruszewski, flux_nonconservative_chandrashekar_isothermal diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index 5a466686..60dc4ce7 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -991,10 +991,17 @@ end # "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal # hydrostatic state under gravitational field" # [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) -@inline function hydrostatic_reconstruction(u_ll, u_rr, equations::CompressibleEulerEquationsWithGravityNoPressure2D) +@inline function hydrostatic_reconstruction(u_ll, u_rr, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) # Unpack left and right states - rho_ll, rho_v_ll, rho_e_ll, phi_ll = u_ll - rho_rr, rho_v_rr, rho_e_rr, phi_rr = u_rr + 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_e_ll = u_ll[4] + rho_e_rr = u_rr[4] + + psi_ll = p_ll / rho_ll * log(rho_ll) + psi_rr = p_rr / rho_rr * log(rho_rr) # Compute equilibrium potential (# TODO: For general case we need to add phi_num - phi_exact)) psi_eq_ll = psi_ll # phi_ll - phi_exact_ll @@ -1006,15 +1013,13 @@ end # Compute equlibrium state rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available rho_eq_rr = exp(psi_eq_rr / (RT0)) - p_eq_ll = RT0*rho_eq_ll - p_eq_rr = RT0*rho_eq_rr - rho_v_eq_ll = 0.0 - rho_v_eq_rr = 0.0 - rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll - rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + p_eq_ll = RT0 * rho_eq_ll + p_eq_rr = RT0 * rho_eq_rr + rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr - u_eq_ll = (rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) - u_eq_rr = (rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr) # Compute residual contribution u_res_ll = u_ll - u_eq_ll @@ -1030,15 +1035,13 @@ end # Compute reconstructed equilibrium state rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available rho_eq_rr = exp(psi_eq_rr / (RT0)) - p_eq_ll = RT0*rho_eq_ll - p_eq_rr = RT0*rho_eq_rr - rho_v_eq_ll = 0.0 - rho_v_eq_rr = 0.0 - rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll - rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr + p_eq_ll = RT0 * rho_eq_ll + p_eq_rr = RT0 * rho_eq_rr + rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr - u_eq_ll = SVector(rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) - u_eq_rr = SVector(rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr) # Compute reconstructed state u_star_ll = u_eq_ll + u_res_ll @@ -1046,5 +1049,4 @@ end return u_star_ll, u_star_rr end - end # @muladd From a01d79ffb208655682d10676b586edf50addbd9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 20 Mar 2026 06:30:04 +0100 Subject: [PATCH 06/18] Added isvalid function --- src/equations/compressible_euler_gravity_nopressure_2d.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index 60dc4ce7..2142227a 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -965,6 +965,14 @@ end return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) end +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) + if u[1] <= 0 || pressure(u, equations) <= 0 + return false + end + return true +end + ####################################################################### # Hydrostatic reconstruction From b46918e4afc4855787bb81db155a55befca57ae1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 20 Mar 2026 06:30:55 +0100 Subject: [PATCH 07/18] Modified test to use pure FV (not well balanced) --- .../dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index 695b5b01..99b508ea 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -35,6 +35,7 @@ volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) #volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +volume_integral = VolumeIntegralPureLGLFiniteVolume(volume_flux_fv = surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) @@ -83,7 +84,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = (SubcellLimiterIDPCorrection(),) +stage_callbacks = () #SubcellLimiterIDPCorrection(), sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback From bcf8a193399efb9dce57ab56ab2b4807e3e22329 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 20 Mar 2026 07:40:25 +0100 Subject: [PATCH 08/18] Implemented get_boundary_outer_state --- ...ompressible_euler_gravity_nopressure_2d.jl | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index 2142227a..25d7e4f9 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -966,7 +966,8 @@ end end # State validation for Newton-bisection method of subcell IDP limiting -@inline function Base.isvalid(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) +@inline function Base.isvalid(u, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) if u[1] <= 0 || pressure(u, equations) <= 0 return false end @@ -1057,4 +1058,27 @@ end return u_star_ll, u_star_rr end + +# Get outer state for the min/max limiters +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_slip_wall), + orientation, boundary_index, + mesh::TreeMesh{2}, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + dg, cache, indices...) + if orientation == 1 + #if boundary_index == 1 # Element is on the right, boundary on the left + return SVector(u_inner[1], + u_inner[2] - 2 * u_inner[2], + u_inner[3], + u_inner[4], + u_inner[5]) + else + return SVector(u_inner[1], + u_inner[2], + u_inner[3] - 2 * u_inner[3], + u_inner[4], + u_inner[5]) + end +end end # @muladd From 45f3e84535b032c2896a1b4bf79a0d86e24463a1 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 24 Mar 2026 13:11:54 +0100 Subject: [PATCH 09/18] fix hydrostatic reconstruction --- ...lixir_euler_gravity_equilibrium_subcell.jl | 19 ++++++++++--------- ...ompressible_euler_gravity_nopressure_2d.jl | 13 +++++++------ 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index 99b508ea..ba425282 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -22,20 +22,21 @@ initial_condition = initial_condition_constant volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) -surface_flux = (FluxHydrostaticReconstruction(flux_kennedy_gruber, +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, DissipationLocalLaxFriedrichs()), hydrostatic_reconstruction), FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, hydrostatic_reconstruction)) -polydeg = 3 +polydeg = 0 basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"],) -volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -#volume_integral = VolumeIntegralFluxDifferencing(volume_flux) -volume_integral = VolumeIntegralPureLGLFiniteVolume(volume_flux_fv = surface_flux) +# volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; +# volume_flux_dg = volume_flux, +# volume_flux_fv = surface_flux) +# +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +#volume_integral = VolumeIntegralPureLGLFiniteVolume(volume_flux_fv = surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) @@ -62,14 +63,14 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 10000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, analysis_polydeg = polydeg, save_analysis = true) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 10, +save_solution = SaveSolutionCallback(dt = 10, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) #extra_node_variables = (:limiting_coefficient,) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index 2142227a..f3581393 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -1000,7 +1000,8 @@ end # hydrostatic state under gravitational field" # [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) @inline function hydrostatic_reconstruction(u_ll, u_rr, - equations::CompressibleEulerEquationsWithGravityNoPressure2D) + equations::Union{CompressibleEulerEquationsWithGravityNoPressure2D, + CompressibleEulerEquationsWithGravity2D}) # Unpack left and right states 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) @@ -1041,15 +1042,15 @@ end psi_eq_rr = psi_rr + phi_rr - phi_star # Compute reconstructed equilibrium state - rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_ll = exp(psi_eq_ll / (RT0)) # Replace RT0 with R*T0 once available rho_eq_rr = exp(psi_eq_rr / (RT0)) p_eq_ll = RT0 * rho_eq_ll p_eq_rr = RT0 * rho_eq_rr - rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll - rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr + rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_star + rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_star - u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll) - u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_star) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_star) # Compute reconstructed state u_star_ll = u_eq_ll + u_res_ll From b8aa60bb2dd14dbc36e618c575d15bb64cf8e81b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 24 Mar 2026 20:03:37 +0100 Subject: [PATCH 10/18] Added blast in atmosphere example with cheap subsonic outflow --- .../elixir_euler_gravity_blast_subcell.jl | 188 ++++++++++++++++++ ...lixir_euler_gravity_equilibrium_subcell.jl | 4 +- ...ompressible_euler_gravity_nopressure_2d.jl | 20 ++ 3 files changed, 210 insertions(+), 2 deletions(-) create mode 100644 examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl new file mode 100644 index 00000000..a05434a2 --- /dev/null +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl @@ -0,0 +1,188 @@ + +using OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiAtmo + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravityNoPressure2D(1.4) + +function initial_condition_constant(x, t, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho = exp(-x[2]) + v1 = 0.0 + v2 = 0.0 + p = exp(-x[2]) + inicenter = SVector(0, 0.1) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + # Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + # r0 = 0.5 # = more reasonable setup + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + + p = r > r0 ? p : p + p0_inner + + prim = SVector(rho, v1, v2, p, x[2]) + return prim2cons(prim, equations) +end + +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_slip_wall), + orientation, boundary_index, + mesh::TreeMesh{2}, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + dg, cache, indices...) + if orientation == 1 + #if boundary_index == 1 # Element is on the right, boundary on the left + return SVector(u_inner[1], + u_inner[2] - 2 * u_inner[2], + u_inner[3], + u_inner[4], + u_inner[5]) + else + return SVector(u_inner[1], + u_inner[2], + u_inner[3] - 2 * u_inner[3], + u_inner[4], + u_inner[5]) + end +end + +# See Section 2.3 of the reference below for a discussion of robust +# subsonic inflow/outflow boundary conditions. +# +# - Jan-Reneé Carlson (2011) +# Inflow/Outflow Boundary Conditions with Application to FUN3D. +# [NASA TM 20110022658](https://ntrs.nasa.gov/citations/20110022658) +@inline function boundary_condition_subsonic(u_inner, orientation::Integer, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_loc, v1_loc, v2_loc, p_loc, phi_loc = cons2prim(u_inner, equations) + + # For subsonic boundary: Take pressure from initial condition + p_out = pressure(initial_condition_constant(x, t, equations), equations) + + prim = SVector(rho_loc, v1_loc, v2_loc, p_out, phi_loc) + u_surface = prim2cons(prim, equations) + + flux_conservative, flux_nonconservative = surface_flux_function + + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner))) + else # orientation == 2 + normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner))) + end + + if isodd(direction) + return -flux_conservative.numerical_flux(u_inner, u_surface, -normal_direction, + equations), + -flux_nonconservative.numerical_flux(u_inner, u_surface, -normal_direction, + equations) + else + return flux_conservative.numerical_flux(u_inner, u_surface, normal_direction, + equations), + flux_nonconservative.numerical_flux(u_inner, u_surface, normal_direction, + equations) + end +end + +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_subsonic), + orientation, boundary_index, + mesh::TreeMesh{2}, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + dg, cache, indices...) + rho_loc, v1_loc, v2_loc, p_loc, phi_loc = cons2prim(u_inner, equations) + x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) + + # For subsonic boundary: Take pressure from initial condition + p_out = pressure(initial_condition_constant(x, t, equations), equations) + + prim = SVector(rho_loc, v1_loc, v2_loc, p_out, phi_loc) + u_surface = prim2cons(prim, equations) + return u_surface +end + +initial_condition = initial_condition_constant + +volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) +surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) + +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction), + FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, + hydrostatic_reconstruction)) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + local_twosided_variables_cons = ["rho"]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +#volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_max = (4.0, 7.7) +coordinates_min = (-4.0, -0.3) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 100_000, + periodicity = (false, false)) + +boundary_conditions = (; + x_neg = boundary_condition_subsonic, + x_pos = boundary_condition_subsonic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_subsonic) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg, + save_analysis = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim, + extra_node_variables = (:limiting_coefficient,)) # + +stepsize_callback = StepsizeCallback(cfl = 0.25) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., + callback = callbacks); diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index ba425282..8eed67a7 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -20,9 +20,9 @@ end initial_condition = initial_condition_constant volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) -surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isothermal) -surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, DissipationLocalLaxFriedrichs()), +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, + DissipationLocalLaxFriedrichs()), hydrostatic_reconstruction), FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, hydrostatic_reconstruction)) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index e6ab8144..a8eae572 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -500,6 +500,26 @@ Trixi.n_nonconservative_terms(::FluxNonConservativeChandrashekarIsothermal) = 2 const flux_nonconservative_chandrashekar_isothermal = FluxNonConservativeChandrashekarIsothermal() +@inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) + _, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + + # TODO: we assume R*T constant... Read from aux vars + RT = 1.0 # p_ll / rho_ll + + e_ll = exp(phi_ll / RT) + e_rr = exp(phi_rr / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = -rho_ll * RT * e_ll * (1 / e_rr - 1 / e_ll) + (p_rr - p_ll) + + f0 = zero(eltype(u_ll)) + return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], f0, + f0) +end + @inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) From a81ff2e7d5aa92f386dc21dfb5c37877be5718c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 28 Apr 2026 13:45:21 +0200 Subject: [PATCH 11/18] Added field for auxiliary variables in equation --- .../elixir_euler_gravity_blast_subcell.jl | 120 ++++++++------ ...lixir_euler_gravity_equilibrium_subcell.jl | 2 +- ...ompressible_euler_gravity_nopressure_2d.jl | 150 ++++++++++-------- 3 files changed, 157 insertions(+), 115 deletions(-) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl index a05434a2..af06a397 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl @@ -26,31 +26,31 @@ function initial_condition_constant(x, t, p = r > r0 ? p : p + p0_inner - prim = SVector(rho, v1, v2, p, x[2]) + prim = SVector(rho, v1, v2, p, x[2], 1.0) return prim2cons(prim, equations) end -@inline function Trixi.get_boundary_outer_state(u_inner, t, - boundary_condition::typeof(boundary_condition_slip_wall), - orientation, boundary_index, - mesh::TreeMesh{2}, - equations::CompressibleEulerEquationsWithGravityNoPressure2D, - dg, cache, indices...) - if orientation == 1 - #if boundary_index == 1 # Element is on the right, boundary on the left - return SVector(u_inner[1], - u_inner[2] - 2 * u_inner[2], - u_inner[3], - u_inner[4], - u_inner[5]) - else - return SVector(u_inner[1], - u_inner[2], - u_inner[3] - 2 * u_inner[3], - u_inner[4], - u_inner[5]) - end -end +# @inline function Trixi.get_boundary_outer_state(u_inner, t, +# boundary_condition::typeof(boundary_condition_slip_wall), +# orientation, boundary_index, +# mesh::TreeMesh{2}, +# equations::CompressibleEulerEquationsWithGravityNoPressure2D, +# dg, cache, indices...) +# if orientation == 1 +# #if boundary_index == 1 # Element is on the right, boundary on the left +# return SVector(u_inner[1], +# u_inner[2] - 2 * u_inner[2], +# u_inner[3], +# u_inner[4], +# u_inner[5]) +# else +# return SVector(u_inner[1], +# u_inner[2], +# u_inner[3] - 2 * u_inner[3], +# u_inner[4], +# u_inner[5]) +# end +# end # See Section 2.3 of the reference below for a discussion of robust # subsonic inflow/outflow boundary conditions. @@ -62,27 +62,14 @@ end direction, x, t, surface_flux_function, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_loc, v1_loc, v2_loc, p_loc, phi_loc = cons2prim(u_inner, equations) - - # For subsonic boundary: Take pressure from initial condition - p_out = pressure(initial_condition_constant(x, t, equations), equations) - - prim = SVector(rho_loc, v1_loc, v2_loc, p_out, phi_loc) - u_surface = prim2cons(prim, equations) + u_surface, normal_direction = get_outer_state_and_normal(u_inner, orientation, + direction, x, t, equations) flux_conservative, flux_nonconservative = surface_flux_function - - # get the appropriate normal vector from the orientation - if orientation == 1 - normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner))) - else # orientation == 2 - normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner))) - end - if isodd(direction) - return -flux_conservative.numerical_flux(u_inner, u_surface, -normal_direction, + return -flux_conservative.numerical_flux(u_inner, u_surface, normal_direction, equations), - -flux_nonconservative.numerical_flux(u_inner, u_surface, -normal_direction, + -flux_nonconservative.numerical_flux(u_inner, u_surface, normal_direction, equations) else return flux_conservative.numerical_flux(u_inner, u_surface, normal_direction, @@ -94,19 +81,62 @@ end @inline function Trixi.get_boundary_outer_state(u_inner, t, boundary_condition::typeof(boundary_condition_subsonic), - orientation, boundary_index, + orientation, direction, mesh::TreeMesh{2}, equations::CompressibleEulerEquationsWithGravityNoPressure2D, dg, cache, indices...) - rho_loc, v1_loc, v2_loc, p_loc, phi_loc = cons2prim(u_inner, equations) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) - # For subsonic boundary: Take pressure from initial condition - p_out = pressure(initial_condition_constant(x, t, equations), equations) + u_surface, normal_direction = get_outer_state_and_normal(u_inner, orientation, + direction, x, t, equations) + return u_surface +end - prim = SVector(rho_loc, v1_loc, v2_loc, p_out, phi_loc) +@inline function get_outer_state_and_normal(u_inner, orientation::Integer, + direction, x, t, + equations::CompressibleEulerEquationsWithGravityNoPressure2D) + rho_loc, v1_loc, v2_loc, p_loc, phi_loc, aux_loc = cons2prim(u_inner, equations) + + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner))) + else # orientation == 2 + normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner))) + end + + # Flip the normal vector if necessary + if isodd(direction) + normal_direction = -normal_direction + end + + # Outer state + rho_out, v1_out, v2_out, p_out, phi_out, aux_out = cons2prim(initial_condition_constant(x, + t, + equations), + equations) + + # Impose outflow... No inflow allowed + v_normal = (v1_loc * normal_direction[1] + v2_loc * normal_direction[2]) / + sqrt(equations.gamma * p_loc / rho_loc) + if v_normal <= -1 + # Supersonic inflow (state from outside).. do nothing + elseif v_normal < 0 + # Subsonic inflow (state from outside, except for pressure) + p_out = p_loc + elseif v_normal < 1 + # Subsonic outflow (state from inside, except for pressure) + rho_out = rho_loc + v1_out = v1_loc + v2_out = v2_loc + else # v_normal >= 1 + # supersonic outflow (state from inside) + return u_inner, normal_direction + end + + # Get outer state + prim = SVector(rho_out, v1_out, v2_out, p_out, phi_loc, aux_loc) u_surface = prim2cons(prim, equations) - return u_surface + return u_surface, normal_direction end initial_condition = initial_condition_constant diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index 8eed67a7..61515ab7 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -13,7 +13,7 @@ function initial_condition_constant(x, t, v1 = 0.0 v2 = 0.0 p = exp(-x[2]) - prim = SVector(rho, v1, v2, p, x[2]) + prim = SVector(rho, v1, v2, p, x[2], 1.0) return prim2cons(prim, equations) end diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index a8eae572..c4e4c54b 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -43,7 +43,7 @@ References: - 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 CompressibleEulerEquationsWithGravityNoPressure2D{RealT <: Real} <: - Trixi.AbstractCompressibleEulerEquations{2, 5} + Trixi.AbstractCompressibleEulerEquations{2, 6} gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications @@ -54,16 +54,20 @@ struct CompressibleEulerEquationsWithGravityNoPressure2D{RealT <: Real} <: end Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravityNoPressure2D) = True() + +# The auxiliary variable is used to store the mean temperature of the element for isothermal-equilibrium-preserving discretizations Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravityNoPressure2D) = ("rho", "rho_v1", "rho_v2", "rho_etot", - "phi") + "phi", + "aux") Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravityNoPressure2D) = ("rho", "v1", "v2", "p", - "phi") + "phi", + "aux") """ boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, @@ -135,11 +139,13 @@ Should be used together with [`UnstructuredMesh2D`](@ref). zero(eltype(u_inner)), zero(eltype(u_inner)), zero(eltype(u_inner)), + zero(eltype(u_inner)), zero(eltype(u_inner))) * norm_, SVector(zero(eltype(u_inner)), noncons * normal[1], noncons * normal[2], zero(eltype(u_inner)), + zero(eltype(u_inner)), zero(eltype(u_inner))) * norm_) end @@ -194,7 +200,7 @@ end # Calculate 2D flux for a single point @inline function Trixi.flux(u, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, rho_v1, rho_v2, rho_etot, phi = u + rho, rho_v1, rho_v2, rho_etot, phi, _ = u v1 = rho_v1 / rho v2 = rho_v2 / rho p = (equations.gamma - 1) * @@ -210,7 +216,7 @@ end f3 = rho_v2 * v2 f4 = (rho_etot + p) * v2 end - return SVector(f1, f2, f3, f4, zero(eltype(u))) + return SVector(f1, f2, f3, f4, zero(eltype(u)), zero(eltype(u))) end # Calculate 2D flux for a single point in the normal direction @@ -226,7 +232,7 @@ end f2 = rho_v_normal * v1 f3 = rho_v_normal * v2 f4 = (rho_etot + p) * v_normal - return SVector(f1, f2, f3, f4, zero(eltype(u))) + return SVector(f1, f2, f3, f4, zero(eltype(u)), zero(eltype(u))) end """ @@ -248,8 +254,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::CompressibleEulerEquationsWithGravityNoPressure2D) # 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, 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 = 0.5f0 * (rho_ll + rho_rr) @@ -276,14 +282,14 @@ The modification is in the energy flux to guarantee pressure equilibrium and was f1 * phi_avg end - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end @inline function Trixi.flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravityNoPressure2D) # 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, 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] @@ -304,7 +310,7 @@ end 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)) + f1 * phi_avg - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end """ @@ -345,7 +351,7 @@ Kinetic energy preserving two-point flux by f4 = (rho_avg * etot_avg + p_avg) * v2_avg end - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end @inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, @@ -370,7 +376,7 @@ end f3 = f1 * v2_avg f4 = f1 * etot_avg + p_avg * v_dot_n_avg - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end """ @@ -391,8 +397,8 @@ See also @inline function Trixi.flux_ranocha(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) # 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, 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) @@ -423,14 +429,14 @@ See also 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + f1 * phi_avg end - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end @inline function Trixi.flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravityNoPressure2D) # 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, 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] @@ -455,14 +461,14 @@ end 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))) + return SVector(f1, f2, f3, f4, zero(eltype(u_ll)), zero(eltype(u_ll))) end function flux_nonconservative_waruszewski_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) - rho_rr, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + rho_rr, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 p_jump = (p_rr - p_ll) @@ -470,13 +476,13 @@ function flux_nonconservative_waruszewski_etal(u_ll, u_rr, f0 = zero(eltype(u_ll)) return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], - f0, f0) + f0, f0, f0) end function flux_nonconservative_waruszewski_etal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) - rho_rr, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + rho_rr, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 p_jump = (p_rr - p_ll) @@ -484,9 +490,9 @@ function flux_nonconservative_waruszewski_etal(u_ll, u_rr, orientation::Integer, f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, noncons, f0, f0, f0) + return SVector(f0, noncons, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, noncons, f0, f0) + return SVector(f0, f0, noncons, f0, f0, f0) end end @@ -503,8 +509,8 @@ const flux_nonconservative_chandrashekar_isothermal = FluxNonConservativeChandra @inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) - _, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) # TODO: we assume R*T constant... Read from aux vars RT = 1.0 # p_ll / rho_ll @@ -517,14 +523,14 @@ const flux_nonconservative_chandrashekar_isothermal = FluxNonConservativeChandra f0 = zero(eltype(u_ll)) return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], f0, - f0) + f0, f0) end @inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll = cons2prim(u_ll, equations) - _, _, _, p_rr, phi_rr = cons2prim(u_rr, equations) + rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) # TODO: we assume R*T constant... Read from aux vars RT = 1.0 # p_ll / rho_ll @@ -537,9 +543,9 @@ end f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, noncons, f0, f0, f0) + return SVector(f0, noncons, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, noncons, f0, f0) + return SVector(f0, f0, noncons, f0, f0, f0) end end @@ -553,25 +559,25 @@ end # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 if nonconservative_term == 1 - rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll, _ = cons2prim(u_ll, equations) e_ll = exp(phi_ll / RT) flux = -rho_ll * RT * e_ll f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, flux, f0, f0, f0) + return SVector(f0, flux, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, flux, f0, f0) + return SVector(f0, f0, flux, f0, f0, f0) end else #if nonconservative_term == 2 flux = 1 f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, flux, f0, f0, f0) + return SVector(f0, flux, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, flux, f0, f0) + return SVector(f0, f0, flux, f0, f0, f0) end end end @@ -581,8 +587,8 @@ end equations::CompressibleEulerEquationsWithGravityNoPressure2D, nonconservative_type::Trixi.NonConservativeJump, nonconservative_term::Integer) - 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, phi_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr, _ = cons2prim(u_rr, equations) # TODO: we assume R*T constant... Read from aux vars RT = 1.0 @@ -596,18 +602,18 @@ end f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, flux, f0, f0, f0) + return SVector(f0, flux, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, flux, f0, f0) + return SVector(f0, f0, flux, f0, f0, f0) end else #if nonconservative_term == 2 flux = p_rr - p_ll f0 = zero(eltype(u_ll)) if orientation == 1 - return SVector(f0, flux, f0, f0, f0) + return SVector(f0, flux, f0, f0, f0, f0) else #if orientation == 2 - return SVector(f0, f0, flux, f0, f0) + return SVector(f0, f0, flux, f0, f0, f0) end end end @@ -663,8 +669,8 @@ References: else # orientation == 2 f3 = f3 + p end - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) + f0 = zero(eltype(u_ll)) + return SVector(f1, f2, f3, f4, f0, f0) end # We add the "upwinding" pressure terms here, but not the average pressure terms as in CompressibleEulerEquationsWithGravity2D @@ -696,11 +702,11 @@ end f1, f2, f3, f4, _ = u_rr * v f4 = f4 + p_rr * v end - + f0 = zero(eltype(u_ll)) return SVector(f1, f2 + p * normal_direction[1], f3 + p * normal_direction[2], - f4, zero(eltype(u_ll))) + f4, f0, f0) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the @@ -798,7 +804,8 @@ end c * u[2] + s * u[3], -s * u[2] + c * u[3], u[4], - u[5]) + u[5], + u[6]) end # Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction @@ -821,7 +828,8 @@ end c * u[2] - s * u[3], s * u[2] + c * u[3], u[4], - u[5]) + u[5], + u[6]) end @inline function Trixi.max_abs_speeds(u, @@ -835,20 +843,20 @@ end # Convert conservative variables to primitive @inline function Trixi.cons2prim(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, rho_v1, rho_v2, rho_etot, phi = u + rho, rho_v1, rho_v2, rho_etot, phi, aux = u v1 = rho_v1 / rho v2 = rho_v2 / rho p = (equations.gamma - 1) * (rho_etot - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) - rho * phi) - return SVector(rho, v1, v2, p, phi) + return SVector(rho, v1, v2, p, phi, aux) end # Convert conservative variables to entropy (see, e.g., Waruszewski et al. (2022)) @inline function Trixi.cons2entropy(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, rho_v1, rho_v2, rho_etot, phi = u + rho, rho_v1, rho_v2, rho_etot, phi, aux = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -863,7 +871,7 @@ end w3 = rho_p * v2 w4 = -rho_p - return SVector(w1, w2, w3, w4, phi) + return SVector(w1, w2, w3, w4, phi, aux) end @inline function Trixi.entropy2cons(w, @@ -875,6 +883,7 @@ end # instead of `-rho * s / (gamma - 1)` V1, V2, V3, V5, _ = w .* (gamma - 1) phi = w[5] + aux = w[6] # s = specific entropy s = gamma - V1 + (V2^2 + V3^2) / (2 * V5) - V5 * phi @@ -886,18 +895,18 @@ end rho_v1 = rho_iota * V2 rho_v2 = rho_iota * V3 rho_etot = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5)) + rho * phi - return SVector(rho, rho_v1, rho_v2, rho_etot, phi) + return SVector(rho, rho_v1, rho_v2, rho_etot, phi, aux) end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, v1, v2, p, phi = prim + rho, v1, v2, p, phi, aux = prim rho_v1 = rho * v1 rho_v2 = rho * 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) + return SVector(rho, rho_v1, rho_v2, rho_etot, phi, aux) end @inline function Trixi.density(u, @@ -908,7 +917,7 @@ end @inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, rho_v1, rho_v2, rho_etot, phi = u + rho, rho_v1, rho_v2, rho_etot, phi, _ = u p = (equations.gamma - 1) * (rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho - rho * phi) return p @@ -916,7 +925,7 @@ end @inline function Trixi.density_pressure(u, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho, rho_v1, rho_v2, rho_etot, phi = u + rho, rho_v1, rho_v2, rho_etot, phi, _ = u rho_times_p = (equations.gamma - 1) * (rho * rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) return rho_times_p @@ -982,7 +991,8 @@ end λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) diss = -0.5f0 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) + f0 = zero(eltype(u_ll)) + return SVector(diss[1], diss[2], diss[3], diss[4], f0, f0) end # State validation for Newton-bisection method of subcell IDP limiting @@ -1024,8 +1034,8 @@ end equations::Union{CompressibleEulerEquationsWithGravityNoPressure2D, CompressibleEulerEquationsWithGravity2D}) # Unpack left and right states - 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, phi_ll, aux_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr, aux_rr = cons2prim(u_rr, equations) rho_e_ll = u_ll[4] rho_e_rr = u_rr[4] @@ -1048,8 +1058,8 @@ end rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr - u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll) - u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll, aux_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr, aux_rr) # Compute residual contribution u_res_ll = u_ll - u_eq_ll @@ -1070,8 +1080,8 @@ end rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_star rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_star - u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_star) - u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_star) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_star, aux_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_star, aux_rr) # Compute reconstructed state u_star_ll = u_eq_ll + u_res_ll @@ -1093,13 +1103,15 @@ end u_inner[2] - 2 * u_inner[2], u_inner[3], u_inner[4], - u_inner[5]) + u_inner[5], + u_inner[6]) else return SVector(u_inner[1], u_inner[2], u_inner[3] - 2 * u_inner[3], u_inner[4], - u_inner[5]) + u_inner[5], + u_inner[6]) end end end # @muladd From eda2e800578731b4e0fc2eedf973bf01136fd06b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 28 Apr 2026 14:59:41 +0200 Subject: [PATCH 12/18] Use auxiliary variable to save the equilibrium R*T (using left value for nonconservative terms) --- ...ompressible_euler_gravity_nopressure_2d.jl | 58 +++++++++++-------- 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index c4e4c54b..cb55d7e5 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -509,11 +509,14 @@ const flux_nonconservative_chandrashekar_isothermal = FluxNonConservativeChandra @inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + rho_ll, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) - # TODO: we assume R*T constant... Read from aux vars - RT = 1.0 # p_ll / rho_ll + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll e_ll = exp(phi_ll / RT) e_rr = exp(phi_rr / RT) @@ -529,11 +532,14 @@ end @inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - rho_ll, _, _, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + rho_ll, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) - # TODO: we assume R*T constant... Read from aux vars - RT = 1.0 # p_ll / rho_ll + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll e_ll = exp(phi_ll / RT) e_rr = exp(phi_rr / RT) @@ -554,15 +560,13 @@ end equations::CompressibleEulerEquationsWithGravityNoPressure2D, nonconservative_type::Trixi.NonConservativeLocal, nonconservative_term::Integer) - # TODO: we assume R*T constant... Read from aux vars - RT = 1.0 # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 if nonconservative_term == 1 - rho_ll, v1_ll, v2_ll, p_ll, phi_ll, _ = cons2prim(u_ll, equations) - e_ll = exp(phi_ll / RT) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + e_ll = exp(phi_ll / RT_ll) - flux = -rho_ll * RT * e_ll + flux = -rho_ll * RT_ll * e_ll f0 = zero(eltype(u_ll)) if orientation == 1 @@ -587,11 +591,14 @@ end equations::CompressibleEulerEquationsWithGravityNoPressure2D, nonconservative_type::Trixi.NonConservativeJump, nonconservative_term::Integer) - rho_ll, v1_ll, v2_ll, p_ll, phi_ll, _ = cons2prim(u_ll, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr, phi_rr, _ = cons2prim(u_rr, equations) - # TODO: we assume R*T constant... Read from aux vars - RT = 1.0 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 if nonconservative_term == 1 @@ -1034,8 +1041,8 @@ end equations::Union{CompressibleEulerEquationsWithGravityNoPressure2D, CompressibleEulerEquationsWithGravity2D}) # Unpack left and right states - rho_ll, v1_ll, v2_ll, p_ll, phi_ll, aux_ll = cons2prim(u_ll, equations) - rho_rr, v1_rr, v2_rr, p_rr, phi_rr, aux_rr = cons2prim(u_rr, equations) + rho_ll, v1_ll, v2_ll, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, phi_rr, RT_rr = cons2prim(u_rr, equations) rho_e_ll = u_ll[4] rho_e_rr = u_rr[4] @@ -1047,19 +1054,22 @@ end psi_eq_ll = psi_ll # phi_ll - phi_exact_ll psi_eq_rr = psi_rr # phi_rr - phi_exact_rr - # Set RT0 to 1 for now. - RT0 = 1 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT0 = RT_ll # Compute equlibrium state - rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available + rho_eq_ll = exp(psi_eq_ll / (RT0)) rho_eq_rr = exp(psi_eq_rr / (RT0)) p_eq_ll = RT0 * rho_eq_ll p_eq_rr = RT0 * rho_eq_rr rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr - u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll, aux_ll) - u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr, aux_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_ll, RT_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_rr, RT_rr) # Compute residual contribution u_res_ll = u_ll - u_eq_ll @@ -1073,15 +1083,15 @@ end psi_eq_rr = psi_rr + phi_rr - phi_star # Compute reconstructed equilibrium state - rho_eq_ll = exp(psi_eq_ll / (RT0)) # Replace RT0 with R*T0 once available + rho_eq_ll = exp(psi_eq_ll / (RT0)) rho_eq_rr = exp(psi_eq_rr / (RT0)) p_eq_ll = RT0 * rho_eq_ll p_eq_rr = RT0 * rho_eq_rr rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_star rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_star - u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_star, aux_ll) - u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_star, aux_rr) + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, rho_e_eq_ll, phi_star, RT_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, rho_e_eq_rr, phi_star, RT_rr) # Compute reconstructed state u_star_ll = u_eq_ll + u_res_ll From cd0336f8a98c3ab8dc15b1e0dcf445c8e23fc5b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 28 Apr 2026 19:22:40 +0200 Subject: [PATCH 13/18] Added Chandrashekar flux for curvilinear solvers and outer state function for P4est meshes --- ...ompressible_euler_gravity_nopressure_2d.jl | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index cb55d7e5..52917ee1 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -586,6 +586,34 @@ end end end +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + nonconservative_type::Trixi.NonConservativeLocal, + nonconservative_term::Integer) + rho_ll, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + + f0 = zero(eltype(u_ll)) + + if nonconservative_term == 1 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll + + e_ll = exp(phi_ll / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = -rho_ll * RT * e_ll + else # nonconservative_term == 2 + noncons = one(eltype(u_ll)) + end + + return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], f0, + f0, f0) +end + @inline function flux_nonconservative_chandrashekar_isothermal(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquationsWithGravityNoPressure2D, @@ -625,6 +653,36 @@ end end end +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + nonconservative_type::Trixi.NonConservativeJump, + nonconservative_term::Integer) + rho_ll, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) + + f0 = zero(eltype(u_ll)) + + if nonconservative_term == 1 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll + + e_ll = exp(phi_ll / RT) + e_rr = exp(phi_rr / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = (1 / e_rr - 1 / e_ll) + + else # nonconservative_term == 2 + noncons = (p_rr - p_ll) + end + + return SVector(f0, noncons, noncons, f0, f0, f0) +end + """ FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsWithGravityNoPressure2D) @@ -1124,4 +1182,21 @@ end u_inner[6]) end end + +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_slip_wall), + normal_direction::AbstractVector, + mesh::P4estMesh{2}, + equations::CompressibleEulerEquationsWithGravityNoPressure2D, + dg, cache, indices...) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction + + return SVector(u_inner[1], + u_inner[2] - 2 * u_normal[1], + u_inner[3] - 2 * u_normal[2], + u_inner[4], + u_inner[5], + u_inner[6]) +end end # @muladd From 615dad682ba64e8a76efad833600e1d0c5a7fcba Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 29 Apr 2026 16:04:50 +0200 Subject: [PATCH 14/18] add callback to compute element mean temperature --- .../elixir_euler_gravity_blast_subcell.jl | 3 + ...lixir_euler_gravity_equilibrium_subcell.jl | 3 + src/TrixiAtmo.jl | 2 + src/callbacks_step/callbacks_step.jl | 1 + src/callbacks_step/mean_temperature.jl | 62 +++++++++++++++++++ 5 files changed, 71 insertions(+) create mode 100644 src/callbacks_step/mean_temperature.jl diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl index af06a397..fe6031c3 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl @@ -194,6 +194,8 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) +mean_temperature_callback = MeanTemperatureCallback() + save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true, @@ -204,6 +206,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.25) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + mean_temperature_callback, save_solution, stepsize_callback) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index 61515ab7..7a2bd2a8 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -70,6 +70,8 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) +mean_temperature_callback = MeanTemperatureCallback() + save_solution = SaveSolutionCallback(dt = 10, save_initial_solution = true, save_final_solution = true, @@ -79,6 +81,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + mean_temperature_callback, save_solution, stepsize_callback) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 41404465..0e0a7df7 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -105,4 +105,6 @@ export AtmosphereLayers, AtmosphereLayersRainyBubble export examples_dir +export MeanTemperatureCallback + end # module TrixiAtmo diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl index 10df9bba..6198454b 100644 --- a/src/callbacks_step/callbacks_step.jl +++ b/src/callbacks_step/callbacks_step.jl @@ -2,3 +2,4 @@ include("analysis_covariant.jl") include("save_solution_covariant.jl") include("stepsize_dg2d.jl") include("save_solution_2d_manifold_in_3d_cartesian.jl") +include("mean_temperature.jl") diff --git a/src/callbacks_step/mean_temperature.jl b/src/callbacks_step/mean_temperature.jl new file mode 100644 index 00000000..854b919a --- /dev/null +++ b/src/callbacks_step/mean_temperature.jl @@ -0,0 +1,62 @@ + +# This method is called to determine whether the callback should be activated +function mean_temperature_callback(u, t, integrator) + # TODO: Here we could do a check and only activate this periodically + return true +end + +# Actual callback function +function mean_temperature_callback(integrator) + semi = integrator.p + mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(semi) + u = Trixi.wrap_array(integrator.u, semi) + + update_mean_temperature!(u, mesh, equations, solver, cache) + + return nothing +end + +# This function computes the mean temperature R*T in each element and stores it in the last variable +# of the solution vector +function update_mean_temperature!(u, mesh::Trixi.AbstractMesh{2}, equations, solver, cache) + @unpack weights = solver.basis + @unpack inverse_jacobian = cache.elements + + Trixi.@threaded for element in eachelement(solver, cache) + # compute the element local mean temperature + RT_mean = zero(eltype(u)) + total_volume = zero(eltype(u)) + for j in eachnode(solver), i in eachnode(solver) + volume_jacobian = abs(inv(Trixi.get_inverse_jacobian(inverse_jacobian, + mesh, + i, j, element))) + u_node = get_node_vars(u, equations, solver, i, j, element) + p_node = pressure(u_node, equations) + rho_node = u_node[1] + RT_node = p_node / rho_node + + RT_mean += RT_node * weights[i] * weights[j] * volume_jacobian + + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + RT_mean = RT_mean / total_volume + + # Update the solution vector + for j in eachnode(solver), i in eachnode(solver) + u[end,i,j,element] = RT_mean + end + end + + return nothing +end + +""" + MeanTemperatureCallback() + + Callback to compute the mean temperature `R*T` in each element and store it in the last variable + of the solution vector. +""" +function MeanTemperatureCallback() + return Trixi.DiscreteCallback(mean_temperature_callback, mean_temperature_callback) +end \ No newline at end of file From 260c2d4dcf98e4f938faa2812cc7638da0ce1298 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 29 Apr 2026 17:49:07 +0200 Subject: [PATCH 15/18] Added 3D version of the Euler equations with gravity and pressure in the noncons term --- src/TrixiAtmo.jl | 1 + ...ompressible_euler_gravity_nopressure_3d.jl | 685 ++++++++++++++++++ src/equations/equations.jl | 1 + 3 files changed, 687 insertions(+) create mode 100644 src/equations/compressible_euler_gravity_nopressure_3d.jl diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 0e0a7df7..0a043c28 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -61,6 +61,7 @@ export CompressibleMoistEulerEquations2D, CompressibleEulerEnergyEquationsWithGravity3D, CompressibleEulerEquationsWithGravity2D, CompressibleEulerEquationsWithGravityNoPressure2D, + CompressibleEulerEquationsWithGravityNoPressure3D, CompressibleEulerInternalEnergyEquationsWithGravity2D export GlobalCartesianCoordinates, GlobalSphericalCoordinates diff --git a/src/equations/compressible_euler_gravity_nopressure_3d.jl b/src/equations/compressible_euler_gravity_nopressure_3d.jl new file mode 100644 index 00000000..95552c7e --- /dev/null +++ b/src/equations/compressible_euler_gravity_nopressure_3d.jl @@ -0,0 +1,685 @@ +# 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""" + CompressibleEulerEquationsWithGravityNoPressure3D(gamma) + +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_{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_{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_{tot} +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_{tot}`` 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_{tot} - \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 CompressibleEulerEquationsWithGravityNoPressure3D{RealT <: Real} <: + Trixi.AbstractCompressibleEulerEquations{3, 7} + 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 CompressibleEulerEquationsWithGravityNoPressure3D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravityNoPressure3D) = True() + +# The auxiliary variable is used to store the mean temperature of the element for isothermal-equilibrium-preserving discretizations +Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravityNoPressure3D) = ("rho", + "rho_v1", + "rho_v2", + "rho_v3", + "rho_etot", + "phi", + "aux") +Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravityNoPressure3D) = ("rho", + "v1", + "v2", + "v3", + "p", + "phi", + "aux") + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + +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::CompressibleEulerEquationsWithGravityNoPressure3D) + 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 + rho_local, v1, v2, v3, p_local, _ = cons2prim(u_inner, equations) + + v_normal = v1 * normal[1] + v2 * normal[2] + v3 * normal[3] + + # 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 Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0 + sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5f0 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0 + A = 2 / ((equations.gamma + 1) * rho_local) + B = p_local * (equations.gamma - 1) / (equations.gamma + 1) + p_star = p_local + + 0.5f0 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # Compute non-conservative term: Since the geopotential is a continuous function, it does not act at the BC + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + # TODO: compute with surface_flux_function + # surface_flux, surface_noncons = surface_flux_function + # surface_noncons(u_ll, u_rr, + # normal_direction::AbstractVector, + # equations::CompressibleEulerEquationsWithGravityNoPressure3D) + noncons = (p_star - p_local) + + # For the slip wall we directly set the flux as the normal velocity is zero + return (SVector(zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_, + SVector(zero(eltype(u_inner)), + noncons * normal[1], + noncons * normal[2], + noncons * normal[3], + 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::CompressibleEulerEquationsWithGravityNoPressure3D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(one(eltype(u_inner)), zero(eltype(u_inner)), + zero(eltype(u_inner))) + elseif orientation == 2 + normal_direction = SVector(zero(eltype(u_inner)), one(eltype(u_inner)), + zero(eltype(u_inner))) + else # orientation == 3 + normal_direction = SVector(zero(eltype(u_inner)), zero(eltype(u_inner)), + one(eltype(u_inner))) + 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::CompressibleEulerEquationsWithGravityNoPressure3D) + +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::CompressibleEulerEquationsWithGravityNoPressure3D) +# # 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 3D 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::CompressibleEulerEquationsWithGravityNoPressure3D) + rho_etot = u[5] + rho, v1, v2, v3, p, _ = cons2prim(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + + v3 * normal_direction[3] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = rho_v_normal * v1 + f3 = rho_v_normal * v2 + f4 = rho_v_normal * v3 + f5 = (rho_etot + p) * v_normal + return SVector(f1, f2, f3, f4, f5, zero(eltype(u)), zero(eltype(u))) +end + +@inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + # Unpack left and right state + rho_etot_ll = u_ll[5] + rho_etot_rr = u_rr[5] + rho_ll, v1_ll, v2_ll, v3_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Average each factor of products in flux + rho_avg = 0.5f0 * (rho_ll + rho_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + + v3_avg * normal_direction[3] + 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 + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * v3_avg + f5 = f1 * etot_avg + p_avg * v_dot_n_avg + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll)), zero(eltype(u_ll))) +end + +@inline function (noncons_flux::FluxNonConservativeChandrashekarIsothermal)(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho_ll, _, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + _, _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) + + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll + + e_ll = exp(phi_ll / RT) + e_rr = exp(phi_rr / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = -rho_ll * RT * e_ll * (1 / e_rr - 1 / e_ll) + (p_rr - p_ll) + + f0 = zero(eltype(u_ll)) + return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], + noncons * normal_direction[3], f0, + f0, f0) +end + +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D, + nonconservative_type::Trixi.NonConservativeLocal, + nonconservative_term::Integer) + rho_ll, _, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + + f0 = zero(eltype(u_ll)) + + if nonconservative_term == 1 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll + + e_ll = exp(phi_ll / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = -rho_ll * RT * e_ll + else # nonconservative_term == 2 + noncons = one(eltype(u_ll)) + end + + return SVector(f0, noncons * normal_direction[1], noncons * normal_direction[2], + noncons * normal_direction[3], f0, + f0, f0) +end + +@inline function flux_nonconservative_chandrashekar_isothermal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D, + nonconservative_type::Trixi.NonConservativeJump, + nonconservative_term::Integer) + rho_ll, _, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + _, _, _, _, p_rr, phi_rr, _ = cons2prim(u_rr, equations) + + f0 = zero(eltype(u_ll)) + + if nonconservative_term == 1 + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT = RT_ll + + e_ll = exp(phi_ll / RT) + e_rr = exp(phi_rr / RT) + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = (1 / e_rr - 1 / e_ll) + + else # nonconservative_term == 2 + noncons = (p_rr - p_ll) + end + + return SVector(f0, noncons, noncons, noncons, f0, f0, f0) +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + +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) +""" +# # We add the "upwinding" pressure terms here, but not the average pressure terms as in CompressibleEulerEquationsWithGravity2D +# @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, +# equations::CompressibleEulerEquationsWithGravityNoPressure3D) +# c = flux_lmars.speed_of_sound + +# # Unpack left and right state +# rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) +# rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + +# v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +# v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + +# # 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 * 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 +# f0 = zero(eltype(u_ll)) +# return SVector(f1, +# f2 + p * normal_direction[1], +# f3 + p * normal_direction[2], +# f4, f0, f0) +# end + +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr, _ = cons2prim(u_rr, equations) + + v_normal_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + v_normal_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + + 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 + +@inline function Trixi.max_abs_speeds(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho, v1, v2, v3, p, _ = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c, abs(v3) + c +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho, rho_v1, rho_v2, rho_v3, rho_etot, phi, aux = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + p = (equations.gamma - 1) * + (rho_etot - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - rho * phi) + + return SVector(rho, v1, v2, v3, p, phi, aux) +end + +# Convert conservative variables to entropy (see, e.g., Waruszewski et al. (2022)) +@inline function Trixi.cons2entropy(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho, v1, v2, v3, p, phi, aux = cons2prim(u, equations) + + v_square = v1^2 + v2^2 + v3^2 + s = log(p) - equations.gamma * log(rho) + rho_p = rho / p + + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + rho_p * (0.5f0 * v_square - phi) + w2 = rho_p * v1 + w3 = rho_p * v2 + w4 = rho_p * v3 + w5 = -rho_p + + return SVector(w1, w2, w3, w4, w5, phi, aux) +end + +@inline function Trixi.entropy2cons(w, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + # See Waruszewski et al. (2022) + @unpack gamma = equations + + # convert to entropy `-rho * s` + # instead of `-rho * s / (gamma - 1)` + V1, V2, V3, V4, V5, _ = w .* (gamma - 1) + phi = w[6] + aux = w[7] + + # s = specific entropy + s = gamma - V1 + (V2^2 + V3^2 + V4^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_v3 = rho_iota * V4 + rho_etot = rho_iota * (1 - (V2^2 + V3^2 + V4^2) / (2 * V5)) + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_v3, rho_etot, phi, aux) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho, v1, v2, v3, p, phi, aux = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_v3 = rho * v3 + rho_etot = p * equations.inv_gamma_minus_one + + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_v3, rho_etot, phi, aux) +end + +@inline function Trixi.density(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho = u[1] + return rho +end + +@inline function Trixi.pressure(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho, rho_v1, rho_v2, rho_v3, rho_etot, phi, _ = u + p = (equations.gamma - 1) * + (rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - rho * phi) + return p +end + +# @inline function Trixi.density_pressure(u, +# equations::CompressibleEulerEquationsWithGravityNoPressure3D) +# rho, rho_v1, rho_v2, rho_etot, phi, _ = u +# rho_times_p = (equations.gamma - 1) * +# (rho * rho_etot - 0.5f0 * (rho_v1^2 + rho_v2^2) - rho^2 * phi) +# return rho_times_p +# end + +# # Calculate thermodynamic entropy for a conservative state `cons` +# @inline function entropy_thermodynamic(cons, +# equations::CompressibleEulerEquationsWithGravityNoPressure3D) +# # Pressure +# p = (equations.gamma - 1) * +# (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]) + +# return s +# end + +# # Calculate mathematical entropy for a conservative state `cons` +# @inline function entropy_math(cons, +# equations::CompressibleEulerEquationsWithGravityNoPressure3D) +# # 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::CompressibleEulerEquationsWithGravityNoPressure3D) = entropy_math(cons, +# equations) + +# Calculate total energy for a conservative state `cons` +@inline Trixi.energy_total(cons, ::CompressibleEulerEquationsWithGravityNoPressure3D) = cons[5] + +# # Calculate kinetic energy for a conservative state `cons` +# @inline function Trixi.energy_kinetic(u, +# equations::CompressibleEulerEquationsWithGravityNoPressure3D) +# rho, rho_v1, rho_v2, rho_etot, _ = 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::CompressibleEulerEquationsWithGravityNoPressure3D) +# return energy_total(cons, equations) - energy_kinetic(cons, equations) - +# cons[1] * cons[5] +# end + +@inline function Trixi.velocity(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the +# gravitational potential +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5f0 * λ * (u_rr - u_ll) + f0 = zero(eltype(u_ll)) + return SVector(diss[1], diss[2], diss[3], diss[4], diss[5], f0, f0) +end + +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + if u[1] <= 0 || pressure(u, equations) <= 0 + return false + end + return true +end + +####################################################################### +# Hydrostatic reconstruction + +# Hydrostatic reconstruction from the paper: +# Ziming Chen, Yingjuan Zhang, Gang Li, Shouguo Qian (2022) +# "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal +# hydrostatic state under gravitational field" +# [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) +@inline function hydrostatic_reconstruction(u_ll, u_rr, + equations::CompressibleEulerEquationsWithGravityNoPressure3D) + # Unpack left and right states + rho_ll, _, _, _, p_ll, phi_ll, RT_ll = cons2prim(u_ll, equations) + rho_rr, _, _, _, p_rr, phi_rr, RT_rr = cons2prim(u_rr, equations) + + rho_e_ll = u_ll[5] + rho_e_rr = u_rr[5] + + psi_ll = p_ll / rho_ll * log(rho_ll) + psi_rr = p_rr / rho_rr * log(rho_rr) + + # Compute equilibrium potential (# TODO: For general case we need to add phi_num - phi_exact)) + psi_eq_ll = psi_ll # + phi_ll - phi_exact_ll + psi_eq_rr = psi_rr # + phi_rr - phi_exact_rr + + # We use the mean element temperature on the left to keep the iso-thermal correction + # completely element-local, as in: + # - Chandrashekar, P., & Zenk, M. (2017). Well-balanced nodal discontinuous Galerkin method for Euler + # equations with gravity. Journal of Scientific Computing, 71(3), 1062-1093. + RT0 = RT_ll + + # Compute equlibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0 * rho_eq_ll + p_eq_rr = RT0 * rho_eq_rr + rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_ll + rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_rr + + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, 0.0, rho_e_eq_ll, phi_ll, RT_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, 0.0, rho_e_eq_rr, phi_rr, RT_rr) + + # Compute residual contribution + u_res_ll = u_ll - u_eq_ll + u_res_rr = u_rr - u_eq_rr + + # Reconstruct phi + phi_star = max(phi_ll, phi_rr) + + # Compute reconstructed equilibrium potential + psi_eq_ll = psi_ll + phi_ll - phi_star + psi_eq_rr = psi_rr + phi_rr - phi_star + + # Compute reconstructed equilibrium state + rho_eq_ll = exp(psi_eq_ll / (RT0)) + rho_eq_rr = exp(psi_eq_rr / (RT0)) + p_eq_ll = RT0 * rho_eq_ll + p_eq_rr = RT0 * rho_eq_rr + rho_e_eq_ll = p_eq_ll / (equations.gamma - 1) + rho_eq_ll * phi_star + rho_e_eq_rr = p_eq_rr / (equations.gamma - 1) + rho_eq_rr * phi_star + + u_eq_ll = SVector(rho_eq_ll, 0.0, 0.0, 0.0, rho_e_eq_ll, phi_star, RT_ll) + u_eq_rr = SVector(rho_eq_rr, 0.0, 0.0, 0.0, rho_e_eq_rr, phi_star, RT_rr) + + # Compute reconstructed state + u_star_ll = u_eq_ll + u_res_ll + u_star_rr = u_eq_rr + u_res_rr + + return u_star_ll, u_star_rr +end + +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_slip_wall), + normal_direction::AbstractVector, + mesh::P4estMesh{3}, + equations::CompressibleEulerEquationsWithGravityNoPressure3D, + dg, cache, indices...) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3] + + normal_direction[3] * u_inner[4]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction + + return SVector(u_inner[1], + u_inner[2] - 2 * u_normal[1], + u_inner[3] - 2 * u_normal[2], + u_inner[4] - 2 * u_normal[3], + u_inner[5], + u_inner[6], + u_inner[7]) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 27743404..7970cd09 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -355,4 +355,5 @@ include("shallow_water_3d.jl") include("reference_data.jl") include("compressible_euler_gravity_2d.jl") include("compressible_euler_gravity_nopressure_2d.jl") +include("compressible_euler_gravity_nopressure_3d.jl") end # @muladd From 4920db5611fcdc7f2751edff91be4ed0ade3599d Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 30 Apr 2026 08:59:37 +0200 Subject: [PATCH 16/18] rename hydrostatic_reconstruction -> hydrostatic_reconstruction_isothermal --- .../elixir_euler_gravity_blast_subcell.jl | 4 +- ...lixir_euler_gravity_equilibrium_subcell.jl | 4 +- hydrostatic_reconstruction.jl | 78 ------------------- src/TrixiAtmo.jl | 2 +- ...ompressible_euler_gravity_nopressure_2d.jl | 4 +- 5 files changed, 7 insertions(+), 85 deletions(-) delete mode 100644 hydrostatic_reconstruction.jl diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl index fe6031c3..038df903 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_blast_subcell.jl @@ -146,9 +146,9 @@ surface_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isotherm surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, DissipationLocalLaxFriedrichs()), - hydrostatic_reconstruction), + hydrostatic_reconstruction_isothermal), FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, - hydrostatic_reconstruction)) + hydrostatic_reconstruction_isothermal)) polydeg = 3 basis = LobattoLegendreBasis(polydeg) diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl index 7a2bd2a8..43ee60e6 100644 --- a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium_subcell.jl @@ -23,9 +23,9 @@ volume_flux = (flux_kennedy_gruber, flux_nonconservative_chandrashekar_isotherma surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_kennedy_gruber, DissipationLocalLaxFriedrichs()), - hydrostatic_reconstruction), + hydrostatic_reconstruction_isothermal), FluxHydrostaticReconstruction(flux_nonconservative_chandrashekar_isothermal, - hydrostatic_reconstruction)) + hydrostatic_reconstruction_isothermal)) polydeg = 0 basis = LobattoLegendreBasis(polydeg) diff --git a/hydrostatic_reconstruction.jl b/hydrostatic_reconstruction.jl deleted file mode 100644 index b476c1a6..00000000 --- a/hydrostatic_reconstruction.jl +++ /dev/null @@ -1,78 +0,0 @@ -struct FluxHydrostaticReconstruction{NumericalFlux, HydrostaticReconstruction} - numerical_flux::NumericalFlux - hydrostatic_reconstruction::HydrostaticReconstruction -end - -@inline function (numflux::FluxHydrostaticReconstruction)(u_ll, u_rr, - orientation_or_normal_direction, - equations::Trixi.AbstractEquations) - @unpack numerical_flux, hydrostatic_reconstruction = numflux - - # Create the reconstructed left/right solution states in conservative form - u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) - - # Use the reconstructed states to compute the numerical surface flux - return numerical_flux(u_ll_star, u_rr_star, orientation_or_normal_direction, - equations) -end - -# Hydrostatic reconstruction from the paper: -# Ziming Chen, Yingjuan Zhang, Gang Li, Shouguo Qian (2022) -# "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal -# hydrostatic state under gravitational field" -# [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) -@inline function hydrostatic_reconstruction(u_ll, u_rr, equations::CompressibleEulerEquationsWithGravityNoPressure2D) - # Unpack left and right states - rho_ll, rho_v_ll, rho_e_ll, phi_ll = u_ll - rho_rr, rho_v_rr, rho_e_rr, phi_rr = u_rr - - # Compute equilibrium potential (# TODO: For general case we need to add phi_num - phi_exact)) - psi_eq_ll = psi_ll # phi_ll - phi_exact_ll - psi_eq_rr = psi_rr # phi_rr - phi_exact_rr - - # Set RT0 to 1 for now. - RT0 = 1 - - # Compute equlibrium state - rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available - rho_eq_rr = exp(psi_eq_rr / (RT0)) - p_eq_ll = RT0*rho_eq_ll - p_eq_rr = RT0*rho_eq_rr - rho_v_eq_ll = 0.0 - rho_v_eq_rr = 0.0 - rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll - rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr - - u_eq_ll = (rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) - u_eq_rr = (rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) - - # Compute residual contribution - u_res_ll = u_ll - u_eq_ll - u_res_rr = u_rr - u_eq_rr - - # Reconstruct phi - phi_star = max(phi_ll, phi_rr) - - # Compute reconstructed equilibrium potential - psi_eq_ll = psi_ll + phi_ll - phi_star - psi_eq_rr = psi_rr + phi_rr - phi_star - - # Compute reconstructed equilibrium state - rho_eq_ll = exp(psi_eq_ll / (RT0)) # How do I get T0? # Replace RT0 with R*T0 once available - rho_eq_rr = exp(psi_eq_rr / (RT0)) - p_eq_ll = RT0*rho_eq_ll - p_eq_rr = RT0*rho_eq_rr - rho_v_eq_ll = 0.0 - rho_v_eq_rr = 0.0 - rho_e_eq_ll = p_eq_ll / (gamma - 1) + rho_eq_ll * phi_ll - rho_e_eq_rr = p_eq_rr / (gamma - 1) + rho_eq_rr * phi_rr - - u_eq_ll = SVector(rho_eq_ll, rho_v_eq_ll, rho_e_eq_ll, phi_ll) - u_eq_rr = SVector(rho_eq_rr, rho_v_eq_rr, rho_e_eq_rr, phi_rr) - - # Compute reconstructed state - u_star_ll = u_eq_ll + u_res_ll - u_star_rr = u_eq_rr + u_res_rr - - return u_star_ll, u_star_rr -end \ No newline at end of file diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 0a043c28..7f842428 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -69,7 +69,7 @@ export GlobalCartesianCoordinates, GlobalSphericalCoordinates export NonlinearSolveDG export flux_chandrashekar, FluxLMARS, FluxHydrostaticReconstruction, - hydrostatic_reconstruction + hydrostatic_reconstruction_isothermal export flux_nonconservative_waruszewski, flux_nonconservative_chandrashekar_isothermal diff --git a/src/equations/compressible_euler_gravity_nopressure_2d.jl b/src/equations/compressible_euler_gravity_nopressure_2d.jl index 52917ee1..5cb1d19f 100644 --- a/src/equations/compressible_euler_gravity_nopressure_2d.jl +++ b/src/equations/compressible_euler_gravity_nopressure_2d.jl @@ -1090,12 +1090,12 @@ end equations) end -# Hydrostatic reconstruction from the paper: +# Hydrostatic reconstruction for the isothermal equilibrium from the paper: # Ziming Chen, Yingjuan Zhang, Gang Li, Shouguo Qian (2022) # "A well-balanced Runge-Kutta discontinuous Galerkin method for the Euler equations in isothermal # hydrostatic state under gravitational field" # [DOI:10.1016/j.camwa.2022.05.025](https://doi.org/10.1016/j.camwa.2022.05.025) -@inline function hydrostatic_reconstruction(u_ll, u_rr, +@inline function hydrostatic_reconstruction_isothermal(u_ll, u_rr, equations::Union{CompressibleEulerEquationsWithGravityNoPressure2D, CompressibleEulerEquationsWithGravity2D}) # Unpack left and right states From b8d78790394e7985a0ac0159a9c4f8e3db31a367 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 4 May 2026 09:54:56 +0200 Subject: [PATCH 17/18] add 3D version of mean temperature callback --- src/callbacks_step/mean_temperature.jl | 36 ++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/callbacks_step/mean_temperature.jl b/src/callbacks_step/mean_temperature.jl index 854b919a..a76ffb24 100644 --- a/src/callbacks_step/mean_temperature.jl +++ b/src/callbacks_step/mean_temperature.jl @@ -51,6 +51,42 @@ function update_mean_temperature!(u, mesh::Trixi.AbstractMesh{2}, equations, sol return nothing end +# This function computes the mean temperature R*T in each element and stores it in the last variable +# of the solution vector +# TODO: Setup 3D test case to validate this function +function update_mean_temperature!(u, mesh::Trixi.AbstractMesh{3}, equations, solver, cache) + @unpack weights = solver.basis + @unpack inverse_jacobian = cache.elements + + Trixi.@threaded for element in eachelement(solver, cache) + # compute the element local mean temperature + RT_mean = zero(eltype(u)) + total_volume = zero(eltype(u)) + for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver) + volume_jacobian = abs(inv(Trixi.get_inverse_jacobian(inverse_jacobian, + mesh, + i, j, k, element))) + u_node = get_node_vars(u, equations, solver, i, j, k, element) + p_node = pressure(u_node, equations) + rho_node = u_node[1] + RT_node = p_node / rho_node + + RT_mean += RT_node * weights[i] * weights[j] * weights[k] * volume_jacobian + + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + end + # normalize with the total volume + RT_mean = RT_mean / total_volume + + # Update the solution vector + for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver) + u[end,i,j,k,element] = RT_mean + end + end + + return nothing +end + """ MeanTemperatureCallback() From 9da83de503c87b71ebc9183f97618f78efa0bf9b Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 6 May 2026 09:19:15 +0200 Subject: [PATCH 18/18] avoid saving intermediate solutions when using MeanTemperatureCallback --- src/callbacks_step/mean_temperature.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/mean_temperature.jl b/src/callbacks_step/mean_temperature.jl index a76ffb24..d313cd84 100644 --- a/src/callbacks_step/mean_temperature.jl +++ b/src/callbacks_step/mean_temperature.jl @@ -94,5 +94,6 @@ end of the solution vector. """ function MeanTemperatureCallback() - return Trixi.DiscreteCallback(mean_temperature_callback, mean_temperature_callback) + return Trixi.DiscreteCallback(mean_temperature_callback, mean_temperature_callback, + save_positions = (false, false)) end \ No newline at end of file