Merge branch 'main' into sc/converters_coupling

Project.toml

@@ -6,6 +6,7 @@ version = "0.6.5-pre"
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
@@ -52,6 +53,7 @@ TrixiMakieExt = "Makie"
 [compat]
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
+DataStructures = "0.18.15"
 DiffEqCallbacks = "2.25"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"
@@ -84,7 +86,7 @@ StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.4.3"
+T8code = "0.4.3, 0.5"
 TimerOutputs = "0.5.7"
 Triangulate = "2.0"
 TriplotBase = "0.1"

examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl

@@ -67,7 +67,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval = 100,
+save_solution = SaveSolutionCallback(dt = 0.1,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)

src/Trixi.jl

@@ -37,7 +37,8 @@ using SciMLBase: CallbackSet, DiscreteCallback,
 import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   AbstractODEIntegrator, init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
-                  terminate!, remake
+                  terminate!, remake, add_tstop!, has_tstop, first_tstop
+
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
@@ -70,6 +71,7 @@ using TriplotBase: TriplotBase
 using TriplotRecipes: DGTriPseudocolor
 @reexport using SimpleUnPack: @unpack
 using SimpleUnPack: @pack!
+using DataStructures: BinaryHeap, FasterForward, extract_all!
 
 # finite difference SBP operators
 using SummationByPartsOperators: AbstractDerivativeOperator,

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)
 

src/time_integration/methods_2N.jl

@@ -204,6 +204,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator2N, dt)
     integrator.dt = dt
 end
 
+# Required e.g. for `glm_speed_callback` 
+function get_proposed_dt(integrator::SimpleIntegrator2N)
+    return integrator.dt
+end
+
 # stop the time integration
 function terminate!(integrator::SimpleIntegrator2N)
     integrator.finalstep = true

src/time_integration/methods_3Sstar.jl

@@ -282,6 +282,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator3Sstar, dt)
     integrator.dt = dt
 end
 
+# Required e.g. for `glm_speed_callback` 
+function get_proposed_dt(integrator::SimpleIntegrator3Sstar)
+    return integrator.dt
+end
+
 # stop the time integration
 function terminate!(integrator::SimpleIntegrator3Sstar)
     integrator.finalstep = true