DGMulti shock capturing part 2: SBP approximation types (#2942)

* changes related to cleanup of existing shock capturing

remove extra end

remove unnecessary comment

* add 3d shock capturing elixir

add 3d test

* formatting

* Adding shock capturing routines for SBP on Line/Quad elements

* add test for non-Gauss SBP shock capturing

* fix naming convention

* updating tests

Removing curved shock capturing test since we haven't implemented low order free-stream metrics

* update CI values with x86 run values

* update initial condition to avoid x86/aarch64 differences with CI vals

* Use weak C0 blast IC for DGMulti Euler shock capturing

Replace the polar-coordinate blast setup with a Cartesian weak C0 initial condition so we avoid floating-point trouble from polar evaluation and from sampling exactly on the discontinuity.

Update the GaussSBP and SBP regression L2/L∞ norms in test_dgmulti_2d.jl for elixir_euler_shockcapturing.jl to match the new IC.

Made-with: Cursor

* Apply suggestions from code review

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>

* comments on 2 * f_ij and sparse operator column entries rows[j]

* fix formatting

* Update src/solvers/dgmulti/shock_capturing.jl

* Update examples/dgmulti_2d/elixir_euler_shockcapturing.jl

* Update src/solvers/dgmulti/shock_capturing.jl

* Update src/solvers/dgmulti/shock_capturing.jl

---------

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
Co-authored-by: Daniel Doehring <daniel.doehring@rwth-aachen.de>

Author: jlchan

examples/dgmulti_2d/elixir_euler_shockcapturing.jl

@@ -6,16 +6,48 @@ using Trixi
 
 equations = CompressibleEulerEquations2D(1.4)
 
-initial_condition = initial_condition_weak_blast_wave
-
-# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
-# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
-# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
-# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
-# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
-# `StepsizeCallback` (CFL-Condition) and less diffusion.
-surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
+# A continuous version of the blast wave initial condition to avoid floating point 
+# issues when evaluating polar coordinates or evaluating at the discontinuity
+function initial_condition_weak_C0_blast_wave(x, t,
+                                              equations::CompressibleEulerEquations2D)
+    RealT = eltype(x)
+    inicenter = SVector(0, 0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    rho_outer = one(RealT)
+    v1_outer = zero(RealT)
+    v2_outer = zero(RealT)
+    p_outer = one(RealT)
+    rho_inner = 1.1691
+    v1_inner = 0.1882
+    v2_inner = 0.1882
+    p_inner = 1.245
+
+    # Calculate primitive variables
+    if r > 0.5f0
+        rho = rho_outer
+        v1 = v1_outer
+        v2 = v2_outer
+        p = p_outer
+    elseif isapprox(r, 0.5f0)
+        rho = 0.5f0 * (rho_outer + rho_inner)
+        v1 = 0.5f0 * (v1_outer + v1_inner)
+        v2 = 0.5f0 * (v2_outer + v2_inner)
+        p = 0.5f0 * (p_outer + p_inner)
+    else
+        rho = rho_inner
+        v1 = v1_inner
+        v2 = v2_inner
+        p = p_inner
+    end
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_weak_C0_blast_wave
+
+surface_flux = FluxLaxFriedrichs()
 volume_flux = flux_ranocha
 
 polydeg = 3

src/solvers/dgmulti/flux_differencing.jl

@@ -129,9 +129,14 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP,
     solution_container = initialize_dgmulti_solution_container(mesh, equations, dg,
                                                                uEltype)
 
+    # this calls the `create_cache` for the shock capturing volume integral                          
+    volume_integral_cache = create_cache(mesh, equations, dg.volume_integral,
+                                         dg, RealT, uEltype)
+
     return (; md, Qrst_skew, dxidxhatj = md.rstxyzJ,
             invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
-            solution_container, du_local_threaded)
+            solution_container, du_local_threaded,
+            volume_integral_cache...)
 end
 
 # most general create_cache: works for `DGMultiFluxDiff{<:Polynomial}`

src/solvers/dgmulti/shock_capturing.jl

@@ -30,6 +30,24 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
             element_to_element_connectivity)
 end
 
+function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
+                      volume_integral::Union{VolumeIntegralShockCapturingHGType,
+                                             VolumeIntegralPureLGLFiniteVolume},
+                      dg::DGMultiFluxDiffSBP, RealT, uEltype) where {NDIMS}
+    element_to_element_connectivity = build_element_to_element_connectivity(mesh, dg)
+
+    # create skew-symmetric parts of sparse hybridized operators for low order scheme.
+    sparse_SBP_operators, _ = StartUpDG.sparse_low_order_SBP_operators(dg.basis)
+    sparse_SBP_operators = map(A -> 0.5f0 * (A - A'), sparse_SBP_operators)
+
+    # Find the joint sparsity pattern of the entire matrix. We store the sparsity pattern as
+    # an adjoint for faster iteration through the rows.
+    sparsity_pattern = sum(map(A -> abs.(A)', sparse_SBP_operators)) .> 100 * eps()
+
+    return (; sparse_SBP_operators, sparsity_pattern,
+            element_to_element_connectivity)
+end
+
 # this method is used when the indicator is constructed as for shock-capturing volume integrals
 function create_cache(::Type{IndicatorHennemannGassner}, equations::AbstractEquations,
                       basis::DGMultiBasis{NDIMS}) where {NDIMS}
@@ -305,6 +323,11 @@ function volume_integral_kernel!(du, u, element, mesh::DGMultiMesh,
         u_i = u_local[i]
         du_i = zero(u_i)
         for id in nzrange(A_base, i)
+            # nonzero column indices for row i of the sparse operator. 
+            # note that because Julia uses SparseMatrixCSC, rows[id] 
+            # are efficient to access. We assume here that `sparsity_pattern`
+            # is symmetric (which is true since A_base is skew-symmetric), 
+            # so nonzero row indices are the same as nonzero column indices.
             j = rows[id]
             u_j = u_local[j]
 
@@ -317,6 +340,13 @@ function volume_integral_kernel!(du, u, element, mesh::DGMultiMesh,
             # note that we do not need to normalize `normal_direction_ij` since
             # it is typically normalized within the flux computation.
             f_ij = volume_flux_fv(u_i, u_j, normal_direction_ij, equations)
+
+            # the factor of 2 is for consistency; for example, if f_ij is the central 
+            # flux, flux differencing with a differentiation matrix should recover the 
+            # flux derivative via
+            #   \sum_j 2 * D_ij * f_ij = \sum_j 2 * D_ij * 0.5 * (f(u_i) + f(u_j))
+            #                          = f(u_i) \sum_j D_ij + \sum_j D_ij f(u_j)
+            #                          = 0 (since \sum_j D_ij = 0) + (D * f(u))_i
             du_i = du_i + 2 * f_ij
         end
         rhs_local[i] = du_i
@@ -325,4 +355,54 @@ function volume_integral_kernel!(du, u, element, mesh::DGMultiMesh,
     # TODO: factor this out to avoid calling it twice during calc_volume_integral!
     return project_rhs_to_gauss_nodes!(du, rhs_local, element, mesh, dg, cache, alpha)
 end
+
+# Calculates the volume integral corresponding to an algebraic low order method for
+# DGMultiFluxDiffSBP (traditional SBP operators with LGL-type nodes).
+# Unlike GaussSBP, the solution lives at nodes that include the face nodes (at positions
+# `rd.Fmask`), so no entropy projection is needed. We build the extended [interior; face]
+# vector in-kernel and project back by scattering face contributions to Fmask positions.
+function volume_integral_kernel!(du, u, element, mesh::DGMultiMesh,
+                                 have_nonconservative_terms::False, equations,
+                                 volume_integral::VolumeIntegralPureLGLFiniteVolume,
+                                 dg::DGMultiFluxDiffSBP, cache, alpha = true)
+    (; volume_flux_fv) = volume_integral
+
+    (; inv_wq, sparsity_pattern) = cache
+    A_base, row_ids, rows, _ = sparse_operator_data(sparsity_pattern)
+    for i in row_ids
+        u_i = u[i, element]
+        du_i = zero(u_i)
+        for id in nzrange(A_base, i)
+            # nonzero column indices for row i of the sparse operator. 
+            # note that because Julia uses SparseMatrixCSC, rows[id] 
+            # are efficient to access. We assume here that `sparsity_pattern`
+            # is symmetric (which is true since A_base is skew-symmetric), 
+            # so nonzero row indices are the same as nonzero column indices.
+            j = rows[id]
+            u_j = u[j, element]
+
+            # compute (Q_1[i,j], Q_2[i,j], ...) where Q_i = ∑_j dxidxhatj * Q̂_j
+            geometric_matrix = get_low_order_geometric_matrix(i, j, element, mesh,
+                                                              cache)
+            reference_operator_entries = get_sparse_operator_entries(i, j, mesh,
+                                                                     cache)
+            normal_direction_ij = geometric_matrix * reference_operator_entries
+
+            # note that we do not need to normalize `normal_direction_ij` since
+            # it is typically normalized within the flux computation.
+            f_ij = volume_flux_fv(u_i, u_j, normal_direction_ij, equations)
+
+            # the factor of 2 is for consistency; for example, if f_ij is the central 
+            # flux, flux differencing with a differentiation matrix should recover the 
+            # flux derivative via
+            #   \sum_j 2 * D_ij * f_ij = \sum_j 2 * D_ij * 0.5 * (f(u_i) + f(u_j))
+            #                          = f(u_i) \sum_j D_ij + \sum_j D_ij f(u_j)
+            #                          = 0 (since \sum_j D_ij = 0) + (D * f(u))_i
+            du_i = du_i + 2 * f_ij
+        end
+        du[i, element] = du[i, element] + alpha * du_i * inv_wq[i]
+    end
+
+    return nothing
+end
 end # @muladd

test/test_dgmulti_2d.jl

@@ -421,16 +421,37 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
                         cells_per_dimension=4, tspan=(0.0, 0.1),
                         l2=[
-                            0.05685180852320552,
-                            0.04308097439005265,
-                            0.04308097439005263,
-                            0.21098250258804
+                            0.04459267863712344,
+                            0.04466070180923881,
+                            0.044660701809239055,
+                            0.16853101962882136
                         ],
                         linf=[
-                            0.2360805191601203,
-                            0.16684117462697776,
-                            0.16684117462697767,
-                            0.8573034682049414
+                            0.18547384414190626,
+                            0.20141544103810083,
+                            0.20141544103810224,
+                            0.6954971316946836
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
+@trixi_testset "elixir_euler_shockcapturing.jl (SBP)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
+                        cells_per_dimension=4, tspan=(0.0, 0.1),
+                        basis=DGMultiBasis(Quad(), 3, approximation_type = SBP()),
+                        l2=[
+                            0.03901629442791245,
+                            0.03830525262399631,
+                            0.03830525262399637,
+                            0.14746814850975215
+                        ],
+                        linf=[
+                            0.16680108480009226,
+                            0.21077037377694813,
+                            0.2107703737769464,
+                            0.6400748796169138
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)