Skip to content

Remove boundary_interpolation matrix for DGSEM #2741

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 38 commits into from
Jan 23, 2026
Merged

Remove boundary_interpolation matrix for DGSEM #2741

Show file tree
Hide file tree
Changes from 35 commits
Commits
Show all changes
38 commits
Select commit Hold shift + click to select a range
bfc56c8
explain dhat
DanielDoehring Jan 5, 2026
2ed77f7
comment boundary matrix
DanielDoehring Jan 5, 2026
ab2f786
comments
DanielDoehring Jan 5, 2026
f86991a
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 6, 2026
8b3fbfb
rev
DanielDoehring Jan 6, 2026
6afc4f4
comments
DanielDoehring Jan 6, 2026
3e17288
rename
DanielDoehring Jan 6, 2026
624fb0a
comment
DanielDoehring Jan 6, 2026
46d7237
comments
DanielDoehring Jan 6, 2026
b55dbe5
update
DanielDoehring Jan 6, 2026
c3acb3e
comment
DanielDoehring Jan 6, 2026
24b9e12
compress boundary matrix
DanielDoehring Jan 6, 2026
30e3960
rev
DanielDoehring Jan 6, 2026
d4d8471
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 6, 2026
83128e1
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 7, 2026
d07cbdf
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 8, 2026
e80f40e
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 10, 2026
4d085b4
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 13, 2026
b1748cf
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 13, 2026
eba129b
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 14, 2026
b9b841a
Merge branch 'main' into ExplainDHat
DanielDoehring Jan 14, 2026
ea1e52e
rm comment
DanielDoehring Jan 16, 2026
262a133
one factor
DanielDoehring Jan 18, 2026
c97574a
Merge branch 'main' into BoundaryInterpolationMatrix
DanielDoehring Jan 18, 2026
66004f4
pcom
DanielDoehring Jan 18, 2026
e45a7dd
bf
DanielDoehring Jan 18, 2026
0f6e40e
fmt
DanielDoehring Jan 18, 2026
b1740cb
bf
DanielDoehring Jan 18, 2026
d754fe8
Merge branch 'main' into BoundaryInterpolationMatrix
DanielDoehring Jan 19, 2026
406a739
rename
DanielDoehring Jan 22, 2026
9344f13
test
DanielDoehring Jan 22, 2026
b0732b9
Merge branch 'main' into BoundaryInterpolationMatrix
DanielDoehring Jan 22, 2026
b3d4e5a
pm
DanielDoehring Jan 22, 2026
73d41ba
comm
DanielDoehring Jan 22, 2026
17c43b1
Merge branch 'main' into BoundaryInterpolationMatrix
DanielDoehring Jan 23, 2026
250eb1c
factor
DanielDoehring Jan 23, 2026
bcc1f60
Merge branch 'BoundaryInterpolationMatrix' of github.com:DanielDoehri…
DanielDoehring Jan 23, 2026
4fe577c
Merge branch 'main' into BoundaryInterpolationMatrix
DanielDoehring Jan 23, 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
3 changes: 0 additions & 3 deletions src/auxiliary/precompile.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,9 +181,6 @@ function _precompile_manual_()
StaticArrays.SVector{nnodes_, RealT},
# InverseVandermondeLegendre
Matrix{RealT},
# BoundaryMatrix
#StaticArrays.SArray{Tuple{nnodes_,2},RealT,2,2*nnodes_},
Matrix{RealT},
# DerivativeMatrix
#StaticArrays.SArray{Tuple{nnodes_,nnodes_},RealT,2,nnodes_^2},
Matrix{RealT}}
Expand Down
13 changes: 2 additions & 11 deletions src/solvers/dgsem/basis_lobatto_legendre.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,15 +17,13 @@ This exceptional case is currently only supported for TreeMesh!
struct LobattoLegendreBasis{RealT <: Real, NNODES,
VectorT <: AbstractVector{RealT},
InverseVandermondeLegendre <: AbstractMatrix{RealT},
BoundaryMatrix <: AbstractMatrix{RealT},
DerivativeMatrix <: AbstractMatrix{RealT}} <:
AbstractBasisSBP{RealT}
nodes::VectorT
weights::VectorT
inverse_weights::VectorT

inverse_vandermonde_legendre::InverseVandermondeLegendre
boundary_interpolation::BoundaryMatrix # Compressed form of "Lhat"

derivative_matrix::DerivativeMatrix # strong form derivative matrix "D"
derivative_split::DerivativeMatrix # strong form derivative matrix minus boundary terms
Expand All @@ -40,19 +38,16 @@ function Adapt.adapt_structure(to, basis::LobattoLegendreBasis)
nodes = SVector{<:Any, RealT}(basis.nodes)
weights = SVector{<:Any, RealT}(basis.weights)
inverse_weights = SVector{<:Any, RealT}(basis.inverse_weights)
boundary_interpolation = adapt(to, basis.boundary_interpolation)
derivative_matrix = adapt(to, basis.derivative_matrix)
derivative_split = adapt(to, basis.derivative_split)
derivative_split_transpose = adapt(to, basis.derivative_split_transpose)
derivative_hat = adapt(to, basis.derivative_hat)
return LobattoLegendreBasis{RealT, nnodes(basis), typeof(nodes),
typeof(inverse_vandermonde_legendre),
typeof(boundary_interpolation),
typeof(derivative_matrix)}(nodes,
weights,
inverse_weights,
inverse_vandermonde_legendre,
boundary_interpolation,
derivative_matrix,
derivative_split,
derivative_split_transpose,
Expand All @@ -67,10 +62,6 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)

_, inverse_vandermonde_legendre = vandermonde_legendre(nodes_, RealT)

boundary_interpolation = zeros(RealT, nnodes_, 2)
boundary_interpolation[:, 1] = calc_Lhat(-one(RealT), nodes_, weights_)
boundary_interpolation[:, 2] = calc_Lhat(one(RealT), nodes_, weights_)

derivative_matrix = polynomial_derivative_matrix(nodes_)
derivative_split = calc_Dsplit(nodes_, weights_)
derivative_split_transpose = Matrix(derivative_split')
Expand All @@ -89,11 +80,9 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)

return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes),
typeof(inverse_vandermonde_legendre),
typeof(boundary_interpolation),
typeof(derivative_matrix)}(nodes, weights,
inverse_weights,
inverse_vandermonde_legendre,
boundary_interpolation,
derivative_matrix,
derivative_split,
derivative_split_transpose,
Expand Down Expand Up @@ -531,6 +520,8 @@ end

# Calculate M^{-1} * L(x), where L(x) is the Lagrange polynomial
# vector at point x.
# Not required for the DGSEM with LGL basis, as boundary evaluations
# collapse to boundary node evaluations.
function calc_Lhat(x, nodes, weights)
n_nodes = length(nodes)
wbary = barycentric_weights(nodes)
Expand Down
13 changes: 6 additions & 7 deletions src/solvers/dgsem_p4est/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -761,7 +761,7 @@ function calc_surface_integral!(du, u,
equations,
surface_integral::SurfaceIntegralWeakForm,
dg::DGSEM, cache)
@unpack boundary_interpolation = dg.basis
@unpack inverse_weights = dg.basis
@unpack surface_flux_values = cache.elements

# Note that all fluxes have been computed with outward-pointing normal vectors.
Expand All @@ -771,30 +771,29 @@ function calc_surface_integral!(du, u,
#
# We also use explicit assignments instead of `+=` to let `@muladd` turn these
# into FMAs (see comment at the top of the file).
factor_1 = boundary_interpolation[1, 1]
factor_2 = boundary_interpolation[nnodes(dg), 2]
boundary_interpolation = inverse_weights[1] # For LGL basis: Identical to boundary interpolation at x = ±1
Comment thread
ranocha marked this conversation as resolved.
Outdated
@threaded for element in eachelement(dg, cache)
for l in eachnode(dg)
for v in eachvariable(equations)
# surface at -x
du[v, 1, l, element] = (du[v, 1, l, element] +
surface_flux_values[v, l, 1, element] *
factor_1)
boundary_interpolation)
Comment thread
DanielDoehring marked this conversation as resolved.
Outdated

# surface at +x
du[v, nnodes(dg), l, element] = (du[v, nnodes(dg), l, element] +
surface_flux_values[v, l, 2, element] *
factor_2)
boundary_interpolation)
Comment thread
DanielDoehring marked this conversation as resolved.
Outdated

# surface at -y
du[v, l, 1, element] = (du[v, l, 1, element] +
surface_flux_values[v, l, 3, element] *
factor_1)
boundary_interpolation)
Comment thread
DanielDoehring marked this conversation as resolved.
Outdated

# surface at +y
du[v, l, nnodes(dg), element] = (du[v, l, nnodes(dg), element] +
surface_flux_values[v, l, 4, element] *
factor_2)
boundary_interpolation)
Comment thread
DanielDoehring marked this conversation as resolved.
Outdated
end
end
end
Expand Down
25 changes: 14 additions & 11 deletions src/solvers/dgsem_p4est/dg_2d_parabolic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -958,16 +958,15 @@ function calc_gradient_surface_integral!(gradients,
mesh::P4estMesh{2}, # for dispatch only
equations_parabolic::AbstractEquationsParabolic,
dg::DGSEM, cache)
@unpack boundary_interpolation = dg.basis
@unpack inverse_weights = dg.basis
@unpack surface_flux_values = cache.elements
@unpack contravariant_vectors = cache.elements

gradients_x, gradients_y = gradients

# We also use explicit assignments instead of `+=` to let `@muladd` turn these
# into FMAs (see comment at the top of the file).
factor_1 = boundary_interpolation[1, 1]
factor_2 = boundary_interpolation[nnodes(dg), 2]
boundary_interpolation = inverse_weights[1] # For LGL basis: Identical to boundary interpolation at x = ±1
@threaded for element in eachelement(dg, cache)
for l in eachnode(dg)
for v in eachvariable(equations_parabolic)
Expand All @@ -985,7 +984,8 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 1,
element] *
factor_1 * normal_direction_x)
boundary_interpolation *
normal_direction_x)

# surface at +x
normal_direction_x, _ = get_normal_direction(2,
Expand All @@ -998,7 +998,7 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 2,
element] *
factor_2 *
boundary_interpolation *
normal_direction_x)

# surface at -y
Expand All @@ -1012,7 +1012,8 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 3,
element] *
factor_1 * normal_direction_x)
boundary_interpolation *
normal_direction_x)

# surface at +y
normal_direction_x, _ = get_normal_direction(4,
Expand All @@ -1025,7 +1026,7 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 4,
element] *
factor_2 *
boundary_interpolation *
normal_direction_x)

# Compute y-component of gradients
Expand All @@ -1041,7 +1042,8 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 1,
element] *
factor_1 * normal_direction_y)
boundary_interpolation *
normal_direction_y)

# surface at +x
_, normal_direction_y = get_normal_direction(2,
Expand All @@ -1054,7 +1056,7 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 2,
element] *
factor_2 *
boundary_interpolation *
normal_direction_y)

# surface at -y
Expand All @@ -1068,7 +1070,8 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 3,
element] *
factor_1 * normal_direction_y)
boundary_interpolation *
normal_direction_y)

# surface at +y
_, normal_direction_y = get_normal_direction(4,
Expand All @@ -1081,7 +1084,7 @@ function calc_gradient_surface_integral!(gradients,
surface_flux_values[v,
l, 4,
element] *
factor_2 *
boundary_interpolation *
normal_direction_y)
end
end
Expand Down
26 changes: 14 additions & 12 deletions src/solvers/dgsem_p4est/dg_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -931,7 +931,7 @@ function calc_surface_integral!(du, u,
equations,
surface_integral::SurfaceIntegralWeakForm,
dg::DGSEM, cache)
@unpack boundary_interpolation = dg.basis
@unpack inverse_weights = dg.basis
@unpack surface_flux_values = cache.elements

# Note that all fluxes have been computed with outward-pointing normal vectors.
Expand All @@ -941,43 +941,45 @@ function calc_surface_integral!(du, u,
#
# We also use explicit assignments instead of `+=` to let `@muladd` turn these
# into FMAs (see comment at the top of the file).
factor_1 = boundary_interpolation[1, 1]
factor_2 = boundary_interpolation[nnodes(dg), 2]
boundary_interpolation = inverse_weights[1] # For LGL basis: Identical to boundary interpolation at x = ±1
@threaded for element in eachelement(dg, cache)
for m in eachnode(dg), l in eachnode(dg)
for v in eachvariable(equations)
# surface at -x
du[v, 1, l, m, element] = (du[v, 1, l, m, element] +
surface_flux_values[v, l, m, 1, element] *
factor_1)
surface_flux_values[v, l, m, 1,
element] *
boundary_interpolation)

# surface at +x
du[v, nnodes(dg), l, m, element] = (du[v, nnodes(dg), l, m, element] +
surface_flux_values[v, l, m, 2,
element] *
factor_2)
boundary_interpolation)

# surface at -y
du[v, l, 1, m, element] = (du[v, l, 1, m, element] +
surface_flux_values[v, l, m, 3, element] *
factor_1)
surface_flux_values[v, l, m, 3,
element] *
boundary_interpolation)

# surface at +y
du[v, l, nnodes(dg), m, element] = (du[v, l, nnodes(dg), m, element] +
surface_flux_values[v, l, m, 4,
element] *
factor_2)
boundary_interpolation)

# surface at -z
du[v, l, m, 1, element] = (du[v, l, m, 1, element] +
surface_flux_values[v, l, m, 5, element] *
factor_1)
surface_flux_values[v, l, m, 5,
element] *
boundary_interpolation)

# surface at +z
du[v, l, m, nnodes(dg), element] = (du[v, l, m, nnodes(dg), element] +
surface_flux_values[v, l, m, 6,
element] *
factor_2)
boundary_interpolation)
end
end
end
Expand Down
Loading
Loading