Intermediate PR for #1930

examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl

@@ -0,0 +1,68 @@
+using Trixi
+using OrdinaryDiffEq
+
+equations = InviscidBurgersEquation1D()
+
+###############################################################################
+# setup the GSBP DG discretization that uses the Gauss operators from 
+# Chan, Del Rey Fernandez, Carpenter (2019). 
+# [https://doi.org/10.1137/18M1209234](https://doi.org/10.1137/18M1209234)
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ec
+
+polydeg = 3
+basis = DGMultiBasis(Line(), polydeg, approximation_type = GaussSBP())
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = first)
+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)
+
+###############################################################################
+#  setup the 1D mesh
+
+cells_per_dimension = (32,)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (-1.0,), coordinates_max = (1.0,),
+                   periodicity = true)
+
+###############################################################################
+#  setup the semidiscretization and ODE problem
+
+semi = SemidiscretizationHyperbolic(mesh,
+                                    equations,
+                                    initial_condition_convergence_test,
+                                    dg)
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+#  setup the callbacks
+
+# prints a summary of the simulation setup and resets the timers
+summary_callback = SummaryCallback()
+
+# analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100, uEltype = real(dg))
+
+# handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+# collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+
+# ###############################################################################
+# # run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, save_everystep = false, callback = callbacks);

src/equations/inviscid_burgers_1d.jl

@@ -168,6 +168,7 @@ end
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::InviscidBurgersEquation1D) = u
+@inline entropy2cons(u, equation::InviscidBurgersEquation1D) = u
 
 # Calculate entropy for a conservative state `cons`
 @inline entropy(u::Real, ::InviscidBurgersEquation1D) = 0.5 * u^2

src/solvers/dgmulti/flux_differencing_gauss_sbp.jl

@@ -51,6 +51,29 @@ struct TensorProductGaussFaceOperator{NDIMS, OperatorType <: AbstractGaussOperat
     nfaces::Int
 end
 
+function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
+                                        dg::DGMulti{1, Line, GaussSBP})
+    rd = dg.basis
+
+    rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg))
+    interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1])
+
+    nnodes_1d = length(rq1D)
+    face_indices_tensor_product = nothing # not needed in 1D; we fall back to mul!
+
+    num_faces = 2
+
+    T_op = typeof(operator)
+    Tm = typeof(interp_matrix_gauss_to_face_1d)
+    Tw = typeof(inv.(wq1D))
+    Tf = typeof(rd.wf)
+    Ti = typeof(face_indices_tensor_product)
+    return TensorProductGaussFaceOperator{1, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d,
+                                                                   inv.(wq1D), rd.wf,
+                                                                   face_indices_tensor_product,
+                                                                   nnodes_1d, num_faces)
+end
+
 # constructor for a 2D operator
 function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
                                         dg::DGMulti{2, Quad, GaussSBP})
@@ -126,6 +149,21 @@ end
     end
 end
 
+@inline function tensor_product_gauss_face_operator!(out::AbstractVector,
+                                                     A::TensorProductGaussFaceOperator{1,
+                                                                                       Interpolation},
+                                                     x::AbstractVector)
+    mul!(out, A.interp_matrix_gauss_to_face_1d, x)
+end
+
+@inline function tensor_product_gauss_face_operator!(out::AbstractVector,
+                                                     A::TensorProductGaussFaceOperator{1,
+                                                                                       <:Projection},
+                                                     x::AbstractVector)
+    mul!(out, A.interp_matrix_gauss_to_face_1d', x)
+    @. out *= A.inv_volume_weights_1d
+end
+
 # By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
 # Since these FMAs can increase the performance of many numerical algorithms,
 # we need to opt-in explicitly.
@@ -352,7 +390,7 @@ end
 # For now, this is mostly the same as `create_cache` for DGMultiFluxDiff{<:Polynomial}.
 # In the future, we may modify it so that we can specialize additional parts of GaussSBP() solvers.
 function create_cache(mesh::DGMultiMesh, equations,
-                      dg::DGMultiFluxDiff{<:GaussSBP, <:Union{Quad, Hex}}, RealT,
+                      dg::DGMultiFluxDiff{<:GaussSBP, <:Union{Line, Quad, Hex}}, RealT,
                       uEltype)
 
     # call general Polynomial flux differencing constructor

src/solvers/dgmulti/types.jl

@@ -347,6 +347,9 @@ function SimpleKronecker(NDIMS, A, eltype_A = eltype(A))
     return SimpleKronecker{NDIMS, typeof(A), typeof(tmp_storage)}(A, tmp_storage)
 end
 
+# fall back to mul! for a 1D Kronecker product
+LinearAlgebra.mul!(b, A_kronecker::SimpleKronecker{1}, x) = mul!(b, A_kronecker.A, x)
+
 # Computes `b = kron(A, A) * x` in an optimized fashion
 function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{2}, x_in)
     @unpack A = A_kronecker

test/test_dgmulti_1d.jl

@@ -29,6 +29,22 @@ isdir(outdir) && rm(outdir, recursive = true)
     end
 end
 
+@trixi_testset "elixir_burgers_gauss_shock_capturing.jl " begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_burgers_gauss_shock_capturing.jl"),
+                        cells_per_dimension=(8,), tspan=(0.0, 0.1),
+                        l2=[0.445804588167854],
+                        linf=[0.74780611426038])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_euler_flux_diff.jl " begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_flux_diff.jl"),
                         cells_per_dimension=(16,),