Remove boundary_interpolation matrix for DGSEM

src/auxiliary/precompile.jl

@@ -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}}

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -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
@@ -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,
@@ -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')
@@ -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,
@@ -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)

src/solvers/dgsem_p4est/dg_2d.jl

@@ -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.
@@ -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]
+    factor = inverse_weights[1] # For LGL basis: Identical to weighted boundary interpolation at x = ±1
     @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)
+                                        factor)
 
                 # surface at +x
                 du[v, nnodes(dg), l, element] = (du[v, nnodes(dg), l, element] +
                                                  surface_flux_values[v, l, 2, element] *
-                                                 factor_2)
+                                                 factor)
 
                 # surface at -y
                 du[v, l, 1, element] = (du[v, l, 1, element] +
                                         surface_flux_values[v, l, 3, element] *
-                                        factor_1)
+                                        factor)
 
                 # surface at +y
                 du[v, l, nnodes(dg), element] = (du[v, l, nnodes(dg), element] +
                                                  surface_flux_values[v, l, 4, element] *
-                                                 factor_2)
+                                                 factor)
             end
         end
     end

src/solvers/dgsem_p4est/dg_2d_parabolic.jl

@@ -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]
+    factor = inverse_weights[1] # For LGL basis: Identical to weighted boundary interpolation at x = ±1
     @threaded for element in eachelement(dg, cache)
         for l in eachnode(dg)
             for v in eachvariable(equations_parabolic)
@@ -985,7 +984,8 @@ function calc_gradient_surface_integral!(gradients,
                                                  surface_flux_values[v,
                                                                      l, 1,
                                                                      element] *
-                                                 factor_1 * normal_direction_x)
+                                                 factor *
+                                                 normal_direction_x)
 
                 # surface at +x
                 normal_direction_x, _ = get_normal_direction(2,
@@ -998,7 +998,7 @@ function calc_gradient_surface_integral!(gradients,
                                                           surface_flux_values[v,
                                                                               l, 2,
                                                                               element] *
-                                                          factor_2 *
+                                                          factor *
                                                           normal_direction_x)
 
                 # surface at -y
@@ -1012,7 +1012,8 @@ function calc_gradient_surface_integral!(gradients,
                                                  surface_flux_values[v,
                                                                      l, 3,
                                                                      element] *
-                                                 factor_1 * normal_direction_x)
+                                                 factor *
+                                                 normal_direction_x)
 
                 # surface at +y
                 normal_direction_x, _ = get_normal_direction(4,
@@ -1025,7 +1026,7 @@ function calc_gradient_surface_integral!(gradients,
                                                           surface_flux_values[v,
                                                                               l, 4,
                                                                               element] *
-                                                          factor_2 *
+                                                          factor *
                                                           normal_direction_x)
 
                 # Compute y-component of gradients
@@ -1041,7 +1042,8 @@ function calc_gradient_surface_integral!(gradients,
                                                  surface_flux_values[v,
                                                                      l, 1,
                                                                      element] *
-                                                 factor_1 * normal_direction_y)
+                                                 factor *
+                                                 normal_direction_y)
 
                 # surface at +x
                 _, normal_direction_y = get_normal_direction(2,
@@ -1054,7 +1056,7 @@ function calc_gradient_surface_integral!(gradients,
                                                           surface_flux_values[v,
                                                                               l, 2,
                                                                               element] *
-                                                          factor_2 *
+                                                          factor *
                                                           normal_direction_y)
 
                 # surface at -y
@@ -1068,7 +1070,8 @@ function calc_gradient_surface_integral!(gradients,
                                                  surface_flux_values[v,
                                                                      l, 3,
                                                                      element] *
-                                                 factor_1 * normal_direction_y)
+                                                 factor *
+                                                 normal_direction_y)
 
                 # surface at +y
                 _, normal_direction_y = get_normal_direction(4,
@@ -1081,7 +1084,7 @@ function calc_gradient_surface_integral!(gradients,
                                                           surface_flux_values[v,
                                                                               l, 4,
                                                                               element] *
-                                                          factor_2 *
+                                                          factor *
                                                           normal_direction_y)
             end
         end

src/solvers/dgsem_p4est/dg_3d.jl

@@ -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.
@@ -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]
+    factor = inverse_weights[1] # For LGL basis: Identical to weighted 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] *
+                                           factor)
 
                 # 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)
+                                                    factor)
 
                 # 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] *
+                                           factor)
 
                 # 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)
+                                                    factor)
 
                 # 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] *
+                                           factor)
 
                 # 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)
+                                                    factor)
             end
         end
     end