Sc/compressible euler multicomponent 2d fluxes (#1681)

* Added coupling converters.

* Added generic converter_function for structured 2d meshes.

* Added fluxes to 2d compressible multicomponent Euler.

* Removed backup file.

* Removed redundant project toml file.

* Removed redundant Project toml file.

* Removed redundant project toml file.

* Removed redundant project toml files.

* Removed rdundant files.

* Removed p4est coupling test elixir.

This is too experimental to be in any branch yet.

* Removed coupling converter functions.

These should be part of a different PR.

* Removed further coupling converter functions.

* Removed redundant swap file.

* Reverted coupling semidiscretization to the main branch's by hand.

This should be part of a different PR.

* Added example elixir for the multicomponent Euler equations on structured meshes.

* Added test for structured 2d Euler multicomponent elixir.

* Reduced resolution of Euler multicomponent elixir to reduce compute time.

* Applied autoformatter on multicomponent Euler.

* Autoformatter on test.

* Update src/equations/compressible_euler_multicomponent_2d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Added checks for type instability for the structured mesh multicomponent Euler equations.

* Autoformatted structured 2d test.

* Added consistency check for new multicomponent Euler fluxes.

* Added rotational invariance test for CompressibleEulerMulticomponentEquations2D.

* Added rotations.

* Added rotated surface for fluxes with only one u-vector.

* Corrected type instability in rotation.

* Applied autoformatter on 2d multicomponent Euler equations.

* Removed redundant speed functions from multicomponent Euler equations.

* Applied autoformatter on unit test.

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

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

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

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,

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!

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),