From a0f52b10090d3b472180ddd66f1553a98a8c6746 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 08:45:06 -0500 Subject: [PATCH 1/7] add 1D Gauss tensor product functionality --- .../dgmulti/flux_differencing_gauss_sbp.jl | 38 ++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 9059caf87f6..3e1a3edf67d 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,19 @@ 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 +388,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 From ea159e9d9d39c12de468a4e2e718a468ef26d604 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 08:45:15 -0500 Subject: [PATCH 2/7] add 1D kronecker product fallback --- src/solvers/dgmulti/types.jl | 3 +++ 1 file changed, 3 insertions(+) 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 From d73b746b2a3929b1dcc6a7b6926c84c6169fe1e8 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 08:54:26 -0500 Subject: [PATCH 3/7] add Burgers' shock capturing example --- .../elixir_burgers_gauss_shock_capturing.jl | 68 +++++++++++++++++++ src/equations/inviscid_burgers_1d.jl | 1 + 2 files changed, 69 insertions(+) create mode 100644 examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl 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..d08ba4cc776 --- /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 et al. + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ec + +polydeg = 3 +basis = DGMultiBasis(Line(), polydeg, approximation_type = GaussSBP()) + +indicator_function(u, ::InviscidBurgersEquation1D) = u[1] + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = indicator_function) +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); \ No newline at end of file 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 From 59626cf353c12ae2dfeae8b56d12dd8484b22ffc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 08:54:31 -0500 Subject: [PATCH 4/7] add test --- test/test_dgmulti_1d.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl index e470de71efb..0eb709f03a9 100644 --- a/test/test_dgmulti_1d.jl +++ b/test/test_dgmulti_1d.jl @@ -29,6 +29,21 @@ 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,), From 3f2da7f6d3f85aaca893a3ddd47dfe3b609b701f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 7 May 2024 16:07:20 +0200 Subject: [PATCH 5/7] format --- .../dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl | 2 +- src/solvers/dgmulti/flux_differencing_gauss_sbp.jl | 6 ++++-- test/test_dgmulti_1d.jl | 7 ++++--- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl index d08ba4cc776..7921d8a5d31 100644 --- a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl +++ b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl @@ -65,4 +65,4 @@ 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); \ No newline at end of file + dt = 1.0, save_everystep = false, callback = callbacks); diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 3e1a3edf67d..63a37f6780b 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -150,13 +150,15 @@ end end @inline function tensor_product_gauss_face_operator!(out::AbstractVector, - A::TensorProductGaussFaceOperator{1, Interpolation}, + 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}, + A::TensorProductGaussFaceOperator{1, + <:Projection}, x::AbstractVector) mul!(out, A.interp_matrix_gauss_to_face_1d', x) @. out *= A.inv_volume_weights_1d diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl index 0eb709f03a9..0d083cf9a72 100644 --- a/test/test_dgmulti_1d.jl +++ b/test/test_dgmulti_1d.jl @@ -30,10 +30,11 @@ isdir(outdir) && rm(outdir, recursive = true) end @trixi_testset "elixir_burgers_gauss_shock_capturing.jl " begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_gauss_shock_capturing.jl"), + @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]) + l2=[0.445804588167854], + linf=[0.74780611426038]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From b61275accd27c75381b217d0194f656fbdda7141 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 7 May 2024 09:22:41 -0500 Subject: [PATCH 6/7] Update examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl Co-authored-by: Hendrik Ranocha --- examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl index 7921d8a5d31..38a95b34743 100644 --- a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl +++ b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl @@ -12,13 +12,11 @@ volume_flux = flux_ec polydeg = 3 basis = DGMultiBasis(Line(), polydeg, approximation_type = GaussSBP()) -indicator_function(u, ::InviscidBurgersEquation1D) = u[1] - indicator_sc = IndicatorHennemannGassner(equations, basis, alpha_max = 0.5, alpha_min = 0.001, alpha_smooth = true, - variable = indicator_function) + variable = first) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) From 02e670245008b57e39d54bf51889db167e99b552 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 09:23:45 -0500 Subject: [PATCH 7/7] add DOI --- examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl index d08ba4cc776..8891267413d 100644 --- a/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl +++ b/examples/dgmulti_1d/elixir_burgers_gauss_shock_capturing.jl @@ -4,7 +4,9 @@ using OrdinaryDiffEq equations = InviscidBurgersEquation1D() ############################################################################### -# setup the GSBP DG discretization that uses the Gauss operators from Chan et al. +# 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