diff --git a/examples/elixir_euler_gravity_held_suarez.jl b/examples/elixir_euler_gravity_held_suarez.jl new file mode 100644 index 00000000..67fb5390 --- /dev/null +++ b/examples/elixir_euler_gravity_held_suarez.jl @@ -0,0 +1,185 @@ +# Held-Suarez test case +# Following Souza et al.: +# The Flux-Differencing Discontinuous Galerkin Method Applied to an Idealized Fully +# Compressible Nonhydrostatic Dry Atmosphere + +using OrdinaryDiffEq +using Trixi, TrixiAtmo +using LinearAlgebra + + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquationsWithGravity3D(gamma) + +function initial_condition_isothermal(x, t, equations::CompressibleEulerEquationsWithGravity3D) + # equation (60) in the paper + temperature = 285 # T_I + gas_constant = 287 # R_d + surface_pressure = 1e5 # p_0 + radius_earth = 6.371229e6 # r_planet + gravitational_acceleration = 9.80616 # g + c_v = 717.5 # Specific heat capacity of dry air at constant volume + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + r = z + radius_earth + + # pressure + # geopotential formulation? + p = surface_pressure * + exp(gravitational_acceleration * + (radius_earth^2 / r - radius_earth) / + (gas_constant * temperature)) + + # density (via ideal gas law) + rho = p / (gas_constant * temperature) + + # geopotential + phi = gravitational_acceleration * (radius_earth - radius_earth^2 / r) + + E = c_v * temperature + phi + + return SVector(rho, 0, 0, 0, rho * E, phi) +end + +@inline function source_terms_coriolis(u, x, t, + equations::CompressibleEulerEquationsWithGravity3D) + radius_earth = 6.371229e6 # r_planet + angular_velocity = 7.29212e-5 # Ω + # 7.27220521664e-05 (Giraldo) + + r = norm(x) + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + r = z + radius_earth + + # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4] + du2 = 2 * angular_velocity * u[3] + du3 = -2 * angular_velocity * u[2] + + return SVector(0, du2, du3, 0, 0, 0) +end + +function cartesian_to_sphere(x) + r = norm(x) + lambda = atan(x[2], x[1]) + if lambda < 0 + lambda += 2 * pi + end + phi = asin(x[3] / r) + + return lambda, phi, r +end + +@inline function source_terms_hs_relaxation(u, x, t, + equations::CompressibleEulerEquationsWithGravity3D) + # equations (55)-(58) in the paper + secondsperday = 60*60*24 + radius_earth = 6.371229e6 # r_planet + k_f = 1/secondsperday # Damping scale for momentum + k_a = 1/(40*secondsperday) # Polar relaxation scale + k_s = 1/(4*secondsperday) # Equatorial relaxation scale + T_min = 200 # Minimum equilibrium temperature + T_equator = 315 # Equatorial equilibrium temperature + surface_pressure = 1e5 # p_0 + deltaT = 60 # Latitudinal temperature difference + deltaTheta = 10 # Vertical temperature difference + sigma_b = 0.7 # Dimensionless damping height + gas_constant = 287 # R_d + c_v = 717.5 # Specific heat capacity of dry air at constant volume + c_p = 1004.5 # Specific heat capacity of dry air at constant pressur + + _, _, _, _, pressure = cons2prim(u, equations) + lon, lat, r = cartesian_to_sphere(x) + temperature = pressure / (u[1] * gas_constant) + + sigma = pressure / surface_pressure # "p_0 instead of instantaneous surface pressure" + delta_sigma = max(0, (sigma-sigma_b)/(1-sigma_b)) # "height factor" + k_v = k_f * delta_sigma + k_T = k_a + (k_s - k_a) * delta_sigma * cos(lat)^4 + + T_equi = max(T_min, + (T_equator - deltaT * sin(lat)^2 - deltaTheta * log(sigma) * cos(lat)^2) * + sigma^(gas_constant/c_p)) + + # project onto r, normalize! @. Yₜ.c.uₕ -= k_v * Y.c.uₕ + # Make sure that r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + r = z + radius_earth + dotprod = (u[2] * x[1] + u[3] * x[2] + u[4] * x[3]) / (r*r) + + du2 = -k_v * (u[2] - dotprod * x[1]) + du3 = -k_v * (u[3] - dotprod * x[2]) + du4 = -k_v * (u[4] - dotprod * x[3]) + + du5 = -k_T * u[1] * c_v * (temperature - T_equi) + + return SVector(0, du2, du3, du4, du5, 0) +end + +@inline function source_terms_held_suarez(u, x, t, + equations::CompressibleEulerEquationsWithGravity3D) + return source_terms_coriolis(u,x,t,equations) + + source_terms_hs_relaxation(u,x,t,equations) +end + +initial_condition = initial_condition_isothermal + +boundary_conditions = Dict(:inside => boundary_condition_slip_wall, + :outside => boundary_condition_slip_wall) + +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) + +# Giraldo: (10,8), polydeg 4 +# baroclinic instability: (16, 8), polydeg 5 +lat_lon_trees_per_dim = 8 +layers = 4 +polydeg = 3 + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +mesh = Trixi.T8codeMeshCubedSphere(lat_lon_trees_per_dim, layers, 6.371229e6, 30000.0, + polydeg = polydeg, initial_refinement_level = 0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_held_suarez, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 365 * 24 * 60 * 60) # time in seconds for 1 year +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 5000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim, + output_directory = "out_heldsuarez") + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution) + +############################################################################### +# run the simulation + +# Use a Runge-Kutta method with automatic (error based) time step size control +# Enable threading of the RK method for better performance on multiple threads +sol = solve(ode, + RDPK3SpFSAL49(thread = Trixi.True()); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., maxiters=1e8, + callback = callbacks); diff --git a/examples/elixir_euler_gravity_warmbubble_3d.jl b/examples/elixir_euler_gravity_warmbubble_3d.jl new file mode 100644 index 00000000..9a4c1e63 --- /dev/null +++ b/examples/elixir_euler_gravity_warmbubble_3d.jl @@ -0,0 +1,136 @@ + +using OrdinaryDiffEqSSPRK, 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::CompressibleEulerEquationsWithGravity3D) + 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, 0.0f0, rho * E, phi) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquationsWithGravity3D(warm_bubble_setup.gamma) + +initial_condition = warm_bubble_setup + +volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski) +surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski) + +polydeg = 4 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +trees_per_dimension = (40, 20, 2) +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0, 2_000.0) +mesh = P4estMesh(trees_per_dimension; + polydeg = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 0, + periodicity = (true, false, true)) + +boundary_conditions = Dict(: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/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index b9311b47..0ea88fb7 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -15,7 +15,7 @@ using Printf: @sprintf using Static: True, False using StrideArrays: PtrArray using StaticArrayInterface: static_size -using LinearAlgebra: cross, norm, dot, det +using LinearAlgebra: cross, norm, normalize, dot, det using Reexport: @reexport using LoopVectorization: @turbo using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace @@ -32,7 +32,8 @@ include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D, - SplitCovariantShallowWaterEquations2D, CompressibleEulerEquationsWithGravity2D + SplitCovariantShallowWaterEquations2D, + CompressibleEulerEquationsWithGravity2D, CompressibleEulerEquationsWithGravity3D export GlobalCartesianCoordinates, GlobalSphericalCoordinates diff --git a/src/equations/compressible_euler_gravity_3d.jl b/src/equations/compressible_euler_gravity_3d.jl new file mode 100644 index 00000000..e14dae79 --- /dev/null +++ b/src/equations/compressible_euler_gravity_3d.jl @@ -0,0 +1,987 @@ +# 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""" + CompressibleEulerEquationsWithGravity3D(gamma) + +The compressible Euler equations +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho e \\ Phi +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ \rho v_1 v_3 \\ ( \rho e +p) v_1 \\ 0 +\end{pmatrix} ++ +\frac{\partial}{\partial y} +\begin{pmatrix} +\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ \rho v_1 v_3 \\ ( \rho e +p) v_2 \\ 0 +\end{pmatrix} ++ +\frac{\partial}{\partial z} +\begin{pmatrix} +\rho v_3 \\ \rho v_1 v_3 \\ \rho v_2 v_3 \\ \rho v_3^2 + p \\ ( \rho e +p) v_3 \\ 0 +\end{pmatrix} += +\begin{pmatrix} +0 \\ - \rho \partial \Phi / \partial x \\ - \rho \partial \Phi / \partial y \\ - \rho \partial \Phi / \partial z \\ 0 \\ 0 +\end{pmatrix} +``` +for an ideal gas with ratio of specific heat `gamma` in three space dimensions. + +Here, ``\rho`` is the density, ``v_1``, ``v_2``, ``v_3`` the velocities, ``e`` the specific total energy, +which includes the internal, kinetik, and geopotential energy, `\Phi` is the gravitational +geopotential, and +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2 + v_2^2 + v_3^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 CompressibleEulerEquationsWithGravity3D{RealT <: Real} <: + Trixi.AbstractCompressibleEulerEquations{3, 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 + + function CompressibleEulerEquationsWithGravity3D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +Trixi.have_nonconservative_terms(::CompressibleEulerEquationsWithGravity3D) = True() +function Trixi.varnames(::typeof(cons2cons), ::CompressibleEulerEquationsWithGravity3D) + ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "phi") +end +function Trixi.varnames(::typeof(cons2prim), ::CompressibleEulerEquationsWithGravity3D) + ("rho", "v1", "v2", "v3", "p", "phi") +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::CompressibleEulerEquationsWithGravity3D) + +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) +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity3D) + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # Some vector that can't be identical to normal_vector (unless normal_vector == 0) + tangent1 = SVector(normal_direction[2], normal_direction[3], -normal_direction[1]) + # Orthogonal projection + tangent1 -= dot(normal, tangent1) * normal + tangent1 = normalize(tangent1) + + # Third orthogonal vector + tangent2 = normalize(cross(normal_direction, tangent1)) + + # rotate the internal solution state + u_local = Trixi.rotate_to_x(u_inner, normal, tangent1, tangent2, equations) + + # compute the primitive variables + rho_local, v_normal, v_tangent1, v_tangent2, 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], + p_star * normal[3], + 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)), + zero(eltype(u_inner))) * norm_) +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquationsWithGravity3D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquationsWithGravity3D) + # get the appropriate normal vector from the orientation + RealT = eltype(u_inner) + if orientation == 1 + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) + elseif orientation == 2 + normal_direction = SVector(zero(RealT), one(RealT), zero(RealT)) + else # orientation == 3 + normal_direction = SVector(zero(RealT), zero(RealT), one(RealT)) + 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::CompressibleEulerEquationsWithGravity3D) + +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::CompressibleEulerEquationsWithGravity3D) + # 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, -y, and -z sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, + equations) + 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 +@inline function Trixi.flux(u, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_e, phi = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + p = (equations.gamma - 1) * + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - rho * phi) + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = rho_v1 * v3 + f5 = (rho_e + p) * v1 + elseif orientation == 2 + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = rho_v2 * v3 + f5 = (rho_e + p) * v2 + else + f1 = rho_v3 + f2 = rho_v3 * v1 + f3 = rho_v3 * v2 + f4 = rho_v3 * v3 + p + f5 = (rho_e + p) * v3 + end + return SVector(f1, f2, f3, f4, f5, 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::CompressibleEulerEquationsWithGravity3D) + rho_e = 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 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + f4 = rho_v_normal * v3 + p * normal_direction[3] + f5 = (rho_e + p) * v_normal + return SVector(f1, f2, f3, f4, f5, zero(eltype(u))) +end + +""" + flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity3D) + +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::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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) + p_avg = 0.5f0 * (p_ll + p_rr) + kin_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_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 = f1 * v3_avg + f5 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg + elseif orientation == 2 + 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 = f1 * v3_avg + f5 = p_avg * v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg + else + pv3_avg = 0.5f0 * (p_ll * v3_rr + p_rr * v3_ll) + f1 = rho_avg * v3_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * v3_avg + p_avg + f5 = p_avg * v3_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv3_avg + end + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + # 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 = 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 + v3_ll * v3_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 * v3_avg + p_avg * normal_direction[3] + f5 = (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)) + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +""" + flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity3D) + +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::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + rho_e_ll = u_ll[5] + rho_e_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.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (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.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_avg * v_dot_n_avg + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = f1 * v3_avg + p_avg * normal_direction[3] + f5 = f1 * e_avg + p_avg * v_dot_n_avg + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +""" + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsWithGravity3D) + +Entropy conserving two-point flux by +- Chandrashekar (2013) + Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes + for Compressible Euler and Navier-Stokes Equations + [DOI: 10.4208/cicp.170712.010313a](https://doi.org/10.4208/cicp.170712.010313a) +""" +@inline function Trixi.flux_chandrashekar(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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) + + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2 + v3_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2 + v3_rr^2) + + # Compute the necessary mean values + rho_avg = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + velocity_square_avg = specific_kin_ll + specific_kin_rr + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + p_mean + f3 = f1 * v2_avg + f4 = f1 * v3_avg + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + elseif orientation == 2 + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_mean + f4 = f1 * v3_avg + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + else + f1 = rho_mean * v3_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * v3_avg + p_mean + f5 = f1 * 0.5f0 * + (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + end + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + specific_kin_ll = 0.5f0 * (v1_ll^2 + v2_ll^2 + v3_ll^2) + specific_kin_rr = 0.5f0 * (v1_rr^2 + v2_rr^2 + v3_rr^2) + + # Compute the necessary mean values + rho_avg = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + velocity_square_avg = specific_kin_ll + specific_kin_rr + + # Multiply with average of normal velocities + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_mean * normal_direction[1] + f3 = f1 * v2_avg + p_mean * normal_direction[2] + f4 = f1 * v3_avg + p_mean * normal_direction[3] + f5 = f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +""" + flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity3D) + +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::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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) + + # 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) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_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 * v3_avg + f5 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5f0 * (p_ll * v1_rr + p_rr * v1_ll) + elseif orientation == 2 + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_avg + f4 = f1 * v3_avg + f5 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5f0 * (p_ll * v2_rr + p_rr * v2_ll) + else # orientation == 3 + f1 = rho_mean * v3_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + f4 = f1 * v3_avg + p_avg + f5 = f1 * + (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5f0 * (p_ll * v3_rr + p_rr * v3_ll) + end + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity3D) + # Unpack left and right state + 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_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + # 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) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_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 * v3_avg + p_avg * normal_direction[3] + f5 = (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)) + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +function flux_nonconservative_waruszewski(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity3D) + rho_ll, _, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, _, phi_rr = u_rr + + 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], + noncons * normal_direction[3], + f0, f0) +end + +function flux_nonconservative_waruszewski(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity3D) + rho_ll, _, _, _, _, phi_ll = u_ll + rho_rr, _, _, _, _, phi_rr = u_rr + + 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, f0) + elseif orientation == 2 + return SVector(f0, f0, noncons, f0, f0, f0) + else #if orientation == 3 + return SVector(f0, f0, f0, noncons, f0, f0) + end +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity3D) + +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) +""" +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity3D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + 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) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + + 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, f5 = v * u_ll + f5 = f5 + p_ll * v + else + f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v + end + + if orientation == 1 + f2 += p + elseif orientation == 2 + f3 += p + else # orientation == 3 + f4 += p + end + + return SVector(f1, f2, f3, f4, f5, zero(eltype(u_ll))) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquationsWithGravity3D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + 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_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + # 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, f5, _ = u_ll * v + f5 = f5 + p_ll * v + else + f1, f2, f3, f4, f5, _ = u_rr * v + f5 = f5 + p_rr * v + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4 + p * normal_direction[3], + f5, 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::CompressibleEulerEquationsWithGravity3D) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll, _ = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr, _ = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + elseif orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + else # orientation == 3 + v_ll = v3_ll + v_rr = v3_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + λ_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::CompressibleEulerEquationsWithGravity3D) + 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 + +# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsWithGravity3D) + 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) + + 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) + elseif orientation == 2 # y-direction + λ_min = v2_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v2_rr + sqrt(equations.gamma * p_rr / rho_rr) + else # z-direction + λ_min = v3_ll - sqrt(equations.gamma * p_ll / rho_ll) + λ_max = v3_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::CompressibleEulerEquationsWithGravity3D) + 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 + +# Rotate normal vector to x-axis; normal, tangent1 and tangent2 need to be orthonormal +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the directions +# has been normalized prior to this rotation of the state vector +@inline function Trixi.rotate_to_x(u, normal_vector, tangent1, tangent2, + equations::CompressibleEulerEquationsWithGravity3D) + # Multiply with [ 1 0 0 0 0; + # 0 ― normal_vector ― 0; + # 0 ― tangent1 ― 0; + # 0 ― tangent2 ― 0; + # 0 0 0 0 1 ] + return SVector(u[1], + normal_vector[1] * u[2] + normal_vector[2] * u[3] + + normal_vector[3] * u[4], + tangent1[1] * u[2] + tangent1[2] * u[3] + tangent1[3] * u[4], + tangent2[1] * u[2] + tangent2[2] * u[3] + tangent2[3] * u[4], + u[5], u[6]) +end + +# Rotate x-axis to normal vector; normal, tangent1 and tangent2 need to be orthonormal +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the directions +# has been normalized prior to this back-rotation of the state vector +@inline function Trixi.rotate_from_x(u, normal_vector, tangent1, tangent2, + equations::CompressibleEulerEquationsWithGravity3D) + # Multiply with [ 1 0 0 0 0; + # 0 | | | 0; + # 0 normal_vector tangent1 tangent2 0; + # 0 | | | 0; + # 0 0 0 0 1 ] + return SVector(u[1], + normal_vector[1] * u[2] + tangent1[1] * u[3] + tangent2[1] * u[4], + normal_vector[2] * u[2] + tangent1[2] * u[3] + tangent2[2] * u[4], + normal_vector[3] * u[2] + tangent1[3] * u[3] + tangent2[3] * u[4], + u[5], u[6]) +end + +@inline function Trixi.max_abs_speeds(u, + equations::CompressibleEulerEquationsWithGravity3D) + 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::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_e, phi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + p = (equations.gamma - 1) * + (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - rho * phi) + + return SVector(rho, v1, v2, v3, p, phi) +end + +# Convert conservative variables to entropy +@inline function Trixi.cons2entropy(u, + equations::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_e, phi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_square = v1^2 + v2^2 + v3^2 + p = (equations.gamma - 1) * (rho_e - 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 * v3 + w5 = -rho_p + + return SVector(w1, w2, w3, w4, w5, phi) +end + +@inline function Trixi.entropy2cons(w, + equations::CompressibleEulerEquationsWithGravity3D) + # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD + # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + @unpack gamma = equations + + # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986) + # instead of `-rho * s / (gamma - 1)` + V1, V2, V3, V4, V5, _ = w .* (gamma - 1) + phi = w[6] + + # s = specific entropy, eq. (53) + V_square = V2^2 + V3^2 + V4^2 + s = gamma - V1 + V_square / (2 * V5) - V5 * phi + + # eq. (52) + rho_iota = ((gamma - 1) / (-V5)^gamma)^(equations.inv_gamma_minus_one) * + exp(-s * equations.inv_gamma_minus_one) + + # eq. (51) + rho = -rho_iota * V5 + rho_v1 = rho_iota * V2 + rho_v2 = rho_iota * V3 + rho_v3 = rho_iota * V4 + rho_e = rho_iota * (1 - V_square / (2 * V5)) + rho * phi + return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, phi) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, + equations::CompressibleEulerEquationsWithGravity3D) + rho, v1, v2, v3, p, phi = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_v3 = rho * v3 + rho_e = 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_e, phi) +end + +@inline function Trixi.density(u, equations::CompressibleEulerEquationsWithGravity3D) + rho = u[1] + return rho +end + +@inline function Trixi.velocity(u, equations::CompressibleEulerEquationsWithGravity3D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +@inline function Trixi.velocity(u, orientation::Int, + equations::CompressibleEulerEquationsWithGravity3D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + +@inline function Trixi.pressure(u, equations::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + p = (equations.gamma - 1) * + (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - rho * phi) + return p +end + +@inline function Trixi.density_pressure(u, + equations::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + rho_times_p = (equations.gamma - 1) * + (rho * rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) - rho^2 * phi) + return rho_times_p +end + +# Calculate thermodynamic entropy for a conservative state `u` +@inline function Trixi.entropy_thermodynamic(u, + equations::CompressibleEulerEquationsWithGravity3D) + rho, _ = u + p = pressure(u, equations) + + # Thermodynamic entropy + s = log(p) - equations.gamma * log(rho) + + return s +end + +# Calculate mathematical entropy for a conservative state `cons` +@inline function Trixi.entropy_math(cons, + equations::CompressibleEulerEquationsWithGravity3D) + S = -entropy_thermodynamic(cons, equations) * cons[1] * + equations.inv_gamma_minus_one + # Mathematical entropy + + return S +end + +# Default entropy is the mathematical entropy +@inline function Trixi.entropy(cons, equations::CompressibleEulerEquationsWithGravity3D) + entropy_math(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline Trixi.energy_total(cons, ::CompressibleEulerEquationsWithGravity3D) = cons[5] + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, + equations::CompressibleEulerEquationsWithGravity3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + return 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho +end + +# Calculate internal energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, + equations::CompressibleEulerEquationsWithGravity3D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) - + cons[1] * cons[6] +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::CompressibleEulerEquationsWithGravity3D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], diss[4], diss[5], zero(eltype(u_ll))) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 7fabd517..7ebe771d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -327,4 +327,5 @@ include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") include("reference_data.jl") include("compressible_euler_gravity_2d.jl") +include("compressible_euler_gravity_3d.jl") end # @muladd diff --git a/test/runtests.jl b/test/runtests.jl index 2f1ef05e..0fe3a2a6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,7 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "euler_gravity" include("test_2d_euler_gravity.jl") + include("test_3d_euler_gravity.jl") end @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" diff --git a/test/test_3d_euler_gravity.jl b/test/test_3d_euler_gravity.jl new file mode 100644 index 00000000..12ffddeb --- /dev/null +++ b/test/test_3d_euler_gravity.jl @@ -0,0 +1,26 @@ +module TestEulerGravity3D + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "Held-Suarez" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_gravity_held_suarez.jl"), + l2=[0.0, 0.0, 0.0, 0.0, 0.0], + linf=[0.0, 0.0, 0.0, 0.0, 0.0], + tspan=(0.0, 1.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 + +end # module