diff --git a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl new file mode 100644 index 00000000000..b0632b1f978 --- /dev/null +++ b/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); diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index f2387f26ba7..130196a4929 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/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 diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 9059caf87f6..63a37f6780b 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/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 diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index 813bc67061e..ef9d7d2bf09 100644 --- a/src/solvers/dgmulti/types.jl +++ b/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 diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl index e470de71efb..0d083cf9a72 100644 --- a/test/test_dgmulti_1d.jl +++ b/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,),