Skip to content

Speed up idp correction function #2644

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
+155 −122
Open

Speed up idp correction function #2644

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
5d5701b
Add alternative idp correction
bennibolm Nov 10, 2025
221de04
Implement suggestions
bennibolm Nov 12, 2025
71e0ebe
Adapt test results
bennibolm Nov 12, 2025
1ef7c5e
Remove old version
bennibolm Nov 12, 2025
850ca81
Merge branch 'main' into speedup-idp-correction
DanielDoehring Nov 12, 2025
e407366
Apply updated 3d version
bennibolm Nov 13, 2025
7a12dd0
Merge branch 'main' into speedup-idp-correction
bennibolm Nov 19, 2025
e1d0170
Merge branch 'main' into speedup-idp-correction
bennibolm Dec 3, 2025
2051bac
Merge branch 'main' into speedup-idp-correction
DanielDoehring Dec 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 55 additions & 43 deletions src/callbacks_stage/subcell_limiter_idp_correction_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,55 +13,67 @@ function perform_idp_correction!(u, dt,
@unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes
@unpack alpha = dg.volume_integral.limiter.cache.subcell_limiter_coefficients

# The following code implements the IDP correction in flux-differencing form:
# u[v, i, j, element] += -dt * inverse_jacobian *
# (inverse_weights[i] *
# ((1 - alpha_1_ip1) * antidiffusive_flux1_ip1[v] - (1 - alpha_1) * antidiffusive_flux1[v]) +
# inverse_weights[j] *
# ((1 - alpha_2_jp1) * antidiffusive_flux2_jp1[v] - (1 - alpha_2) * antidiffusive_flux2[v]))

# For LGL nodes, the high-order and low-order fluxes at element interfaces are equal
# and therefore, the antidiffusive fluxes are zero there.
# To avoid adding zeros and speed up the simulation, we directly loop over the subcell
# interfaces.

@threaded for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
# Perform correction in 1st/x-direction
for j in eachnode(dg), i in 2:nnodes(dg)
# Apply to right node
alpha1 = max(alpha[i - 1, j, element], alpha[i, j, element])

# Sign switch as in apply_jacobian!
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
mesh, i, j, element)
flux1 = get_node_vars(antidiffusive_flux1_R, equations, dg,
i, j, element)
dg_factor = -dt * inverse_jacobian * inverse_weights[i] * (1 - alpha1)
multiply_add_to_node_vars!(u, dg_factor, flux1,
equations, dg, i, j, element)

# Apply to left node
# Sign switch as in apply_jacobian!
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
mesh, i - 1, j, element)
flux1_ip1 = get_node_vars(antidiffusive_flux1_L, equations, dg,
i, j, element)
dg_factor = dt * inverse_jacobian * inverse_weights[i - 1] * (1 - alpha1)
multiply_add_to_node_vars!(u, dg_factor, flux1_ip1,
equations, dg, i - 1, j, element)
end

# Note: For LGL nodes, the high-order and low-order fluxes at element interfaces are equal.
# Therefore, the antidiffusive fluxes are zero.
# To avoid accessing zero entries, we directly use zero vectors instead.
if i > 1 # Not at "left" boundary node
alpha1 = max(alpha[i - 1, j, element], alpha[i, j, element])
alpha_flux1 = (1 - alpha1) *
get_node_vars(antidiffusive_flux1_R, equations, dg,
i, j, element)
else # At "left" boundary node
alpha_flux1 = zero(SVector{nvariables(equations), eltype(u)})
end
if i < nnodes(dg) # Not at "right" boundary node
alpha1_ip1 = max(alpha[i, j, element], alpha[i + 1, j, element])
alpha_flux1_ip1 = (1 - alpha1_ip1) *
get_node_vars(antidiffusive_flux1_L, equations, dg,
i + 1, j, element)
else # At "right" boundary node
alpha_flux1_ip1 = zero(SVector{nvariables(equations), eltype(u)})
end
if j > 1 # Not at "bottom" boundary node
alpha2 = max(alpha[i, j - 1, element], alpha[i, j, element])
alpha_flux2 = (1 - alpha2) *
get_node_vars(antidiffusive_flux2_R, equations, dg,
i, j, element)
else # At "bottom" boundary node
alpha_flux2 = zero(SVector{nvariables(equations), eltype(u)})
end
if j < nnodes(dg) # Not at "top" boundary node
alpha2_jp1 = max(alpha[i, j, element], alpha[i, j + 1, element])
alpha_flux2_jp1 = (1 - alpha2_jp1) *
get_node_vars(antidiffusive_flux2_L, equations, dg,
i, j + 1, element)
else # At "top" boundary node
alpha_flux2_jp1 = zero(SVector{nvariables(equations), eltype(u)})
end
# Perform correction in 2nd/y-direction
for j in 2:nnodes(dg), i in eachnode(dg)
# Apply to right node
alpha2 = max(alpha[i, j - 1, element], alpha[i, j, element])

for v in eachvariable(equations)
u[v, i, j, element] += dt * inverse_jacobian *
(inverse_weights[i] *
(alpha_flux1_ip1[v] - alpha_flux1[v]) +
inverse_weights[j] *
(alpha_flux2_jp1[v] - alpha_flux2[v]))
end
# Sign switch as in apply_jacobian!
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
mesh, i, j, element)
flux2 = get_node_vars(antidiffusive_flux2_R, equations, dg,
i, j, element)
dg_factor = -dt * inverse_jacobian * inverse_weights[j] * (1 - alpha2)
multiply_add_to_node_vars!(u, dg_factor, flux2,
equations, dg, i, j, element)

# Apply to left node
# Sign switch as in apply_jacobian!
inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
mesh, i, j - 1, element)
flux2_jp1 = get_node_vars(antidiffusive_flux2_L, equations, dg,
i, j, element)
dg_factor = dt * inverse_jacobian * inverse_weights[j - 1] * (1 - alpha2)
multiply_add_to_node_vars!(u, dg_factor, flux2_jp1,
equations, dg, i, j - 1, element)
end
end

Expand Down
16 changes: 8 additions & 8 deletions test/test_tree_2d_euler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -360,16 +360,16 @@ end
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_euler_sedov_blast_wave_sc_subcell.jl"),
l2=[
0.4227549115123529,
0.14825759222652649,
0.14825759222682933,
0.6164668313131949
0.4227191130908862,
0.14825292449073538,
0.14822591031295396,
0.6164645445036752
],
linf=[
1.6391908143728386,
0.8344433355906021,
0.8344433355966195,
6.450305752671201
1.6394237885082292,
0.8374761298606049,
0.8322520901940953,
6.4503170484248855
],
tspan=(0.0, 1.0),
initial_refinement_level=4,
Expand Down
18 changes: 9 additions & 9 deletions test/test_tree_2d_eulermulti.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,18 +95,18 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem")
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"),
l2=[
73.41054363926742,
1.5072038797716156,
57405.58964098063,
0.17877099207437225,
0.010085388785440972
73.4069631318976,
1.5072402964611578,
57402.22856847899,
0.1787601960286661,
0.010085224169684756
],
linf=[
213.59140793740318,
24.57625853486584,
152498.21319871658,
213.591502983954,
24.576256265492944,
152497.46588262887,
0.5911106543157919,
0.09936092838440383
0.09936092397887711
],
initial_refinement_level=3,
tspan=(0.0, 0.001))
Expand Down
Loading