DGMulti shock capturing part 1: cleanup (#2940)

* changes related to cleanup of existing shock capturing

remove extra end

remove unnecessary comment

* add 3d shock capturing elixir

add 3d test

* formatting

* fix naming convention

* Update src/solvers/dgmulti/shock_capturing.jl

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

* Update src/solvers/dgmulti/shock_capturing.jl

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

* Apply suggestions from code review

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

* formatting

* num_modes -> nmodes

---------

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Author: 3 people

examples/dgmulti_3d/elixir_euler_shockcapturing.jl

@@ -0,0 +1,52 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+initial_condition = initial_condition_weak_blast_wave
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+
+polydeg = 3
+basis = DGMultiBasis(Hex(), polydeg, approximation_type = GaussSBP())
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+dg = DGMulti(basis,
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = volume_integral)
+
+cells_per_dimension = (4, 4, 4)
+mesh = DGMultiMesh(dg, cells_per_dimension;
+                   coordinates_min = (-2.0, -2.0, -2.0),
+                   coordinates_max = (2.0, 2.0, 2.0),
+                   periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    boundary_conditions = boundary_condition_periodic)
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval = 10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     solution_variables = cons2prim)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-6, reltol = 1.0e-6,
+            ode_default_options()..., callback = callbacks);

src/solvers/dgmulti/dg.jl

@@ -143,6 +143,13 @@ In particular, not the face nodes themselves are returned.
     return Base.OneTo(dg.basis.Nfq * mesh.md.num_elements)
 end
 
+# The `nmodes` functions returns the number of polynomial modes for a degree N
+# approximation on a specific type of element. For tensor product elements (Quad, Hex),
+# this is the total number of modes in all dimensions, i.e., (N+1)^NDIMS.
+@inline nmodes(N, ::Line) = N + 1
+@inline nmodes(N, ::Quad) = (N + 1)^2
+@inline nmodes(N, ::Hex) = (N + 1)^3
+
 # interface with semidiscretization_hyperbolic
 wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
 wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode

src/solvers/dgmulti/flux_differencing.jl

@@ -73,12 +73,28 @@ end
 
 # Build element-to-element connectivity from face-to-face connectivity.
 # Used for smoothing of shock capturing blending parameters (see `apply_smoothing!`).
+#
+# Here, `mesh.md.FToF` is a `(num_faces_per_element × num_elements)` array where
+# `FToF[f, e]` stores the global face index of the neighbor of local face `f` on
+# element `e`. 
+#
+# Global face indices are laid out as
+#
+#   global_face_index = (element_index - 1) * num_faces + local_face_index,
+# 
+# so that the element index can be recovered by integer division:
+#
+#   element_index = (global_face_index - 1) ÷ num_faces + 1.
+# 
+# For a non-periodic boundary face, `FToF[f, e]` points back to face `f` of element 
+# `e` itself, so boundary elements are listed as their own neighbor.
 function build_element_to_element_connectivity(mesh::DGMultiMesh, dg::DGMulti)
     face_to_face_connectivity = mesh.md.FToF
     element_to_element_connectivity = similar(face_to_face_connectivity)
     for e in axes(face_to_face_connectivity, 2)
         for f in axes(face_to_face_connectivity, 1)
             neighbor_face_index = face_to_face_connectivity[f, e]
+
             # Reverse-engineer element index from face index. Assumes all elements
             # have the same number of faces.
             neighbor_element_index = ((neighbor_face_index - 1) ÷ dg.basis.num_faces) +

src/solvers/dgmulti/shock_capturing.jl

@@ -1,3 +1,9 @@
+# Since `@muladd` can fuse multiply-add operations and thus improve performance in
+# the flux differencing loops, we opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
 # by default, return an empty tuple for volume integral caches
 function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
                       volume_integral::VolumeIntegralShockCapturingHGType,
@@ -12,40 +18,52 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
     # create sparse hybridized operators for low order scheme
     Qrst, E = StartUpDG.sparse_low_order_SBP_operators(dg.basis)
     Brst = map(n -> Diagonal(n .* dg.basis.wf), dg.basis.nrstJ)
-    sparse_hybridized_SBP_operators = map((Q, B) -> 0.5 * [Q-Q' E'*B; -B*E zeros(size(B))],
-                                          Qrst, Brst)
+    sparse_SBP_operators = map((Q, B) -> 0.5f0 * [Q-Q' E'*B; -B*E zeros(size(B))],
+                               Qrst, Brst)
 
     # 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_hybridized_SBP_operators)) .>
+    sparsity_pattern = sum(map(A -> abs.(A)', sparse_SBP_operators)) .>
                        100 * eps()
 
-    return (; sparse_hybridized_SBP_operators, sparsity_pattern,
+    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::RefElemData{NDIMS}) where {NDIMS}
+                      basis::DGMultiBasis{NDIMS}) where {NDIMS}
     uEltype = real(basis)
     alpha = Vector{uEltype}()
     alpha_tmp = similar(alpha)
 
-    MVec = MVector{nnodes(basis), uEltype}
-    indicator_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
-    modal_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
+    MVec_nodes = MVector{nnodes(basis), uEltype}
+    indicator_threaded = MVec_nodes[MVec_nodes(undef) for _ in 1:Threads.maxthreadid()]
+    MVec_modes = MVector{nmodes(basis.N, basis.element_type), uEltype}
+    modal_threaded = MVec_modes[MVec_modes(undef) for _ in 1:Threads.maxthreadid()]
+
+    inverse_vandermonde = calc_inverse_vandermonde(basis)
 
+    return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, inverse_vandermonde)
+end
+
+# calculates the inverse of the Vandermonde matrix for shock capturing purposes.
+# This version is for tensor product elements (Line, Quad, Hex)
+function calc_inverse_vandermonde(basis::DGMultiBasis{NDIMS, <:Union{Line, Quad, Hex}}) where {NDIMS}
     # initialize inverse Vandermonde matrices at Gauss-Legendre nodes
     (; N) = basis
     lobatto_node_coordinates_1D, _ = StartUpDG.gauss_lobatto_quad(0, 0, N)
     VDM_1D = StartUpDG.vandermonde(Line(), N, lobatto_node_coordinates_1D)
     inverse_vandermonde = SimpleKronecker(NDIMS, inv(VDM_1D))
 
-    return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, inverse_vandermonde)
+    return inverse_vandermonde
 end
 
 function (indicator_hg::IndicatorHennemannGassner)(u, mesh::DGMultiMesh,
-                                                   equations, dg::DGMulti{NDIMS}, cache;
+                                                   equations,
+                                                   dg::DGMulti{NDIMS,
+                                                               <:Union{Line, Quad, Hex}},
+                                                   cache;
                                                    kwargs...) where {NDIMS}
     (; alpha_max, alpha_min, alpha_smooth, variable) = indicator_hg
     (; alpha, alpha_tmp, indicator_threaded, modal_threaded, inverse_vandermonde) = indicator_hg.cache
@@ -56,7 +74,7 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::DGMultiMesh,
     end
 
     # magic parameters
-    threshold = 0.5 * 10^(-1.8 * (dg.basis.N + 1)^0.25)
+    threshold = 0.5f0 * 10^(-1.8 * (dg.basis.N + 1)^0.25)
     parameter_s = log((1 - 0.0001) / 0.0001)
 
     @threaded for element in eachelement(mesh, dg)
@@ -101,7 +119,8 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::DGMultiMesh,
             energy_frac_1 = zero(total_energy)
         end
         if !(iszero(total_energy_clip1))
-            energy_frac_2 = (total_energy_clip1 - total_energy_clip2) / total_energy_clip1
+            energy_frac_2 = (total_energy_clip1 - total_energy_clip2) /
+                            total_energy_clip1
         else
             energy_frac_2 = zero(total_energy_clip1)
         end
@@ -134,15 +153,15 @@ end
 # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha
 function apply_smoothing!(mesh::DGMultiMesh, alpha, alpha_tmp, dg::DGMulti, cache)
 
-    # Copy alpha values such that smoothing is indpedenent of the element access order
+    # Copy alpha values such that smoothing is independent of the element access order
     alpha_tmp .= alpha
 
     # smooth alpha with its neighboring value
-    for element in eachelement(mesh, dg)
+    @threaded for element in eachelement(mesh, dg)
         for face in Base.OneTo(StartUpDG.num_faces(dg.basis.element_type))
             neighboring_element = cache.element_to_element_connectivity[face, element]
             alpha_neighbor = alpha_tmp[neighboring_element]
-            alpha[element] = max(alpha[element], 0.5 * alpha_neighbor)
+            alpha[element] = max(alpha[element], 0.5f0 * alpha_neighbor)
         end
     end
 
@@ -194,16 +213,16 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
 end
 
 function get_sparse_operator_entries(i, j, mesh::DGMultiMesh{1}, cache)
-    return SVector(cache.sparse_hybridized_SBP_operators[1][i, j])
+    return SVector(cache.sparse_SBP_operators[1][i, j])
 end
 
 function get_sparse_operator_entries(i, j, mesh::DGMultiMesh{2}, cache)
-    Qr, Qs = cache.sparse_hybridized_SBP_operators
+    Qr, Qs = cache.sparse_SBP_operators
     return SVector(Qr[i, j], Qs[i, j])
 end
 
 function get_sparse_operator_entries(i, j, mesh::DGMultiMesh{3}, cache)
-    Qr, Qs, Qt = cache.sparse_hybridized_SBP_operators
+    Qr, Qs, Qt = cache.sparse_SBP_operators
     return SVector(Qr[i, j], Qs[i, j], Qt[i, j])
 end
 
@@ -244,7 +263,7 @@ function get_contravariant_matrix(i, element, mesh::DGMultiMesh{3}, cache)
 end
 
 function get_avg_contravariant_matrix(i, j, element, mesh::DGMultiMesh, cache)
-    return 0.5 * (get_contravariant_matrix(i, element, mesh, cache) +
+    return 0.5f0 * (get_contravariant_matrix(i, element, mesh, cache) +
             get_contravariant_matrix(j, element, mesh, cache))
 end
 
@@ -290,7 +309,8 @@ function volume_integral_kernel!(du, u, element, mesh::DGMultiMesh,
             u_j = u_local[j]
 
             # 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)
+            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
 
@@ -305,3 +325,4 @@ 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
+end # @muladd

src/solvers/dgmulti/types.jl

@@ -23,19 +23,21 @@ const DGMultiFluxDiff{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxTyp
                                                       <:SurfaceIntegralWeakForm,
                                                       <:Union{VolumeIntegralFluxDifferencing,
                                                               VolumeIntegralShockCapturingHGType,
-                                                              VolumeIntegralAdaptiveEC_WF_DG}} where {
-                                                                                                      NDIMS
-                                                                                                      }
+                                                              VolumeIntegralAdaptiveEC_WF_DG,
+                                                              VolumeIntegralPureLGLFiniteVolume}} where {
+                                                                                                         NDIMS
+                                                                                                         }
 
 const DGMultiFluxDiffSBP{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxType,
                                                          <:SurfaceIntegralWeakForm,
                                                          <:Union{VolumeIntegralFluxDifferencing,
-                                                                 VolumeIntegralShockCapturingHGType}} where {
-                                                                                                             NDIMS,
-                                                                                                             ApproxType <:
-                                                                                                             Union{SBP,
-                                                                                                                   AbstractDerivativeOperator}
-                                                                                                             }
+                                                                 VolumeIntegralShockCapturingHGType,
+                                                                 VolumeIntegralPureLGLFiniteVolume}} where {
+                                                                                                            NDIMS,
+                                                                                                            ApproxType <:
+                                                                                                            Union{SBP,
+                                                                                                                  AbstractDerivativeOperator}
+                                                                                                            }
 
 const DGMultiSBP{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxType,
                                                  SurfaceIntegral,
@@ -169,6 +171,10 @@ Constructs a basis for DGMulti solvers. Returns a "StartUpDG.RefElemData" object
   For more info, see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/).
 
 """
+const DGMultiBasis{NDIMS, element_type, approximation_type} = StartUpDG.RefElemData{NDIMS,
+                                                                                    element_type,
+                                                                                    approximation_type}
+
 function DGMultiBasis(element_type, polydeg; approximation_type = Polynomial(),
                       kwargs...)
     return RefElemData(element_type, approximation_type, polydeg; kwargs...)

test/test_dgmulti_3d.jl

@@ -344,6 +344,18 @@ end
                                          n_cells_max = 0,
                                          RealT = Float64)
 end
+
+@trixi_testset "elixir_euler_shockcapturing.jl (Hex, GaussSBP)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
+                        cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1),
+                        l2=[1.33894396e-02, 7.62751979e-03, 7.62751979e-03,
+                            7.62751979e-03, 4.93582917e-02],
+                        linf=[2.18776976e-01, 1.04524872e-01, 1.04524872e-01,
+                            1.04524872e-01, 8.01212059e-01])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory