Skip to content

Commit ff6ba9b

Browse files
SimonCanDanielDoehringranocha
authored
Added abstract vector for max_abs_speed_naive and CompressibleEulerMulticomponentEquations2D. (#2607)
* Added n::AbstractVector version of max_abs_speed_naive for CompressibleEulerMulticomponentEquations2D. * Renamed normal vector to be more consistent with the rest of the code. * Added test for this PR. * Added elixir for 2d multicomponent Euler using p4est meshes. * Update examples/p4est_2d_dgsem/elixir_eulermulti_convergence_ec.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Changed the surface flux to flux_lax_friedrichs. This fixed the bug that the new max_abs_speed was not being used. * Replaced example elixir with shock. * Replaced test for p4est euler multi. * Update src/equations/compressible_euler_multicomponent_2d.jl * Apply suggestions from code review * Update src/equations/compressible_euler_multicomponent_2d.jl * Apply suggestions from code review * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl Co-authored-by: Hendrik Ranocha <[email protected]> * More explicit types to avoid tpye instabilities. * Corrected test values. * Increased the memory for the multicomponent Euler test. * Applied auto formatter. * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/equations/compressible_euler_multicomponent_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Added unit test for min_max_speed_naive. --------- Co-authored-by: Daniel Doehring <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent 2564ef7 commit ff6ba9b

File tree

4 files changed

+188
-0
lines changed

4 files changed

+188
-0
lines changed

‎examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
# 1) Dry Air 2) SF_6
5+
equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
6+
gas_constants = (0.287, 1.578))
7+
8+
"""
9+
initial_condition_shock(coordinates, t, equations::CompressibleEulerEquations2D)
10+
11+
Shock traveling from left to right where it interacts with a Perturbed interface.
12+
"""
13+
@inline function initial_condition_shock(x, t,
14+
equations::CompressibleEulerMulticomponentEquations2D)
15+
rho_0 = 1.25 # kg/m^3
16+
p_0 = 101325 # Pa
17+
T_0 = 293 # K
18+
u_0 = 352 # m/s
19+
d = 20
20+
w = 40
21+
22+
if x[1] < 25
23+
# Shock region.
24+
v1 = 0.35
25+
v2 = 0.0
26+
rho1 = 1.72
27+
rho2 = 0.03
28+
p = 1.57
29+
elseif (x[1] <= 30) || (x[1] <= 70 && abs(125 - x[2]) > w / 2)
30+
# Intermediate region.
31+
v1 = 0.0
32+
v2 = 0.0
33+
rho1 = 1.25
34+
rho2 = 0.03
35+
p = 1.0
36+
else
37+
(x[1] <= 70 + d)
38+
# SF_6 region.
39+
v1 = 0.0
40+
v2 = 0.0
41+
rho1 = 0.03
42+
rho2 = 6.03 #SF_6
43+
p = 1.0
44+
end
45+
46+
return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
47+
end
48+
49+
# Define the simulation.
50+
initial_condition = initial_condition_shock
51+
52+
surface_flux = flux_lax_friedrichs
53+
volume_flux = flux_ranocha
54+
basis = LobattoLegendreBasis(3)
55+
56+
limiter_idp = SubcellLimiterIDP(equations, basis;
57+
positivity_variables_cons = ["rho" * string(i)
58+
for i in eachcomponent(equations)])
59+
60+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
61+
volume_flux_dg = volume_flux,
62+
volume_flux_fv = surface_flux)
63+
64+
solver = DGSEM(basis, surface_flux, volume_integral)
65+
66+
coordinates_min = (0.0, 0.0)
67+
coordinates_max = (250.0, 250.0)
68+
mesh = P4estMesh((32, 32), polydeg = 3, coordinates_min = (0.0, 0.0),
69+
coordinates_max = (250.0, 250.0), periodicity = (false, true),
70+
initial_refinement_level = 0)
71+
72+
# Completely free outflow
73+
function boundary_condition_outflow(u_inner, normal_direction::AbstractVector,
74+
x, t,
75+
surface_flux_function,
76+
equations)
77+
# Calculate the boundary flux entirely from the internal solution state
78+
return flux(u_inner, normal_direction, equations)
79+
end
80+
81+
boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
82+
:x_pos => boundary_condition_outflow)
83+
84+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
85+
boundary_conditions = boundary_conditions)
86+
87+
###############################################################################
88+
# ODE solvers, callbacks etc.
89+
90+
tspan = (0.0, 5.0)
91+
ode = semidiscretize(semi, tspan)
92+
93+
summary_callback = SummaryCallback()
94+
95+
analysis_interval = 10
96+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
97+
save_analysis = true)
98+
99+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
100+
101+
save_solution = SaveSolutionCallback(dt = 2.0,
102+
save_initial_solution = true,
103+
save_final_solution = true,
104+
solution_variables = cons2prim)
105+
106+
stepsize_callback = StepsizeCallback(cfl = 0.2)
107+
108+
callbacks = CallbackSet(summary_callback,
109+
analysis_callback,
110+
alive_callback,
111+
save_solution,
112+
stepsize_callback)
113+
114+
###############################################################################
115+
# run the simulation
116+
117+
stage_callbacks = (SubcellLimiterIDPCorrection(),
118+
BoundsCheckCallback(save_errors = false, interval = 100))
119+
# `interval` is used when calling this elixir in the tests with `save_errors=true`.
120+
121+
@time begin
122+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
123+
dt = 0.1, # solve needs some value here but it will be overwritten by the stepsize_callback
124+
ode_default_options()..., callback = callbacks)
125+
end

‎src/equations/compressible_euler_multicomponent_2d.jl‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -689,6 +689,43 @@ end
689689
return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
690690
end
691691

692+
# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
693+
@inline function max_abs_speed_naive(u_ll::SVector{N, T}, u_rr::SVector{N, T},
694+
normal_direction::SVector{2, T},
695+
equations::CompressibleEulerMulticomponentEquations2D) where {
696+
N,
697+
T <:
698+
AbstractFloat
699+
}
700+
# Unpack conservative variables
701+
rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
702+
rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
703+
704+
# Get densities and gammas
705+
rho_ll = T(density(u_ll, equations))
706+
rho_rr = T(density(u_rr, equations))
707+
gamma_ll = T(totalgamma(u_ll, equations))
708+
gamma_rr = T(totalgamma(u_rr, equations))
709+
710+
# Velocity components
711+
v_ll_vec = SVector(rho_v1_ll / rho_ll, rho_v2_ll / rho_ll)
712+
v_rr_vec = SVector(rho_v1_rr / rho_rr, rho_v2_rr / rho_rr)
713+
714+
# Project velocities onto the direction normal_direction.
715+
v_ll = dot(v_ll_vec, normal_direction)
716+
v_rr = dot(v_rr_vec, normal_direction)
717+
718+
# Compute pressures
719+
p_ll = (gamma_ll - one(T)) * (rho_e_ll - T(0.5) * dot(v_ll_vec, v_ll_vec) * rho_ll)
720+
p_rr = (gamma_rr - one(T)) * (rho_e_rr - T(0.5) * dot(v_rr_vec, v_rr_vec) * rho_rr)
721+
722+
# Sound speeds
723+
c_ll = sqrt(gamma_ll * p_ll / rho_ll)
724+
c_rr = sqrt(gamma_rr * p_rr / rho_rr)
725+
726+
return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
727+
end
728+
692729
# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
693730
@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
694731
equations::CompressibleEulerMulticomponentEquations2D)

‎test/test_p4est_2d.jl‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -560,6 +560,27 @@ end
560560
@test_allocations(Trixi.rhs!, semi, sol, 1000)
561561
end
562562

563+
@trixi_testset "elixir_eulermulti_shock.jl" begin
564+
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock.jl"),
565+
l2=[
566+
0.09238286359859513,
567+
0.0006429842672568709,
568+
0.22741647050108965,
569+
0.09129657317724366,
570+
0.2714793219290268
571+
],
572+
linf=[
573+
0.6949550920191041,
574+
0.027123073012782554,
575+
1.7445063388384803,
576+
1.1616284674760593,
577+
4.340140043930651
578+
])
579+
# Ensure that we do not have excessive memory allocations
580+
# (e.g., from type instabilities)
581+
@test_allocations(Trixi.rhs!, semi, sol, 15000)
582+
end
583+
563584
@trixi_testset "elixir_mhd_alfven_wave.jl" begin
564585
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
565586
l2=[1.0513414461545583e-5, 1.0517900957166411e-6,

‎test/test_unit.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1978,6 +1978,11 @@ end
19781978
@test max_abs_speed_naive(u_ll, u_rr, orientation, equations) ≈
19791979
max_abs_speed(u_ll, u_rr, orientation, equations)
19801980
end
1981+
1982+
@test max_abs_speed_naive(u_ll, u_rr, 1, equations) ≈
1983+
max_abs_speed_naive(u_ll, u_rr, SVector(1.0, 0.0), equations)
1984+
@test max_abs_speed_naive(u_ll, u_rr, 2, equations) ≈
1985+
max_abs_speed_naive(u_ll, u_rr, SVector(0.0, 1.0), equations)
19811986
end
19821987
end
19831988

0 commit comments

Comments
 (0)