diff --git a/.typos.toml b/.typos.toml new file mode 100644 index 00000000..97e44980 --- /dev/null +++ b/.typos.toml @@ -0,0 +1,2 @@ +[default.extend-words] +claus = "claus" diff --git a/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl new file mode 100644 index 00000000..3e211abb --- /dev/null +++ b/examples/euler/dry_air/buoyancy/elixir_euler_gravity_warmbubble.jl @@ -0,0 +1,136 @@ + +using OrdinaryDiffEqLowStorageRK +using Trixi, TrixiAtmo + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquationsWithGravity2D) + RealT = eltype(x) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000 + center_z = 2000 + # radius of perturbation + radius = 2000 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300 + potential_temperature_perturbation = zero(RealT) + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5f0 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + # T = potential_temperature - g / (c_p) * x[2] + + # density + rho = p / (R * T) + + # Geopotential + phi = g * x[2] + + v1 = 20 + v2 = 0 + E = c_v * T + 0.5f0 * (v1^2 + v2^2) + phi + return SVector(rho, rho * v1, rho * v2, rho * E, phi) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma) + +initial_condition = warm_bubble_setup + +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski_etal) +surface_flux = (FluxLMARS(340.0), flux_zero) + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, -5000.0) +coordinates_max = (20_000.0, 15_000.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000, + periodicity = (true, false)) + +boundary_conditions = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 10.0, #interval = 1, #dt = 10.0, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl new file mode 100644 index 00000000..38fb7398 --- /dev/null +++ b/examples/euler/dry_air/mountain_flow/elixir_euler_gravity_schaer_mountains.jl @@ -0,0 +1,134 @@ +using OrdinaryDiffEqLowStorageRK +using TrixiAtmo +using Trixi + +#schaer mountain test case + +# Initial condition +function initial_condition_schaer_mountain(x, t, + equations::CompressibleEulerEquationsWithGravity2D) + g = 9.81 + c_p = 1004.0 + c_v = 717.0 + gamma = c_p / c_v + + # Exner pressure from hydrostatic balance for x[2] + potential_temperature_int = 280.0 #constant of integration + bvfrequency = 0.01 #Brunt-Väisälä frequency + + exner = 1 + + g^2 / (c_p * potential_temperature_int * bvfrequency^2) * + (exp(-bvfrequency^2 / g * x[2]) - 1) + + # mean potential temperature + potential_temperature = potential_temperature_int * exp(bvfrequency^2 / g * x[2]) + + # temperature + T = potential_temperature * exner + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # density + rho = p / (R * T) + + #Geopotential + phi = g * x[2] + + v1 = 10.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + phi + return SVector(rho, rho * v1, rho * v2, rho * E, phi) +end + +############################################################################### + +function mapping(xi_, eta_) + xi = xi_ * 25_000 #xi_ * 10_000.0 + eta = eta_ * 10_500 + 10_500# eta_ * 500.0 + 500.0 + #upper boundary + H = 21_000.0 + + #topography + h_c = 250.0 + lambda_c = 4000.0 + a_c = 5000.0 + + topo = -h_c * exp(-(xi / a_c)^2) * cos(pi * xi / lambda_c)^2 + + x = xi + y = H * (eta - topo) / (H - topo) + return SVector(x, y) +end + +# Create curved mesh with 200 x 100 elements +polydeg = 3 +cells_per_dimension = (60, 30) +mesh = P4estMesh(cells_per_dimension; polydeg = polydeg, mapping = mapping, + periodicity = false) + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravity2D(1004.0 / 717.0) + +initial_condition = initial_condition_schaer_mountain + +boundary_conditions_dirichlet = (x_neg = BoundaryConditionDirichlet(initial_condition_schaer_mountain), + x_pos = BoundaryConditionDirichlet(initial_condition_schaer_mountain), + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +basis = LobattoLegendreBasis(polydeg) + +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski_etal) +surface_flux = (FluxLMARS(340.0), flux_zero) + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-25_000.0, 0.0) +coordinates_max = (25_000.0, 21_000.0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions_dirichlet) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 60 * 60)# * 10) # 10h = 36000 s + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +solution_variables = cons2prim + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = solution_variables) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl new file mode 100644 index 00000000..491d3c05 --- /dev/null +++ b/examples/euler/dry_air/tests/elixir_euler_gravity_equilibrium.jl @@ -0,0 +1,74 @@ +using OrdinaryDiffEqLowStorageRK +using TrixiAtmo +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquationsWithGravity2D(1.4) + +function initial_condition_constant(x, t, + equations::CompressibleEulerEquationsWithGravity2D) + rho = exp(-x[2]) + v1 = 0.0 + v2 = 0.0 + p = exp(-x[2]) + prim = SVector(rho, v1, v2, p, x[2]) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_constant + +volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski_etal) +surface_flux = (flux_lax_friedrichs, flux_zero) + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10_000, + periodicity = false) + +boundary_conditions = (; x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 9d72ed0c..021c5226 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -59,6 +59,7 @@ export CompressibleMoistEulerEquations2D, CompressibleEulerPotentialTemperatureEquationsWithGravity3D, CompressibleEulerEnergyEquationsWithGravity2D, CompressibleEulerEnergyEquationsWithGravity3D, + CompressibleEulerEquationsWithGravity2D CompressibleEulerInternalEnergyEquationsWithGravity2D export GlobalCartesianCoordinates, GlobalSphericalCoordinates @@ -67,6 +68,8 @@ export NonlinearSolveDG export flux_chandrashekar, FluxLMARS +export flux_nonconservative_waruszewski + export flux_nonconservative_zeros, flux_nonconservative_ec, flux_nonconservative_surface_simplified, source_terms_geometric_coriolis, source_terms_coriolis, source_terms_coriolis_lagrange_multiplier, diff --git a/src/equations/compressible_euler_gravity_2d.jl b/src/equations/compressible_euler_gravity_2d.jl new file mode 100644 index 00000000..dbd4c343 --- /dev/null +++ b/src/equations/compressible_euler_gravity_2d.jl @@ -0,0 +1,850 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent +@doc raw""" + CompressibleEulerEquationsWithGravity2D(gamma) + +The compressible Euler equations with gravity 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 CompressibleEulerEquationsWithGravity2D{RealT <: Real} <: + Trixi.AbstractCompressibleEulerEquations{2, 5} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + function CompressibleEulerEquationsWithGravity2D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravity2D) = True() +Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravity2D) = ("rho", + "rho_v1", + "rho_v2", + "rho_etot", + "phi") +Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravity2D) = ("rho", + "v1", + "v2", + "p", + "phi") + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + +Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book +- Eleuterio F. Toro (2009) + Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + 3rd edition + [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + +Should be used together with [`UnstructuredMesh2D`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_local = Trixi.rotate_to_x(u_inner, normal, equations) + + # compute the primitive variables + rho_local, v_normal, v_tangent, p_local, _ = cons2prim(u_local, equations) + + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A 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 + + # For the slip wall we directly set the flux as the normal velocity is zero + return (SVector(zero(eltype(u_inner)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_, + SVector(zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_) +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquationsWithGravity2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + # get the appropriate normal vector from the orientation + if orientation == 1 + normal_direction = SVector(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::CompressibleEulerEquationsWithGravity2D) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity2D) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + fluxes = boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, equations) + boundary_flux = (-fluxes[1], -fluxes[2]) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux +end + +# Calculate 2D flux for a single point +@inline function Trixi.flux(u, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_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 + p + f3 = rho_v1 * v2 + f4 = (rho_etot + p) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + 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::CompressibleEulerEquationsWithGravity2D) + 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 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + 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::CompressibleEulerEquationsWithGravity2D) + +This flux is is a modification of the original kinetic energy preserving two-point flux by +- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) + Kinetic energy and entropy preserving schemes for compressible flows + by split convective forms + [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058) + +The modification is in the energy flux to guarantee pressure equilibrium and was developed by +- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) + Preventing spurious pressure oscillations in split convective form discretizations for + compressible flows + [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060) +""" +@inline function Trixi.flux_shima_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, 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 + p_avg + f3 = f1 * v2_avg + f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + + f1 * phi_avg + else + pv2_avg = 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + f1 = rho_avg * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + 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::CompressibleEulerEquationsWithGravity2D) + # 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 + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = (f1 * velocity_square_avg + + p_avg * v_dot_n_avg * equations.inv_gamma_minus_one + + 0.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::CompressibleEulerEquationsWithGravity2D) + +Kinetic energy preserving two-point flux by +- Kennedy and Gruber (2008) + Reduced aliasing formulations of the convective terms within the + Navier-Stokes equations for a compressible fluid + [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020) +""" +@inline function Trixi.flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_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 + p_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 + p_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::CompressibleEulerEquationsWithGravity2D) + # 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 + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + 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::CompressibleEulerEquationsWithGravity2D) + +Entropy conserving and kinetic energy preserving two-point flux by +- Hendrik Ranocha (2018) + Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods + for Hyperbolic Balance Laws + [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) +See also +- Hendrik Ranocha (2020) + Entropy Conserving and Kinetic Energy Preserving Numerical Methods for + the Euler Equations Using Summation-by-Parts Operators + [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) +""" +@inline function Trixi.flux_ranocha(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, 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) + 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 orientation + if orientation == 1 + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + p_avg + f3 = f1 * v2_avg + f4 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + f1 * phi_avg + else + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + f4 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + 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::CompressibleEulerEquationsWithGravity2D) + # 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) + 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_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + + 0.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::CompressibleEulerEquationsWithGravity2D) + rho_ll, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, phi_rr = u_rr + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) + + 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::CompressibleEulerEquationsWithGravity2D) + rho_ll, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, phi_rr = u_rr + + # We omit the 0.5 in the density average since Trixi.jl always multiplies the non-conservative flux with 0.5 + noncons = ln_mean(rho_ll, rho_rr) * (phi_rr - phi_ll) + + f0 = zero(eltype(u_ll)) + if orientation == 1 + return SVector(f0, noncons, f0, f0, f0) + else #if orientation == 2 + return SVector(f0, f0, noncons, f0, f0) + end +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" +# The struct is already defined in CompressibleEulerEquations2D + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = 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 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4, _ = v * u_ll + f4 = f4 + p_ll * v + else + f1, f2, f3, f4, _ = v * u_rr + f4 = f4 + p_rr * v + end + + if orientation == 1 + f2 = f2 + p + else # orientation == 2 + f3 = f3 + p + end + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll, _ = 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 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4, _ = u_ll * v + f4 = f4 + p_ll * v + else + f1, f2, f3, f4, _ = u_rr * v + f4 = f4 + p_rr * v + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4, zero(eltype(u_ll))) +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Calculate minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + if orientation == 1 # x-direction + λ_min = v1_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v1_rr + sqrt(equations.gamma * p_rr / rho_rr) + else # y-direction + λ_min = v2_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v2_rr + sqrt(equations.gamma * p_rr / rho_rr) + end + + return λ_min, λ_max +end + +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity2D) + rho_ll, v1_ll, v2_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, _ = cons2prim(u_rr, equations) + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + # The v_normals are already scaled by the norm + λ_min = v_normal_ll - sqrt(equations.gamma * p_ll / rho_ll) * norm_ + λ_max = v_normal_rr + sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + return λ_min, λ_max +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this rotation of the state vector +@inline function Trixi.rotate_to_x(u, normal_vector, + equations::CompressibleEulerEquationsWithGravity2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4], + u[5]) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function Trixi.rotate_from_x(u, normal_vector, + equations::CompressibleEulerEquationsWithGravity2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 t_1 0; + # 0 n_2 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] - s * u[3], + s * u[2] + c * u[3], + u[4], + u[5]) +end + +@inline function Trixi.max_abs_speeds(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, v1, v2, p, _ = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_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::CompressibleEulerEquationsWithGravity2D) + 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::CompressibleEulerEquationsWithGravity2D) + # 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::CompressibleEulerEquationsWithGravity2D) + 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::CompressibleEulerEquationsWithGravity2D) + rho = u[1] + return rho +end + +@inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravity2D) + 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::CompressibleEulerEquationsWithGravity2D) + 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::CompressibleEulerEquationsWithGravity2D) + # 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::CompressibleEulerEquationsWithGravity2D) + # Mathematical entropy + S = -entropy_thermodynamic(cons, equations) * cons[1] * + equations.inv_gamma_minus_one + + return S +end + +# Default entropy is the mathematical entropy +@inline Trixi.entropy(cons, equations::CompressibleEulerEquationsWithGravity2D) = entropy_math(cons, + equations) + +# Calculate total energy for a conservative state `cons` +@inline Trixi.energy_total(cons, ::CompressibleEulerEquationsWithGravity2D) = cons[4] + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, + equations::CompressibleEulerEquationsWithGravity2D) + rho, rho_v1, rho_v2, rho_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::CompressibleEulerEquationsWithGravity2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) - + cons[1] * cons[5] +end + +@inline function Trixi.velocity(u, equations::CompressibleEulerEquationsWithGravity2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + return SVector(v1, v2) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the +# gravitational potential +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity2D) + λ = 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 diff --git a/src/equations/equations.jl b/src/equations/equations.jl index e045f68c..623ee48c 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -353,4 +353,5 @@ include("compressible_euler_energy_with_gravity_3d.jl") include("compressible_euler_internal_energy_with_gravity_2d.jl") include("shallow_water_3d.jl") include("reference_data.jl") +include("compressible_euler_gravity_2d.jl") end # @muladd diff --git a/test/runtests.jl b/test/runtests.jl index 4a3a57d9..4d2c35ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -69,6 +69,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) end end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_gravity" + include("test_2d_euler_gravity.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "spherical_advection" @testset verbose=true showtiming=true "Spherical advection tests" begin include("test_spherical_advection.jl") diff --git a/test/test_2d_euler_gravity.jl b/test/test_2d_euler_gravity.jl new file mode 100644 index 00000000..8dc9bc53 --- /dev/null +++ b/test/test_2d_euler_gravity.jl @@ -0,0 +1,58 @@ +module Test2DEulerGravity + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = joinpath(EXAMPLES_DIR, "euler/dry_air") + +@trixi_testset "elixir_euler_gravity_equilibrium (Shima + LLF)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixi_testset "elixir_euler_gravity_equilibrium (Kennedy-Gruber)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0], + volume_flux=(flux_kennedy_gruber, + flux_nonconservative_waruszewski_etal), + surface_flux=(flux_kennedy_gruber, + flux_nonconservative_waruszewski_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixi_testset "elixir_euler_gravity_equilibrium (Ranocha)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tests", + "elixir_euler_gravity_equilibrium.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0], + volume_flux=(flux_ranocha, flux_nonconservative_waruszewski_etal), + surface_flux=(flux_ranocha, flux_nonconservative_waruszewski_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +end # module diff --git a/test/test_type.jl b/test/test_type.jl index 81bba1bd..974f3519 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -198,6 +198,82 @@ end end end +@timed_testset "Compressible Euler Total Energy With Gravity 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquationsWithGravity2D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), zero(RealT), zero(RealT), + one(RealT), one(RealT)) + + normal_direction = SVector(one(RealT), one(RealT)) + surface_flux_function = (flux_lax_friedrichs, flux_zero) + orientations = [1, 2] + directions = [1, 2] + flux_lmars = FluxLMARS(RealT(340.0)) + + for direction in directions, orientation in orientations + @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation, + direction, + x, t, + surface_flux_function, + equations)) == + SVector{5, RealT} + end + for orientation in orientations + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux(u_ll, orientation, + equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_lmars(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_waruszewski_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred flux_lmars(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_nonconservative_waruszewski_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred boundary_condition_slip_wall(u, normal_direction, + x, t, + surface_flux_function, + equations)) == SVector{5, RealT} + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction, + equations)) == RealT + end +end + @timed_testset "Compressible Euler Potential Temperature 3D" begin for RealT in (Float32, Float64) equations = @inferred CompressibleEulerPotentialTemperatureEquations3D(c_p = RealT(1004),