Fixed cons2entropy files and implemented entropy2cons for CompressibleEulerMulticomponent1D and CompressibleEulerMulticomponent2D

src/equations/compressible_euler_multicomponent_1d.jl

@@ -461,6 +461,7 @@ end
 # Convert conservative variables to entropy
 @inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations1D)
     @unpack cv, gammas, gas_constants = equations
+
     rho_v1, rho_e = u
 
     rho = density(u, equations)
@@ -480,21 +481,86 @@ end
     s = log(p) - gamma * log(rho) - log(gas_constant)
     rho_p = rho / p
     T = (rho_e - 0.5 * rho * v_square) / (help1)
-    entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant *
-                                                                  ((gamma - s) /
-                                                                   (gamma - 1.0) -
-                                                                   (0.5 * v_square *
-                                                                    rho_p))
+
+    entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] *
+                                                                   (1 - log(T)) +
+                                                                   gas_constants[i] *
+                                                                   (1 + log(u[i + 2])) -
+                                                                   v1^2 / (2 * T))
                                                                   for i in eachcomponent(equations))
 
     w1 = gas_constant * v1 * rho_p
-    w2 = gas_constant * (-1.0 * rho_p)
+    w2 = gas_constant * (-rho_p)
 
     entrop_other = SVector{2, real(equations)}(w1, w2)
 
     return vcat(entrop_other, entrop_rho)
 end
 
+# Convert entropy variables to conservative variables
+@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations1D)
+    @unpack gammas, gas_constants, cv, cp = equations
+    T = -1 / w[2]
+    v1 = w[1] * T
+    cons_rho = SVector{ncomponents(equations), real(equations)}(exp(1 /
+                                                                    gas_constants[i] *
+                                                                    (-cv[i] *
+                                                                     log(-w[2]) -
+                                                                     cp[i] + w[i + 2] -
+                                                                     0.5 * w[1]^2 /
+                                                                     w[2]))
+                                                                for i in eachcomponent(equations))
+
+    rho = zero(cons_rho[1])
+    help1 = zero(cons_rho[1])
+    help2 = zero(cons_rho[1])
+    p = zero(cons_rho[1])
+    for i in eachcomponent(equations)
+        rho += cons_rho[i]
+        help1 += cons_rho[i] * cv[i] * gammas[i]
+        help2 += cons_rho[i] * cv[i]
+        p += cons_rho[i] * gas_constants[i] * T
+    end
+    u1 = rho * v1
+    gamma = help1 / help2
+    u2 = p / (gamma - 1) + 0.5 * rho * v1^2
+    cons_other = SVector{2, real(equations)}(u1, u2)
+    return vcat(cons_other, cons_rho)
+end
+
+@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations1D)
+    @unpack cv, gammas, gas_constants = equations
+    rho_v1, rho_e = u
+    rho = density(u, equations)
+    T = temperature(u, equations)
+
+    total_entropy = zero(u[1])
+    for i in eachcomponent(equations)
+        total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2]))
+    end
+
+    return total_entropy
+end
+
+@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations1D)
+    @unpack cv, gammas, gas_constants = equations
+
+    rho_v1, rho_e = u
+
+    rho = density(u, equations)
+    help1 = zero(rho)
+
+    for i in eachcomponent(equations)
+        help1 += u[i + 2] * cv[i]
+    end
+
+    v1 = rho_v1 / rho
+    v_square = v1^2
+    T = (rho_e - 0.5 * rho * v_square) / help1
+
+    return T
+end
+
 """
     totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D)
 

src/equations/compressible_euler_multicomponent_2d.jl

@@ -665,22 +665,57 @@ end
     s = log(p) - gamma * log(rho) - log(gas_constant)
     rho_p = rho / p
     T = (rho_e - 0.5 * rho * v_square) / (help1)
-    entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant *
-                                                                  ((gamma - s) /
-                                                                   (gamma - 1.0) -
-                                                                   (0.5 * v_square *
-                                                                    rho_p))
+
+    entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] *
+                                                                   (1 - log(T)) +
+                                                                   gas_constants[i] *
+                                                                   (1 + log(u[i + 3])) -
+                                                                   v_square / (2 * T))
                                                                   for i in eachcomponent(equations))
 
     w1 = gas_constant * v1 * rho_p
     w2 = gas_constant * v2 * rho_p
-    w3 = gas_constant * rho_p * (-1)
+    w3 = gas_constant * (-rho_p)
 
     entrop_other = SVector{3, real(equations)}(w1, w2, w3)
 
     return vcat(entrop_other, entrop_rho)
 end
 
+# Convert entropy variables to conservative variables
+@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D)
+    @unpack gammas, gas_constants, cp, cv = equations
+    T = -1 / w[3]
+    v1 = w[1] * T
+    v2 = w[2] * T
+    v_squared = v1^2 + v2^2
+    cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] -
+                                                                     cv[i] *
+                                                                     (1 - log(T)) +
+                                                                     v_squared /
+                                                                     (2 * T)) /
+                                                                    gas_constants[i] -
+                                                                    1)
+                                                                for i in eachcomponent(equations))
+
+    rho = zero(cons_rho[1])
+    help1 = zero(cons_rho[1])
+    help2 = zero(cons_rho[1])
+    p = zero(cons_rho[1])
+    for i in eachcomponent(equations)
+        rho += cons_rho[i]
+        help1 += cons_rho[i] * cv[i] * gammas[i]
+        help2 += cons_rho[i] * cv[i]
+        p += cons_rho[i] * gas_constants[i] * T
+    end
+    u1 = rho * v1
+    u2 = rho * v2
+    gamma = help1 / help2
+    u3 = p / (gamma - 1) + 0.5 * rho * v_squared
+    cons_other = SVector{3, real(equations)}(u1, u2, u3)
+    return vcat(cons_other, cons_rho)
+end
+
 # Convert primitive to conservative variables
 @inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D)
     @unpack cv, gammas = equations
@@ -700,6 +735,39 @@ end
     return vcat(cons_other, cons_rho)
 end
 
+@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D)
+    @unpack cv, gammas, gas_constants = equations
+    rho = density(u, equations)
+    T = temperature(u, equations)
+
+    total_entropy = zero(u[1])
+    for i in eachcomponent(equations)
+        total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3]))
+    end
+
+    return total_entropy
+end
+
+@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D)
+    @unpack cv, gammas, gas_constants = equations
+
+    rho_v1, rho_v2, rho_e = u
+
+    rho = density(u, equations)
+    help1 = zero(rho)
+
+    for i in eachcomponent(equations)
+        help1 += u[i + 3] * cv[i]
+    end
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v_square = v1^2 + v2^2
+    T = (rho_e - 0.5 * rho * v_square) / help1
+
+    return T
+end
+
 """
     totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D)
 

test/test_tree_1d_eulermulti.jl

@@ -2,6 +2,7 @@ module TestExamples1DEulerMulti
 
 using Test
 using Trixi
+using ForwardDiff
 
 include("test_trixi.jl")
 
@@ -10,6 +11,18 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 @testset "Compressible Euler Multicomponent" begin
 #! format: noindent
 
+@trixi_testset "Testing entropy2cons and cons2entropy" begin
+    using ForwardDiff
+    gammas = (1.3272378792562836, 1.5269959187969864, 1.8362285750521512, 1.0409061360276926, 1.4652015053812224, 1.3626493264184423)
+    gas_constants = (1.817636851910076, 6.760820475922636, 5.588953939749113, 6.31574782981543, 3.362932038038397, 3.212779569399733)
+    equations = CompressibleEulerMulticomponentEquations1D(gammas=SVector{length(gammas)}(gammas...), 
+                                                                gas_constants=SVector{length(gas_constants)}(gas_constants...))
+    u = [-1.4632513788889214, 0.9908786980927811, 0.2909066990257628, 0.6256623915420473, 0.4905882754313441, 0.14481800501749112, 1.0333532872771651, 0.6805599818745411]
+    w = cons2entropy(u, equations)
+    @test w ≈ ForwardDiff.gradient(u -> total_entropy(u, equations), u)
+    @test entropy2cons(w, equations) ≈ u
+end
+
 @trixi_testset "elixir_eulermulti_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
                         l2=[0.15330089521538684, 0.4417674632047301,

test/test_tree_2d_eulermulti.jl

@@ -10,6 +10,18 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
 @testset "Compressible Euler Multicomponent" begin
 #! format: noindent
 
+@trixi_testset "Testing entropy2cons and cons2entropy" begin
+    using ForwardDiff
+    gammas = (1.1546412974182538, 1.1171560258914812, 1.097107661471476, 1.0587601652669245, 1.6209889683979308, 1.6732209755396386, 1.2954303574165822)
+    gas_constants = (5.969461071171914, 3.6660802003290183, 6.639008614675539, 8.116604827140456, 6.190706056680031, 1.6795013743693712, 2.197737590916966)
+    equations = CompressibleEulerMulticomponentEquations2D(gammas=SVector{length(gammas)}(gammas...), 
+                                                                gas_constants=SVector{length(gas_constants)}(gas_constants...))
+    u = [-1.7433292819144075, 0.8844413258376495, 0.6050737175812364, 0.8261998359817043, 1.0801186290896465, 0.505654488367698, 0.6364415555805734, 0.851669392285058, 0.31219606420306223, 1.0930477805612038]
+    w = cons2entropy(u, equations)
+    @test w ≈ ForwardDiff.gradient(u -> total_entropy(u, equations), u)
+    @test entropy2cons(w, equations) ≈ u  
+end
+
 # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the
 #       simulations are set up with dimensional states. For example, the reference pressure in SI
 #       units is 101325 Pa, i.e., pressure has values of O(10^5)