diff --git a/examples/dgmulti_2d/elixir_euler_shockcapturing.jl b/examples/dgmulti_2d/elixir_euler_shockcapturing.jl index 124e53cd6ae..11a8e1b035b 100644 --- a/examples/dgmulti_2d/elixir_euler_shockcapturing.jl +++ b/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 diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl index 5b3e6c68955..0abe70ffd09 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/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}` diff --git a/src/solvers/dgmulti/shock_capturing.jl b/src/solvers/dgmulti/shock_capturing.jl index 25e1150c5bc..5927f76db2d 100644 --- a/src/solvers/dgmulti/shock_capturing.jl +++ b/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 diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index e99eaf54bba..4a561f1c73b 100644 --- a/test/test_dgmulti_2d.jl +++ b/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)