-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathelixir_euler_kelvin_helmholtz_instability_amr.jl
109 lines (85 loc) · 4.12 KB
/
elixir_euler_kelvin_helmholtz_instability_amr.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
using OrdinaryDiffEq
using Trixi
###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)
"""
initial_condition(x, t, equations::CompressibleEulerEquations2D)
A version of the classical Kelvin-Helmholtz instability based on
- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
of the Euler Equations
[arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
"""
function initial_condition(x, t, equations::CompressibleEulerEquations2D)
# change discontinuity to tanh
# typical resolution 128^2, 256^2
# domain size is [-1,+1]^2
slope = 15
amplitude = 0.02
B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
rho = 0.5 + 0.75 * B
v1 = 0.5 * (B - 1)
v2 = 0.1 * sin(2 * pi * x[1])
p = 1.0
return prim2cons(SVector(rho, v1, v2, p), equations)
end
surface_flux = flux_lax_friedrichs
volume_flux = flux_chandrashekar
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max=0.002,
alpha_min=0.0001,
alpha_smooth=true,
variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg=volume_flux,
volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)
coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level=5,
n_cells_max=100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
###############################################################################
# ODE solvers, callbacks etc.
tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval=analysis_interval)
amr_indicator = IndicatorHennemannGassner(semi,
alpha_max=1.0,
alpha_min=0.0001,
alpha_smooth=false,
variable=Trixi.density)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
base_level=4,
med_level=0, med_threshold=0.0003, # med_level = current level
max_level=6, max_threshold=0.003)
amr_callback = AMRCallback(semi, amr_controller,
interval=1,
adapt_initial_condition=true,
adapt_initial_condition_only_refine=true)
callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
amr_callback)
###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(),
save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary
###############################################################################
# plot the numerical result at the final time
using Plots
figdir = joinpath(dirname(@__DIR__), "figures")
fontsizes = (xtickfontsize=18, ytickfontsize=18, xguidefontsize=20, yguidefontsize=20, legendfontsize=18)
pd = PlotData2D(sol)
plot(pd["rho"], title="", colorbar_title="", size=(650, 500), clim=(0.47, 2.3); fontsizes...)
savefig(joinpath(figdir, "kelvin_helmholtz_density_t" * string(round(Int, last(tspan))) * ".pdf"))
plot(getmesh(pd), xlabel="\$x\$", ylabel="\$y\$", size=(580, 500), linewidth=2, grid=false; fontsizes...)
savefig(joinpath(figdir, "kelvin_helmholtz_mesh_t" * string(round(Int, last(tspan))) * ".pdf"))