Merge branch 'main' into sc/p4est-view-enhanced

trixi-framework · Jan 20, 2025 · 663a435 · 663a435
2 parents 369ef54 + 6149104
commit 663a435
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 43 deletions.
diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl
@@ -40,31 +40,27 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)
     nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT)
     inverse_weights_ = inv.(weights_)
 
-    _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_, RealT)
+    _, 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 = 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_')
-    derivative_dhat_ = calc_dhat(nodes_, weights_)
+    derivative_matrix = polynomial_derivative_matrix(nodes_)
+    derivative_split = calc_dsplit(nodes_, weights_)
+    derivative_split_transpose = Matrix(derivative_split')
+    derivative_dhat = calc_dhat(nodes_, weights_)
 
-    # type conversions to get the requested real type and enable possible
-    # optimizations of runtime performance and latency
+    # Type conversions to enable possible optimizations of runtime performance 
+    # and latency
     nodes = SVector{nnodes_, RealT}(nodes_)
     weights = SVector{nnodes_, RealT}(weights_)
     inverse_weights = SVector{nnodes_, RealT}(inverse_weights_)
 
-    inverse_vandermonde_legendre = convert.(RealT, inverse_vandermonde_legendre_)
-    boundary_interpolation = convert.(RealT, boundary_interpolation_)
-
-    # Usually as fast as `SMatrix` (when using `let` in the volume integral/`@threaded`)
-    derivative_matrix = Matrix{RealT}(derivative_matrix_)
-    derivative_split = Matrix{RealT}(derivative_split_)
-    derivative_split_transpose = Matrix{RealT}(derivative_split_transpose_)
-    derivative_dhat = Matrix{RealT}(derivative_dhat_)
+    # We keep the matrices above stored using the standard `Matrix` type 
+    # since this is usually as fast as `SMatrix`
+    # (when using `let` in the volume integral/`@threaded`)
+    # and reduces latency
 
     return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes),
                                 typeof(inverse_vandermonde_legendre),
@@ -163,17 +159,15 @@ function MortarL2(basis::LobattoLegendreBasis)
     RealT = real(basis)
     nnodes_ = nnodes(basis)
 
-    forward_upper_ = calc_forward_upper(nnodes_, RealT)
-    forward_lower_ = calc_forward_lower(nnodes_, RealT)
-    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
-    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
-
-    # type conversions to get the requested real type and enable possible
-    # optimizations of runtime performance and latency
+    forward_upper = calc_forward_upper(nnodes_, RealT)
+    forward_lower = calc_forward_lower(nnodes_, RealT)
+    reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
+    reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
 
-    # Usually as fast as `SMatrix` but better for latency
-    forward_upper = Matrix{RealT}(forward_upper_)
-    forward_lower = Matrix{RealT}(forward_lower_)
+    # We keep the matrices above stored using the standard `Matrix` type 
+    # since this is usually as fast as `SMatrix`
+    # (when using `let` in the volume integral/`@threaded`)
+    # and reduces latency
 
     # TODO: Taal performance
     #       Check the performance of different implementations of `mortar_fluxes_to_elements!`
@@ -183,8 +177,6 @@ function MortarL2(basis::LobattoLegendreBasis)
     #       `@tullio` when the matrix sizes are not necessarily static.
     # reverse_upper = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_upper_)
     # reverse_lower = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_lower_)
-    reverse_upper = Matrix{RealT}(reverse_upper_)
-    reverse_lower = Matrix{RealT}(reverse_lower_)
 
     LobattoLegendreMortarL2{RealT, nnodes_, typeof(forward_upper),
                             typeof(reverse_upper)}(forward_upper, forward_lower,

diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl
@@ -22,7 +22,7 @@ Create a discontinuous Galerkin spectral element method (DGSEM) using a
 """
 const DGSEM = DG{Basis} where {Basis <: LobattoLegendreBasis}
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(basis::LobattoLegendreBasis,
                surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm(),
@@ -32,7 +32,7 @@ function DGSEM(basis::LobattoLegendreBasis,
               typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral)
 end
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(basis::LobattoLegendreBasis,
                surface_integral::AbstractSurfaceIntegral,
                volume_integral = VolumeIntegralWeakForm(),
@@ -41,7 +41,7 @@ function DGSEM(basis::LobattoLegendreBasis,
               typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral)
 end
 
-# TODO: Deprecated in v0.3 (no longer documented)
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(RealT, polydeg::Integer,
                surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm(),
@@ -51,6 +51,7 @@ function DGSEM(RealT, polydeg::Integer,
     return DGSEM(basis, surface_flux, volume_integral, mortar)
 end
 
+# This API is no longer documented, and we recommend avoiding its public use.
 function DGSEM(polydeg, surface_flux = flux_central,
                volume_integral = VolumeIntegralWeakForm())
     DGSEM(Float64, polydeg, surface_flux, volume_integral)

diff --git a/src/solvers/dgsem/l2projection.jl b/src/solvers/dgsem/l2projection.jl
@@ -29,9 +29,9 @@ function calc_forward_upper(n_nodes, RealT = Float64)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] + 1), nodes, wbary)
+        poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] + 1), nodes, wbary)
         for i in 1:n_nodes
             operator[j, i] = poly[i]
         end
@@ -49,9 +49,9 @@ function calc_forward_lower(n_nodes, RealT = Float64)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] - 1), nodes, wbary)
+        poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] - 1), nodes, wbary)
         for i in 1:n_nodes
             operator[j, i] = poly[i]
         end
@@ -70,12 +70,12 @@ function calc_reverse_upper(n_nodes, ::Val{:gauss}, RealT = Float64)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] + 1),
+        poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] + 1),
                                                   gauss_nodes, gauss_wbary)
         for i in 1:n_nodes
-            operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i]
+            operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i]
         end
     end
 
@@ -97,12 +97,12 @@ function calc_reverse_lower(n_nodes, ::Val{:gauss}, RealT = Float64)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
-    operator = zeros(n_nodes, n_nodes)
+    operator = zeros(RealT, n_nodes, n_nodes)
     for j in 1:n_nodes
-        poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] - 1),
+        poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] - 1),
                                                   gauss_nodes, gauss_wbary)
         for i in 1:n_nodes
-            operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i]
+            operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i]
         end
     end