Skip to content

Commit b9376e8

Browse files
DanielDoehringDanielDoehring
authored
Close-up Mach 3 Cylinder Bow-Shock Elixir (#2503)
* Close-up mach 3 cylinder bow-shock elixir * Update examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl * comments * comments * commnets
1 parent 13fd168 commit b9376e8

File tree

2 files changed

+215
-0
lines changed

2 files changed

+215
-0
lines changed

‎examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl‎

Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
1+
using Trixi
2+
using OrdinaryDiffEqSSPRK
3+
4+
###############################################################################
5+
# Geometry & boundary conditions
6+
7+
# Mapping to create a "close-up" mesh around the second quadrant of a cylinder,
8+
# implemented by Georgii Oblapenko. If you use this in your own work, please cite:
9+
#
10+
# - G. Oblapenko and A. Tarnovskiy (2024)
11+
# Reproducibility Repository for the paper:
12+
# Entropy-stable fluxes for high-order Discontinuous Galerkin simulations of high-enthalpy flows.
13+
# [DOI: 10.5281/zenodo.13981615](https://doi.org/10.5281/zenodo.13981615)
14+
# [GitHub](https://github.com/knstmrd/paper_ec_trixi_chem)
15+
#
16+
# as well as the corresponding paper:
17+
# - G. Oblapenko and M. Torrilhon (2025)
18+
# Entropy-conservative high-order methods for high-enthalpy gas flows.
19+
# Computers & Fluids, 2025.
20+
# [DOI: 10.1016/j.compfluid.2025.106640](https://doi.org/10.1016/j.compfluid.2025.106640)
21+
#
22+
# The mapping produces the following geometry & shock (indicated by the asterisks `* `):
23+
# ____x_neg____
24+
# | |
25+
# | |
26+
# | |
27+
# | * |
28+
# | * y
29+
# | Inflow * _
30+
# | state * p
31+
# x * o
32+
# _ * s
33+
# n * |
34+
# e * |
35+
# g Shock .
36+
# | * .
37+
# | * . <- x_pos
38+
# | * .
39+
# | * . (Cylinder)
40+
# |_______y_neg_______.
41+
function mapping_cylinder_shock_fitted(xi_, eta_,
42+
cylinder_radius, spline_points)
43+
shock_shape = [
44+
(spline_points[1], 0.0), # Shock position on the stagnation line (`y_neg`, y = 0)
45+
(spline_points[2], spline_points[2]), # Shock position at -45° angle
46+
(0.0, spline_points[3]) # Shock position at outflow (`y_pos`, x = x_max)
47+
] # 3 points that define the geometry of the mesh which follows the shape of the shock (known a-priori)
48+
R = [sqrt(shock_shape[i][1]^2 + shock_shape[i][2]^2) for i in 1:3] # 3 radii
49+
50+
# Construct spline with form R[1] + c2 * eta_01^2 + c3 * eta_01^3,
51+
# chosen such that derivative w.r.t eta_01 is 0 at eta_01 = 0 such that
52+
# we have symmetry along the stagnation line (`y_neg`, y = 0).
53+
#
54+
# A single cubic spline doesn't fit the shock perfectly,
55+
# but is the simplest curve that does a reasonable job and it also can be easily computed analytically.
56+
# The choice of points on the stagnation line and outflow region is somewhat self-evident
57+
# (capture the minimum and maximum extent of the shock stand-off),
58+
# and the point at the 45 degree angle seemed the most logical to add
59+
# since it only requires one additional value (and not two),
60+
# simplifies the math a bit, and the angle lies exactly in between the other angles.
61+
spline_matrix = [1.0 1.0; 0.25 0.125]
62+
spline_RHS = [R[3] - R[1], R[2] - R[1]]
63+
spline_coeffs = spline_matrix \ spline_RHS # c2, c3
64+
65+
eta_01 = (eta_ + 1) / 2 # Transform `eta_` in [-1, 1] to `eta_01` in [0, 1]
66+
# "Flip" `xi_` in [-1, 1] to `xi_01` in [0, 1] since
67+
# shock positions where originally for first quadrant, here we use second quadrant
68+
xi_01 = (-xi_ + 1) / 2
69+
70+
R_outer = R[1] + spline_coeffs[1] * eta_01^2 + spline_coeffs[2] * eta_01^3
71+
72+
angle = -π / 4 + eta_ * π / 4 # Angle runs from -90° to 0°
73+
r = (cylinder_radius + xi_01 * (R_outer - cylinder_radius))
74+
75+
return SVector(round(r * sin(angle); digits = 8), round(r * cos(angle); digits = 8))
76+
end
77+
78+
@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D)
79+
# set the freestream flow parameters
80+
rho_freestream = equations.gamma
81+
v1 = 3.0 # => Mach 3 for unity speed of sound
82+
v2 = 0
83+
p_freestream = 1
84+
prim = SVector(rho_freestream, v1, v2, p_freestream)
85+
return prim2cons(prim, equations)
86+
end
87+
88+
@inline function boundary_condition_supersonic_inflow(u_inner,
89+
normal_direction::AbstractVector,
90+
x, t,
91+
surface_flux_function,
92+
equations::CompressibleEulerEquations2D)
93+
u_boundary = initial_condition_mach3_flow(x, t, equations)
94+
flux = Trixi.flux(u_boundary, normal_direction, equations)
95+
96+
return flux
97+
end
98+
99+
# For physical significance of boundary conditions, see sketch at `mapping_cylinder_shock_fitted`
100+
boundary_conditions = Dict(:x_neg => boundary_condition_supersonic_inflow, # Supersonic inflow
101+
:y_neg => boundary_condition_slip_wall, # Induce symmetry by slip wall
102+
:y_pos => boundary_condition_do_nothing, # Free outflow
103+
:x_pos => boundary_condition_slip_wall) # Cylinder
104+
105+
###############################################################################
106+
# Equations, mesh and solver
107+
108+
gamma = 1.4
109+
equations = CompressibleEulerEquations2D(gamma)
110+
111+
polydeg = 3
112+
basis = LobattoLegendreBasis(polydeg)
113+
114+
surface_flux = flux_hllc
115+
volume_flux = flux_ranocha
116+
117+
shock_indicator = IndicatorHennemannGassner(equations, basis,
118+
alpha_max = 0.5,
119+
alpha_min = 0.001,
120+
alpha_smooth = true,
121+
variable = density_pressure)
122+
volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
123+
volume_flux_dg = volume_flux,
124+
volume_flux_fv = surface_flux)
125+
126+
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
127+
volume_integral = volume_integral)
128+
129+
trees_per_dimension = (25, 25)
130+
131+
cylinder_radius = 0.5
132+
# Follow from a-priori known shock shape, originally for first qaudrant,
133+
# here transformed to second quadrant, see `mapping_cylinder_shock_fitted`.
134+
spline_points = [1.32, 1.05, 2.25]
135+
cylinder_mapping = (xi, eta) -> mapping_cylinder_shock_fitted(xi, eta,
136+
cylinder_radius,
137+
spline_points)
138+
139+
mesh = P4estMesh(trees_per_dimension,
140+
polydeg = polydeg,
141+
mapping = cylinder_mapping,
142+
periodicity = false)
143+
144+
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
145+
volume_integral = volume_integral)
146+
147+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_mach3_flow,
148+
solver,
149+
boundary_conditions = boundary_conditions)
150+
151+
###############################################################################
152+
# Semidiscretization & callbacks
153+
154+
tspan = (0.0, 5.0) # More or less stationary shock position reached
155+
ode = semidiscretize(semi, tspan)
156+
157+
summary_callback = SummaryCallback()
158+
159+
analysis_callback = AnalysisCallback(semi, interval = 5000)
160+
alive_callback = AliveCallback(alive_interval = 200)
161+
162+
save_solution = SaveSolutionCallback(dt = 0.25,
163+
save_initial_solution = true,
164+
save_final_solution = true,
165+
solution_variables = cons2prim)
166+
167+
amr_controller = ControllerThreeLevel(semi, shock_indicator;
168+
base_level = 0,
169+
med_level = 1, med_threshold = 0.175,
170+
max_level = 3, max_threshold = 0.35)
171+
172+
amr_callback = AMRCallback(semi, amr_controller,
173+
interval = 25,
174+
adapt_initial_condition = true,
175+
adapt_initial_condition_only_refine = true)
176+
177+
callbacks = CallbackSet(summary_callback,
178+
analysis_callback, alive_callback,
179+
save_solution, amr_callback)
180+
181+
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
182+
variables = (Trixi.density, pressure))
183+
184+
###############################################################################
185+
# Run the simulation
186+
187+
sol = solve(ode, SSPRK33(stage_limiter! = stage_limiter!, thread = Trixi.True());
188+
dt = 1.6e-5, # Fixed timestep works decent here
189+
ode_default_options()..., callback = callbacks);

‎test/test_p4est_2d.jl‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1016,6 +1016,32 @@ end
10161016
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
10171017
end
10181018
end
1019+
1020+
@trixi_testset "elixir_euler_cylinder_bowshock_mach3.jl" begin
1021+
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
1022+
"elixir_euler_cylinder_bowshock_mach3.jl"),
1023+
tspan=(0.0, 1e-3),
1024+
l2=[
1025+
0.03787745781612722,
1026+
0.03339276348608649,
1027+
0.05301001151898993,
1028+
0.2868802674001281
1029+
],
1030+
linf=[
1031+
2.5347156069842978,
1032+
2.6657123832452414,
1033+
3.786891603220761,
1034+
21.305497055838977
1035+
])
1036+
# Ensure that we do not have excessive memory allocations
1037+
# (e.g., from type instabilities)
1038+
let
1039+
t = sol.t[end]
1040+
u_ode = sol.u[end]
1041+
du_ode = similar(u_ode)
1042+
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
1043+
end
1044+
end
10191045
end
10201046

10211047
# Clean up afterwards: delete Trixi.jl output directory

0 commit comments

Comments
 (0)