Rename derivative_dhat => derivative_hat, explain minus signs

src/auxiliary/precompile.jl

@@ -181,9 +181,8 @@ function _precompile_manual_()
                                     StaticArrays.SVector{nnodes_, RealT},
                                     # InverseVandermondeLegendre
                                     Matrix{RealT},
-                                    # BoundaryMatrix
-                                    #StaticArrays.SArray{Tuple{nnodes_,2},RealT,2,2*nnodes_},
-                                    Matrix{RealT},
+                                    # VectorT2, EntriesBoundaryMatrix
+                                    StaticArrays.SVector{2, RealT},
                                     # DerivativeMatrix
                                     #StaticArrays.SArray{Tuple{nnodes_,nnodes_},RealT,2,nnodes_^2},
                                     Matrix{RealT}}
@@ -309,16 +308,16 @@ function _precompile_manual_()
     Base.precompile(Tuple{Type{LobattoLegendreBasis}, Int})
     for RealT in (Float64,)
         Base.precompile(Tuple{Type{LobattoLegendreBasis}, RealT, Int})
-        @assert Base.precompile(Tuple{typeof(Trixi.calc_dhat), Vector{RealT},
+        @assert Base.precompile(Tuple{typeof(Trixi.calc_Dhat), Vector{RealT},
                                       Vector{RealT}})
-        @assert Base.precompile(Tuple{typeof(Trixi.calc_dsplit), Vector{RealT},
+        @assert Base.precompile(Tuple{typeof(Trixi.calc_Dsplit), Vector{RealT},
                                       Vector{RealT}})
         @assert Base.precompile(Tuple{typeof(Trixi.polynomial_derivative_matrix),
                                       Vector{RealT}})
         @assert Base.precompile(Tuple{typeof(Trixi.polynomial_interpolation_matrix),
                                       Vector{RealT}, Vector{RealT}})
         @assert Base.precompile(Tuple{typeof(Trixi.barycentric_weights), Vector{RealT}})
-        @assert Base.precompile(Tuple{typeof(Trixi.calc_lhat), RealT, Vector{RealT},
+        @assert Base.precompile(Tuple{typeof(Trixi.calc_Lhat), RealT, Vector{RealT},
                                       Vector{RealT}})
         @assert Base.precompile(Tuple{typeof(Trixi.lagrange_interpolating_polynomials),
                                       RealT, Vector{RealT}, Vector{RealT}})

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -17,21 +17,20 @@ This exceptional case is currently only supported for TreeMesh!
 struct LobattoLegendreBasis{RealT <: Real, NNODES,
                             VectorT <: AbstractVector{RealT},
                             InverseVandermondeLegendre <: AbstractMatrix{RealT},
-                            BoundaryMatrix <: AbstractMatrix{RealT},
+                            VectorT2 <: AbstractVector{RealT},
                             DerivativeMatrix <: AbstractMatrix{RealT}} <:
        AbstractBasisSBP{RealT}
     nodes::VectorT
     weights::VectorT
     inverse_weights::VectorT
 
     inverse_vandermonde_legendre::InverseVandermondeLegendre
-    boundary_interpolation::BoundaryMatrix # lhat
+    boundary_interpolation_entries::VectorT2 # Compressed form of "Lhat"
 
-    derivative_matrix::DerivativeMatrix # strong form derivative matrix
+    derivative_matrix::DerivativeMatrix # strong form derivative matrix "D"
     derivative_split::DerivativeMatrix # strong form derivative matrix minus boundary terms
     derivative_split_transpose::DerivativeMatrix # transpose of `derivative_split`
-    derivative_dhat::DerivativeMatrix # weak form matrix "dhat",
-    # negative adjoint wrt the SBP dot product
+    derivative_hat::DerivativeMatrix # weak form matrix "Dhat", negative adjoint wrt the SBP dot product
 end
 
 function Adapt.adapt_structure(to, basis::LobattoLegendreBasis)
@@ -41,23 +40,23 @@ 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)
+    boundary_interpolation_entries = adapt(to, basis.boundary_interpolation_entries)
     derivative_matrix = adapt(to, basis.derivative_matrix)
     derivative_split = adapt(to, basis.derivative_split)
     derivative_split_transpose = adapt(to, basis.derivative_split_transpose)
-    derivative_dhat = adapt(to, basis.derivative_dhat)
+    derivative_hat = adapt(to, basis.derivative_hat)
     return LobattoLegendreBasis{RealT, nnodes(basis), typeof(nodes),
                                 typeof(inverse_vandermonde_legendre),
-                                typeof(boundary_interpolation),
+                                typeof(boundary_interpolation_entries),
                                 typeof(derivative_matrix)}(nodes,
                                                            weights,
                                                            inverse_weights,
                                                            inverse_vandermonde_legendre,
-                                                           boundary_interpolation,
+                                                           boundary_interpolation_entries,
                                                            derivative_matrix,
                                                            derivative_split,
                                                            derivative_split_transpose,
-                                                           derivative_dhat)
+                                                           derivative_hat)
 end
 
 function LobattoLegendreBasis(RealT, polydeg::Integer)
@@ -69,13 +68,17 @@ 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_)
+    boundary_interpolation[:, 1] = calc_Lhat(-one(RealT), nodes_, weights_)
+    boundary_interpolation[:, 2] = calc_Lhat(one(RealT), nodes_, weights_)
+    # Currently only used parts in the code
+    boundary_interpolation_entries = zeros(RealT, 2)
+    boundary_interpolation_entries[1] = boundary_interpolation[1, 1]
+    boundary_interpolation_entries[2] = boundary_interpolation[end, 2]
 
     derivative_matrix = polynomial_derivative_matrix(nodes_)
-    derivative_split = calc_dsplit(nodes_, weights_)
+    derivative_split = calc_Dsplit(nodes_, weights_)
     derivative_split_transpose = Matrix(derivative_split')
-    derivative_dhat = calc_dhat(nodes_, weights_)
+    derivative_hat = calc_Dhat(nodes_, weights_)
 
     # Type conversions to enable possible optimizations of runtime performance
     # and latency
@@ -90,15 +93,15 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)
 
     return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes),
                                 typeof(inverse_vandermonde_legendre),
-                                typeof(boundary_interpolation),
+                                typeof(boundary_interpolation_entries),
                                 typeof(derivative_matrix)}(nodes, weights,
                                                            inverse_weights,
                                                            inverse_vandermonde_legendre,
-                                                           boundary_interpolation,
+                                                           boundary_interpolation_entries,
                                                            derivative_matrix,
                                                            derivative_split,
                                                            derivative_split_transpose,
-                                                           derivative_dhat)
+                                                           derivative_hat)
 end
 LobattoLegendreBasis(polydeg::Integer) = LobattoLegendreBasis(Float64, polydeg)
 
@@ -410,45 +413,50 @@ end
 
 # TODO: Taal refactor, allow other RealT below and adapt constructors above accordingly
 
-# Calculate the Dhat matrix
-function calc_dhat(nodes, weights)
+# Calculate the Dhat matrix = -M^{-1} D^T M for weak form differentiation.
+# Note that this is the negated version of the matrix that shows up on the RHS of the
+# DG update multiplying the physical flux evaluations.
+function calc_Dhat(nodes, weights)
     n_nodes = length(nodes)
-    dhat = Matrix(polynomial_derivative_matrix(nodes)')
+    Dhat = Matrix(polynomial_derivative_matrix(nodes)')
 
+    # Perform M matrix multplicaitons and negate
     for n in 1:n_nodes, j in 1:n_nodes
-        dhat[j, n] *= -weights[n] / weights[j]
+        Dhat[j, n] *= -weights[n] / weights[j]
     end
 
-    return dhat
+    return Dhat
 end
 
-# Calculate the Dsplit matrix for split-form differentiation: dplit = 2D - M⁻¹B
-function calc_dsplit(nodes, weights)
+# Calculate the Dsplit matrix for split-form differentiation: Dsplit = 2D - M⁻¹B
+# Note that this is the negated version of the matrix that shows up on the RHS of the 
+# DG update multiplying the two-point numerical volume flux evaluations.
+function calc_Dsplit(nodes, weights)
     # Start with 2 x the normal D matrix
-    dsplit = 2 .* polynomial_derivative_matrix(nodes)
+    Dsplit = 2 .* polynomial_derivative_matrix(nodes)
 
-    # Modify to account for
-    dsplit[1, 1] += 1 / weights[1]
-    dsplit[end, end] -= 1 / weights[end]
+    # Modify to account for the weighted boundary terms
+    Dsplit[1, 1] += 1 / weights[1] # B[1, 1] = -1
+    Dsplit[end, end] -= 1 / weights[end] # B[end, end] = 1
 
-    return dsplit
+    return Dsplit
 end
 
 # Calculate the polynomial derivative matrix D.
 # This implements algorithm 37 "PolynomialDerivativeMatrix" from Kopriva's book.
 function polynomial_derivative_matrix(nodes)
     n_nodes = length(nodes)
-    d = zeros(eltype(nodes), n_nodes, n_nodes)
+    D = zeros(eltype(nodes), n_nodes, n_nodes)
     wbary = barycentric_weights(nodes)
 
     for i in 1:n_nodes, j in 1:n_nodes
         if j != i
-            d[i, j] = (wbary[j] / wbary[i]) * 1 / (nodes[i] - nodes[j])
-            d[i, i] -= d[i, j]
+            D[i, j] = (wbary[j] / wbary[i]) * 1 / (nodes[i] - nodes[j])
+            D[i, i] -= D[i, j]
         end
     end
 
-    return d
+    return D
 end
 
 # Calculate and interpolation matrix (Vandermonde matrix) between two given sets of nodes
@@ -525,18 +533,20 @@ function barycentric_weights(nodes)
     return weights
 end
 
-# Calculate Lhat.
-function calc_lhat(x, nodes, weights)
+# Calculate Lhat = M^{-1} * B.
+# Note that this is the negated version of the matrix that shows up on the RHS of the
+# DG update multiplying the two-point numerical surface flux evaluations.
+function calc_Lhat(x, nodes, weights)
     n_nodes = length(nodes)
     wbary = barycentric_weights(nodes)
 
-    lhat = lagrange_interpolating_polynomials(x, nodes, wbary)
+    Lhat = lagrange_interpolating_polynomials(x, nodes, wbary)
 
     for i in 1:n_nodes
-        lhat[i] /= weights[i]
+        Lhat[i] /= weights[i]
     end
 
-    return lhat
+    return Lhat
 end
 
 """

src/solvers/dgsem_p4est/dg_2d.jl

@@ -761,15 +761,15 @@ function calc_surface_integral!(du, u,
                                 equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGSEM, cache)
-    @unpack boundary_interpolation = dg.basis
+    @unpack boundary_interpolation_entries = dg.basis
     @unpack surface_flux_values = cache.elements
 
     # Note that all fluxes have been computed with outward-pointing normal vectors.
     # Access the factors only once before beginning the loop to increase performance.
     # 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_1 = boundary_interpolation_entries[1]
+    factor_2 = boundary_interpolation_entries[2]
     @threaded for element in eachelement(dg, cache)
         for l in eachnode(dg)
             for v in eachvariable(equations)

src/solvers/dgsem_p4est/dg_2d_parabolic.jl

@@ -367,7 +367,7 @@ function calc_volume_integral!(du, flux_viscous,
                                mesh::P4estMesh{2},
                                equations_parabolic::AbstractEquationsParabolic,
                                dg::DGSEM, cache)
-    (; derivative_dhat) = dg.basis
+    (; derivative_hat) = dg.basis
     (; contravariant_vectors) = cache.elements
     flux_viscous_x, flux_viscous_y = flux_viscous
 
@@ -385,7 +385,7 @@ function calc_volume_integral!(du, flux_viscous,
                                                   i, j, element)
             contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
             for ii in eachnode(dg)
-                multiply_add_to_node_vars!(du, derivative_dhat[ii, i],
+                multiply_add_to_node_vars!(du, derivative_hat[ii, i],
                                            contravariant_flux1,
                                            equations_parabolic, dg, ii, j, element)
             end
@@ -396,7 +396,7 @@ function calc_volume_integral!(du, flux_viscous,
                                                   i, j, element)
             contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2
             for jj in eachnode(dg)
-                multiply_add_to_node_vars!(du, derivative_dhat[jj, j],
+                multiply_add_to_node_vars!(du, derivative_hat[jj, j],
                                            contravariant_flux2,
                                            equations_parabolic, dg, i, jj, element)
             end
@@ -772,7 +772,7 @@ function calc_gradient_volume_integral!(gradients, u_transformed,
                                         mesh::P4estMesh{2}, # for dispatch only
                                         equations_parabolic::AbstractEquationsParabolic,
                                         dg::DG, cache)
-    @unpack derivative_dhat = dg.basis
+    @unpack derivative_hat = dg.basis
     @unpack contravariant_vectors = cache.elements
     gradients_x, gradients_y = gradients
 
@@ -784,13 +784,13 @@ function calc_gradient_volume_integral!(gradients, u_transformed,
                                    i, j, element)
 
             for ii in eachnode(dg)
-                multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i],
+                multiply_add_to_node_vars!(gradients_x, derivative_hat[ii, i],
                                            u_node, equations_parabolic, dg,
                                            ii, j, element)
             end
 
             for jj in eachnode(dg)
-                multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j],
+                multiply_add_to_node_vars!(gradients_y, derivative_hat[jj, j],
                                            u_node, equations_parabolic, dg,
                                            i, jj, element)
             end
@@ -958,7 +958,7 @@ function calc_gradient_surface_integral!(gradients,
                                          mesh::P4estMesh{2}, # for dispatch only
                                          equations_parabolic::AbstractEquationsParabolic,
                                          dg::DGSEM, cache)
-    @unpack boundary_interpolation = dg.basis
+    @unpack boundary_interpolation_entries = dg.basis
     @unpack surface_flux_values = cache.elements
     @unpack contravariant_vectors = cache.elements
 
@@ -967,8 +967,8 @@ function calc_gradient_surface_integral!(gradients,
     # Access the factors only once before beginning the loop to increase performance.
     # 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_1 = boundary_interpolation_entries[1]
+    factor_2 = boundary_interpolation_entries[2]
     @threaded for element in eachelement(dg, cache)
         for l in eachnode(dg)
             for v in eachvariable(equations_parabolic)

src/solvers/dgsem_p4est/dg_3d.jl

@@ -931,15 +931,15 @@ function calc_surface_integral!(du, u,
                                 equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGSEM, cache)
-    @unpack boundary_interpolation = dg.basis
+    @unpack boundary_interpolation_entries = dg.basis
     @unpack surface_flux_values = cache.elements
 
     # Note that all fluxes have been computed with outward-pointing normal vectors.
     # Access the factors only once before beginning the loop to increase performance.
     # 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_1 = boundary_interpolation_entries[1]
+    factor_2 = boundary_interpolation_entries[2]
     @threaded for element in eachelement(dg, cache)
         for m in eachnode(dg), l in eachnode(dg)
             for v in eachvariable(equations)