Skip to content

Commit e407366

Browse files
bennibolmbennibolm
committed
Apply updated 3d version
1 parent 850ca81 commit e407366

File tree

1 file changed

+83
-62
lines changed

1 file changed

+83
-62
lines changed

‎src/callbacks_stage/subcell_limiter_idp_correction_3d.jl‎

Lines changed: 83 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -8,77 +8,98 @@
88
function perform_idp_correction!(u, dt,
99
mesh::P4estMesh{3},
1010
equations, dg, cache)
11-
@unpack inverse_weights = dg.basis # Plays role of DG subcell sizes
11+
@unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes
1212
@unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes
1313
@unpack alpha = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
1414

15+
# The following code implements the IDP correction in flux-differencing form:
16+
# u[v, i, j, k, element] += -dt * inverse_jacobian *
17+
# (inverse_weights[i] *
18+
# ((1 - alpha_1_ip1) * antidiffusive_flux1_ip1[v] - (1 - alpha_1) * antidiffusive_flux1[v]) +
19+
# inverse_weights[j] *
20+
# ((1 - alpha_2_jp1) * antidiffusive_flux2_jp1[v] - (1 - alpha_2) * antidiffusive_flux2[v]) +
21+
# inverse_weights[k] *
22+
# ((1 - alpha_3_kp1) * antidiffusive_flux3_kp1[v] - (1 - alpha_3) * antidiffusive_flux3[v]))
23+
24+
# For LGL nodes, the high-order and low-order fluxes at element interfaces are equal
25+
# and therefore, the antidiffusive fluxes are zero there.
26+
# To avoid adding zeros and speed up the simulation, we directly loop over the subcell
27+
# interfaces.
28+
1529
@threaded for element in eachelement(dg, cache)
16-
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
30+
# Perform correction in 1st/x-direction
31+
for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
32+
# Apply to right node
33+
alpha1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element])
34+
35+
# Sign switch as in apply_jacobian!
36+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
37+
mesh, i, j, k, element)
38+
flux1 = get_node_vars(antidiffusive_flux1_R, equations, dg,
39+
i, j, k, element)
40+
dg_factor = -dt * inverse_jacobian * inverse_weights[i] * (1 - alpha1)
41+
multiply_add_to_node_vars!(u, dg_factor, flux1,
42+
equations, dg, i, j, k, element)
43+
44+
# Apply to left node
45+
# Sign switch as in apply_jacobian!
46+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
47+
mesh, i - 1, j, k, element)
48+
flux1_ip1 = get_node_vars(antidiffusive_flux1_L, equations, dg,
49+
i, j, k, element)
50+
dg_factor = dt * inverse_jacobian * inverse_weights[i - 1] * (1 - alpha1)
51+
multiply_add_to_node_vars!(u, dg_factor, flux1_ip1,
52+
equations, dg, i - 1, j, k, element)
53+
end
54+
55+
# Perform correction in 2nd/y-direction
56+
for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
57+
# Apply to right node
58+
alpha2 = max(alpha[i, j - 1, k, element], alpha[i, j, k, element])
59+
1760
# Sign switch as in apply_jacobian!
1861
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
1962
mesh, i, j, k, element)
63+
flux2 = get_node_vars(antidiffusive_flux2_R, equations, dg,
64+
i, j, k, element)
65+
dg_factor = -dt * inverse_jacobian * inverse_weights[j] * (1 - alpha2)
66+
multiply_add_to_node_vars!(u, dg_factor, flux2,
67+
equations, dg, i, j, k, element)
68+
69+
# Apply to left node
70+
# Sign switch as in apply_jacobian!
71+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
72+
mesh, i, j - 1, k, element)
73+
flux2_jp1 = get_node_vars(antidiffusive_flux2_L, equations, dg,
74+
i, j, k, element)
75+
dg_factor = dt * inverse_jacobian * inverse_weights[j - 1] * (1 - alpha2)
76+
multiply_add_to_node_vars!(u, dg_factor, flux2_jp1,
77+
equations, dg, i, j - 1, k, element)
78+
end
2079

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
80+
# Perform correction in 3rd/z-direction
81+
for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
82+
# Apply to right node
83+
alpha3 = max(alpha[i, j, k - 1, element], alpha[i, j, k, element])
7284

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
85+
# Sign switch as in apply_jacobian!
86+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
87+
mesh, i, j, k, element)
88+
flux3 = get_node_vars(antidiffusive_flux3_R, equations, dg,
89+
i, j, k, element)
90+
dg_factor = -dt * inverse_jacobian * inverse_weights[k] * (1 - alpha3)
91+
multiply_add_to_node_vars!(u, dg_factor, flux3,
92+
equations, dg, i, j, k, element)
93+
94+
# Apply to left node
95+
# Sign switch as in apply_jacobian!
96+
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
97+
mesh, i, j, k - 1, element)
98+
flux3_kp1 = get_node_vars(antidiffusive_flux3_L, equations, dg,
99+
i, j, k, element)
100+
dg_factor = dt * inverse_jacobian * inverse_weights[k - 1] * (1 - alpha3)
101+
multiply_add_to_node_vars!(u, dg_factor, flux3_kp1,
102+
equations, dg, i, j, k - 1, element)
82103
end
83104
end
84105

0 commit comments

Comments
 (0)