From b6cacf8f7892d68a6fb0770e5d685f5b2cedaf97 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 2 Mar 2026 22:08:45 +0100 Subject: [PATCH] Gauss-Legendre DGSEM Solver: 1D Parabolic `TreeMesh` --- src/solvers/dgsem_tree/dg_1d.jl | 7 +++-- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 34 +++++++++++++++++++++++ src/solvers/dgsem_tree/dg_2d.jl | 13 ++++++--- test/test_parabolic_1d.jl | 34 +++++++++++++++++++++++ 4 files changed, 82 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index b4c231295b2..1e81efc05a7 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -714,15 +714,18 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 # into FMAs (see comment at the top of the file). @threaded for element in eachelement(dg, cache) for v in eachvariable(equations) + # Aliases for repeatedly accessed variables + surface_flux_minus = surface_flux_values[v, 1, element] + surface_flux_plus = surface_flux_values[v, 2, element] for ii in eachnode(dg) # surface at -x du[v, ii, element] = (du[v, ii, element] - - surface_flux_values[v, 1, element] * + surface_flux_minus * boundary_interpolation_inverse_weights[ii, 1]) # surface at +x du[v, ii, element] = (du[v, ii, element] + - surface_flux_values[v, 2, element] * + surface_flux_plus * boundary_interpolation_inverse_weights[ii, 2]) end end diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index ac06ec4dfac..a2d4fc445fd 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -470,6 +470,40 @@ function calc_surface_integral_gradient!(gradients, return nothing end +function calc_surface_integral_gradient!(gradients, + mesh::TreeMesh{1}, # for dispatch only + equations_parabolic::AbstractEquationsParabolic, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for v in eachvariable(equations_parabolic) + # Aliases for repeatedly accessed variables + surface_flux_minus = surface_flux_values[v, 1, element] + surface_flux_plus = surface_flux_values[v, 2, element] + for ii in eachnode(dg) + # surface at -x + gradients[v, ii, element] = (gradients[v, ii, element] - + surface_flux_minus * + boundary_interpolation_inverse_weights[ii, + 1]) + + # surface at +x + gradients[v, ii, element] = (gradients[v, ii, element] + + surface_flux_plus * + boundary_interpolation_inverse_weights[ii, + 2]) + end + end + end + + return nothing +end + # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7d2efbb7238..49861ee5946 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1197,30 +1197,35 @@ function calc_surface_integral!(du, u, @threaded for element in eachelement(dg, cache) for l in eachnode(dg) for v in eachvariable(equations) + # Aliases for repeatedly accessed variables + surface_flux_minus = surface_flux_values[v, l, 1, element] + surface_flux_plus = surface_flux_values[v, l, 2, element] for ii in eachnode(dg) # surface at -x du[v, ii, l, element] = (du[v, ii, l, element] - - surface_flux_values[v, l, 1, element] * + surface_flux_minus * boundary_interpolation_inverse_weights[ii, 1]) # surface at +x du[v, ii, l, element] = (du[v, ii, l, element] + - surface_flux_values[v, l, 2, element] * + surface_flux_plus * boundary_interpolation_inverse_weights[ii, 2]) end + surface_flux_minus = surface_flux_values[v, l, 3, element] + surface_flux_plus = surface_flux_values[v, l, 4, element] for jj in eachnode(dg) # surface at -y du[v, l, jj, element] = (du[v, l, jj, element] - - surface_flux_values[v, l, 3, element] * + surface_flux_minus * boundary_interpolation_inverse_weights[jj, 1]) # surface at +y du[v, l, jj, element] = (du[v, l, jj, element] + - surface_flux_values[v, l, 4, element] * + surface_flux_plus * boundary_interpolation_inverse_weights[jj, 2]) end diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index dac1035f8a0..75af0c07c60 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -37,6 +37,19 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end +@trixi_testset "TreeMesh1D: elixir_advection_diffusion_ldg.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_advection_diffusion_ldg.jl"), + solver=DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis), + tspan=(0.0, 0.4), + l2=[4.126471023759558e-6], linf=[1.4470099431229677e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) +end + @trixi_testset "TreeMesh1D: elixir_advection_diffusion_gradient_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", "elixir_advection_diffusion_gradient_source_terms.jl"), @@ -342,6 +355,27 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end +@trixi_testset "TreeMesh1D: elixir_navierstokes_viscous_shock.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_navierstokes_viscous_shock.jl"), + solver=DGSEM(polydeg = 3, surface_flux = flux_hlle, + basis_type = GaussLegendreBasis), + l2=[ + 0.00010415910094963455, + 7.569570282227496e-5, + 8.643799824895884e-5 + ], + linf=[ + 0.0004795456761867989, + 0.0003525509032139551, + 0.0004044657250887873 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) +end + @trixi_testset "TreeMesh1D: elixir_navierstokes_viscous_shock_imex.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", "elixir_navierstokes_viscous_shock_imex.jl"),