Use PtrArray for contravariant_vectors

amrueda · Dec 6, 2023 · 245b20d · 245b20d
1 parent d2b5294
commit 245b20d
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 26 deletions.
diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -48,7 +48,7 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
 using IfElse: ifelse
 using LinearMaps: LinearMap
 using LoopVectorization: LoopVectorization, @turbo, indices
-using StaticArrayInterface: static_length # used by LoopVectorization
+using StaticArrayInterface: static_length, static_size # used by LoopVectorization
 using MuladdMacro: @muladd
 using Octavian: Octavian, matmul!
 using Polyester: Polyester, @batch # You know, the cheapest threads you can find...

diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl
@@ -13,7 +13,7 @@ mutable struct P4estElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIM
     # [jacobian_i, jacobian_j, node_i, node_j, node_k, element] where jacobian_i is the first index of the Jacobian matrix,...
     jacobian_matrix::Array{RealT, NDIMSP3}
     # Contravariant vectors, scaled by J, in Kopriva's blue book called Ja^i_n (i index, n dimension)
-    contravariant_vectors::Array{RealT, NDIMSP3}   # [dimension, index, node_i, node_j, node_k, element]
+    contravariant_vectors::PtrArray{RealT, NDIMSP3}   # [dimension, index, node_i, node_j, node_k, element]
     # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix)
     inverse_jacobian::Array{RealT, NDIMSP1}   # [node_i, node_j, node_k, element]
     # Buffer for calculated surface flux
@@ -106,10 +106,10 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re
     _contravariant_vectors = Vector{RealT}(undef,
                                            ndims_spa^2 * nnodes(basis)^NDIMS *
                                            nelements)
-    contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors),
-                                        (ndims_spa, ndims_spa,
-                                         ntuple(_ -> nnodes(basis), NDIMS)...,
-                                         nelements))
+    contravariant_vectors = PtrArray(pointer(_contravariant_vectors),
+                                     (StaticInt(ndims_spa), ndims_spa,
+                                      ntuple(_ -> nnodes(basis), NDIMS)...,
+                                      nelements))
 
     _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements)
     inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian),

diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl
@@ -22,12 +22,12 @@ function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG
 end
 
 # Extract contravariant vector Ja^i (i = index) as SVector
-# TODO: Here, size() causes a lot of allocations! Find an alternative to improve performance
 @inline function get_contravariant_vector(index, contravariant_vectors, indices...)
     SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]),
-                   Val(size(contravariant_vectors, 1))))
+                   Val(static_size(contravariant_vectors, 1))))
 end
 
+
 @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
                                                   orientation,
                                                   boundary_condition::BoundaryConditionPeriodic,

diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl
@@ -75,9 +75,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17
         flux2 = flux(u_node, 2, equations)
         # Compute the contravariant flux by taking the scalar product of the
         # first contravariant vector Ja^1 and the flux vector
-        #Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
-        Ja11 = contravariant_vectors[1, 1, i, j, element]
-        Ja12 = contravariant_vectors[2, 1, i, j, element]
+        Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
         contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
         for ii in eachnode(dg)
             multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i],
@@ -87,9 +85,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17
 
         # Compute the contravariant flux by taking the scalar product of the
         # second contravariant vector Ja^2 and the flux vector
-        #Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
-        Ja21 = contravariant_vectors[1, 2, i, j, element]
-        Ja22 = contravariant_vectors[2, 2, i, j, element]
+        Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
         contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2
         for jj in eachnode(dg)
             multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j],
@@ -122,10 +118,7 @@ end
 
         # Compute the contravariant flux by taking the scalar product of the
         # first contravariant vector Ja^1 and the flux vector
-        #Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
-        Ja11 = contravariant_vectors[1, 1, i, j, element]
-        Ja12 = contravariant_vectors[2, 1, i, j, element]
-        Ja13 = contravariant_vectors[3, 1, i, j, element]
+        Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
         contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3
         for ii in eachnode(dg)
             multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i],
@@ -135,21 +128,15 @@ end
 
         # Compute the contravariant flux by taking the scalar product of the
         # second contravariant vector Ja^2 and the flux vector
-        #Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
-        Ja21 = contravariant_vectors[1, 2, i, j, element]
-        Ja22 = contravariant_vectors[2, 2, i, j, element]
-        Ja23 = contravariant_vectors[3, 2, i, j, element]
+        Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
         contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3
         for jj in eachnode(dg)
             multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j],
                                        contravariant_flux2, equations, dg, i, jj,
                                        element)
         end
 
-        #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element)
-        Ja31 = contravariant_vectors[1, 3, i, j, element]
-        Ja32 = contravariant_vectors[2, 3, i, j, element]
-        Ja33 = contravariant_vectors[3, 3, i, j, element]
+        Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element)
         for v in eachvariable(equations)
             du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] +
                                      Ja33 * flux3[v])