Skip to content

Commit c4c6f41

Browse files
bennibolmbennibolm
committed
Initial implementation of 3d idp mortars
1 parent 9bffb3f commit c4c6f41

File tree

7 files changed

+1271
-4
lines changed

7 files changed

+1271
-4
lines changed

‎examples/tree_3d_dgsem/elixir_euler_amr_sc_subcell.jl‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
using Trixi
2+
3+
###############################################################################
4+
# semidiscretization of the compressible Euler equations
5+
6+
equations = CompressibleEulerEquations3D(1.4)
7+
8+
"""
9+
initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D)
10+
11+
A Gaussian pulse in the density with constant velocity and pressure; reduces the
12+
compressible Euler equations to the linear advection equations.
13+
"""
14+
function initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D)
15+
# rho = 1 + exp(-(x[1]^2 + x[2]^2 + x[3]^2)) / 2
16+
rho = 1 + exp(-(x[1]^2 + x[2]^2)) / 2
17+
v1 = 1
18+
v2 = 1
19+
v3 = 0.0#1
20+
rho_v1 = rho * v1
21+
rho_v2 = rho * v2
22+
rho_v3 = 0.0#rho * v3
23+
p = 1
24+
rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2 + v3^2)
25+
return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e)
26+
end
27+
initial_condition = initial_condition_density_pulse
28+
29+
surface_flux = flux_lax_friedrichs
30+
volume_flux = flux_ranocha
31+
polydeg = 3
32+
basis = LobattoLegendreBasis(polydeg)
33+
limiter_idp = SubcellLimiterIDP(equations, basis;
34+
positivity_variables_cons = [],
35+
positivity_variables_nonlinear = [],
36+
local_twosided_variables_cons = []) # required for testing
37+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
38+
volume_flux_dg = volume_flux,
39+
volume_flux_fv = surface_flux)
40+
mortar = MortarIDP(equations, basis;
41+
# basis_function = :piecewise_constant,
42+
positivity_variables_cons = ["rho"],
43+
positivity_variables_nonlinear = [pressure],
44+
pure_low_order = false)
45+
solver = DGSEM(basis, surface_flux, volume_integral, mortar)
46+
47+
coordinates_min = (-5.0, -5.0, -5.0)
48+
coordinates_max = (5.0, 5.0, 5.0)
49+
refinement_patches = ((type = "box", coordinates_min = (0.0, -1.0, -1.0),
50+
coordinates_max = (1.0, 1.0, 1.0)),
51+
(type = "box", coordinates_min = (0.0, -0.5, -0.5),
52+
coordinates_max = (0.5, 0.5, 0.5)))
53+
mesh = TreeMesh(coordinates_min, coordinates_max,
54+
initial_refinement_level = 4,
55+
refinement_patches = refinement_patches,
56+
n_cells_max = 100_000)
57+
58+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
59+
60+
###############################################################################
61+
# ODE solvers, callbacks etc.
62+
63+
tspan = (0.0, 1.0)
64+
ode = semidiscretize(semi, tspan)
65+
66+
summary_callback = SummaryCallback()
67+
68+
analysis_interval = 100
69+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
70+
extra_analysis_integrals = (entropy,))
71+
72+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
73+
74+
save_restart = SaveRestartCallback(interval = 100,
75+
save_final_restart = true)
76+
77+
save_solution = SaveSolutionCallback(interval = 100,
78+
save_initial_solution = true,
79+
save_final_solution = true,
80+
solution_variables = cons2prim)
81+
82+
amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
83+
base_level = 3,
84+
med_level = 4, med_threshold = 1.05,
85+
max_level = 5, max_threshold = 1.3)
86+
amr_callback = AMRCallback(semi, amr_controller,
87+
interval = 5,
88+
adapt_initial_condition = true,
89+
adapt_initial_condition_only_refine = true)
90+
91+
stepsize_callback = StepsizeCallback(cfl = 0.9)
92+
93+
callbacks = CallbackSet(summary_callback,
94+
analysis_callback, alive_callback,
95+
# amr_callback,
96+
save_solution,
97+
stepsize_callback);
98+
99+
###############################################################################
100+
# run the simulation
101+
102+
stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())
103+
104+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
105+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
106+
callback = callbacks);

‎src/callbacks_stage/subcell_limiter_idp_correction_3d.jl‎

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,4 +84,128 @@ function perform_idp_correction!(u, dt,
8484

8585
return nothing
8686
end
87+
88+
function perform_idp_mortar_correction(u, dt, mesh::TreeMesh{3}, equations, dg, cache)
89+
(; orientations, limiting_factor) = cache.mortars
90+
91+
(; surface_flux_values) = cache.elements
92+
(; surface_flux_values_high_order) = cache.antidiffusive_fluxes
93+
(; boundary_interpolation) = dg.basis
94+
95+
for mortar in eachmortar(dg, cache)
96+
if isapprox(limiting_factor[mortar], one(eltype(limiting_factor)))
97+
continue
98+
end
99+
100+
if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
101+
if orientations[mortar] == 1
102+
direction_small = 1
103+
direction_large = 2
104+
elseif orientations[mortar] == 2
105+
direction_small = 3
106+
direction_large = 4
107+
else # orientations[mortar] == 3
108+
direction_small = 5
109+
direction_large = 6
110+
end
111+
# In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
112+
# This sign switch is directly applied to the boundary interpolation factors here.
113+
factor_small = boundary_interpolation[1, 1]
114+
factor_large = -boundary_interpolation[nnodes(dg), 2]
115+
else # large_sides[mortar] == 2 -> small elements on left side
116+
if orientations[mortar] == 1
117+
direction_small = 2
118+
direction_large = 1
119+
elseif orientations[mortar] == 2
120+
direction_small = 4
121+
direction_large = 3
122+
else # orientations[mortar] == 3
123+
direction_small = 6
124+
direction_large = 5
125+
end
126+
# In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
127+
# This sign switch is directly applied to the boundary interpolation factors here.
128+
factor_large = boundary_interpolation[1, 1]
129+
factor_small = -boundary_interpolation[nnodes(dg), 2]
130+
end
131+
132+
for j in eachnode(dg), i in eachnode(dg)
133+
if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
134+
if orientations[mortar] == 1
135+
# L2 mortars in x-direction
136+
indices_small = (1, i, j)
137+
indices_large = (nnodes(dg), i, j)
138+
elseif orientations[mortar] == 2
139+
# L2 mortars in y-direction
140+
indices_small = (i, 1, j)
141+
indices_large = (i, nnodes(dg), j)
142+
else # orientations[mortar] == 3
143+
# L2 mortars in z-direction
144+
indices_small = (i, j, 1)
145+
indices_large = (i, j, nnodes(dg))
146+
end
147+
else # large_sides[mortar] == 2 -> small elements on left side
148+
if orientations[mortar] == 1
149+
# L2 mortars in x-direction
150+
indices_small = (nnodes(dg), i, j)
151+
indices_large = (1, i, j)
152+
elseif orientations[mortar] == 2
153+
# L2 mortars in y-direction
154+
indices_small = (i, nnodes(dg), j)
155+
indices_large = (i, 1, j)
156+
else # orientations[mortar] == 3
157+
# L2 mortars in z-direction
158+
indices_small = (i, j, nnodes(dg))
159+
indices_large = (i, j, 1)
160+
end
161+
end
162+
163+
# small elements
164+
for small_element_index in 1:4
165+
small_element = cache.mortars.neighbor_ids[small_element_index, mortar]
166+
inverse_jacobian_small = get_inverse_jacobian(cache.elements.inverse_jacobian,
167+
mesh, indices_small...,
168+
small_element)
169+
170+
flux_small_high_order = get_node_vars(surface_flux_values_high_order,
171+
equations, dg,
172+
i, j, direction_small,
173+
small_element)
174+
flux_small_low_order = get_node_vars(surface_flux_values, equations, dg,
175+
i, j, direction_small,
176+
small_element)
177+
flux_difference_small = factor_small *
178+
(flux_small_high_order .- flux_small_low_order)
179+
180+
multiply_add_to_node_vars!(u,
181+
dt * inverse_jacobian_small *
182+
(1 - limiting_factor[mortar]),
183+
flux_difference_small, equations, dg,
184+
indices_small..., small_element)
185+
end
186+
187+
# large element
188+
large_element = cache.mortars.neighbor_ids[5, mortar]
189+
inverse_jacobian_large = get_inverse_jacobian(cache.elements.inverse_jacobian,
190+
mesh, indices_large...,
191+
large_element)
192+
193+
flux_large_high_order = get_node_vars(surface_flux_values_high_order,
194+
equations, dg, i, j, direction_large,
195+
large_element)
196+
flux_large_low_order = get_node_vars(surface_flux_values, equations, dg,
197+
i, j, direction_large, large_element)
198+
flux_difference_large = factor_large *
199+
(flux_large_high_order .- flux_large_low_order)
200+
201+
multiply_add_to_node_vars!(u,
202+
dt * inverse_jacobian_large *
203+
(1 - limiting_factor[mortar]),
204+
flux_difference_large, equations, dg,
205+
indices_large..., large_element)
206+
end
207+
end
208+
209+
return nothing
210+
end
87211
end # @muladd

0 commit comments

Comments
 (0)