Speed up idp correction function

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -13,55 +13,73 @@ 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[i, j, element] *
+    #    (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]))
+    # with
+    # alpha_1 = max(alpha[i - 1, j, element], alpha[i, j, element]),
+    # alpha_1_ip1 = max(alpha[i, j, element], alpha[i + 1, j, element])
+    # and equivalently for alpha_2 and alpha_2_jp1.
+
+    # 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)
+            # Subcell interface between nodes (i - 1, j) and (i, j)
+            alpha1 = max(alpha[i - 1, j, element], alpha[i, j, element])
+
+            # Apply to right node (i, j)
             # 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 (i - 1, j)
+            # 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)
+            # Subcell interface between nodes (i, j - 1) and (i, j)
+            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
+            # Apply to right node (i, j)
+            # 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 (i, j - 1)
+            # 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
 

src/callbacks_stage/subcell_limiter_idp_correction_3d.jl

@@ -8,77 +8,105 @@
 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[i, j, k, element] *
+    #    (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]))
+    # with
+    # alpha_1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element]),
+    # alpha_1_ip1 = max(alpha[i, j, k, element], alpha[i + 1, j, k, element])
+    # and equivalently for alpha_2, alpha_2_jp1, alpha_3, alpha_3_kp1.
+
+    # 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)
+            # Subcell interface between nodes (i - 1, j, k) and (i, j, k)
+            alpha1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element])
+
+            # Apply to right node (i, j, k)
+            # 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 (i - 1, j, k)
+            # 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)
+            # Subcell interface between nodes (i, j - 1, k) and (i, j, k)
+            alpha2 = max(alpha[i, j - 1, k, element], alpha[i, j, k, element])
+
+            # Apply to right node (i, j, k)
             # 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 (i, j - 1, k)
+            # 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)
+            # Subcell interface between nodes (i, j, k - 1) and (i, j, k)
+            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
+            # Apply to right node (i, j, k)
+            # 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 (i, j, k - 1)
+            # 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
 

test/test_tree_2d_euler.jl

@@ -402,16 +402,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,

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))