More abstract containers (#2632)

* More `AbstractContainers`

* fmt

* fmt

* fix

* bb

* bf

* bf

* bf

* rm

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: DanielDoehring

src/solvers/dgsem/containers.jl

@@ -0,0 +1,23 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+abstract type AbstractElementContainer <: AbstractContainer end
+function nelements end
+
+abstract type AbstractInterfaceContainer <: AbstractContainer end
+function ninterfaces end
+abstract type AbstractMPIInterfaceContainer <: AbstractContainer end
+function nmpiinterfaces end
+
+abstract type AbstractBoundaryContainer <: AbstractContainer end
+function nboundaries end
+
+abstract type AbstractMortarContainer <: AbstractContainer end
+function nmortars end
+abstract type AbstractMPIMortarContainer <: AbstractContainer end
+function nmpimortars end
+end # @muladd

src/solvers/dgsem/dgsem.jl

@@ -132,5 +132,6 @@ end
     return u_mean / total_volume # normalize with the total volume
 end
 
+include("containers.jl")
 include("calc_volume_integral.jl")
 end # @muladd

src/solvers/dgsem_p4est/containers.jl

@@ -14,7 +14,7 @@ mutable struct P4estElementContainer{NDIMS, RealT <: Real, uEltype <: Real,
                                      ArrayuEltypeNDIMSP2 <:
                                      DenseArray{uEltype, NDIMSP2},
                                      VectoruEltype <: DenseVector{uEltype}} <:
-               AbstractContainer
+               AbstractElementContainer
     # Physical coordinates at each node
     node_coordinates::ArrayRealTNDIMSP2 # [orientation, node_i, node_j, node_k, element]
 
@@ -222,7 +222,7 @@ mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, NDIMSP2,
                                        IdsVector <: DenseVector{Int},
                                        IndicesVector <:
                                        DenseVector{NTuple{NDIMS, Symbol}}} <:
-               AbstractContainer
+               AbstractInterfaceContainer
     u::uArray                   # [primary/secondary, variable, i, j, interface]
     neighbor_ids::IdsMatrix     # [primary/secondary, interface]
     node_indices::IndicesMatrix # [primary/secondary, interface]
@@ -344,7 +344,7 @@ mutable struct P4estBoundaryContainer{NDIMS, uEltype <: Real, NDIMSP1,
                                       IndicesVector <:
                                       DenseVector{NTuple{NDIMS, Symbol}},
                                       uVector <: DenseVector{uEltype}} <:
-               AbstractContainer
+               AbstractBoundaryContainer
     u::uArray                   # [variables, i, j, boundary]
     neighbor_ids::IdsVector     # [boundary]
     node_indices::IndicesVector # [boundary]
@@ -510,7 +510,7 @@ mutable struct P4estMortarContainer{NDIMS, uEltype <: Real, NDIMSP1, NDIMSP3,
                                     IdsVector <: DenseVector{Int},
                                     IndicesVector <:
                                     DenseVector{NTuple{NDIMS, Symbol}}} <:
-               AbstractContainer
+               AbstractMortarContainer
     u::uArray # [small/large side, variable, position, i, j, mortar]
     neighbor_ids::IdsMatrix # [position, mortar]
     node_indices::IndicesMatrix # [small/large, mortar]

src/solvers/dgsem_p4est/containers_parallel.jl

@@ -11,7 +11,7 @@ mutable struct P4estMPIInterfaceContainer{NDIMS, uEltype <: Real, NDIMSP2,
                                           IndicesVector <:
                                           DenseVector{NTuple{NDIMS, Symbol}},
                                           uVector <: DenseVector{uEltype}} <:
-               AbstractContainer
+               AbstractMPIInterfaceContainer
     u::uArray                   # [primary/secondary, variable, i, j, interface]
     local_neighbor_ids::VecInt  # [interface]
     node_indices::IndicesVector # [interface]
@@ -124,7 +124,7 @@ mutable struct P4estMPIMortarContainer{NDIMS, uEltype <: Real, RealT <: Real, ND
                                        NDIMSP2, NDIMSP3,
                                        uArray <: DenseArray{uEltype, NDIMSP3},
                                        uVector <: DenseVector{uEltype}} <:
-               AbstractContainer
+               AbstractMPIMortarContainer
     u::uArray                                      # [small/large side, variable, position, i, j, mortar]
     local_neighbor_ids::Vector{Vector{Int}}        # [mortar][ids]
     local_neighbor_positions::Vector{Vector{Int}}  # [mortar][positions]

src/solvers/dgsem_structured/containers.jl

@@ -6,7 +6,7 @@
 #! format: noindent
 
 struct StructuredElementContainer{NDIMS, RealT <: Real, uEltype <: Real,
-                                  NDIMSP1, NDIMSP2, NDIMSP3}
+                                  NDIMSP1, NDIMSP2, NDIMSP3} <: AbstractElementContainer
     # Physical coordinates at each node
     node_coordinates::Array{RealT, NDIMSP2} # [orientation, node_i, node_j, node_k, element]
     # ID of neighbor element in negative direction in orientation

src/solvers/dgsem_tree/containers.jl

@@ -8,7 +8,7 @@
 # Dimension independent code related to containers of the DG solver
 # with the mesh type TreeMesh
 
-abstract type AbstractTreeElementContainer <: AbstractContainer end
+abstract type AbstractTreeElementContainer <: AbstractElementContainer end
 
 # Return number of elements
 @inline nelements(elements::AbstractTreeElementContainer) = length(elements.cell_ids)
@@ -30,7 +30,7 @@ In particular, not the elements themselves are returned.
 @inline Base.real(elements::AbstractTreeElementContainer) = eltype(elements.node_coordinates)
 @inline Base.eltype(elements::AbstractTreeElementContainer) = eltype(elements.surface_flux_values)
 
-abstract type AbstractTreeInterfaceContainer <: AbstractContainer end
+abstract type AbstractTreeInterfaceContainer <: AbstractInterfaceContainer end
 
 # Return number of interfaces
 @inline ninterfaces(interfaces::AbstractTreeInterfaceContainer) = length(interfaces.orientations)
@@ -39,7 +39,19 @@ abstract type AbstractTreeInterfaceContainer <: AbstractContainer end
 # Return number of equation variables
 @inline nvariables(interfaces::AbstractTreeInterfaceContainer) = size(interfaces.u, 2)
 
-abstract type AbstractTreeBoundaryContainer <: AbstractContainer end
+abstract type AbstractTreeMPIInterfaceContainer <: AbstractMPIInterfaceContainer end
+
+# Return number of interfaces
+@inline function nmpiinterfaces(mpi_interfaces::AbstractTreeMPIInterfaceContainer)
+    return length(mpi_interfaces.orientations)
+end
+# Return number of interface nodes for 2D and 3D. For 1D hard-coded to 1 interface node.
+@inline nnodes(interfaces::AbstractTreeMPIInterfaceContainer) = size(interfaces.u, 3)
+# Return number of equation variables
+@inline nvariables(interfaces::AbstractTreeMPIInterfaceContainer) = size(interfaces.u,
+                                                                         2)
+
+abstract type AbstractTreeBoundaryContainer <: AbstractBoundaryContainer end
 
 # Return number of boundaries
 @inline nboundaries(boundaries::AbstractTreeBoundaryContainer) = length(boundaries.orientations)
@@ -48,11 +60,18 @@ abstract type AbstractTreeBoundaryContainer <: AbstractContainer end
 # Return number of equation variables
 @inline nvariables(boundaries::AbstractTreeBoundaryContainer) = size(boundaries.u, 2)
 
-abstract type AbstractTreeL2MortarContainer <: AbstractContainer end
+abstract type AbstractTreeL2MortarContainer <: AbstractMortarContainer end
 
 # Return number of L2 mortars
 @inline nmortars(l2mortars::AbstractTreeL2MortarContainer) = length(l2mortars.orientations)
 
+abstract type AbstractTreeL2MPIMortarContainer <: AbstractMPIMortarContainer end
+
+# Return number of L2 mortars
+@inline function nmpimortars(mpi_l2mortars::AbstractTreeL2MPIMortarContainer)
+    length(mpi_l2mortars.orientations)
+end
+
 function reinitialize_containers!(mesh::TreeMesh, equations, dg::DGSEM, cache)
     # Get new list of leaf cells
     leaf_cell_ids = local_leaf_cells(mesh.tree)

src/solvers/dgsem_tree/containers_2d.jl

@@ -471,13 +471,6 @@ function init_boundaries!(boundaries, elements, mesh::TreeMesh2D)
     return boundaries.n_boundaries_per_direction
 end
 
-abstract type AbstractTreeL2MortarContainer2D <: AbstractTreeL2MortarContainer end
-
-# Return number of mortar nodes (L2 mortars are only h-adaptive, not p-adaptive)
-@inline nnodes(mortars::AbstractTreeL2MortarContainer2D) = size(mortars.u_upper, 3)
-# Return number of equation variables
-@inline nvariables(mortars::AbstractTreeL2MortarContainer2D) = size(mortars.u_upper, 2)
-
 # Container data structure (structure-of-arrays style) for DG L2 mortars
 # Positions/directions for orientations = 1, large_sides = 2:
 # mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
@@ -489,7 +482,7 @@ abstract type AbstractTreeL2MortarContainer2D <: AbstractTreeL2MortarContainer e
 # lower = 1 |    |
 #           |    |
 mutable struct TreeL2MortarContainer2D{uEltype <: Real} <:
-               AbstractTreeL2MortarContainer2D
+               AbstractTreeL2MortarContainer
     u_upper::Array{uEltype, 4}  # [leftright, variables, i, mortars]
     u_lower::Array{uEltype, 4}  # [leftright, variables, i, mortars]
     neighbor_ids::Array{Int, 2} # [position, mortars]
@@ -502,6 +495,11 @@ mutable struct TreeL2MortarContainer2D{uEltype <: Real} <:
     _neighbor_ids::Vector{Int}
 end
 
+# Return number of mortar nodes (L2 mortars are only h-adaptive, not p-adaptive)
+@inline nnodes(mortars::TreeL2MortarContainer2D) = size(mortars.u_upper, 3)
+# Return number of equation variables
+@inline nvariables(mortars::TreeL2MortarContainer2D) = size(mortars.u_upper, 2)
+
 # See explanation of Base.resize! for the element container
 function Base.resize!(mortars::TreeL2MortarContainer2D, capacity)
     n_nodes = nnodes(mortars)
@@ -741,7 +739,7 @@ end
 
 # Container data structure (structure-of-arrays style) for DG MPI interfaces
 mutable struct TreeMPIInterfaceContainer2D{uEltype <: Real} <:
-               AbstractTreeInterfaceContainer
+               AbstractTreeMPIInterfaceContainer
     u::Array{uEltype, 4}            # [leftright, variables, i, interfaces]
     # Note: `local_neighbor_ids` stores the MPI-local neighbors, but with globally valid index!
     local_neighbor_ids::Vector{Int} # [interfaces]
@@ -789,12 +787,6 @@ function TreeMPIInterfaceContainer2D{uEltype}(capacity::Integer, n_variables,
                                                 remote_sides, _u)
 end
 
-# TODO: Taal, rename to ninterfaces?
-# Return number of interfaces
-@inline function nmpiinterfaces(mpi_interfaces::TreeMPIInterfaceContainer2D)
-    length(mpi_interfaces.orientations)
-end
-
 # Create MPI interface container and initialize MPI interface data in `elements`.
 function init_mpi_interfaces(cell_ids, mesh::TreeMesh2D,
                              elements::TreeElementContainer2D)
@@ -913,7 +905,7 @@ end
 # lower = 1 |    |
 #           |    |
 mutable struct TreeMPIL2MortarContainer2D{uEltype <: Real} <:
-               AbstractTreeL2MortarContainer2D
+               AbstractTreeL2MPIMortarContainer
     u_upper::Array{uEltype, 4} # [leftright, variables, i, mortars]
     u_lower::Array{uEltype, 4} # [leftright, variables, i, mortars]
     # Note: `local_neighbor_ids` stores the MPI-local neighbors, but with globally valid index!
@@ -927,6 +919,11 @@ mutable struct TreeMPIL2MortarContainer2D{uEltype <: Real} <:
     _u_lower::Vector{uEltype}
 end
 
+# Return number of mortar nodes (L2 mortars are only h-adaptive, not p-adaptive)
+@inline nnodes(mortars::TreeMPIL2MortarContainer2D) = size(mortars.u_upper, 3)
+# Return number of equation variables
+@inline nvariables(mortars::TreeMPIL2MortarContainer2D) = size(mortars.u_upper, 2)
+
 # See explanation of Base.resize! for the element container
 function Base.resize!(mpi_mortars::TreeMPIL2MortarContainer2D, capacity)
     n_nodes = nnodes(mpi_mortars)
@@ -979,11 +976,6 @@ function TreeMPIL2MortarContainer2D{uEltype}(capacity::Integer, n_variables,
                                                _u_upper, _u_lower)
 end
 
-# Return number of L2 mortars
-@inline function nmpimortars(mpi_l2mortars::TreeMPIL2MortarContainer2D)
-    length(mpi_l2mortars.orientations)
-end
-
 # Create MPI mortar container and initialize MPI mortar data in `elements`.
 function init_mpi_mortars(cell_ids, mesh::TreeMesh2D,
                           elements::TreeElementContainer2D,

src/solvers/dgsem_unstructured/containers_2d.jl

@@ -6,7 +6,8 @@
 #! format: noindent
 
 # Container data structure (structure-of-arrays style) for DG elements on curved unstructured mesh
-struct UnstructuredElementContainer2D{RealT <: Real, uEltype <: Real}
+struct UnstructuredElementContainer2D{RealT <: Real, uEltype <: Real} <:
+       AbstractElementContainer
     node_coordinates::Array{RealT, 4}      # [ndims, nnodes, nnodes, nelement]
     jacobian_matrix::Array{RealT, 5}       # [ndims, ndims, nnodes, nnodes, nelement]
     inverse_jacobian::Array{RealT, 3}      # [nnodes, nnodes, nelement]
@@ -117,12 +118,13 @@ function init_element!(elements, element, basis::LobattoLegendreBasis,
 end
 
 # generic container for the interior interfaces of an unstructured mesh
-struct UnstructuredInterfaceContainer2D{uEltype <: Real}
-    u::Array{uEltype, 4} # [primary/secondary, variables, i, interfaces]
-    start_index::Vector{Int}       # [interfaces]
-    index_increment::Vector{Int}       # [interfaces]
-    element_ids::Array{Int, 2}     # [primary/secondary, interfaces]
-    element_side_ids::Array{Int, 2}     # [primary/secondary, interfaces]
+struct UnstructuredInterfaceContainer2D{uEltype <: Real} <:
+       AbstractInterfaceContainer
+    u::Array{uEltype, 4}            # [primary/secondary, variables, i, interfaces]
+    start_index::Vector{Int}        # [interfaces]
+    index_increment::Vector{Int}    # [interfaces]
+    element_ids::Array{Int, 2}      # [primary/secondary, interfaces]
+    element_side_ids::Array{Int, 2} # [primary/secondary, interfaces]
 end
 
 # Construct an empty curved interface container to be filled later with neighbour
@@ -255,12 +257,13 @@ end
 
 # TODO: Clean-up meshes. Find a better name since it's also used for other meshes
 # generic container for the boundary interfaces of an unstructured mesh
-struct UnstructuredBoundaryContainer2D{RealT <: Real, uEltype <: Real}
-    u::Array{uEltype, 3} # [variables, i, boundaries]
-    element_id::Vector{Int}       # [boundaries]
-    element_side_id::Vector{Int}       # [boundaries]
-    node_coordinates::Array{RealT, 3}   # [ndims, nnodes, boundaries]
-    name::Vector{Symbol}    # [boundaries]
+struct UnstructuredBoundaryContainer2D{RealT <: Real, uEltype <: Real} <:
+       AbstractBoundaryContainer
+    u::Array{uEltype, 3}              # [variables, i, boundaries]
+    element_id::Vector{Int}           # [boundaries]
+    element_side_id::Vector{Int}      # [boundaries]
+    node_coordinates::Array{RealT, 3} # [ndims, nnodes, boundaries]
+    name::Vector{Symbol}              # [boundaries]
 end
 
 # construct an empty curved boundary container to be filled later with neighbour

test/test_unstructured_2d.jl

@@ -28,12 +28,6 @@ isdir(outdir) && rm(outdir, recursive = true)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
-    semi32 = Trixi.trixi_adapt(Array, Float32, semi)
-    @test real(semi32.solver) == Float32
-    @test real(semi32.solver.basis) == Float32
-    @test real(semi32.solver.mortar) == Float32
-    # TODO: remake ignores the mesh as well
-    @test real(semi32.mesh) == Float64
 end
 
 @trixi_testset "elixir_euler_free_stream.jl" begin