diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl new file mode 100644 index 0000000..33216b2 --- /dev/null +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -0,0 +1,112 @@ +using OrdinaryDiffEq +using Trixi, TrixiAtmo +using TrixiAtmo: source_terms_gravity + +function initial_condition_warm_bubble(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + g = equations.g + c_p = equations.c_p + c_v = equations.c_v + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # 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.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * 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_ref) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + T = potential_temperature * exner + + # density + rho = p / (R * T) + v1 = 20.0 + v2 = 0.0 + + return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerPotentialTemperatureEquations2D() + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_theta +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +cells_per_dimension = (64, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) +initial_condition = initial_condition_warm_bubble + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_gravity, + 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 = 1000 + +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 = 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), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl new file mode 100644 index 0000000..c3faf15 --- /dev/null +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -0,0 +1,77 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerPotentialTemperatureEquations2D() + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = SVector(0, 0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) + + return TrixiAtmo.prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +volume_flux = flux_theta +solver = DGSEM(polydeg = 3, surface_flux = flux_theta, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (32, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, true)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# 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) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + 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); +summary_callback() # print the timer summary diff --git a/examples/elixir_euler_potential_temperature_gravity_wave.jl b/examples/elixir_euler_potential_temperature_gravity_wave.jl new file mode 100644 index 0000000..880e41c --- /dev/null +++ b/examples/elixir_euler_potential_temperature_gravity_wave.jl @@ -0,0 +1,102 @@ +using OrdinaryDiffEq +using Trixi, TrixiAtmo +using TrixiAtmo: source_terms_gravity + +function initial_condition_gravity_wave(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + g = equations.g + c_p = equations.c_p + c_v = equations.c_v + # center of perturbation + x_c = 100_000.0 + a = 5_000 + H = 10_000 + # perturbation in potential temperature + potential_temperature_ref = 300.0 * exp(0.01^2 / g * x[2]) + potential_temperature_perturbation = 0.01 * sinpi(x[2] / H) / (1 + (x[1] - x_c)^2 / a^2) + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 + g^2 / (c_p * 300.0 * 0.01^2) * (exp(-0.01^2 / g * x[2]) - 1) + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + T = potential_temperature * exner + + # density + rho = p / (R * T) + v1 = 20.0 + v2 = 0.0 + + return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerPotentialTemperatureEquations2D() + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(340.0) + +solver = DGSEM(basis, surface_flux) + +coordinates_min = (0.0, 0.0) +coordinates_max = (300_000.0, 10_000.0) + +cells_per_dimension = (600, 20) # Delta x = Delta z = 1 km +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) +initial_condition = initial_condition_gravity_wave + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_gravity, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3000.0) # 1000 seconds final time + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +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 = 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), + maxiters = 1.0e7, + dt = 1e-1, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index d34c63d..8571896 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_LMARS, source_terms_geopotential, cons2drypot +using TrixiAtmo: source_terms_geopotential, cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equations @@ -61,7 +61,7 @@ source_term = source_terms_geopotential polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index bdf7b4e..d406bab 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo -using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble, - flux_LMARS +using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble using NLsolve: nlsolve ############################################################################### @@ -242,7 +241,7 @@ source_term = source_terms_moist_bubble polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index 02ab73e..4646a8b 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -3,7 +3,7 @@ using Trixi, TrixiAtmo using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, source_terms_phase_change, source_terms_nonhydrostatic_rayleigh_sponge, - cons2drypot, flux_LMARS + cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equation @@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e6c4608..fcfb4c0 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -24,8 +24,10 @@ include("solvers/solvers.jl") include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D +export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, flux_LMARS +export flux_chandrashekar, flux_theta, flux_theta_rhoAM, flux_theta_physentropy, + flux_theta_physentropy2 export examples_dir diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl new file mode 100644 index 0000000..926d7bc --- /dev/null +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -0,0 +1,380 @@ +using Trixi +using Trixi: ln_mean, stolarsky_mean +import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy, FluxLMARS + +@muladd begin +#! format: noindent +struct CompressibleEulerPotentialTemperatureEquations2D{RealT <: Real} <: + AbstractCompressibleEulerPotentialTemperatureEquations{2, 4} + p_0::RealT + c_p::RealT + c_v::RealT + g::RealT + R::RealT + gamma::RealT +end + +function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Float64) + p_0 = 100_000.0 + c_p = 1004.0 + c_v = 717.0 + R = c_p - c_v + gamma = c_p / c_v + return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R, + gamma) +end + +function varnames(::typeof(cons2cons), + ::CompressibleEulerPotentialTemperatureEquations2D) + ("rho", "rho_v1", "rho_v2", "rho_theta") +end + +varnames(::typeof(cons2prim), ::CompressibleEulerPotentialTemperatureEquations2D) = ("rho", + "v1", + "v2", + "p1") + +# Calculate 1D flux for a single point in the normal direction. +# Note, this directional vector is not normalized. +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = (rho_v_normal) * v1 + p * normal_direction[1] + f3 = (rho_v_normal) * v2 + p * normal_direction[2] + f4 = (rho_theta) * v_normal + return SVector(f1, f2, f3, f4) +end + +@inline function flux(u, orientation::Integer, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = (rho_theta) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = (rho_theta) * v2 + end + + return SVector(f1, f2, f3, f4) +end + +# Slip-wall boundary condition +# Determine the boundary numerical surface flux for a slip wall condition. +# Imposes a zero normal velocity at the wall. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleEulerPotentialTemperatureEquations2D) + @unpack gamma = equations + 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 = 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.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2.0 / ((gamma + 1) * rho_local) + B = p_local * (gamma - 1) / (gamma + 1) + p_star = p_local + + 0.5f0 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4.0 * 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))) * norm_ +end + +# Fix sign for structured mesh. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # 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) + 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 + +@inline function source_terms_gravity(u, x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, _, _, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, + zero(eltype(u))) +end + +# Rotate momentum flux. The same as in compressible Euler. +@inline function rotate_to_x(u, normal_vector::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # 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]) +end + +# Low Mach number approximate Riemann solver (LMARS) from +# X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S. +# Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian +# Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, +# https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + a = 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] + + norm_ = norm(normal_direction) + + rho = 0.5f0 * (rho_ll + rho_rr) + + p_interface = 0.5f0 * (p_ll + p_rr) - 0.5f0 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f1, f2, f3, f4 = u_ll * v_interface + else + f1, f2, f3, f4 = u_rr * v_interface + end + + return SVector(f1, + f2 + p_interface * normal_direction[1], + f3 + p_interface * normal_direction[2], + f4) +end + +## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation + +@inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + _, _, _, rho_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_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 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + return SVector(f1, f2, f3, f4) +end + +@inline function flux_theta_rhoAM(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + _, _, _, rho_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = 0.5f0 * (rho_ll + rho_rr) + + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_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 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + return SVector(f1, f2, f3, f4) +end + +@inline function flux_theta_physentropy(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + _, _, _, rho_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + #rho_mean = ln_mean(rho_ll, rho_rr) + #rho_mean = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + #f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f1 = f4 * rho_mean + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + + return SVector(f1, f2, f3, f4) +end + +@inline function flux_theta_physentropy2(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + _, _, _, rho_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + #rho_mean = 0.5f0 * (rho_ll + rho_rr) + #rho_mean = ln_mean(rho_ll/rho_theta_ll, rho_rr/rho_theta_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_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 = inv_ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr) * f1 + return SVector(f1, f2, f3, f4) +end + +@inline function prim2cons(prim, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, v1, v2, p = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_theta = (p / equations.p_0)^(1 / equations.gamma) * equations.p_0 / equations.R + return SVector(rho, rho_v1, rho_v2, rho_theta) +end + +@inline function cons2prim(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + return SVector(rho, v1, v2, p) +end + +@inline function cons2cons(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) + return u +end + +@inline function cons2entropy(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + + k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma + w1 = -0.5f0 * rho_v1^2 / (rho)^2 - 0.5f0 * rho_v2^2 / (rho)^2 + w2 = rho_v1 / rho + w3 = rho_v2 / rho + w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1) + + return SVector(w1, w2, w3, w4) +end + +@inline function entropy_math(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # Mathematical entropy + p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma + + U = (p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1])) + + return U +end + +# Default entropy is the mathematical entropy +@inline function entropy(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) + entropy_math(cons, equations) +end + +@inline function energy_total(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) + entropy(cons, equations) +end + +@inline function energy_kinetic(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) + return 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1]) +end + +@inline function max_abs_speeds(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) + rho, v1, v2, p = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c +end +end # @muladd diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 638cc39..53557af 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -5,7 +5,7 @@ using Trixi using Trixi: ln_mean, inv_ln_mean import Trixi: varnames, flux_chandrashekar, boundary_condition_slip_wall, cons2prim, cons2entropy, max_abs_speed_naive, max_abs_speeds, - entropy, energy_total, flux + entropy, energy_total, flux, FluxLMARS @muladd begin #! format: noindent @@ -23,7 +23,6 @@ struct CompressibleMoistEulerEquations2D{RealT <: Real} <: kappa::RealT # ratio of the gas constant R_d gamma::RealT # = inv(kappa- 1); can be used to write slow divisions as fast multiplications L_00::RealT # latent heat of evaporation at 0 K - a::RealT end function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) @@ -38,9 +37,8 @@ function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) gamma = c_pd / c_vd # = 1/(1 - kappa) kappa = 1 - inv(gamma) L_00 = 3147620.0 - a = 360.0 return CompressibleMoistEulerEquations2D{RealT}(p_0, c_pd, c_vd, R_d, c_pv, c_vv, - R_v, c_pl, g, kappa, gamma, L_00, a) + R_v, c_pl, g, kappa, gamma, L_00) end # Calculate 1D flux for a single point. @@ -120,13 +118,13 @@ end if v_normal <= 0.0 sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * - inv(gamma - 1)) + (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + inv(gamma - 1)) else # v_normal > 0.0 A = 2.0 / ((gamma + 1) * rho_local) B = p_local * (gamma - 1) / (gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B))) end @@ -370,11 +368,11 @@ end # positive even power with default value 2 gamma = 2.0 # relaxation coefficient > 0 - alpha = 0.5 + alpha = 0.5f0 tau_s = zero(eltype(u)) if z > z_s - tau_s = alpha * sin(0.5 * (z - z_s) * inv(z_top - z_s))^(gamma) + tau_s = alpha * sin(0.5f0 * (z - z_s) * inv(z_top - z_s))^(gamma) end return SVector(zero(eltype(u)), @@ -400,7 +398,7 @@ end v2 = rho_v2 / rho # Inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -452,9 +450,10 @@ end # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleMoistEulerEquations2D) - @unpack a = equations + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, rho_qv_rr, rho_ql_rr = u_rr @@ -474,9 +473,10 @@ end # Compute the necessary interface flux components norm_ = norm(normal_direction) - rho = 0.5 * (rho_ll + rho_rr) - p_interface = 0.5 * (p_ll + p_rr) - beta * 0.5 * a * rho * (v_rr - v_ll) / norm_ - v_interface = 0.5 * (v_ll + v_rr) - beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + rho = 0.5f0 * (rho_ll + rho_rr) + p_interface = 0.5f0 * (p_ll + p_rr) - beta * 0.5f0 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - + beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ if (v_interface > 0) f1, f2, f3, f4, f5, f6 = u_ll * v_interface @@ -549,7 +549,7 @@ end g_v = L_00 + (c_pv - s_v) * T g_l = (c_pl - s_l) * T - w1 = g_d - 0.5 * v_square + w1 = g_d - 0.5f0 * v_square w2 = v1 w3 = v2 w4 = -1 @@ -566,7 +566,7 @@ end rho_v2 = rho * v2 rho_qv = rho * qv rho_ql = rho * ql - rho_E = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + rho_E = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end @@ -670,7 +670,7 @@ end v2 = rho_v2 / rho # inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute Temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -716,7 +716,7 @@ end v2 = rho_v2 / rho # inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute Temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -950,8 +950,8 @@ end v2_rr = rho_v2_rr / rho_rr # inner energy - rho_e_ll = (rho_E_ll - 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll)) - rho_e_rr = (rho_E_rr - 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr)) + rho_e_ll = (rho_E_ll - 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll)) + rho_e_rr = (rho_E_rr - 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr)) # Absolute Temperature T_ll = (rho_e_ll - L_00 * rho_qv_ll) / @@ -977,18 +977,18 @@ end inv_T_mean = inv_ln_mean(inv(T_ll), inv(T_rr)) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v1_square_avg = 0.5 * (v1_ll^2 + v1_rr^2) - v2_square_avg = 0.5 * (v2_ll^2 + v2_rr^2) - rho_qd_avg = 0.5 * (rho_qd_ll + rho_qd_rr) - rho_qv_avg = 0.5 * (rho_qv_ll + rho_qv_rr) - rho_ql_avg = 0.5 * (rho_ql_ll + rho_ql_rr) - inv_T_avg = 0.5 * (inv(T_ll) + inv(T_rr)) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v1_square_avg = 0.5f0 * (v1_ll^2 + v1_rr^2) + v2_square_avg = 0.5f0 * (v2_ll^2 + v2_rr^2) + rho_qd_avg = 0.5f0 * (rho_qd_ll + rho_qd_rr) + rho_qv_avg = 0.5f0 * (rho_qv_ll + rho_qv_rr) + rho_ql_avg = 0.5f0 * (rho_ql_ll + rho_ql_rr) + inv_T_avg = 0.5f0 * (inv(T_ll) + inv(T_rr)) v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg p_int = inv(inv_T_avg) * (R_d * rho_qd_avg + R_v * rho_qv_avg + R_q * rho_ql_avg) - K_avg = 0.5 * (v1_square_avg + v2_square_avg) + K_avg = 0.5f0 * (v1_square_avg + v2_square_avg) f_1d = rho_qd_mean * v_dot_n_avg f_1v = rho_qv_mean * v_dot_n_avg diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a833960..6fe1256 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -2,4 +2,8 @@ using Trixi: AbstractEquations abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end + include("compressible_moist_euler_2d_lucas.jl") +include("compressible_euler_potential_temperature_2d.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 541cfd8..9635f2e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,4 +20,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" include("test_spherical_advection.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "euler_potential_temperature" + include("test_2d_euler_potential_temperature.jl") + end end diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl new file mode 100644 index 0000000..863b3e7 --- /dev/null +++ b/test/test_2d_euler_potential_temperature.jl @@ -0,0 +1,89 @@ +module TestExamples2DEulerPotentialTemperature + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") # TODO - This is a repetition from Trixi.jl + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory for examples? + +@trixiatmo_testset "elixir_euler_potential_temperature_dry_bubble" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_dry_bubble.jl"), + l2=[ + 1.300427456987984e-6, + 2.6010873806861595e-5, + 0.0006660120005093007, + 9.51074191163579e-6, + ], + linf=[ + 1.031131432183141e-5, + 0.00020623855042956052, + 0.006392070867001616, + 7.841424786647622e-5, + ], + polydeg=3, + tspan=(0.0, 0.1)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixiatmo_testset "elixir_euler_potential_temperature_gravity_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_gravity_wave.jl"), + l2=[ + 8.92434879991671e-9, + 1.8131676850378482e-7, + 2.6374650049135436e-5, + 1.4620631924879524e-6, + ], + linf=[ + 7.544232794032268e-8, + 1.654911628179434e-6, + 0.00023724858495745153, + 1.8884994915424613e-5, + ], + polydeg=3, + 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 + +@trixiatmo_testset "elixir_euler_potential_temperature_ec" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_ec.jl"), + l2=[ + 0.06174365254988615, + 0.05018008812040643, + 0.05018774429827882, + 0.0057820937341387935, + ], + linf=[ + 0.2932352942992503, + 0.3108954213591686, + 0.31082193857647905, + 0.027465217019157606, + ]) + # 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_trixi_consistency.jl b/test/test_trixi_consistency.jl index b13b5d3..477a3b8 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -49,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true) trixi_include(trixi_elixir, equations = equations_moist, volume_flux = flux_chandrashekar, - surface_flux = flux_LMARS, + surface_flux = Trixi.FluxLMARS(360.0), maxiters = maxiters) errors_atmo = Main.analysis_callback(Main.sol)