diff --git a/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl b/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl new file mode 100644 index 00000000000..95f71f38130 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulermulti_convergence_ec.jl @@ -0,0 +1,55 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler multicomponent equations +equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, 0.4)) + +initial_condition = initial_condition_convergence_test + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +cells_per_dimension = (16, 16) +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# 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 = 0.5) + +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/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 7b437f4a1b4..ecd3bc80c0a 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -270,6 +270,29 @@ end return vcat(f_other, f_rho) end +# Calculate 1D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + gamma = totalgamma(u, equations) + p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + + f_rho = densities(u, v_normal, equations) + f1 = rho_v1 * v_normal + p * normal_direction[1] + f2 = rho_v2 * v_normal + p * normal_direction[2] + f3 = (rho_e + p) * v_normal + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + """ flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) @@ -446,6 +469,76 @@ See also return vcat(f_other, f_rho) end +@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, cv = equations + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + 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] + + # helpful variables + help1_ll = zero(v1_ll) + help1_rr = zero(v1_rr) + enth_ll = zero(v1_ll) + enth_rr = zero(v1_rr) + for i in eachcomponent(equations) + enth_ll += u_ll[i + 3] * gas_constants[i] + enth_rr += u_rr[i + 3] * gas_constants[i] + help1_ll += u_ll[i + 3] * cv[i] + help1_rr += u_rr[i + 3] * cv[i] + end + + # temperature and pressure + T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + (v_dot_n_ll + v_dot_n_rr) + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] + f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) @@ -491,6 +584,50 @@ end return (abs(v1) + c, abs(v2) + c) end +@inline function rotate_to_x(u, normal_vector, + equations::CompressibleEulerMulticomponentEquations2D) + # 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 + # [ n_1 n_2 0 0; + # t_1 t_2 0 0; + # 0 0 1 0 + # 0 0 0 1] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] + s * u[2], + -s * u[1] + c * u[2], + u[3], + densities...) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function rotate_from_x(u, normal_vector, + equations::CompressibleEulerMulticomponentEquations2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ n_1 t_1 0 0; + # n_2 t_2 0 0; + # 0 0 1 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] - s * u[2], + s * u[1] + c * u[2], + u[3], + densities...) +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 71782644b17..43be04f745d 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -61,6 +61,23 @@ struct FluxRotated{NumericalFlux} numerical_flux::NumericalFlux end +# Rotated surface flux computation (2D version) +@inline function (flux_rotated::FluxRotated)(u, + normal_direction::AbstractVector, + equations::AbstractEquations{2}) + @unpack numerical_flux = flux_rotated + + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal_vector = normal_direction / norm_ + + u_rotated = rotate_to_x(u, normal_vector, equations) + + f = numerical_flux(u_rotated, 1, equations) + + return rotate_from_x(f, normal_vector, equations) * norm_ +end + # Rotated surface flux computation (2D version) @inline function (flux_rotated::FluxRotated)(u_ll, u_rr, normal_direction::AbstractVector, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 96202e00f58..1addc29e3e6 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -247,6 +247,32 @@ end end end +@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 1.5123651627525257e-5, + 1.51236516273878e-5, + 2.4544918394022538e-5, + 5.904791661362391e-6, + 1.1809583322724782e-5, + ], + linf=[ + 8.393471747591974e-5, + 8.393471748258108e-5, + 0.00015028562494778797, + 3.504466610437795e-5, + 7.00893322087559e-5, + ]) + # 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 Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), # Expected errors are exactly the same as with TreeMesh! diff --git a/test/test_unit.jl b/test/test_unit.jl index 7f774e11d3b..d2e744da62f 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -296,9 +296,11 @@ end end Trixi.move_connectivity!(c::MyContainer, first, last, destination) = c Trixi.delete_connectivity!(c::MyContainer, first, last) = c - Trixi.reset_data_structures!(c::MyContainer) = (c.data = Vector{Int}(undef, - c.capacity + 1); - c) + function Trixi.reset_data_structures!(c::MyContainer) + (c.data = Vector{Int}(undef, + c.capacity + 1); + c) + end function Base.:(==)(c1::MyContainer, c2::MyContainer) return (c1.capacity == c2.capacity && c1.length == c2.length && @@ -611,6 +613,18 @@ end @test_throws ArgumentError TimeSeriesCallback(semi, [1.0 1.0 1.0; 2.0 2.0 2.0]) end +@timed_testset "Consistency check for single point flux: CEMCE" begin + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, 0.4)) + u = SVector(0.1, -0.5, 1.0, 1.0, 2.0) + + orientations = [1, 2] + for orientation in orientations + @test flux(u, orientation, equations) ≈ + flux_ranocha(u, u, orientation, equations) + end +end + @timed_testset "Consistency check for HLL flux (naive): CEE" begin flux_hll = FluxHLL(min_max_speed_naive) @@ -1221,6 +1235,29 @@ end end @testset "FluxRotated vs. direct implementation" begin + @timed_testset "CompressibleEulerMulticomponentEquations2D" begin + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + gas_constants = (0.4, + 0.4)) + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0), + SVector(-0.1, -0.3, 1.2, 1.3, 1.4)] + + f_std = flux + f_rot = FluxRotated(f_std) + println(typeof(f_std)) + println(typeof(f_rot)) + for u in u_values, + normal_direction in normal_directions + + @test f_rot(u, normal_direction, equations) ≈ + f_std(u, normal_direction, equations) + end + end + @timed_testset "CompressibleEulerEquations2D" begin equations = CompressibleEulerEquations2D(1.4) normal_directions = [SVector(1.0, 0.0),