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

[jlchan]

Adding the stationary contact-preserving matrix dissipation operator of Winters, Derigs, Gassner, and Walch (2017) from https://doi.org/10.1016/j.jcp.2016.12.006. Only CompressibleEulerEquations1D, CompressibleEulerEquations2D are implemented at the moment.

The implementation is directly adapted from Atum.jl

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.87%. Comparing base (31e3c8f) to head (e40b99b).

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2291      +/-   ##
==========================================
+ Coverage   96.86%   96.87%   +0.01%     
==========================================
  Files         490      490              
  Lines       39518    39593      +75     
==========================================
+ Hits        38279    38354      +75     
  Misses       1239     1239              

Flag	Coverage Δ
unittests	96.87% <100.00%> (+0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[andrewwinters5000]

Thanks for adding this @jlchan ! My only additional thought would be to add the 3D version for the Euler equations while we are at it.

[andrewwinters5000]

@jlchan I took the liberty of preparing what I think is a working version of the 3D variant. It would need an additional test or so but it is in the file below

wdgw-matrix-diss.jl

@inline function (dissipation::DissipationMatrixWintersEtal)(u_ll, u_rr,
                                                             normal_direction::AbstractVector,
                                                             equations::CompressibleEulerEquations3D)
    (; gamma) = equations

    norm_ = norm(normal_direction)
    unit_normal_direction = normal_direction / norm_

    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)

    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)
    v3_avg = avg(v3_ll, v3_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 + 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)

    v_avg_normal = dot(SVector(v1_avg, v2_avg, v3_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 # scaled repeated eigenvalue in the tangential direction

    v1_minus_c = v1_avg - c_bar * unit_normal_direction[1]
    v2_minus_c = v2_avg - c_bar * unit_normal_direction[2]
    v3_minus_c = v3_avg - c_bar * unit_normal_direction[3]
    v1_plus_c = v1_avg + c_bar * unit_normal_direction[1]
    v2_plus_c = v2_avg + c_bar * unit_normal_direction[2]
    v3_plus_c = v3_avg + c_bar * unit_normal_direction[3]
    v1_tangential = v1_avg - v_avg_normal * unit_normal_direction[1]
    v2_tangential = v2_avg - v_avg_normal * unit_normal_direction[2]
    v3_tangential = v3_avg - v_avg_normal * unit_normal_direction[3]

    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_v3_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 +
                                             v3_minus_c * entropy_var_rho_v3_jump
    velocity_plus_c_dot_entropy_vars_jump = v1_plus_c * entropy_var_rho_v1_jump +
                                            v2_plus_c * entropy_var_rho_v2_jump +
                                            v3_plus_c * entropy_var_rho_v3_jump
    velocity_avg_dot_vjump = v1_avg * entropy_var_rho_v1_jump +
                             v2_avg * entropy_var_rho_v2_jump +
                             v3_avg * entropy_var_rho_v3_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,
                                            entropy_var_rho_v3_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))

    dissipation_rho_v3 = (w1 * v3_minus_c +
                          w2 * v3_avg +
                          w3 * v3_plus_c +
                          lambda_4 * (entropy_var_rho_v3_jump -
                           unit_normal_direction[3] * entropy_var_v_normal_jump +
                           entropy_var_rho_e_jump * v3_tangential))

    v_tangential_dot_entropy_vars_jump = v1_tangential * entropy_var_rho_v1_jump +
                                         v2_tangential * entropy_var_rho_v2_jump +
                                         v3_tangential * entropy_var_rho_v3_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 + v3_avg^2 - v_avg_normal^2)))

    return -0.5f0 *
           SVector(dissipation_rho, dissipation_rho_v1, dissipation_rho_v2,
                   dissipation_rho_v2, dissipation_rhoe) * norm_
end

[jlchan]

@jlchan I took the liberty of preparing what I think is a working version of the 3D variant. It would need an additional test or so but it is in the file below

Awesome, thanks! I'll add it into the PR but it probably will be next week (some travel this weekend).

[andrewwinters5000]

@jlchan I took the liberty of preparing what I think is a working version of the 3D variant. It would need an additional test or so but it is in the file below

Awesome, thanks! I'll add it into the PR but it probably will be next week (some travel this weekend).

Apologies, ignore the previous code as there were mistakes in the way the tangential directions were handled. I got confused with the heavily extracted way of computing the matrix dissipation term. Below is an implementation reminiscent of how it is done in FLUXO where we use rotational invariance to rotate the solution, compute the dissipation term, and then back-rotate the result.

diss-matrix-3d.jl

@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 = avg(v1_ll, v1_rr)
    v2_avg = avg(v2_ll, v2_rr)
    v3_avg = avg(v3_ll, v3_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 + 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

    # 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
    for j in 1:5
       diss[j] *= D[j]
    end
    # Finally apply the remaining eigenvector matrix
    diss = R * diss

    # Step 4:
    # Do not forget to backrotate my young padawan and then return with proper normalization scaling
    return -0.5f0 * rotate_from_x(diss, normal_vector, tangent1, tangent2, equations) * norm_
end

[andrewwinters5000]

Unsure if we will include the 3D version in this PR but I just left a comment about it for completeness.

[andrewwinters5000]

If we decide to include the 3D we might also want to mention that the 1D/2D came from Atum.jl and the 3D was helped by FLUXO (https://github.com/project-fluxo/fluxo/blob/master/src/equation/navierstokes/riemann.f90#L901)