diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 4caaff8fc17..c9abba46219 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -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 diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl index 09269a89593..6622b58ac38 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl @@ -8,77 +8,98 @@ function perform_idp_correction!(u, dt, mesh::P4estMesh{3}, equations, dg, cache) - @unpack inverse_weights = dg.basis # Plays role of DG subcell sizes + @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_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, k, 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]) + + # inverse_weights[k] * + # ((1 - alpha_3_kp1) * antidiffusive_flux3_kp1[v] - (1 - alpha_3) * antidiffusive_flux3[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 k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + # Perform correction in 1st/x-direction + for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg) + # Apply to right node + alpha1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element]) + + # Sign switch as in apply_jacobian! + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, k, element) + flux1 = get_node_vars(antidiffusive_flux1_R, equations, dg, + i, j, k, element) + dg_factor = -dt * inverse_jacobian * inverse_weights[i] * (1 - alpha1) + multiply_add_to_node_vars!(u, dg_factor, flux1, + equations, dg, i, j, k, element) + + # Apply to left node + # Sign switch as in apply_jacobian! + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i - 1, j, k, element) + flux1_ip1 = get_node_vars(antidiffusive_flux1_L, equations, dg, + i, j, k, 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, k, element) + end + + # Perform correction in 2nd/y-direction + for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg) + # Apply to right node + alpha2 = max(alpha[i, j - 1, k, element], alpha[i, j, k, element]) + # Sign switch as in apply_jacobian! inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, mesh, i, j, k, element) + flux2 = get_node_vars(antidiffusive_flux2_R, equations, dg, + i, j, k, element) + dg_factor = -dt * inverse_jacobian * inverse_weights[j] * (1 - alpha2) + multiply_add_to_node_vars!(u, dg_factor, flux2, + equations, dg, i, j, k, element) + + # Apply to left node + # Sign switch as in apply_jacobian! + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j - 1, k, element) + flux2_jp1 = get_node_vars(antidiffusive_flux2_L, equations, dg, + i, j, k, 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, k, 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, k, element], alpha[i, j, k, element]) - alpha_flux1 = (1 - alpha1) * - get_node_vars(antidiffusive_flux1_R, equations, dg, - i, j, k, 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, k, element], alpha[i + 1, j, k, element]) - alpha_flux1_ip1 = (1 - alpha1_ip1) * - get_node_vars(antidiffusive_flux1_L, equations, dg, - i + 1, j, k, 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, k, element], alpha[i, j, k, element]) - alpha_flux2 = (1 - alpha2) * - get_node_vars(antidiffusive_flux2_R, equations, dg, - i, j, k, 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, k, element], alpha[i, j + 1, k, element]) - alpha_flux2_jp1 = (1 - alpha2_jp1) * - get_node_vars(antidiffusive_flux2_L, equations, dg, - i, j + 1, k, element) - else # At "top" boundary node - alpha_flux2_jp1 = zero(SVector{nvariables(equations), eltype(u)}) - end - if k > 1 # Not at "front" boundary node - alpha3 = max(alpha[i, j, k - 1, element], alpha[i, j, k, element]) - alpha_flux3 = (1 - alpha3) * - get_node_vars(antidiffusive_flux3_R, equations, dg, - i, j, k, element) - else # At "front" boundary node - alpha_flux3 = zero(SVector{nvariables(equations), eltype(u)}) - end - if k < nnodes(dg) # Not at "back" boundary node - alpha3_kp1 = max(alpha[i, j, k, element], alpha[i, j, k + 1, element]) - alpha_flux3_kp1 = (1 - alpha3_kp1) * - get_node_vars(antidiffusive_flux3_L, equations, dg, - i, j, k + 1, element) - else # At "back" boundary node - alpha_flux3_kp1 = zero(SVector{nvariables(equations), eltype(u)}) - end + # Perform correction in 3rd/z-direction + for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg) + # Apply to right node + alpha3 = max(alpha[i, j, k - 1, element], alpha[i, j, k, element]) - for v in eachvariable(equations) - u[v, i, j, k, element] += dt * inverse_jacobian * - (inverse_weights[i] * - (alpha_flux1_ip1[v] - alpha_flux1[v]) + - inverse_weights[j] * - (alpha_flux2_jp1[v] - alpha_flux2[v]) + - inverse_weights[k] * - (alpha_flux3_kp1[v] - alpha_flux3[v])) - end + # Sign switch as in apply_jacobian! + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, k, element) + flux3 = get_node_vars(antidiffusive_flux3_R, equations, dg, + i, j, k, element) + dg_factor = -dt * inverse_jacobian * inverse_weights[k] * (1 - alpha3) + multiply_add_to_node_vars!(u, dg_factor, flux3, + equations, dg, i, j, k, element) + + # Apply to left node + # Sign switch as in apply_jacobian! + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, k - 1, element) + flux3_kp1 = get_node_vars(antidiffusive_flux3_L, equations, dg, + i, j, k, element) + dg_factor = dt * inverse_jacobian * inverse_weights[k - 1] * (1 - alpha3) + multiply_add_to_node_vars!(u, dg_factor, flux3_kp1, + equations, dg, i, j, k - 1, element) end end diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index da5305f0217..e395a77f40c 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -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, diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 9254ccaaeaf..0dfea90561e 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -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))