diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 0ec083a..1491ee8 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: FluxLMARS, source_terms_gravity, flux_theta function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl index 2570783..c3faf15 100644 --- a/examples/elixir_euler_potential_temperature_ec.jl +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_theta + ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerPotentialTemperatureEquations2D() 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..64900ef --- /dev/null +++ b/examples/elixir_euler_potential_temperature_gravity_wave.jl @@ -0,0 +1,102 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +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 1256a9c..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: FluxLMARS, source_terms_geopotential, cons2drypot +using TrixiAtmo: source_terms_geopotential, cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equations diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 956c4a2..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, - FluxLMARS +using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble using NLsolve: nlsolve ############################################################################### diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index dd8de4f..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, FluxLMARS + cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equation diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e4f64c3..fcfb4c0 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,8 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, FluxLMARS, flux_theta +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 index 1b9aee8..926d7bc 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -225,6 +225,88 @@ end 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 @@ -267,7 +349,7 @@ end # Mathematical entropy p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma - U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1])) + U = (p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1])) return U end diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl index 737e6e5..863b3e7 100644 --- a/test/test_2d_euler_potential_temperature.jl +++ b/test/test_2d_euler_potential_temperature.jl @@ -34,6 +34,33 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 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"),