diff --git a/src/solvers/dgsem/containers.jl b/src/solvers/dgsem/containers.jl new file mode 100644 index 00000000000..a2c05b80b0b --- /dev/null +++ b/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 diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 979f2d0577d..546fffde840 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/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 diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 3f74f699f19..2f660bd4bca 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/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] diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 4b6fc813703..ba46c8f7c20 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/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] diff --git a/src/solvers/dgsem_structured/containers.jl b/src/solvers/dgsem_structured/containers.jl index 839c917749b..032e2286a48 100644 --- a/src/solvers/dgsem_structured/containers.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/containers.jl b/src/solvers/dgsem_tree/containers.jl index dd0d10c26e6..648bfbc9cb6 100644 --- a/src/solvers/dgsem_tree/containers.jl +++ b/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) diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 0c477ae396f..54b0263cb67 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/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, diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index eec5df5f9de..b12c4ed0039 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/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 diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 61e35cdf563..7f06ee3e75e 100644 --- a/test/test_unstructured_2d.jl +++ b/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