Add Winters et al matrix dissipation for 1D, 2D Euler

src/Trixi.jl

@@ -190,7 +190,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
-       DissipationLaxFriedrichsEntropyVariables,
+       DissipationLaxFriedrichsEntropyVariables, DissipationMatrixWintersEtal,
        FluxLaxFriedrichs, max_abs_speed_naive,
        FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt,
        FluxLMARS,

src/equations/compressible_euler_1d.jl

@@ -708,6 +708,64 @@ end
     λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
 end
 
+@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             normal_direction::AbstractVector,
+                                                             equations::CompressibleEulerEquations1D)
+    gamma = equations.gamma
+
+    norm_ = norm(normal_direction)
+    unit_normal_direction = normal_direction / norm_
+
+    rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
+
+    b_ll = rho_ll / (2 * p_ll)
+    b_rr = rho_rr / (2 * p_rr)
+
+    rho_log = ln_mean(rho_ll, rho_rr)
+    b_log = ln_mean(b_ll, b_rr)
+    v1_avg = avg(v1_ll, v1_rr)
+    p_avg = avg(rho_ll, rho_rr) / (2 * avg(b_ll, b_rr))
+    v1_squared_bar = v1_ll * v1_rr
+    h_bar = gamma / (2 * b_log * (gamma - 1)) + 0.5f0 * v1_squared_bar
+    c_bar = sqrt(gamma * p_avg / rho_log)
+
+    v1_minus_c = v1_avg - c_bar * unit_normal_direction[1]
+    v1_plus_c = v1_avg + c_bar * unit_normal_direction[1]
+    v_avg_normal = v1_avg * unit_normal_direction[1]
+
+    entropy_vars_jump = cons2entropy(u_rr, equations) - cons2entropy(u_ll, equations)
+
+    lambda_1 = abs(v_avg_normal - c_bar) * rho_log / (2 * gamma)
+    lambda_2 = abs(v_avg_normal) * rho_log * (gamma - 1) / gamma
+    lambda_3 = abs(v_avg_normal + c_bar) * rho_log / (2 * gamma)
+    lambda_4 = abs(v_avg_normal) * p_avg
+
+    entropy_var_rho_jump, entropy_var_rho_v1_jump, entropy_var_rho_e_jump = entropy_vars_jump
+
+    w1 = lambda_1 * (entropy_var_rho_jump + v1_minus_c * entropy_var_rho_v1_jump +
+          (h_bar - c_bar * v_avg_normal) * entropy_var_rho_e_jump)
+    w2 = lambda_2 * (entropy_var_rho_jump + v1_avg * entropy_var_rho_v1_jump +
+          v1_squared_bar / 2 * entropy_var_rho_e_jump)
+    w3 = lambda_3 * (entropy_var_rho_jump + v1_plus_c * entropy_var_rho_v1_jump +
+          (h_bar + c_bar * v_avg_normal) * entropy_var_rho_e_jump)
+
+    dissipation_rho = w1 + w2 + w3
+
+    dissipation_rho_v1 = (w1 * v1_minus_c +
+                          w2 * v1_avg +
+                          w3 * v1_plus_c)
+
+    dissipation_rhoe = (w1 * (h_bar - c_bar * v_avg_normal) +
+                        w2 * 0.5f0 * v1_squared_bar +
+                        w3 * (h_bar + c_bar * v_avg_normal) +
+                        lambda_4 *
+                        (entropy_var_rho_e_jump * (v1_avg * v1_avg - v_avg_normal^2)))
+
+    return -0.5f0 * SVector(dissipation_rho, dissipation_rho_v1, dissipation_rhoe) *
+           norm_
+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::CompressibleEulerEquations1D)

src/equations/compressible_euler_2d.jl

@@ -1507,6 +1507,95 @@ end
     return λ_min, λ_max
 end
 
+@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             normal_direction::AbstractVector,
+                                                             equations::CompressibleEulerEquations2D)
+    (; gamma) = equations
+
+    norm_ = norm(normal_direction)
+    unit_normal_direction = normal_direction / norm_
+
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    b_ll = rho_ll / (2 * p_ll)
+    b_rr = rho_rr / (2 * p_rr)
+
+    rho_log = ln_mean(rho_ll, rho_rr)
+    b_log = ln_mean(b_ll, b_rr)
+    v1_avg = avg(v1_ll, v1_rr)
+    v2_avg = avg(v2_ll, v2_rr)
+    p_avg = avg(rho_ll, rho_rr) / (2 * avg(b_ll, b_rr))
+    v_squared_bar = v1_ll * v1_rr + v2_ll * v2_rr
+    h_bar = gamma / (2 * b_log * (gamma - 1)) + 0.5f0 * v_squared_bar
+    c_bar = sqrt(gamma * p_avg / rho_log)
+
+    v_avg_normal = dot(SVector(v1_avg, v2_avg), unit_normal_direction)
+
+    lambda_1 = abs(v_avg_normal - c_bar) * rho_log / (2 * gamma)
+    lambda_2 = abs(v_avg_normal) * rho_log * (gamma - 1) / gamma
+    lambda_3 = abs(v_avg_normal + c_bar) * rho_log / (2 * gamma)
+    lambda_4 = abs(v_avg_normal) * p_avg
+
+    v1_minus_c = v1_avg - c_bar * unit_normal_direction[1]
+    v2_minus_c = v2_avg - c_bar * unit_normal_direction[2]
+    v1_plus_c = v1_avg + c_bar * unit_normal_direction[1]
+    v2_plus_c = v2_avg + c_bar * unit_normal_direction[2]
+    v1_tangential = v1_avg - v_avg_normal * unit_normal_direction[1]
+    v2_tangential = v2_avg - v_avg_normal * unit_normal_direction[2]
+
+    entropy_vars_jump = cons2entropy(u_rr, equations) - cons2entropy(u_ll, equations)
+    entropy_var_rho_jump, entropy_var_rho_v1_jump,
+    entropy_var_rho_v2_jump, entropy_var_rho_e_jump = entropy_vars_jump
+
+    velocity_minus_c_dot_entropy_vars_jump = v1_minus_c * entropy_var_rho_v1_jump +
+                                             v2_minus_c * entropy_var_rho_v2_jump
+    velocity_plus_c_dot_entropy_vars_jump = v1_plus_c * entropy_var_rho_v1_jump +
+                                            v2_plus_c * entropy_var_rho_v2_jump
+    velocity_avg_dot_vjump = v1_avg * entropy_var_rho_v1_jump +
+                             v2_avg * entropy_var_rho_v2_jump
+    w1 = lambda_1 * (entropy_var_rho_jump + velocity_minus_c_dot_entropy_vars_jump +
+          (h_bar - c_bar * v_avg_normal) * entropy_var_rho_e_jump)
+    w2 = lambda_2 * (entropy_var_rho_jump + velocity_avg_dot_vjump +
+          v_squared_bar / 2 * entropy_var_rho_e_jump)
+    w3 = lambda_3 * (entropy_var_rho_jump + velocity_plus_c_dot_entropy_vars_jump +
+          (h_bar + c_bar * v_avg_normal) * entropy_var_rho_e_jump)
+
+    entropy_var_v_normal_jump = dot(SVector(entropy_var_rho_v1_jump,
+                                            entropy_var_rho_v2_jump),
+                                    unit_normal_direction)
+
+    dissipation_rho = w1 + w2 + w3
+
+    dissipation_rho_v1 = (w1 * v1_minus_c +
+                          w2 * v1_avg +
+                          w3 * v1_plus_c +
+                          lambda_4 * (entropy_var_rho_v1_jump -
+                           unit_normal_direction[1] * entropy_var_v_normal_jump +
+                           entropy_var_rho_e_jump * v1_tangential))
+
+    dissipation_rho_v2 = (w1 * v2_minus_c +
+                          w2 * v2_avg +
+                          w3 * v2_plus_c +
+                          lambda_4 * (entropy_var_rho_v2_jump -
+                           unit_normal_direction[2] * entropy_var_v_normal_jump +
+                           entropy_var_rho_e_jump * v2_tangential))
+
+    v_tangential_dot_entropy_vars_jump = v1_tangential * entropy_var_rho_v1_jump +
+                                         v2_tangential * entropy_var_rho_v2_jump
+
+    dissipation_rhoe = (w1 * (h_bar - c_bar * v_avg_normal) +
+                        w2 * 0.5f0 * v_squared_bar +
+                        w3 * (h_bar + c_bar * v_avg_normal) +
+                        lambda_4 * (v_tangential_dot_entropy_vars_jump +
+                         entropy_var_rho_e_jump *
+                         (v1_avg^2 + v2_avg^2 - v_avg_normal^2)))
+
+    return -0.5f0 *
+           SVector(dissipation_rho, dissipation_rho_v1, dissipation_rho_v2,
+                   dissipation_rhoe) * norm_
+end
+
 # Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction
 # has been normalized prior to this rotation of the state vector
 @inline function rotate_to_x(u, normal_vector, equations::CompressibleEulerEquations2D)

src/equations/numerical_fluxes.jl

@@ -261,6 +261,25 @@ function Base.show(io::IO, d::DissipationLaxFriedrichsEntropyVariables)
     print(io, "DissipationLaxFriedrichsEntropyVariables(", d.max_abs_speed, ")")
 end
 
+@doc raw"""
+    DissipationMatrixWintersEtal()
+
+Creates the Roe-like entropy stable matrix dissipation operator from Winters et al. (2017). This operator
+must be used together with an entropy-conservative two-point flux function 
+(e.g., [`flux_ranocha`](@ref) or [`flux_chandrashekar`](@ref)) to yield 
+an entropy-stable surface flux. The surface flux function can be initialized as:
+```julia
+flux_es = FluxPlusDissipation(flux_ec, DissipationMatrixWintersEtal())
+```
+This implementation is adapted from the [Atum.jl library](https://github.com/mwarusz/Atum.jl/blob/c7ed44f2b7972ac726ef345da7b98b0bda60e2a3/src/balancelaws/euler.jl#L198).
+For the derivation of the matrix dissipation operator, see:
+- A. R. Winters, D. Derigs, G. Gassner, S. Walch, A uniquely defined entropy stable matrix dissipation operator 
+  for high Mach number ideal MHD and compressible Euler simulations (2017). Journal of Computational Physics.
+  [DOI: 10.1016/j.jcp.2016.12.006](https://doi.org/10.1016/j.jcp.2016.12.006).
+"""
+struct DissipationMatrixWintersEtal end
+@inline avg(u_ll, u_rr) = 0.5f0 * (u_ll + u_rr)
+
 """
     FluxHLL(min_max_speed=min_max_speed_davis)
 

test/test_unit.jl

@@ -3,6 +3,7 @@ module TestUnit
 using Test
 using Trixi
 
+using LinearAlgebra: norm, dot
 using DelimitedFiles: readdlm
 
 # Use Convex and ECOS to load the extension that extends functions for testing
@@ -1637,6 +1638,61 @@ end
     end
 end
 
+@timed_testset "DissipationMatrixWintersEtal entropy dissipation and consistency tests" begin
+    equations = CompressibleEulerEquations1D(1.4)
+    dissipation_matrix_winters_etal = DissipationMatrixWintersEtal()
+
+    # test constant preservation and entropy dissipation vector
+    u_ll = prim2cons(SVector(1, 0, 2.0), equations)
+    u_rr = prim2cons(SVector(1.1, 0, 2.0), equations)
+    v_ll = cons2entropy(u_ll, equations)
+    v_rr = cons2entropy(u_rr, equations)
+    @test norm(dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) <
+          100 * eps()
+    @test dot(v_ll - v_rr,
+              dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) ≥ 0
+
+    # test non-unit vector
+    u_ll = prim2cons(SVector(rand(), randn(), rand()), equations)
+    u_rr = prim2cons(SVector(rand(), randn(), rand()), equations)
+    v_ll = cons2entropy(u_ll, equations)
+    v_rr = cons2entropy(u_rr, equations)
+    @test dot(v_ll - v_rr,
+              dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)) ≥ 0
+    @test dissipation_matrix_winters_etal(u_ll, u_rr, SVector(0.1), equations) ≈
+          0.1 * dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0), equations)
+
+    equations = CompressibleEulerEquations2D(1.4)
+
+    # test that 2D flux is consistent with 1D matrix flux
+    u_ll = prim2cons(SVector(1, 0.1, 0, 1.0), equations)
+    u_rr = prim2cons(SVector(1.1, -0.2, 0, 2.0), equations)
+    v_ll = cons2entropy(u_ll, equations)
+    v_rr = cons2entropy(u_rr, equations)
+    normal = SVector(1.0, 0.0)
+    ids = [1, 2, 4] # indices of 1D variables/fluxes within the 2D solution
+    @test dissipation_matrix_winters_etal(u_ll, u_rr, normal, equations)[ids] ≈
+          dissipation_matrix_winters_etal(u_ll[ids], u_rr[ids],
+                                          SVector(1.0),
+                                          CompressibleEulerEquations1D(1.4))
+
+    # test 2D entropy dissipation
+    u_ll = prim2cons(SVector(1, 1, -3, 100.0), equations)
+    u_rr = prim2cons(SVector(100, -2, 4, 1.0), equations)
+    v_ll = cons2entropy(u_ll, equations)
+    v_rr = cons2entropy(u_rr, equations)
+    dissipation = dissipation_matrix_winters_etal(u_ll, u_rr, SVector(1.0, 1.0),
+                                                  equations)
+    @test dot(v_ll - v_rr, dissipation) ≥ 0
+
+    # test non-unit vector
+    normal_direction = SVector(1.0, 2.0)
+    @test dissipation_matrix_winters_etal(u_ll, u_rr, normal_direction, equations) ≈
+          norm(normal_direction) * dissipation_matrix_winters_etal(u_ll, u_rr,
+                                          normal_direction / norm(normal_direction),
+                                          equations)
+end
+
 @testset "Equivalent Wave Speed Estimates" begin
     @timed_testset "Linearized Euler 3D" begin
         equations = LinearizedEulerEquations3D(v_mean_global = (0.42, 0.37, 0.7),