-
Notifications
You must be signed in to change notification settings - Fork 0
/
elixir_euler_astro_jet.jl
129 lines (103 loc) · 5.29 KB
/
elixir_euler_astro_jet.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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
using OrdinaryDiffEq
using Trixi
###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations = CompressibleEulerEquations2D(gamma)
# Initial condition adopted from
# - Yong Liu, Jianfang Lu, and Chi-Wang Shu
# An oscillation free discontinuous Galerkin method for hyperbolic systems
# Preprint available at https://tinyurl.com/c76fjtx4
# Mach = 2000 jet
function initial_condition_astro_jet(x, t, equations::CompressibleEulerEquations2D)
@unpack gamma = equations
rho = 0.5
v1 = 0.0
v2 = 0.0
p = 0.4127
# add inflow for t>0 at x=-0.5
# domain size is [-0.5,+0.5]^2
if (t > 0) && (x[1] ≈ -0.5) && (abs(x[2]) < 0.05)
rho = 5.0
v1 = 800 # about Mach number Ma = 2000
v2 = 0.0
p = 0.4127
end
return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_astro_jet
boundary_conditions = (x_neg=BoundaryConditionDirichlet(initial_condition_astro_jet),
x_pos=BoundaryConditionDirichlet(initial_condition_astro_jet),
y_neg=boundary_condition_periodic,
y_pos=boundary_condition_periodic)
surface_flux = flux_lax_friedrichs # HLLC needs more shock capturing (alpha_max)
volume_flux = flux_ranocha # works with Chandrashekar flux as well
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
# shock capturing necessary for this tough example
indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max=0.3,
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 = (-0.5, -0.5)
coordinates_max = ( 0.5, 0.5)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level=6,
periodicity=(false,true),
n_cells_max=100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
###############################################################################
# ODE solvers, callbacks etc.
tspan = (0.0, 0.001)
ode = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()
analysis_interval = 5000
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 = ControllerThreeLevelCombined(semi, amr_indicator, indicator_sc,
base_level=2,
med_level =0, med_threshold=0.0003, # med_level = current level
max_level =8, max_threshold=0.003,
max_threshold_secondary=indicator_sc.alpha_max)
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)
# positivity limiter necessary for this tough example
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
variables=(Trixi.density, pressure))
###############################################################################
# run the simulation
# use adaptive time stepping based on error estimates, time step roughly dt = 1e-7
sol = solve(ode, SSPRK43(stage_limiter!),
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)
# `colorbar_scale=:log10` does not seem to be supported by the default backend GR.
# It is even not properly supported by the PyPlot backend via Plots. Thus, we
# just plot the `log10` of the data and change it manually later...
# The first index of `pd.data` is the density rho.
heatmap(pd.x, pd.y, log10.(pd.data[1]), xguide="\$x\$", yguide="\$y\$",
xlim=extrema(pd.x), ylims=extrema(pd.y), legend=:none, aspect_ratio=:equal,
colorbar=true, size=(700, 500); fontsizes...)
savefig(joinpath(figdir, "jet_density.pdf"))
plot(getmesh(pd), xlabel="\$x\$", ylabel="\$y\$", size=(600, 500), linewidth=2, grid=false; fontsizes...)
savefig(joinpath(figdir, "jet_mesh.pdf"))