Fix type conversion and type instability in MortarL2 (trixi-framework#2242)

* Fix

* Fix

Author: huiyuxie

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -159,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)
+    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
-
-    # 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!`
@@ -179,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,

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