Merge branch 'main' into jc/remove_cache_parabolic_interfaces

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/[email protected].23
+        uses: crate-ci/[email protected].26

NEWS.md

@@ -8,6 +8,7 @@ for human readability.
 
 #### Added
 - AMR for hyperbolic-parabolic equations on 3D `P4estMesh`
+- `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`
 
 ## Changes when updating to v0.6 from v0.5.x
 

src/callbacks_step/summary.jl

@@ -152,7 +152,13 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator;
         Polyester.reset_threads!()
     end
 
-    mpi_isroot() || return nothing
+    # The summary callback should only print information on the root process.
+    # However, all other MPI processes should also reset the timer so that
+    # it can be used to diagnose performance.
+    if !mpi_isroot()
+        reset_timer!(timer())
+        return nothing
+    end
 
     print_startup_message()
 

src/equations/compressible_euler_2d.jl

@@ -1150,7 +1150,7 @@ end
 end
 
 """
-    flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)
+    flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D)
 
 Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
 [Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
@@ -1185,18 +1185,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     if orientation == 1 # x-direction
         vel_L = v1_ll
         vel_R = v1_rr
-        ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2
     elseif orientation == 2 # y-direction
         vel_L = v2_ll
         vel_R = v2_rr
-        ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2
     end
     vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
-    ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2)
+    v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr
+    v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr
+    vel_roe_mag = (v1_roe^2 + v2_roe^2) / sum_sqrt_rho^2
     H_ll = (rho_e_ll + p_ll) / rho_ll
     H_rr = (rho_e_rr + p_rr) / rho_rr
     H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
-    c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag))
     Ssl = min(vel_L - c_ll, vel_roe - c_roe)
     Ssr = max(vel_R + c_rr, vel_roe + c_roe)
     sMu_L = Ssl - vel_L
@@ -1252,6 +1252,98 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     return SVector(f1, f2, f3, f4)
 end
 
+function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector,
+                   equations::CompressibleEulerEquations2D)
+    # Calculate primitive variables and speed of sound
+    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]
+
+    norm_ = norm(normal_direction)
+    norm_sq = norm_ * norm_
+    inv_norm_sq = inv(norm_sq)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    # Obtain left and right fluxes
+    f_ll = flux(u_ll, normal_direction, equations)
+    f_rr = flux(u_rr, normal_direction, equations)
+
+    # Compute Roe averages
+    sqrt_rho_ll = sqrt(rho_ll)
+    sqrt_rho_rr = sqrt(rho_rr)
+    sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr
+
+    v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho
+    v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho
+    vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2]
+    vel_roe_mag = v1_roe^2 + v2_roe^2
+
+    e_ll = u_ll[4] / rho_ll
+    e_rr = u_rr[4] / rho_rr
+
+    H_ll = (u_ll[4] + p_ll) / rho_ll
+    H_rr = (u_rr[4] + p_rr) / rho_rr
+
+    H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_
+
+    Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe)
+    Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe)
+    sMu_L = Ssl - v_dot_n_ll
+    sMu_R = Ssr - v_dot_n_rr
+
+    if Ssl >= 0.0
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+    elseif Ssr <= 0.0
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+    else
+        SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R +
+                 (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R)
+        if Ssl <= 0.0 <= SStar
+            densStar = rho_ll * sMu_L / (Ssl - SStar)
+            enerStar = e_ll +
+                       (SStar - v_dot_n_ll) *
+                       (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L))
+            UStar1 = densStar
+            UStar2 = densStar *
+                     (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq)
+            UStar3 = densStar *
+                     (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq)
+            UStar4 = densStar * enerStar
+            f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1])
+            f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2])
+            f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3])
+            f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4])
+        else
+            densStar = rho_rr * sMu_R / (Ssr - SStar)
+            enerStar = e_rr +
+                       (SStar - v_dot_n_rr) *
+                       (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R))
+            UStar1 = densStar
+            UStar2 = densStar *
+                     (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq)
+            UStar3 = densStar *
+                     (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq)
+            UStar4 = densStar * enerStar
+            f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1])
+            f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2])
+            f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3])
+            f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4])
+        end
+    end
+    return SVector(f1, f2, f3, f4)
+end
+
 """
     min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)
 

src/equations/compressible_euler_3d.jl

@@ -1192,7 +1192,7 @@ end
 end
 
 """
-    flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
+    flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D)
 
 Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
 [Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
@@ -1231,25 +1231,22 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     if orientation == 1 # x-direction
         vel_L = v1_ll
         vel_R = v1_rr
-        ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 +
-                   (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2
     elseif orientation == 2 # y-direction
         vel_L = v2_ll
         vel_R = v2_rr
-        ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 +
-                   (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2
     else # z-direction
         vel_L = v3_ll
         vel_R = v3_rr
-        ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 +
-                   (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2
     end
     vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
-    ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2)
+    v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr
+    v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr
+    v3_roe = sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr
+    vel_roe_mag = (v1_roe^2 + v2_roe^2 + v3_roe^2) / sum_sqrt_rho^2
     H_ll = (rho_e_ll + p_ll) / rho_ll
     H_rr = (rho_e_rr + p_rr) / rho_rr
     H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
-    c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag))
     Ssl = min(vel_L - c_ll, vel_roe - c_roe)
     Ssr = max(vel_R + c_rr, vel_roe + c_roe)
     sMu_L = Ssl - vel_L
@@ -1321,6 +1318,110 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     return SVector(f1, f2, f3, f4, f5)
 end
 
+function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector,
+                   equations::CompressibleEulerEquations3D)
+    # Calculate primitive variables and speed of sound
+    rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations)
+
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
+                 v3_ll * normal_direction[3]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] +
+                 v3_rr * normal_direction[3]
+
+    norm_ = norm(normal_direction)
+    norm_sq = norm_ * norm_
+    inv_norm_sq = inv(norm_sq)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    # Obtain left and right fluxes
+    f_ll = flux(u_ll, normal_direction, equations)
+    f_rr = flux(u_rr, normal_direction, equations)
+
+    # Compute Roe averages
+    sqrt_rho_ll = sqrt(rho_ll)
+    sqrt_rho_rr = sqrt(rho_rr)
+    sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr
+
+    v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho
+    v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho
+    v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) / sum_sqrt_rho
+    vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] +
+              v3_roe * normal_direction[3]
+    vel_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2
+
+    e_ll = u_ll[5] / rho_ll
+    e_rr = u_rr[5] / rho_rr
+
+    H_ll = (u_ll[5] + p_ll) / rho_ll
+    H_rr = (u_rr[5] + p_rr) / rho_rr
+
+    H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_
+
+    Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe)
+    Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe)
+    sMu_L = Ssl - v_dot_n_ll
+    sMu_R = Ssr - v_dot_n_rr
+
+    if Ssl >= 0.0
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+        f5 = f_ll[5]
+    elseif Ssr <= 0.0
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+        f5 = f_rr[5]
+    else
+        SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R +
+                 (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R)
+        if Ssl <= 0.0 <= SStar
+            densStar = rho_ll * sMu_L / (Ssl - SStar)
+            enerStar = e_ll +
+                       (SStar - v_dot_n_ll) *
+                       (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L))
+            UStar1 = densStar
+            UStar2 = densStar *
+                     (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq)
+            UStar3 = densStar *
+                     (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq)
+            UStar4 = densStar *
+                     (v3_ll + (SStar - v_dot_n_ll) * normal_direction[3] * inv_norm_sq)
+            UStar5 = densStar * enerStar
+            f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1])
+            f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2])
+            f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3])
+            f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4])
+            f5 = f_ll[5] + Ssl * (UStar5 - u_ll[5])
+        else
+            densStar = rho_rr * sMu_R / (Ssr - SStar)
+            enerStar = e_rr +
+                       (SStar - v_dot_n_rr) *
+                       (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R))
+            UStar1 = densStar
+            UStar2 = densStar *
+                     (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq)
+            UStar3 = densStar *
+                     (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq)
+            UStar4 = densStar *
+                     (v3_rr + (SStar - v_dot_n_rr) * normal_direction[3] * inv_norm_sq)
+            UStar5 = densStar * enerStar
+            f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1])
+            f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2])
+            f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3])
+            f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4])
+            f5 = f_rr[5] + Ssr * (UStar5 - u_rr[5])
+        end
+    end
+    return SVector(f1, f2, f3, f4, f5)
+end
+
 """
     min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
 

test/test_p4est_2d.jl

@@ -203,6 +203,32 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_sedov.jl with HLLC Flux" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
+                        l2=[
+                            0.4229948321239887,
+                            0.2559038337457483,
+                            0.2559038337457484,
+                            1.2990046683564136,
+                        ],
+                        linf=[
+                            1.4989357969730492,
+                            1.325456585141623,
+                            1.3254565851416251,
+                            6.331283015053501,
+                        ],
+                        surface_flux=flux_hllc,
+                        tspan=(0.0, 0.3))
+    # 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_sedov.jl (HLLE)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
                         l2=[

test/test_p4est_3d.jl

@@ -234,6 +234,34 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
+                        l2=[
+                            8.444868392439035e-16,
+                            4.889826056731442e-15,
+                            2.2921260987087585e-15,
+                            4.268460455702414e-15,
+                            1.1356712092620279e-14,
+                        ],
+                        linf=[
+                            7.749356711883593e-14,
+                            4.513472928735496e-13,
+                            2.9790059308254513e-13,
+                            1.057154364048074e-12,
+                            1.6271428648906294e-12,
+                        ],
+                        tspan=(0.0, 0.1),
+                        surface_flux=flux_hllc)
+    # 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_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
                         l2=[