Add Winters et al matrix dissipation for 1D, 2D Euler (#2291)

* adding Winters et al matrix dissipation in 1D and 2D

* exporting dissipation operator

* adding tests

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* renaming

* fix naming in tests

* define avg

* fix formatting

* Apply suggestions from code review

Co-authored-by: Andrew Winters <[email protected]>

* more formatting

* adding LinAlg: norm

* add LinAlg: dot

* formatting

* adding 3D matrix dissipation operator

* adding 3D matrix dissipation tests

* add modified sod example and test

* credit fluxo

* Apply suggestions from code review

Co-authored-by: Andrew Winters <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>

* remove avg(...) calls

fixing 0.5f -> 0.5f0

* add specialization for orientation::Integer

* adding 2D and 3D tests

* format modified Sod elixir

* add DOI to comment

* fix tree_2d_dgsem test

* Update test/test_tree_2d_euler.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* formatting

---------

Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Andrew Winters <[email protected]>

Author: 3 people

examples/dgmulti_1d/elixir_euler_modified_sod.jl

@@ -0,0 +1,44 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+surface_flux = FluxPlusDissipation(flux_ranocha, DissipationMatrixWintersEtal())
+volume_flux = flux_ranocha
+dg = DGMulti(polydeg = 3, element_type = Line(), approximation_type = SBP(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+equations = CompressibleEulerEquations1D(1.4)
+
+function initial_condition_modified_sod(x, t, ::CompressibleEulerEquations1D)
+    if x[1] < 0.3
+        return prim2cons(SVector(1, 0.75, 1), equations)
+    else
+        return prim2cons(SVector(0.125, 0.0, 0.1), equations)
+    end
+end
+
+initial_condition = initial_condition_modified_sod
+
+cells_per_dimension = (50,)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (0.0,), coordinates_max = (1.0,), periodicity = false)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval = 100)
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(); adaptive = false,
+            dt = 0.5 * estimate_dt(mesh, dg),
+            ode_default_options()...,
+            callback = callbacks);

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 = 0.5f0 * (v1_ll + v1_rr)
+    p_avg = 0.5f0 * (rho_ll + rho_rr) / (b_ll + b_rr) # 2 * b_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 = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (rho_ll + rho_rr) / (b_ll + b_rr) # 2 * b_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/compressible_euler_3d.jl

@@ -1237,6 +1237,100 @@ end
                    u[5])
 end
 
+@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             normal_direction::AbstractVector,
+                                                             equations::CompressibleEulerEquations3D)
+    (; gamma) = equations
+
+    # Step 1:
+    # Rotate solution into the appropriate direction
+
+    norm_ = norm(normal_direction)
+    # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
+    normal_vector = normal_direction / norm_
+
+    # Some vector that can't be identical to normal_vector (unless normal_vector == 0)
+    tangent1 = SVector(normal_direction[2], normal_direction[3], -normal_direction[1])
+    # Orthogonal projection
+    tangent1 -= dot(normal_vector, tangent1) * normal_vector
+    tangent1 = normalize(tangent1)
+
+    # Third orthogonal vector
+    tangent2 = normalize(cross(normal_direction, tangent1))
+
+    u_ll_rotated = rotate_to_x(u_ll, normal_vector, tangent1, tangent2, equations)
+    u_rr_rotated = rotate_to_x(u_rr, normal_vector, tangent1, tangent2, equations)
+
+    # Step 2:
+    # Compute the averages using the rotated variables
+    rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll_rotated, equations)
+    rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr_rotated, 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 = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    v3_avg = 0.5f0 * (v3_ll + v3_rr)
+    p_avg = 0.5f0 * (rho_ll + rho_rr) / (b_ll + b_rr)
+    v_squared_bar = v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr
+    h_bar = gamma / (2 * b_log * (gamma - 1)) + 0.5f0 * v_squared_bar
+    c_bar = sqrt(gamma * p_avg / rho_log)
+
+    # Step 3:
+    # Build the dissipation term as given in Appendix A of the paper 
+    # - 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).
+
+    # Get entropy variables jump in the rotated variables
+    w_jump = cons2entropy(u_rr_rotated, equations) -
+             cons2entropy(u_ll_rotated, equations)
+
+    # Entries of the diagonal scaling matrix where D = ABS(\Lambda)T
+    lambda_1 = abs(v1_avg - c_bar) * rho_log / (2 * gamma)
+    lambda_2 = abs(v1_avg) * rho_log * (gamma - 1) / gamma
+    lambda_3 = abs(v1_avg) * p_avg # scaled repeated eigenvalue in the tangential direction
+    lambda_5 = abs(v1_avg + c_bar) * rho_log / (2 * gamma)
+    D = SVector(lambda_1, lambda_2, lambda_3, lambda_3, lambda_5)
+
+    # Entries of the right eigenvector matrix (others have already been precomputed)
+    r21 = v1_avg - c_bar
+    r25 = v1_avg + c_bar
+    r51 = h_bar - v1_avg * c_bar
+    r52 = 0.5f0 * v_squared_bar
+    r55 = h_bar + v1_avg * c_bar
+
+    # Build R and transpose of R matrices
+    R = @SMatrix [[1;; 1;; 0;; 0;; 1];
+                  [r21;; v1_avg;; 0;; 0;; r25];
+                  [v2_avg;; v2_avg;; 1;; 0;; v2_avg];
+                  [v3_avg;; v3_avg;; 0;; 1;; v3_avg];
+                  [r51;; r52;; v2_avg;; v3_avg;; r55]]
+
+    RT = @SMatrix [[1;; r21;; v2_avg;; v3_avg;; r51];
+                   [1;; v1_avg;; v2_avg;; v3_avg;; r52];
+                   [0;; 0;; 1;; 0;; v2_avg];
+                   [0;; 0;; 0;; 1;; v3_avg];
+                   [1;; r25;; v2_avg;; v3_avg;; r55]]
+
+    # Compute the dissipation term R * D * R^T * [[w]] from right-to-left
+
+    # First comes R^T * [[w]]
+    diss = RT * w_jump
+    # Next multiply with the eigenvalues and Barth scaling
+    diss = D .* diss
+    # Finally apply the remaining eigenvector matrix
+    diss = R * diss
+
+    # Step 4:
+    # Do not forget to backrotate and then return with proper normalization scaling
+    return -0.5f0 * rotate_from_x(diss, normal_vector, tangent1, tangent2, equations) *
+           norm_
+end
+
 # Rotate x-axis to normal vector; normal, tangent1 and tangent2 need to be orthonormal
 # Called inside `FluxRotated` in `numerical_fluxes.jl` so the directions
 # has been normalized prior to this back-rotation of the state vector

src/equations/numerical_fluxes.jl

@@ -261,6 +261,54 @@ 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())
+```
+The 1D and 2D implementations are adapted from the [Atum.jl library](https://github.com/mwarusz/Atum.jl/blob/c7ed44f2b7972ac726ef345da7b98b0bda60e2a3/src/balancelaws/euler.jl#L198).
+The 3D implementation is adapted from the [FLUXO library](https://github.com/project-fluxo/fluxo)
+
+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 function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             orientation::Integer,
+                                                             equations::AbstractEquations{1})
+    return dissipation(u_ll, u_rr, SVector(1), equations)
+end
+
+@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             orientation::Integer,
+                                                             equations::AbstractEquations{2})
+    if orientation == 1
+        return dissipation(u_ll, u_rr, SVector(1, 0), equations)
+    else # orientation == 2
+        return dissipation(u_ll, u_rr, SVector(0, 1), equations)
+    end
+end
+
+@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
+                                                             orientation::Integer,
+                                                             equations::AbstractEquations{3})
+    if orientation == 1
+        return dissipation(u_ll, u_rr, SVector(1, 0, 0), equations)
+    elseif orientation == 2
+        return dissipation(u_ll, u_rr, SVector(0, 1, 0), equations)
+    else # orientation == 3
+        return dissipation(u_ll, u_rr, SVector(0, 0, 1), equations)
+    end
+end
+
 """
     FluxHLL(min_max_speed=min_max_speed_davis)
 

test/test_dgmulti_1d.jl

@@ -165,6 +165,26 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_modified_sod.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "dgmulti_1d",
+                                 "elixir_euler_modified_sod.jl"),
+                        cells_per_dimension=(16,),
+                        l2=[0.26352391505659767, 0.4528974787813885, 0.9310255091126164],
+                        linf=[
+                            0.6268146194274395,
+                            0.8214003799995101,
+                            1.8606901431409795
+                        ])
+    # 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_fdsbp_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         l2=[