From a33ad14b9ae232deda6d118cd5d359d27199c2db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 24 Jan 2024 16:28:14 +0100 Subject: [PATCH] Added new ES numerical fluxes --- src/Trixi.jl | 1 + src/equations/ideal_mhd_multiion_2d.jl | 268 +++++++++++++++++++++++++ 2 files changed, 269 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index 663c07a29e5..eb65828692f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -178,6 +178,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, + DissipationEntropyStableLump, DissipationEntropyStable, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, min_max_speed_chen_noelle, diff --git a/src/equations/ideal_mhd_multiion_2d.jl b/src/equations/ideal_mhd_multiion_2d.jl index 985c260cbe7..c508b6a9fe2 100644 --- a/src/equations/ideal_mhd_multiion_2d.jl +++ b/src/equations/ideal_mhd_multiion_2d.jl @@ -1043,4 +1043,272 @@ Computes the sum of the densities times the sum of the pressures end return rho_total * p_total end + +""" + DissipationEntropyStableLump(max_abs_speed=max_abs_speed_naive) + +Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed +is estimated as +`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, +defaulting to [`max_abs_speed_naive`](@ref). +""" +struct DissipationEntropyStableLump{MaxAbsSpeed} + max_abs_speed::MaxAbsSpeed +end + +DissipationEntropyStableLump() = DissipationEntropyStableLump(max_abs_speed_naive) + +@inline function (dissipation::DissipationEntropyStableLump)(u_ll, u_rr, + orientation_or_normal_direction, + equations::IdealMhdMultiIonEquations2D) + @unpack gammas = equations + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + + w_ll = cons2entropy(u_ll, equations) + w_rr = cons2entropy(u_rr, equations) + prim_ll = cons2prim(u_ll, equations) + prim_rr = cons2prim(u_rr, equations) + B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + # Some global averages + B1_avg = 0.5 * (B1_ll + B1_rr) + B2_avg = 0.5 * (B2_ll + B2_rr) + B3_avg = 0.5 * (B3_ll + B3_rr) + psi_avg = 0.5 * (psi_ll + psi_rr) + + dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) + + beta_plus_ll = 0 + beta_plus_rr = 0 + # Get the lumped dissipation for all components + for k in eachcomponent(equations) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll = get_component(k, prim_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr = get_component(k, prim_rr, equations) + + w1_ll, w2_ll, w3_ll, w4_ll, w5_ll = get_component(k, w_ll, equations) + w1_rr, w2_rr, w3_rr, w4_rr, w5_rr = get_component(k, w_rr, equations) + + # Auxiliary variables + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 + vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 + + # Mean variables + rho_ln = ln_mean(rho_ll, rho_rr) + beta_ln = ln_mean(beta_ll, beta_rr) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) + p_mean = 0.5 * rho_avg / beta_avg + p_star = 0.5 * rho_ln / beta_ln + vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) + vel_avg_norm = v1_avg^2 + v2_avg^2 + v3_avg^2 + E_bar = p_star / (gammas[k] - 1) + 0.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) + + h11 = rho_ln + h22 = rho_ln * v1_avg^2 + p_mean + h33 = rho_ln * v2_avg^2 + p_mean + h44 = rho_ln * v3_avg^2 + p_mean + h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln + + vel_avg_norm * p_mean + + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) + + d1 = -0.5 * λ * h11 * (w1_rr - w1_ll) + d2 = -0.5 * λ * h22 * (w2_rr - w2_ll) + d3 = -0.5 * λ * h33 * (w3_rr - w3_ll) + d4 = -0.5 * λ * h44 * (w4_rr - w4_ll) + d5 = -0.5 * λ * h55 * (w5_rr - w5_ll) + + beta_plus_ll += beta_ll + beta_plus_rr += beta_rr + + set_component!(dissipation, k, d1, d2, d3, d4, d5, equations) + end + + # Set the magnetic field and psi terms + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) + dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2]) + dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3]) + dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end]) + + return dissipation +end + +function Base.show(io::IO, d::DissipationEntropyStableLump) + print(io, "DissipationEntropyStableLump(", d.max_abs_speed, ")") +end + +""" +DissipationEntropyStable(max_abs_speed=max_abs_speed_naive) + +Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed +is estimated as +`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`, +defaulting to [`max_abs_speed_naive`](@ref). +""" +struct DissipationEntropyStable{MaxAbsSpeed} + max_abs_speed::MaxAbsSpeed +end + +DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive) + +@inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr, + orientation_or_normal_direction, + equations::IdealMhdMultiIonEquations2D) + @unpack gammas = equations + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + + w_ll = cons2entropy(u_ll, equations) + w_rr = cons2entropy(u_rr, equations) + prim_ll = cons2prim(u_ll, equations) + prim_rr = cons2prim(u_rr, equations) + B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + # Some global averages + B1_avg = 0.5 * (B1_ll + B1_rr) + B2_avg = 0.5 * (B2_ll + B2_rr) + B3_avg = 0.5 * (B3_ll + B3_rr) + psi_avg = 0.5 * (psi_ll + psi_rr) + + dissipation = zero(MVector{nvariables(equations), eltype(u_ll)}) + + beta_plus_ll = 0 + beta_plus_rr = 0 + # Get the lumped dissipation for all components + for k in eachcomponent(equations) + rho_ll, v1_ll, v2_ll, v3_ll, p_ll = get_component(k, prim_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr = get_component(k, prim_rr, equations) + + w1_ll, w2_ll, w3_ll, w4_ll, w5_ll = get_component(k, w_ll, equations) + w1_rr, w2_rr, w3_rr, w4_rr, w5_rr = get_component(k, w_rr, equations) + + # Auxiliary variables + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 + vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 + + # Mean variables + rho_ln = ln_mean(rho_ll, rho_rr) + beta_ln = ln_mean(beta_ll, beta_rr) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + tau = 1 / (beta_ll + beta_rr) + p_mean = 0.5 * rho_avg / beta_avg + p_star = 0.5 * rho_ln / beta_ln + vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) + vel_avg_norm = v1_avg^2 + v2_avg^2 + v3_avg^2 + E_bar = p_star / (gammas[k] - 1) + 0.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg) + + h11 = rho_ln + h12 = rho_ln * v1_avg + h13 = rho_ln * v2_avg + h14 = rho_ln * v3_avg + h15 = E_bar + d1 = -0.5 * λ * (h11 * (w1_rr - w1_ll) + + h12 * (w2_rr - w2_ll) + + h13 * (w3_rr - w3_ll) + + h14 * (w4_rr - w4_ll) + + h15 * (w5_rr - w5_ll) ) + + h21 = h12 + h22 = rho_ln * v1_avg^2 + p_mean + h23 = h21 * v2_avg + h24 = h21 * v3_avg + h25 = (E_bar + p_mean) * v1_avg + d2 = -0.5 * λ * (h21 * (w1_rr - w1_ll) + + h22 * (w2_rr - w2_ll) + + h23 * (w3_rr - w3_ll) + + h24 * (w4_rr - w4_ll) + + h25 * (w5_rr - w5_ll) ) + + h31 = h13 + h32 = h23 + h33 = rho_ln * v2_avg^2 + p_mean + h34 = h31 * v3_avg + h35 = (E_bar + p_mean) * v2_avg + d3 = -0.5 * λ * (h31 * (w1_rr - w1_ll) + + h32 * (w2_rr - w2_ll) + + h33 * (w3_rr - w3_ll) + + h34 * (w4_rr - w4_ll) + + h35 * (w5_rr - w5_ll) ) + + h41 = h14 + h42 = h24 + h43 = h34 + h44 = rho_ln * v3_avg^2 + p_mean + h45 = (E_bar + p_mean) * v3_avg + d4 = -0.5 * λ * (h41 * (w1_rr - w1_ll) + + h42 * (w2_rr - w2_ll) + + h43 * (w3_rr - w3_ll) + + h44 * (w4_rr - w4_ll) + + h45 * (w5_rr - w5_ll) ) + + h51 = h15 + h52 = h25 + h53 = h35 + h54 = h45 + h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln + + vel_avg_norm * p_mean + + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 )) + d5 = -0.5 * λ * (h51 * (w1_rr - w1_ll) + + h52 * (w2_rr - w2_ll) + + h53 * (w3_rr - w3_ll) + + h54 * (w4_rr - w4_ll) + + h55 * (w5_rr - w5_ll) ) + + beta_plus_ll += beta_ll + beta_plus_rr += beta_rr + + set_component!(dissipation, k, d1, d2, d3, d4, d5, equations) + end + + # Set the magnetic field and psi terms + h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) + + # diagonal entries + dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1]) + dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2]) + dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3]) + dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end]) + # Off-diagonal entries + for k in eachcomponent(equations) + _, _, _, _, w5_ll = get_component(k, w_ll, equations) + _, _, _, _, w5_rr = get_component(k, w_rr, equations) + + dissipation[1] -= 0.5 * λ * h_B_psi * B1_avg * (w5_rr - w5_ll) + dissipation[2] -= 0.5 * λ * h_B_psi * B2_avg * (w5_rr - w5_ll) + dissipation[3] -= 0.5 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll) + dissipation[end] -= 0.5 * λ * h_B_psi * psi_avg * (w5_rr - w5_ll) + + ind_E = 3 + (k - 1) * 5 + 5 + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3]) + dissipation[ind_E] -= 0.5 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end]) + end + + return dissipation +end + +function Base.show(io::IO, d::DissipationEntropyStable) + print(io, "DissipationEntropyStable(", d.max_abs_speed, ")") +end + end # @muladd