Skip to content

Gauss-Legendre DGSEM Solver: 1D Parabolic TreeMesh #2842

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ranocha merged 2 commits into from
Mar 3, 2026
Merged

Gauss-Legendre DGSEM Solver: 1D Parabolic TreeMesh #2842

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b6cacf8
Gauss-Legendre DGSEM Solver: 1D Parabolic `TreeMesh`
DanielDoehring Mar 2, 2026
6c3ca3b
Merge branch 'main' into GL_Parabolic_TreeMesh_1D_Periodic
ranocha Mar 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions src/solvers/dgsem_tree/dg_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Comment thread
ranocha marked this conversation as resolved.
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
Expand Down
34 changes: 34 additions & 0 deletions src/solvers/dgsem_tree/dg_1d_parabolic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
13 changes: 9 additions & 4 deletions src/solvers/dgsem_tree/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
34 changes: 34 additions & 0 deletions test/test_parabolic_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"),
Expand Down Expand Up @@ -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"),
Expand Down
Loading