HLL MHD Breaking changes for Consistency

examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl

@@ -12,7 +12,9 @@ initial_condition = initial_condition_convergence_test
 
 # Get the DG approximation space
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 4,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0.0, 0.0)

examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl

@@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-1.0, -1.0, -1.0)

examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl

@@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 5, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 5,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 # Create the mesh

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl

@@ -32,7 +32,7 @@ initial_condition = initial_condition_briowu_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = flux_hlle
 volume_flux = flux_derigs_etal
 basis = LobattoLegendreBasis(4)
 

examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl

@@ -39,7 +39,7 @@ initial_condition = initial_condition_ryujones_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = flux_hlle
 volume_flux = flux_hindenlang_gassner
 basis = LobattoLegendreBasis(3)
 

examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl

@@ -61,7 +61,7 @@ initial_condition = initial_condition_shu_osher_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = flux_hlle
 volume_flux = flux_hindenlang_gassner
 basis = LobattoLegendreBasis(4)
 

examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl

@@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations2D(gamma)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0.0, 0.0)

examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl

@@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-1.0, -1.0, -1.0)

examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl

@@ -11,7 +11,9 @@ equations = IdealGlmMhdEquations2D(gamma)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg = 7, surface_flux = (flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg = 7,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 # Get the unstructured quad mesh from a file (downloads the file if not available locally)

src/equations/ideal_glm_mhd_1d.jl

@@ -277,6 +277,22 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::IdealGlmMhdEquations1D)
+    rho_ll, rho_v1_ll, _ = u_ll
+    rho_rr, rho_v1_rr, _ = u_rr
+
+    # Calculate primitive variables
+    v1_ll = rho_v1_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+
+    λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+    λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+
+    return λ_min, λ_max
+end
+
 # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations1D)
@@ -298,15 +314,15 @@ end
 end
 
 """
-    min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
 
 Calculate minimum and maximum wave speeds for HLL-type fluxes as in
 - Li (2005)
   An HLLC Riemann solver for magneto-hydrodynamics
   [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
 """
-@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
-                                     equations::IdealGlmMhdEquations1D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        equations::IdealGlmMhdEquations1D)
     rho_ll, rho_v1_ll, _ = u_ll
     rho_rr, rho_v1_rr, _ = u_rr
 

src/equations/ideal_glm_mhd_2d.jl

@@ -775,6 +775,56 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate primitive velocity variables
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+
+    # Approximate the left-most and right-most eigenvalues in the Riemann fan
+    if orientation == 1 # x-direction
+        λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+        λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+    else # y-direction
+        λ_min = v2_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+        λ_max = v2_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+    end
+
+    return λ_min, λ_max
+end
+
+@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate primitive velocity variables
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+
+    v_normal_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2])
+    v_normal_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2])
+
+    c_f_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
+    c_f_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
+
+    # Estimate the min/max eigenvalues in the normal direction
+    λ_min = min(v_normal_ll - c_f_ll, v_normal_rr - c_f_rr)
+    λ_max = max(v_normal_rr + c_f_rr, v_normal_rr + c_f_rr)
+
+    return λ_min, λ_max
+end
+
 # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations2D)
@@ -833,15 +883,15 @@ end
 end
 
 """
-    min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D)
 
 Calculate minimum and maximum wave speeds for HLL-type fluxes as in
 - Li (2005)
   An HLLC Riemann solver for magneto-hydrodynamics
   [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
 """
-@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
-                                     equations::IdealGlmMhdEquations2D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        equations::IdealGlmMhdEquations2D)
     rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
 
@@ -870,8 +920,8 @@ Calculate minimum and maximum wave speeds for HLL-type fluxes as in
     return λ_min, λ_max
 end
 
-@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
-                                     equations::IdealGlmMhdEquations2D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector,
+                                        equations::IdealGlmMhdEquations2D)
     rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
 

src/equations/ideal_glm_mhd_3d.jl

@@ -670,6 +670,64 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::IdealGlmMhdEquations3D)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr
+
+    # Calculate primitive variables and speed of sound
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v3_rr = rho_v3_rr / rho_rr
+
+    # Approximate the left-most and right-most eigenvalues in the Riemann fan
+    if orientation == 1 # x-direction
+        λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+        λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+    elseif orientation == 2 # y-direction
+        λ_min = v2_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+        λ_max = v2_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+    else # z-direction
+        λ_min = v3_ll - calc_fast_wavespeed(u_ll, orientation, equations)
+        λ_max = v3_rr + calc_fast_wavespeed(u_rr, orientation, equations)
+    end
+
+    return λ_min, λ_max
+end
+
+@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::IdealGlmMhdEquations3D)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr
+
+    # Calculate primitive velocity variables
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+    v3_rr = rho_v3_rr / rho_rr
+
+    v_normal_ll = (v1_ll * normal_direction[1] +
+                   v2_ll * normal_direction[2] +
+                   v3_ll * normal_direction[3])
+    v_normal_rr = (v1_rr * normal_direction[1] +
+                   v2_rr * normal_direction[2] +
+                   v3_rr * normal_direction[3])
+
+    # Estimate the min/max eigenvalues in the normal direction
+    λ_min = v_normal_ll - calc_fast_wavespeed(u_ll, normal_direction, equations)
+    λ_max = v_normal_rr + calc_fast_wavespeed(u_rr, normal_direction, equations)
+
+    return λ_min, λ_max
+end
+
 # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations3D)
@@ -742,15 +800,15 @@ end
 end
 
 """
-    min_max_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations3D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations3D)
 
 Calculate minimum and maximum wave speeds for HLL-type fluxes as in
 - Li (2005)
   An HLLC Riemann solver for magneto-hydrodynamics
   [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020)
 """
-@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
-                                     equations::IdealGlmMhdEquations3D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        equations::IdealGlmMhdEquations3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr
 
@@ -787,8 +845,8 @@ Calculate minimum and maximum wave speeds for HLL-type fluxes as in
     return λ_min, λ_max
 end
 
-@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
-                                     equations::IdealGlmMhdEquations3D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector,
+                                        equations::IdealGlmMhdEquations3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr
 

test/test_tree_2d_mhd.jl

@@ -208,7 +208,8 @@ end
                             0.04879208429337193,
                         ],
                         tspan=(0.0, 0.06),
-                        surface_flux=(flux_hll, flux_nonconservative_powell))
+                        surface_flux=(flux_hlle,
+                                      flux_nonconservative_powell))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let

test/test_tree_3d_mhd.jl

@@ -240,7 +240,8 @@ end
                                       1000
                             end
                         end,
-                        surface_flux=(flux_hll, flux_nonconservative_powell),
+                        surface_flux=(flux_hlle,
+                                      flux_nonconservative_powell),
                         volume_flux=(flux_central, flux_nonconservative_powell),
                         coordinates_min=(0.0, 0.0, 0.0),
                         coordinates_max=(1.0, 1.0, 1.0),

test/test_unit.jl

@@ -960,14 +960,12 @@ end
 end
 
 @timed_testset "Consistency check for HLLE flux: MHD" begin
-    # Note: min_max_speed_naive for MHD is essentially min_max_speed_einfeldt
-
     equations = IdealGlmMhdEquations1D(1.4)
     u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1),
         SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2)]
 
     for u in u_values
-        @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)
+        @test flux_hlle(u, u, 1, equations) ≈ flux(u, 1, equations)
     end
 
     equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#
@@ -981,11 +979,11 @@ end
         SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]
 
     for u in u_values, orientation in orientations
-        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
     end
 
     for u in u_values, normal_direction in normal_directions
-        @test flux_hll(u, u, normal_direction, equations) ≈
+        @test flux_hlle(u, u, normal_direction, equations) ≈
               flux(u, normal_direction, equations)
     end
 
@@ -1001,11 +999,11 @@ end
         SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)]
 
     for u in u_values, orientation in orientations
-        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+        @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations)
     end
 
     for u in u_values, normal_direction in normal_directions
-        @test flux_hll(u, u, normal_direction, equations) ≈
+        @test flux_hlle(u, u, normal_direction, equations) ≈
               flux(u, normal_direction, equations)
     end
 end