Skip to content

Commit 850ca81

Browse files
DanielDoehringDanielDoehring
authored
Merge branch 'main' into speedup-idp-correction
2 parents 1ef7c5e + 2e61748 commit 850ca81

14 files changed

+783
-7
lines changed

‎NEWS.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,15 @@ used in the Julia ecosystem. Notable changes will be documented in this file
66
for human readability.
77

88

9+
## Changes in the v0.13 lifecycle
10+
11+
#### Added
12+
- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582]).
13+
In the new version, IDP positivity limiting for conservative variables (using
14+
the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) is supported.
15+
`BoundsCheckCallback` is not supported in 3D yet.
16+
17+
918
## Changes when updating to v0.13 from v0.12.x
1019

1120
#### Changed

‎examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
using Trixi
2+
3+
###############################################################################
4+
# semidiscretization of the compressible Euler equations
5+
6+
equations = CompressibleEulerEquations3D(1.4)
7+
8+
"""
9+
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
10+
11+
The Sedov blast wave setup based on Flash
12+
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node191.html#SECTION010114000000000000000
13+
with smaller strength of the initial discontinuity.
14+
"""
15+
function initial_condition_sedov_blast_wave(x, t,
16+
equations::CompressibleEulerEquations3D)
17+
# Set up polar coordinates
18+
inicenter = SVector(0.0, 0.0, 0.0)
19+
x_norm = x[1] - inicenter[1]
20+
y_norm = x[2] - inicenter[2]
21+
z_norm = x[3] - inicenter[3]
22+
r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
23+
24+
# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node191.html#SECTION010114000000000000000
25+
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
26+
E = 1.0
27+
p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2)
28+
p0_outer = 1.0e-1 # "simpler" setup since positivity limiter for pressure is not yet supported in 3D
29+
30+
# Calculate primitive variables
31+
rho = 1.0
32+
v1 = 0.0
33+
v2 = 0.0
34+
v3 = 0.0
35+
p = r > r0 ? p0_outer : p0_inner
36+
37+
return prim2cons(SVector(rho, v1, v2, v3, p), equations)
38+
end
39+
40+
initial_condition = initial_condition_sedov_blast_wave
41+
42+
surface_flux = flux_lax_friedrichs
43+
volume_flux = flux_ranocha
44+
polydeg = 3
45+
basis = LobattoLegendreBasis(polydeg)
46+
limiter_idp = SubcellLimiterIDP(equations, basis;
47+
positivity_variables_cons = ["rho"])
48+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
49+
volume_flux_dg = volume_flux,
50+
volume_flux_fv = surface_flux)
51+
solver = DGSEM(basis, surface_flux, volume_integral)
52+
53+
coordinates_min = (-1.0, -1.0, -1.0)
54+
coordinates_max = (1.0, 1.0, 1.0)
55+
56+
trees_per_dimension = (8, 8, 8)
57+
mesh = P4estMesh(trees_per_dimension,
58+
polydeg = 1, initial_refinement_level = 0,
59+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
60+
periodicity = true)
61+
62+
# create the semi discretization object
63+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
64+
65+
###############################################################################
66+
# ODE solvers, callbacks etc.
67+
68+
tspan = (0.0, 3.0)
69+
ode = semidiscretize(semi, tspan)
70+
71+
summary_callback = SummaryCallback()
72+
73+
analysis_interval = 100
74+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
75+
76+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
77+
78+
save_solution = SaveSolutionCallback(interval = 10,
79+
save_initial_solution = true,
80+
save_final_solution = true,
81+
extra_node_variables = (:limiting_coefficient,))
82+
83+
stepsize_callback = StepsizeCallback(cfl = 0.5)
84+
85+
callbacks = CallbackSet(summary_callback,
86+
analysis_callback,
87+
alive_callback,
88+
save_solution,
89+
stepsize_callback)
90+
91+
###############################################################################
92+
# run the simulation
93+
94+
stage_callbacks = (SubcellLimiterIDPCorrection(),)
95+
96+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
97+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
98+
callback = callbacks);

‎src/callbacks_stage/subcell_limiter_idp_correction.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,4 +64,5 @@ init_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
6464
finalize_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
6565

6666
include("subcell_limiter_idp_correction_2d.jl")
67+
include("subcell_limiter_idp_correction_3d.jl")
6768
end # @muladd

‎src/callbacks_stage/subcell_limiter_idp_correction_3d.jl‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
2+
# Since these FMAs can increase the performance of many numerical algorithms,
3+
# we need to opt-in explicitly.
4+
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
5+
@muladd begin
6+
#! format: noindent
7+
8+
function perform_idp_correction!(u, dt,
9+
mesh::P4estMesh{3},
10+
equations, dg, cache)
11+
@unpack inverse_weights = dg.basis # Plays role of DG subcell sizes
12+
@unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes
13+
@unpack alpha = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
14+
15+
@threaded for element in eachelement(dg, cache)
16+
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
17+
# Sign switch as in apply_jacobian!
18+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
19+
mesh, i, j, k, element)
20+
21+
# Note: For LGL nodes, the high-order and low-order fluxes at element interfaces are equal.
22+
# Therefore, the antidiffusive fluxes are zero.
23+
# To avoid accessing zero entries, we directly use zero vectors instead.
24+
if i > 1 # Not at "left" boundary node
25+
alpha1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element])
26+
alpha_flux1 = (1 - alpha1) *
27+
get_node_vars(antidiffusive_flux1_R, equations, dg,
28+
i, j, k, element)
29+
else # At "left" boundary node
30+
alpha_flux1 = zero(SVector{nvariables(equations), eltype(u)})
31+
end
32+
if i < nnodes(dg) # Not at "right" boundary node
33+
alpha1_ip1 = max(alpha[i, j, k, element], alpha[i + 1, j, k, element])
34+
alpha_flux1_ip1 = (1 - alpha1_ip1) *
35+
get_node_vars(antidiffusive_flux1_L, equations, dg,
36+
i + 1, j, k, element)
37+
else # At "right" boundary node
38+
alpha_flux1_ip1 = zero(SVector{nvariables(equations), eltype(u)})
39+
end
40+
if j > 1 # Not at "bottom" boundary node
41+
alpha2 = max(alpha[i, j - 1, k, element], alpha[i, j, k, element])
42+
alpha_flux2 = (1 - alpha2) *
43+
get_node_vars(antidiffusive_flux2_R, equations, dg,
44+
i, j, k, element)
45+
else # At "bottom" boundary node
46+
alpha_flux2 = zero(SVector{nvariables(equations), eltype(u)})
47+
end
48+
if j < nnodes(dg) # Not at "top" boundary node
49+
alpha2_jp1 = max(alpha[i, j, k, element], alpha[i, j + 1, k, element])
50+
alpha_flux2_jp1 = (1 - alpha2_jp1) *
51+
get_node_vars(antidiffusive_flux2_L, equations, dg,
52+
i, j + 1, k, element)
53+
else # At "top" boundary node
54+
alpha_flux2_jp1 = zero(SVector{nvariables(equations), eltype(u)})
55+
end
56+
if k > 1 # Not at "front" boundary node
57+
alpha3 = max(alpha[i, j, k - 1, element], alpha[i, j, k, element])
58+
alpha_flux3 = (1 - alpha3) *
59+
get_node_vars(antidiffusive_flux3_R, equations, dg,
60+
i, j, k, element)
61+
else # At "front" boundary node
62+
alpha_flux3 = zero(SVector{nvariables(equations), eltype(u)})
63+
end
64+
if k < nnodes(dg) # Not at "back" boundary node
65+
alpha3_kp1 = max(alpha[i, j, k, element], alpha[i, j, k + 1, element])
66+
alpha_flux3_kp1 = (1 - alpha3_kp1) *
67+
get_node_vars(antidiffusive_flux3_L, equations, dg,
68+
i, j, k + 1, element)
69+
else # At "back" boundary node
70+
alpha_flux3_kp1 = zero(SVector{nvariables(equations), eltype(u)})
71+
end
72+
73+
for v in eachvariable(equations)
74+
u[v, i, j, k, element] += dt * inverse_jacobian *
75+
(inverse_weights[i] *
76+
(alpha_flux1_ip1[v] - alpha_flux1[v]) +
77+
inverse_weights[j] *
78+
(alpha_flux2_jp1[v] - alpha_flux2[v]) +
79+
inverse_weights[k] *
80+
(alpha_flux3_kp1[v] - alpha_flux3[v]))
81+
end
82+
end
83+
end
84+
85+
return nothing
86+
end
87+
end # @muladd

‎src/solvers/dgsem_p4est/dg.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,4 +91,6 @@ include("dg_3d_parabolic.jl")
9191
include("dg_parallel.jl")
9292

9393
include("subcell_limiters_2d.jl")
94+
include("subcell_limiters_3d.jl")
95+
include("dg_3d_subcell_limiters.jl")
9496
end # @muladd

0 commit comments

Comments
 (0)